Diyar, Francesco C. (for, perhaps unknowingly, alleviating my postdoc anxiety), Francesco M. (for being the extravagant event organizer that he is), Greivin, Jordan, Jules, Léo, Mathis, Maxl (for the late night conversations), Pierre (for the entertaining stories), Simon, Tom, Vincent (for keeping the conversations going when Andrei is not there), Wenki, Yann, and Yehudi (for somehow always knowing how to make sense of the French administration). A special mention goes to Thorsten, for being everyone's favorite postdoc. Outside this little world of LPTHE, I am glad for the many movie nights we organized with Erik and Anna since our ICFP days. Last, but not least, this moral support group concludes with Pulkit, Ro, and Bianca. I never imagined at the start of the PhD that this list could become so long. I am

R E S U M É

Le sujet principal de cette thèse est la dynamique périodique en temps des jonctions supraconductrices. En particulier, on se concentre sur l'étude des propriétés spectrales, ainsi que de transport, de boîtes quantiques (QD) couplées à de multiples réservoirs supraconducteurs (S). De tels systèmes supraconducteurs hybrides ont pour caractéristique que la périodicité temporelle peut être simplement réalisée par l'application de tensions constantes (et commensurables) sur les réservoirs. En conséquence de cette périodicité, les états liés d'Andreev, localisés sur la boîte à l'équilibre (en l'absence de tension appliquée), se transforment en des résonances de Floquet-Andreev. Ces résonances ont une largeur en énergie finie, à cause de la présence de réflexions d'Andreev multiples et forment des échelles de Wannier-Stark. À partir des équations de Bogoliubov-de Gennes, périodiques en temps, on utilise la théorie de Floquet pour obtenir un modèle de liaisons fortes dans l'espace de Fourier. Une méthode de fraction continue est ensuite utilisée pour obtenir la résolvante du système, à partir de laquelle les observables dans l'état stationnaire peuvent être calculées. Nous calculons numériquement les spectres de Floquet, ainsi que la composante continue du courant de la jonction S-QD-S et de la bijonction S-QD-S-QD-S, également connue sous le nom de "molécule d'Andreev". Nous montrons que la fonction spectrale peut être mesurée par spectroscopie à effet tunnel. Nous examinons ensuite la bijonction en dehors du régime moléculaire. Nous constatons que le pilotage périodique induit un couplage à longue portée entre les boîtes quantiques, entraînant un comportement interférométrique. Enfin, nous présentons une méthode de dérivation des équations maîtresses markoviennes pour des systèmes fermioniques sans interactions, reposant sur une approximation de résonances étroites. Nous appliquons la méthode à une boîte quantique à niveaux multiples couplée à des réservoirs métalliques, ainsi qu'à une boîte quantique soumise à un pilotage périodique. Pour une boîte quantique supraconductrice, nous dérivons une expression pour les populations des états de Floquet-Andreev.

A B S T R A C T

This thesis explores the time-periodic dynamics of superconducting junctions.

In particular, we study the spectral and transport properties of quantum dot (QD) systems connected to multiple superconducting (S) reservoirs. The unique feature of such systems is that time-periodicity can be realized by the application of constant and commensurate voltage bias. As a result of the periodic driving, the equilibrium Andreev bound states, localized on the dot, turn into nonequilibrium Floquet-Andreev resonances. Using Floquet theory, we map the time-periodic Bogoliubov-de Gennes equations to a tight-binding chain in Fourier space. A continued fraction method is then used to obtain the resolvent of the system, from which observables in the steady state can be calculated. We consider two observables: the dc component of the current and the time-averaged spectral function, which we show can be probed by tunneling spectroscopy. We numerically calculate the Floquet spectra, as well as the dc component of the current, of the S-QD-S junction and the S-QD-S-QD-S bijunction, also known as the "Andreev molecule". We then consider the bijunction away from the molecular regime. We find that periodic driving induces a longrange coupling between the dots, resulting in interferometric behavior. Finally, we present a method for the derivation of Markovian master equations for noninteracting fermionic systems, relying on a narrow resonance approximation. We apply the method to a multilevel QD coupled to metallic reservoirs, as well as to a periodically driven QD. For a superconducting QD, we derive an expression for the populations of the Floquet-Andreev states.

A C K N O W L E D G M E N T S

There are many people without whom this thesis could neither have started, nor kept going, nor finished. First and foremost, I must thank my supervisor, Benoît Douçot, who belongs in all three of these categories, for his patience and his approach to physics, which is, to say the least, inspiring. I owe a special thanks to my two rapporteurs, Fabio Pistolesi and Landry Bretheau, for being incredibly kind and taking the time to carefully read this manuscript, as well as to the members of the jury, Tristan Cren and Cristina Bena. I am especially grateful to Gregory Schehr and Ça glar Girit for making useful suggestions at every meeting of the comité de suivi. In particular, I am grateful to Ça glar Girit for proposing that we look at the MAR current. Going further back in time, I would like to thank the professors at the physics department of the university of Crete, where I first started my studies. In particular, I would like to thank Gregory Psaltakis and the late Theodore Tomaras, for giving me encouragement when I needed it the most. The past three years could have hardly been any fun if it were not for the fellow CDD members of the LPTHE. I consider myself lucky to have been a part of this wonderful group of people: Andrei (for his sense of humor), Anthony, Floquet theory is concerned with systems which are periodic in time. It was introduced towards the end of the 19th century by G. Floquet for the study of linear systems of differential equations with time-periodic coefficients [1]. At that time, the main motivation was the study of classical systems, particularly in the context of celestial mechanics and the problem of the stability of periodic orbits [2,[START_REF] Poincaré | Les méthodes nouvelles de la mécanique céleste[END_REF]. In a quantum mechanical context, the time-dependent Schrödinger equation is always linear, and a natural application of Floquet theory arises whenever the Hamiltonian of the system is time periodic

ih d dt |ψ(t)⟩ = H(t) |ψ(t)⟩ , with H(t + T) = H(t), (1.1)
where T = 2π ω 0 is the period.

floquet engineering is currently one of the principal motivations for studying Floquet systems. The term refers to the fact that it is possible to use periodic driving, such as irradiation by a laser, to control a quantum system and fabricate its properties on demand [START_REF] Holthaus | Floquet engineering with quasienergy bands of periodically driven optical lattices[END_REF][START_REF] Marin Bukov | Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to Floquet engineering[END_REF][START_REF] Basov | Towards properties on demand in quantum materials[END_REF]. Essentially, the driving enlarges the parameter space of the system, providing an additional "experimental knob" which can be used to realize states that are not accessible in equilibrium. A paradigmatic example that demonstrates the principle of Floquet engineering in the classical world is Kapitza's pendulum [START_REF] Kapitza | Dynamic stability of the pendulum with vibrating suspension point[END_REF]: a rigid pendulum whose suspension point oscillates vertically. At a high enough vibration frequency, the upright position of the pendulum is dynamically stabilized, illustrating that the states of a driven system can become very different from those in equilibrium.

There is renewed interest in the study of periodically driven quantum systems, owing both to theoretical progress and to developments in experimental techniques [START_REF] Oka | Floquet engineering of quantum materials[END_REF]. For example, periodic driving has been used to open band gaps in graphene [START_REF] Oka | Photovoltaic Hall effect in graphene[END_REF][START_REF] Broers | Observing light-induced Floquet band gaps in the longitudinal conductivity of graphene[END_REF], to induce topological properties in non-topological materials [START_REF] Lindner | Floquet topological insulator in semiconductor quantum wells[END_REF][START_REF] Mciver | Light-induced anomalous Hall effect in graphene[END_REF], or to realize discrete-time crystals [START_REF] Else | Floquet time crystals[END_REF][START_REF] Zhang | Observation of a discrete time crystal[END_REF]. Moreover, advances in ultrafast spectroscopy have allowed the spectacular imaging of the Floquet bandstructure of a topological insulator in real time [START_REF] Wang | Observation of Floquet-Bloch states on the surface of a topological insulator[END_REF][START_REF] Mahmood | Selective scattering between Floquet-Bloch and Volkov states in a topological insulator[END_REF].

Basics of the Floquet method

The Floquet theorem guarantees that the solutions of Eq. (1.1) have the following form ψ j (t) = e -iε j t/h u j (t) , with u j (t + T) = u j (t) (1.2) i.e., they can be written as a product of a phase factor and a time-periodic function with the same period as the Hamiltonian. Therefore, the Floquet theorem is the temporal analog of the Bloch theorem for systems with spatial periodicity.

Inserting the solution (1.2) into the time-dependent Schrödinger equation,

H(t) -ih d dt u j (t) = ε j u j (t) , (1.3)
one finds that the Floquet states u j (t) are eigenstates of the operator H(t)ih d dt with eigenvalue ε j . However, it is straightforward to verify that one can equivalently write H(t)ih d dt e imω 0 t u j (t) = (ε j + mhω 0 )e imω 0 t u j (t) , (1.4) with ω 0 = 2π/T being the driving frequency. As a result, the quasienergies ε j are only defined up to integer multiples of the frequency ω 0 [START_REF] Zel'dovich | The quasienergy of a quantum-mechanical system subjected to a periodic action[END_REF]. The Shirley approach to the time-dependent problem [START_REF] Shirley | Solution of the Schrödinger equation with a Hamiltonian periodic in time[END_REF] is to expand the time-periodic quantities in Fourier series in an extended Hilbert space H ⊗ T, consisting of a tensor product of the Hilbert space of the system H with the vector space T of time-periodic functions with period T [START_REF] Sambe | Steady states and quasienergies of a quantum-mechanical system in an oscillating field[END_REF]. Written in a matrix representation, the effective Floquet Hamiltonian is a matrix of infinite dimensions

H(t) = ∑
        . . . H 0 -hω 0 H 1 H 2 H -1 H 0 H 1 H -2 H -1 H 0 + hω 0 . . .                  . . . u 1 j u 0 j u -1 j . . .          = ε j          . . . u 1 j u 0 j u -1 j . . .          (1.7)
1.1 floquet systems where each element of the matrix is a block with a dimension equal to the dimension of H. Equation (1.6) means that the initial time-dependent problem has been mapped to a static tight-binding model where time has been traded for a fictional dimension. In this picture, the Fourier modes m, n can be interpreted as positions on the tight-binding lattice, with off-diagonal terms H n̸ =0 describing the hopping between different sites by absorption or emission of n "photons", and H 0 describing an on-site energy. Moreover, the term mhω 0 is analogous to the existence of a fictitious electric field [START_REF] Oka | Floquet engineering of quantum materials[END_REF]. The presence of this term implies that a high frequency produces energetically separated states. This permits a simplification of the problem, since the infinite matrix can be drastically truncated.

Open and closed Floquet systems

Floquet systems are a special class of non-equilibrium systems and, as such, not described by the usual equilibrium thermodynamic rules [START_REF] Kohn | Periodic thermodynamics[END_REF][START_REF] Hone | Statistical mechanics of Floquet systems: the pervasive problem of near degeneracies[END_REF]. For example, the notion of a Floquet ground-state is ill-defined, since the quasienergies play the role of phases and, consequently, cannot be uniquely ordered. Another point that is not trivial and has recently attracted attention is how to determine the populations of Floquet states [START_REF] Iadecola | Occupation of topological Floquet bands in open systems[END_REF][START_REF] Seetharam | Controlled population of Floquet-Bloch states via coupling to Bose and Fermi baths[END_REF][START_REF] Matsyshyn | Fermi-Dirac staircase occupation of Floquet bands and current rectification inside the optical gap of metals: an exact approach[END_REF][START_REF] Shi | The Floquet Fermi liquid[END_REF].

At this point, one should make the distinction between open and closed Floquet systems: In a generic closed system, periodic driving supplies energy that is absorbed until the system thermalizes (heats to an infinite temperature) [START_REF] Lazarides | Equilibrium states of generic quantum systems subject to periodic driving[END_REF]. That is not to say that closed Floquet systems are without interest. Indeed, heating can be exponentially slow, particularly in high-frequency regimes, resulting in long-lived transient states [START_REF] Abanin | Effective Hamiltonians, prethermalization, and slow energy absorption in periodically driven many-body systems[END_REF]. Moreover, the presence of disorder can result in many-body localization which can also lead to an absence of thermalization [START_REF] Khemani | Phase structure of driven quantum systems[END_REF][START_REF] Abanin | Theory of many-body localization in periodically driven systems[END_REF].

A completely different situation arises if the system is open, meaning that it is in contact with a reservoir. Then, the presence of dissipation to the environment acts as a cooling mechanism. It can therefore provide a balance to the periodic drive and, in many cases, the system reaches a non-equilibrium steady state [START_REF] Iadecola | Occupation of topological Floquet bands in open systems[END_REF][START_REF] Seetharam | Controlled population of Floquet-Bloch states via coupling to Bose and Fermi baths[END_REF][START_REF] Matsyshyn | Fermi-Dirac staircase occupation of Floquet bands and current rectification inside the optical gap of metals: an exact approach[END_REF][START_REF] Shi | The Floquet Fermi liquid[END_REF][START_REF] Mori | Floquet states in open quantum systems[END_REF][START_REF] Tsuji | Floquet states[END_REF][START_REF] Park | Steady Floquet-Andreev states in graphene Josephson junctions[END_REF]. Open systems naturally arise in platforms for the study of quantum transport, where a small system is connected to two or more electronic reservoirs. Non-equilibrium states then occur by, for example, keeping the reservoirs at different chemical potentials. In this thesis, we will focus on the case where the reservoirs are superconducting. Due to the Josephson relation (to be defined shortly), periodic driving can then be realized by application of a dc voltage bias.

superconducting devices

In recent years there has been a surge of interest in superconducting devices and circuits. Interest is both theoretical and practical, since, for example, such devices are potentially useful for creating superconducting qubits [START_REF] Clarke | Superconducting quantum bits[END_REF][START_REF] Leijnse | Coupling spin qubits via superconductors[END_REF][START_REF] Kjaergaard | Superconducting qubits: current state of play[END_REF]. Their elementary building block is the Josephson junction, made by two superconductors separated by an insulating material (Fig. 1.1a), and their operation is based on the Josephson effect.

In the words of Anderson [START_REF] Anderson | More is different[END_REF], it was understanding of the importance of broken symmetry in superconductivity that led to Josephson's discovery. A superconductor spontaneously breaks a global U(1) symmetry in the sense that its wavefunction has a preferred macroscopic phase. A consequence of this symmetry breaking is that charge is not well-defined and one needs to make a superposition of electron and hole states to build an excitation. Then, spatial variations of the pair potential lead to a new scattering process, where an electron-like excitation is scattered into a hole-like excitation. This process was independently identified by Saint-James and de Gennes [START_REF] De Gennes | Elementary excitations in the vicinity of a normal metal-superconducting metal contact[END_REF][START_REF] Saint-James | Excitations élémentaires au voisinage de la surface de séparation d'un métal normal et d'un métal supraconducteur[END_REF] and by Andreev [START_REF] Andreev | Thermal conductivity of the intermediate state in superconductors[END_REF][START_REF] Andreev | Thermal conductivity of the intermediate state of superconductors II[END_REF], and goes by the name of Andreev reflection. This process is of particular importance in superconducting junctions, where it is used for the microscopic description of the Josephson effect. Most notably, it leads to the formation of bound states in the weak-link region of the junction, with energies below the superconducting gap. These are called Andreev bound states.

Josephson effect

In 1962, Josephson predicted the existence of a supercurrent due to the coherent tunneling of Cooper pairs between two superconductors separated by a thin insulating layer [START_REF] Josephson | Possible new effects in superconductive tunnelling[END_REF]. The phenomenon can be summarized by the following two equations:

I(t) = I c sin(ϕ(t)) (1.8a) dϕ dt = 2e h V(t) (1.8b)
where ϕ is the phase difference between the two superconductors. In the absence of any applied voltage, the phase difference does not change with time and there is a constant dissipationless current (a supercurrent) flowing through the junction whose amplitude depends on the phase difference

I = I c sin(ϕ), (1.9)
where I c is the critical current. This is the dc Josephson effect. The effect depends on the fact that the macroscopic wavefunctions of the superconductors overlap in the middle region [START_REF] Feynman | The Schrödinger equation in a classical context: a seminar on superconductivity[END_REF].

If a dc voltage bias V is applied through the junction, then the phase difference evolves with time ϕ(t) = ϕ 0 + 2e h Vt. The result is that the supercurrent becomes alternating with a frequency ω J = 2eV h :

I(t) = I c sin(ϕ 0 + ω J t). (1.10)
This is the ac Josephson effect.

The dc and ac Josephson effect are not specific to the case of the tunnel junction originally considered by Josephson. Different hybrid superconducting systems made of superconductors coupled through some sort of weak link can exhibit the effects. We denote such Josephson junctions as S-X-S where S are the superconductors and X is the non-superconducting material. X can be an insulating barrier, in which case one talks of a "tunnel junction", but it can also be a point contact, a short constriction, a normal metal, or another type of material in the place of the insulating layer. Equation (1.8b) is a general quantum-mechanical relation that always holds, while the current-phase relation does not always have the sinusoidal form of Eq. (1.8a). For a review of different types of Josephson junctions and the corresponding current-phase relation, see [START_REF] Bakurskiy | Basic properties of the Josephson effect[END_REF].

In this thesis, we will focus on the case where X is a quantum dot (QD) [START_REF] Franceschi | Hybrid superconductor-quantum dot devices[END_REF][START_REF] Martín-Rodero | Josephson and Andreev transport through quantum dots[END_REF][START_REF] Meden | The Anderson-Josephson quantum dot-a theory perspective[END_REF]. Quantum dots are artificial atoms [START_REF] Kastner | Artificial atoms[END_REF][START_REF] Kouwenhoven | Few-electron quantum dots[END_REF] where electrons are trapped in a small region, resulting in a discrete spectrum. They are usually made of a metal, a semiconductor, a carbon nanotube or small molecules. Crucially, the number of trapped carriers can be controlled with a gate voltage, and this results in a variety of possible regimes for transport through the dot; for example, see [START_REF] Nazarov | Quantum transport: introduction to nanoscience[END_REF]. We will in particular consider the weakly-interacting regime where the dot can be modeled as a single discrete level [START_REF] Vecino | Josephson current through a correlated quantum level: Andreev states and π junction behavior[END_REF].

Andreev reflection

The size of a Cooper pair ξ 0 1 can be estimated from Heisenberg's uncertainty relation. A typical value is a few hundred nanometers, which is much larger than the average distance between electrons in matter ξ 0 ∼ 10 3 λ F , where λ F is the Fermi wavelength. This means that in a superconductor electrons are correlated at mesoscopic distances that are of the order of ξ 0 , which is therefore called the superconducting coherence length. Another direct result is the proximity effect: superconducting correlations are induced in non-superconducting materials that are brought into contact with a superconductor. The proximity effect is intimately related to the Andreev reflection mechanism [START_REF] Pannetier | Andreev reflection and proximity effect[END_REF]. As an example, consider a normal-superconducting interface (N-S). Naively, an electron coming from the normal metal with energy around the Fermi energy cannot tunnel through to the superconductor since there are no available states at energies inside the superconducting gap |E| < ∆. However, since the proximity effect creates non-zero superconducting correlations in the normal region, this means that there is a finite probability that the electron in the state k ↑ will pair up with another electron with opposite momentum and spin -k ↓, forming a Cooper pair. The Cooper pair can then tunnel to the S region. The result is that a charge of 2e is transferred from the N region to the S region [START_REF] Beenakker | Why does a metal-superconductor junction have a resistance?[END_REF]. An alternative way to view this process is to think of the ingoing electron as being reflected on the N-S interface as an outgoing hole of opposite spin. This reflection of an electron into a hole (and vice versa) is called the Andreev reflection and plays an important role in describing the microscopic physics of Josephson junctions.

Andreev bound states

Consider now the possibility of an S-N-S structure. If a right-moving electron in the N region is Andreev reflected as a left-moving hole on the right interface, the hole can then be reflected as a right-moving electron on the left interface. The electron will then be Andreev reflected into a hole on the right interface, and so on. Such multiple Andreev reflection processes between the two superconducting reservoirs will result in the formation of bound states, like a standing wave made of a superposition of electrons and holes. These are the Andreev bound states (ABS) [START_REF] Sauls | Andreev bound states and their signatures[END_REF]. The ABS have energies which come in pairs around zero inside the superconducting gap. In the limit where the length of the N region is short compared to the superconducting coherence length, the ABS energies are given by [START_REF] Martín-Rodero | Josephson and Andreev transport through quantum dots[END_REF] E A = ±∆ 1τ sin 2 (ϕ/2)

1 In BCS theory ξ 0 = hv F ∆ , where v F is the Fermi velocity and ∆ is the size of the superconducting gap. where τ is the transmission probability. As illustrated in Fig. 1.2, two consecutive Andreev reflections are equivalent to the transfer of a Cooper pair from the left superconductor to the right superconductor. This process corresponds to the negative energy bound state. The opposite process transfers Cooper pairs from the right superconductor to the left and corresponds to a positive energy. Therefore, each Andreev bound state results in a supercurrent of opposite direction I(ϕ) = 2e h ∂E A ∂ϕ and contributes to the dc Josephson effect. In recent years, enormous progress has been made in the fabrication and experimental manipulation of superconducting devices. Andreev bound states have been measured with the help of tunneling [START_REF] Pillet | Andreev bound states in supercurrent-carrying carbon nanotubes revealed[END_REF] and microwave spectroscopy [START_REF] Bretheau | Exciting Andreev pairs in a superconducting atomic contact[END_REF][START_REF] Bretheau | Theory of microwave spectroscopy of Andreev bound states with a Josephson junction[END_REF][START_REF] Van Woerkom | Microwave spectroscopy of spinful Andreev bound states in ballistic semiconductor Josephson junctions[END_REF][START_REF] Tosi | Spin-orbit splitting of Andreev states revealed by microwave spectroscopy[END_REF]. Progress towards the realization of an Andreev qubit [START_REF] Despósito | Controlled dephasing of Andreev states in superconducting quantum point contacts[END_REF][START_REF] Zazunov | Andreev level qubit[END_REF][START_REF] Chtchelkatchev | Andreev quantum dots for spin manipulation[END_REF] has also been achieved with the coherent manipulation of the Andreev states [START_REF] Janvier | Coherent manipulation of Andreev states in superconducting atomic contacts[END_REF].

Voltage-biased Josephson junctions

When a voltage difference V is applied across a junction, the current oscillates as a sine function of the Josephson frequency. Taking its time-average would then result in an average of zero Cooper pairs being transmitted through the junction. Therefore, one naively expects that the average current is zero if the voltage difference is smaller than the gap eV < ∆ (if eV exceeds the value of the gap, then quasiparticles can directly tunnel through the junction, producing a dissipative current). In reality, however, one observes a non-zero current even in the subgap region eV < ∆. This is possible because of higher-order Andreev reflection processes. The literature on this topic goes back to the classic paper by Klapwijk, Blonder, and Tinkham [START_REF] Klapwijk | Explanation of subharmonic energy gap structure in superconducting contacts[END_REF]. By traversing a normal region with a potential difference of V, quasiparticles can change their energy by an amount equal to eV. The quasiparticles can then undergo a number n of successive multiple Andreev reflections (MARs), thus changing their energy by neV. Then, whenever eV = 2∆/n, the quasiparticle has gained enough energy to overcome the superconducting gap. As an example, a third-order MAR process, occurring when eV > 2∆/3, is shown in Fig. 1.3a. MAR processes lead to a dc dissipative current, which exhibits jumps when the voltage is an integer subdivision of the gap. This is called the subharmonic gap structure or subgap structure of the current-voltage characteristics [START_REF] Octavio | Subharmonic energy-gap structure in superconducting constrictions[END_REF][START_REF] Averin | Ac Josephson effect in a single quantum channel[END_REF][START_REF] Bratus | Theory of subharmonic gap structure in superconducting mesoscopic tunnel contacts[END_REF][START_REF] Cuevas | Hamiltonian approach to the transport properties of superconducting quantum point contacts[END_REF]. It follows that at small voltage values one needs to take into account an increasing number of higher-order processes in order to calculate the current. This makes the inclusion of MAR processes a theoretical challenge.

In the limit of resonant tunneling through the junction, the subgap structure is greatly modified [START_REF] Martín-Rodero | Josephson and Andreev transport through quantum dots[END_REF][START_REF] Levy Yeyati | Resonant tunneling through a small quantum dot coupled to superconducting leads[END_REF][START_REF] Jonckheere | Nonequilibrium supercurrent through a quantum dot: Current harmonics and proximity effect due to a normal-metal lead[END_REF]: features corresponding to an odd number of MAR are enhanced, while even trajectories are suppressed. Moreover, the subgap current is sensitive to the energy position of the resonant level [START_REF] Johansson | Resonant multiple Andreev reflections in mesoscopic superconducting junctions[END_REF]. An example, calculated with the method presented in Chapter 3, is shown in Fig. 1.3b.

The Andreev molecule

Josephson junctions are nowadays routinely used in experiments as parts of superconducting circuits, and various setups have been developed in which Josephson junctions are used to realize superconducting qubits [START_REF] Kjaergaard | Superconducting qubits: current state of play[END_REF][START_REF] Devoret | Superconducting circuits for quantum information: an outlook[END_REF][START_REF] Kockum | Quantum bits with Josephson junctions[END_REF]. One of the proposed configurations is the Andreev level qubit, which uses the ABS of a Josephson junction as a two-level system [START_REF] Despósito | Controlled dephasing of Andreev states in superconducting quantum point contacts[END_REF][START_REF] Zazunov | Andreev level qubit[END_REF][START_REF] Chtchelkatchev | Andreev quantum dots for spin manipulation[END_REF][START_REF] Janvier | Coherent manipulation of Andreev states in superconducting atomic contacts[END_REF]. At the same time, as experimental control and techniques become more refined, there is increasing interest in the properties of more complex setups, such as transport properties of multiterminal Josephson junctions, where the normal region is coupled to more than two superconducting leads (Fig. 1.1b). In particular, a multiterminal Josephson junction can be used to engineer topological states, even when the junction is made from topologically trivial materials [START_REF] Riwar | Multi-terminal Josephson junctions as topological matter[END_REF]. Moreover, multiterminal configurations are being used to study nonlocal transport effects. For example, there is the possibility to create correlations between Cooper pairs, the so-called quartets [START_REF] Freyn | Production of nonlocal quartets and phase-sensitive entanglement in a superconducting beam splitter[END_REF][START_REF] Jonckheere | Multipair dc Josephson resonances in a biased all-superconducting bijunction[END_REF][START_REF] Pfeffer | Subgap structure in the conductance of a three-terminal Josephson junction[END_REF][START_REF] Cohen | Nonlocal supercurrent of quartets in a three-terminal Josephson junction[END_REF]. . This could be useful for realizing qubits whose coupling, for example, can be tuned by changing their phase difference, but one should have a distance of the junctions which remains comparable to the superconducting coherence length ξ 0 . In a typical superconductor, such as aluminum, ξ 0 ∼ 100 nm. The Andreev molecule and its signatures have been the recent subject of both theoretical [START_REF] Pillet | Nonlocal Josephson effect in Andreev molecules[END_REF][START_REF] Pillet | Scattering description of Andreev molecules[END_REF][START_REF] Kornich | Fine energy splitting of overlapping Andreev bound states in multiterminal superconducting nanostructures[END_REF][START_REF] Kornich | Overlapping Andreev states in semiconducting nanowires: competition of one-dimensional and three-dimensional propagation[END_REF][START_REF] Kocsis | Strong nonlocal tuning of the current-phase relation of a quantum dot based Andreev molecule[END_REF][START_REF] Scherübl | Transport signatures of an Andreev molecule in a quantum dot-superconductor-quantum dot setup[END_REF], as well as experimental studies [START_REF] Kürtössy | Andreev molecule in parallel InAs nanowires[END_REF][START_REF] Matsuo | Observation of nonlocal Josephson effect on double InAs nanowires[END_REF][START_REF] Coraiola | Phase-engineering the Andreev band structure of a three-terminal Josephson junction[END_REF][START_REF] Matsuo | Engineering of anomalous Josephson effect in coherently coupled Josephson junctions[END_REF][START_REF] Matsuo | Phase-dependent Andreev molecules and superconducting gap closing in coherently coupled Josephson junctions[END_REF]. The nonlocal Josephson effect was recently demonstrated [START_REF] Haxell | Demonstration of the nonlocal Josephson effect in Andreev molecules[END_REF]. Recent works have also focused on the occurrence of the Josephson diode effect-an asymmetry of the critical value of the supercurrent depending on its direction-in Andreev molecules [START_REF] Kocsis | Strong nonlocal tuning of the current-phase relation of a quantum dot based Andreev molecule[END_REF][START_REF] Pillet | Josephson diode effect in Andreev molecules[END_REF][START_REF] Hodt | On-off switch and sign change for a nonlocal Josephson diode in spin-valve Andreev molecules[END_REF][START_REF] Matsuo | Josephson diode effect derived from short-range coherent coupling[END_REF].

outline

In Chapter 2, a brief description of BCS theory is given, followed by a study of the S-QD-S system, as well as of the Andreev molecule in equilibrium (in the absence of voltage bias).

In Chapter 3, we present a formalism that can be used to study the Floquet spectrum of multiterminal Josephson junctions biased with commensurate voltages.

In particular, we first study a quantum dot coupled to two superconducting leads (S-QD-S system) before turning our attention to the Andreev molecule geometry (S-QD-S-QD-S system). The focus is on two observables in the steady state: the spectral function, which can be probed with tunneling spectroscopy, and the dc component of the current. The main results are the spectra of the driven Andreev molecule, as well as the corresponding subgap current. Most of the results of this chapter have been published in [START_REF] Keliri | Driven Andreev molecule[END_REF].

In Chapter 4, we study the long-range coupling that develops between the quantum dots of the S-QD-S-QD-S system when the middle superconductor is longer than the superconducting coherence length. We focus on the influence of this Floquet-Tomasch effect on the subgap region of the spectrum. The results of this chapter have been published in [START_REF] Keliri | Long-range coupling between superconducting quantum dots induced by periodic driving[END_REF]. Finally, in Chapter 5, we present preliminary work on the derivation of a Floquet master equation for the S-QD-S system. A general method is presented for the derivation of master equations for non-interacting fermionic systems. The method relies on a "quasiparticle approximation" and gives the populations of the Floquet states at the limit of weak coupling to the reservoirs, when the Floquet-Andreev resonances are sharply peaked.

introduction: historical

At the start of the 20th century, a debate was ongoing: what happens to the electrical resistance of a metal wire when its temperature is lowered to absolute zero? On one side there was Lord Kelvin's argument that electrons would freeze around the atoms, therefore causing the resistance to rise infinitely, while Drude's theory of the free electron gas suggested that the resistance would smoothly disappear, making the metal a perfect conductor at zero temperature [START_REF] Van Delft | The discovery of superconductivity[END_REF][START_REF] Van Delft | History and significance of the discovery of superconductivity by Kamerlingh Onnes in 1911[END_REF]. In 1911, the coldest place on earth was in Kamerlingh Onnes' lab in Leiden. Having managed to liquefy helium three years prior, 1 Onnes was attempting to give a definite answer to the resistance riddle by cooling down pure metals and measuring their electrical conductivity. However, what he found was yet another puzzle, one that would remain unsolved for almost 50 years. In Onnes' own words [START_REF] Onnes | Investigations into the properties of substances at low temperatures, which have led, amongst other things, to the preparation of liquid helium[END_REF]:

"As has been said, the experiment left no doubt that, as far as accuracy of measurement went, the resistance disappeared. At the same time, however, something unexpected occurred. The disappearance did not take place gradually but abruptly. [...] Thus the mercury at 4.2K has entered a new state, which, owing to its particular electrical properties, can be called the state of superconductivity."

A superconductor has many amazing properties (the most notable being zero electrical resistance and consequently possibility to transport current without dissipation-as discovered by Onnes-and perfect diamagnetism and expulsion of magnetic fields-as discovered by Meissner in 1933), but, most importantly, it is a many-body quantum phenomenon expressed macroscopically. An important step towards theoretical understanding was F. London's idea that the superconducting state possess a rigidity under perturbations. Bardeen later proposed that the rigidity of the wavefunction could be explained by the existence of an energy gap that separates the superconducting ground state from the lowest excitation energy.

The complete microscopic theory of superconductivity was finally proposed in 1957 by Bardeen, Cooper and Schrieffer [START_REF] Bardeen | Theory of superconductivity[END_REF]. The BCS theory explained the superconducting state as a condensate of Cooper pairs, pairs of electrons interacting through the exchange of phonons. Shortly after, in 1962, Josephson understood that a supercurrent ought to flow between two superconductors separated by a thin insulating layer (tunnel junction). If, moreover, the superconductors are at different voltages, the supercurrent is alternating, and the voltage to frequency conversion is so accurate, that it provides the standard for the definition of the volt. These two phenomena, now called the dc and ac Josephson effects, are routinely exploited by superconducting devices.

In this chapter I try to give a brief description of BCS theory, mostly following the development used in the excellent book by Schrieffer [START_REF] Schrieffer | Theory of superconductivity[END_REF]. The Bogoliubovde Gennes equation [START_REF] Gennes | Superconductivity of metals and alloys[END_REF] which is useful for inhomogeneous superconductors is presented, and care is taken in discussing the dynamics of Bogoliubov-de Gennes states. The chapter concludes with an application to the S-QD-S junction and the three-terminal S-QD-S-QD-S junction, otherwise known as the Andreev molecule. Throughout this thesis superconductivity always refers to conventional superconductivity, as described by BCS theory.

description of the condensed state

Cooper studied the problem of a pair of electrons which interact above a noninteracting Fermi sea [START_REF] Cooper | Bound electron pairs in a degenerate Fermi gas[END_REF]. He showed that if there is a net attractive interaction, no matter how weak, there exist bound states of electron pairs with opposite momenta and spins. Since the bound states have an energy below the Fermi energy, it is energetically favorable for electrons to pair up, resulting in the formation of a macroscopic number of these Cooper pairs, and an opening of a gap at the Fermi level. The Cooper pairs have some resemblance to bosons, but they are still subject to the Pauli exclusion principle. One can check that the pairs do not truly obey the Bose-Einstein statistics since

c -k↓ c k↑ , c † k↑ c † -k↓ = 1 -( n-k↓ + nk↑ ) c -k↓ c k↑ , c † k ′ ↑ c † -k ′ ↓ = 0, k ̸ = k ′ c -k↓ c k↑ , c -k ′ ↓ c k ′ ↑ = c † k↑ c † -k↓ , c † k ′ ↑ c † -k ′ ↓ = 0. (2.1)
The starting point for BCS theory is really Landau's Fermi liquid theory which describes the normal state of a metal. This is a reasonable starting point since the normal and the superconducting state have a small energy difference (around 10 -8 eV per electron). The wavefunction of the superconducting state is therefore constructed by considering pairing correlations between otherwise noninteracting quasiparticles [START_REF] Schrieffer | Theory of superconductivity[END_REF]. Throughout the thesis we sometimes abusively speak of electrons and holes, although in reality one should call them "quasi-electrons" and "quasi-holes". We will use the language of second quantization and denote the creation and annihilation operators for quasiparticles of momentum k and spin σ (↑ or ↓) as c † kσ and c kσ , respectively. These operators obey the usual fermionic anticommutation rules, {c † kσ ,

c k ′ σ ′ } = δ kk ′ δ σσ ′ , and {c † kσ , c † k ′ σ ′ } = {c kσ , c k ′ σ ′ } = 0.
|ψ 0 ⟩ = ∏ k (u k + v k c † -k↓ c † k↑ ) |0⟩ . (2.2)
The ground state is made of a superposition of occupied and unoccupied pair states, with v k being the probability amplitude to find the pair (k ↑, -k ↓) occupied, while u k is the probability amplitude to find the pair unoccupied. Since the ground state mixes states with different particle numbers, one can have nonzero expectation values of the type

∆ ∝ c -k↓ c k↑ ̸ = 0.
Since this quantity is zero in the normal state, it can serve as an order parameter for the superconducting phase transition. The order parameter is a complex quantity ∆ = |∆|e iϕ , and it is not gauge-invariant, meaning that it is not an observable. What is observable is the magnitude |∆|, which represents the size of the superconducting gap that opens due to the creation of Cooper pairs. Its existence was first verified in 1960 by I. Giaever [START_REF] Giaever | Energy gap in superconductors measured by electron tunneling[END_REF].

Excitations

Once the ground state is constructed, excited states can be produced by adding an electron in the state p ↑, and leaving its pair -p ↓ empty. Due to the Pauli principle, the pair state (p↑, -p↓) is blocked from participating in the pairing interaction, and therefore the energy of the system is increased. The singleparticle excitation spectrum in the superconducting phase is drastically different from the normal phase in the vicinity of the Fermi surface. This is due to the elements of bcs theory existence of a gap which needs to be overcome in order to create an excitation from the ground state. Adding an electron of momentum p and spin ↑

c † p↑ |ψ 0 ⟩ = c † p↑ ∏ k (u k + v k c † -k↓ c † k↑ ) |0⟩ = u p c † p↑ ∏ k̸ =p (u k + v k c † -k↓ c † k↑ ) |0⟩ ≡ u p ψ p↑ (2.3)
or removing an electron of momentum -p and spin ↓

c -p↓ |ψ 0 ⟩ = c -p↓ ∏ k (u k + v k c † -k↓ c † k↑ ) |0⟩ = v p c † p↑ ∏ k̸ =p (u k + v k c † -k↓ c † k↑ ) |0⟩ ≡ v p ψ p↑ .
(2.4)

creates the same state ψ p↑ (aside from a normalization factor). In a similar manner we find that removing an electron from the state p ↑ or adding an electron to the state -p ↓ produces the same state, since

c p↑ |ψ 0 ⟩ = -v p ψ -p↓ , c † -p↓ |ψ 0 ⟩ = u p ψ -p↓ . (2.5)
It is therefore reasonable to define composite quasiparticle creation operators γ † kσ which are a linear combination of bare creation and annihilation operators,

γ † kσ = u k c † kσ + σv k c -k-σ ,
where σ = +1 for spin ↑ and σ = -1 for spin ↓ . The fact that the amplitudes u k , v k coincide with the amplitudes in Eq. (2.2) is a consequence of imposing two conditions on the operators γ † and their adjoints γ :

1. the quasiparticle operators γ † and γ obey the usual anticommutation rules of Fermi-Dirac statistics 2. the ground state is destroyed by application of the annihilation operators

γ kσ |ψ 0 ⟩ = 0, ∀k, σ.
The last condition means that the ground state is the vacuum for the γ quasiparticles. From the above, it follows that the γ † kσ operator can be identified as the creation operator of a quasiparticle in the kσ state

γ † k↑ |ψ 0 ⟩ = ψ k↑ γ † -k↓ |ψ 0 ⟩ = ψ -k↓ , (2.6)
while its hermitian conjugate can be identified as the annihilation operator of that state.

Bogoliubov canonical transformation and diagonalization

The mean-field Hamiltonian for the BCS theory can be written as

H BCS = ∑ kσ ϵ k c † kσ c kσ + ∑ k ∆c † k↑ c † -k↓ + ∆ * c -k↓ c k↑ , (2.7)
where ϵ k is the kinetic energy measured from the Fermi level. Electron number is not conserved in H BCS due to the unusual terms c † k↑ c † -k↓ and c -k↓ c k↑ . Since the electron number is not well defined one expects that the eigenstates of H BCS will be a mixture of electrons and holes. One therefore introduces new operators

γ † kσ = u k c † kσ + σv k c -k-σ . (2.8)
This is the Bogoliubov transformation, or Bogoliubov-Valatin transformation [START_REF] Bogoljubov | A new method in the theory of superconductivity[END_REF][START_REF] Valatin | Comments on the theory of superconductivity[END_REF]. More than a mathematical trick to diagonalize the BCS Hamiltonian, the quasiparticle operators γ † kσ produce the excited states, as we saw previously. Demanding that the operators satisfy the Fermi-Dirac statistics {γ † kσ ,

γ k ′ σ ′ } = δ kk ′ δ σσ ′ , one gets the condition |u k | 2 + |v k | 2 = 1 for the BCS coherence factors. With this constraint, the inverse transformation is c † kσ = u * k γ † kσ -σv k γ -k-σ (2.9a) c kσ = u k γ kσ -σv * k γ † -k-σ (2.9b)
This can be compactly expressed using Nambu spinors

γ † k↑ γ -k↓ = u k v k -v * k u * k c † k↑ c -k↓ ⇔ c † k↑ c -k↓ = u * k -v k v * k u k γ † k↑ γ -k↓ .
(2.10)

By using the Bogoliubov transformation on H BCS we obtain

H BCS = ∑ k ϵ k |u k | 2 -ϵ k |v k | 2 + ∆u * k v k + ∆ * u k v * k γ † k↑ γ k↑ + -ϵ k |u k | 2 + ϵ k |v k | 2 -∆u * k v k -∆ * u k v * k γ -k↓ γ † -k↓ + -2ϵ k u * k v * k + ∆(u * k ) 2 -∆ * (v * k ) 2 γ † k↑ γ † -k↓ + -2ϵ k u k v k + ∆ * u 2 k -∆v 2 k γ -k↓ γ † k↑ .
(2.11)

The BCS coherence factors are complex numbers which we can write as

u k = x k e iϕ/2
(2.12a)

v k = y k e -iϕ/
-2ϵ k x k y k + |∆|(x 2 k -y 2 k ) = 0. (2.13)
Together with the constraint x 2 k + y 2 k = 1, we find that

x 2 k = 1 2   1 + ϵ k ϵ 2 k + |∆| 2   ≡ 1 2 1 + ϵ k E k (2.14a) y 2 k = 1 2   1 - ϵ k ϵ 2 k + |∆| 2   ≡ 1 2 1 - ϵ k E k . (2.14b)
We have introduced the quantity

E k = ϵ 2 k + |∆| 2 .
We will immediately see that it will acquire the meaning of a quasiparticle excitation energy. Indeed, the BCS Hamiltonian is now diagonalized

H BCS = ∑ k ϵ k (x 2 k -y 2 k ) + 2|∆|x k y k γ † k↑ γ k↑ -ϵ k (x 2 k -y 2 k ) + 2|∆|x k y k γ -k↓ γ † -k↓ = ∑ kσ E k γ † kσ γ kσ -∑ k (E k -ϵ k ), (2.15)
and the excitations are gapped since a minimum energy equal to |∆| is required for their creation.

H BCS |ψ 0 ⟩ = E g |ψ 0 ⟩ H BCS ψ pσ = (E g + E p ) ψ pσ .
(2.16)

the bogoliubov-de gennes equations

The Bogoliubov-de Gennes (BdG) equations [START_REF] Gennes | Superconductivity of metals and alloys[END_REF] are an alternative formulation of Eq. (2.15), leading to an eigenvalue equation for the coherence factors. This formalism is often used for describing inhomogeneous superconductors. Equation (2.15) can be rewritten as

H BCS , γ † kσ = E k γ † kσ H BCS , γ kσ = -E k γ kσ (2.17)
At the same time, using the definition of H BCS in Eq. (2.7), one can calculate the commutators 3

H BCS , c † k↑ = ϵ k c † k↑ + ∆ * c -k↓ H BCS , c -k↓ = -ϵ k c -k↓ + ∆c † k↑ (2.18)
Using the Bogoliubov transformation, the γ can be expressed in terms of creation and annihilation operators. Then, comparing coefficients on either side of Eq. (2.17), one arrives at an eigenvalue equation for the amplitudes u, v

ϵ k ∆ ∆ * -ϵ k u k v k = E k u k v k . (2.19)
The excitation energies are now the eigenvalues of the above equation. In the absence of an energy gap, the BdG equations would simply reduce to two decoupled Schrödinger equations: one for electrons, and a time-reversed one for holes. In the superconducting state, though, the two are mixed. It is easy to verify that Eq. (2.19) has a second eigenstate

v * k -u * k T with eigenvalue -E k .
This redundancy is also apparent in the fact that taking E → -E switches the role of the creation and annihilation operators γ † , γ, since creating a quasiparticle with energy E is identical to destroying a quasiparticle with energy -E. This redundancy was introduced by the doubling of the degrees of freedom which was used to bring the BCS Hamiltonian into a quadratic form.

Dynamics of Bogoliubov-de Gennes states

So far we have not involved any time dependence in our problem. We will need to be a bit careful with the time-evolution of the γ quasiparticles. As a reminder, if we are working in the Heisenberg picture, with U(t) = e -iHt being the time evolution operator of the system described by a Hamiltonian H(t), then operators transform as

A H (t) = U † (t)A(t)U(t). (2.20)
It follows that in the most general case their time evolution is given by

i d dt A H (t) = U † (t) A(t), H(t) + i d dt A(t) U(t).
(2.21)

3 Useful commutator relations: We now want to consider the time-evolution of states which evolve under a time-dependent Hamiltonian of BCS type

c † kσ c kσ , c † kσ = c † kσ c † kσ c kσ , c kσ = -c kσ c -k↓ c k↑ , c † k↑ = c -k↓ c † k↑ c † -k↓ , c -k↓ = c † k↑ c † k↑ c † -k↓ , c k↑ = -c † -k↓ c -k↓ c k↑ , c † -k↓ = -c k↑
H(t) = ∑ kσ ϵ k (t)c † kσ c kσ + ∑ k ∆(t)c † k↑ c † -k↓ + ∆ * (t)c -k↓ c k↑ . (2.23)
We make the ansatz that the time-evolved BCS ground state is always the vacuum for the quasiparticle operators

γ(t) |ψ(t)⟩ = 0, ∀t. (2.24)
In other words, the assumption is that |ψ(t)⟩ remains a BCS state at all times.

Taking the time-derivative, the immediate result is that

i d dt γ(t) |ψ(t)⟩ = 0 i d dt γ(t) |ψ(t)⟩ + γ(t)H(t) |ψ(t)⟩ = 0 i d dt γ(t) |ψ(t)⟩ + (γ(t)H(t) -H(t)γ(t)) |ψ(t)⟩ = 0. (2.25)
This means that the time evolution of the γ operators has a minus sign with respect to the usual Heisenberg picture since

i d dt γ(t) = H(t), γ(t) . (2.26)
This is formally equivalent to4 γ(t) = U(t)γ(0)U † (t).

(2.27)

In the case of a time-independent Hamiltonian and, due to Eq. (2.17), the evolution of the quasiparticle operators follows directly We can therefore write γ H (t) = γ(0) = γ. It will often be useful to switch between the γ † , γ quasiparticle operators and the 'bare' c † , c fermionic operators expressed in the Heisenberg picture. Starting from Eq. (2.9), the bare operators c † and c can be expressed in the Heisenberg picture as

γ kσ (t) = e iE k t γ kσ (0) (2.28a) γ † kσ (t) = e -iE k t γ † kσ (0). ( 2
c † σH (t) = u * (t)γ † σ -σv(t)γ -σ (2.30a) c -σH (t) = u(t)γ -σ + σv * (t)γ † σ . (2.30b)
The Heisenberg equation of motion for the c † operator is

i d dt c † σH (t) = U † c † σ , H(t) U(t) = -ϵ(t)c † σH (t) -∆ * (t)c -σ (t) i d dt u * (t)γ † σ -σv(t)γ -σ = -ϵ(t) u * (t)γ † σ -σv(t)γ -σ -∆ * (t) u(t)γ -σ + σv * (t)γ † σ . (2.31)
Matching coefficients on the two sides gives us the time-dependent version of the BdG equation:

i d dt u(t) = ϵ(t)u(t) + ∆(t)v(t) (2.32a) i d dt v(t) = ∆ * (t)u(t) -ϵ(t)v(t), (2.32b)
which is a natural generalisation of Eq. (2.19). In the time-independent case the evolution of the BdG states reduces to

u(t) v(t) = e -iEt u v (2.33)
and Eq. (2.30) gives us

c † σH (t) = e iEt u * γ † σ -σe -iEt vγ -σ (2.34a) c σH (t) = e -iEt uγ σ -σe iEt v * γ † -σ . (2.34b)
The utility of the above equations will be more evident in the next chapters, where we will deal with states that are time-dependent.

elements of bcs theory 2.4 applications 2.4.1 The S-QD-S system in equilibrium

We will now use the BdG formalism in order to derive some basic results about Josephson junctions. Throughout this thesis, we will model the Josephson junction as two superconducting reservoirs coupled to a single-level quantum dot. By solving the Bogoliubov-de Gennes equations, we can derive the subgap spectrum (the Andreev bound states). The superconducting reservoirs are labeled as j = L, R and are described by the BCS Hamiltonian

H r = ∑ jkσ ϵ k c † jkσ c jkσ + ∑ jk ∆ j c † jk↑ c † j-k↓ + ∆ * j c j-k↓ c jk↑ . (2.35)
Without any coupling to the dot, the eigenstates of the leads are the operators

γ † jkσ = u jk c † jkσ + σv jk c j-k-σ (2.36)
which are solutions of the BdG equations with eigenvalues E. The leads are then coupled to the single-level quantum dot through a tunneling term

V = ∑ jkσ J j d † σ c jkσ + c † jkσ d σ , (2.37)
where d † σ is the creation operator on the dot. We will use the Lippmann-Schwinger (LS) method in order to derive the spectrum on the dot. The LS equation is a rewriting of the Schrödinger equation, taking into account boundary conditions of the scattering problem. More specifically, consider an incoming beam of particles which are scattered off some potential V. Then, if in the absence of the potential, the state |ψ 0 ⟩ is an eigenstate of the unperturbed Hamil- tonian H 0 with energy E, a solution for the scattered wave |ψ⟩ is [105]:

|ψ⟩ = |ψ 0 ⟩ + 1 E -H 0 + iη V |ψ⟩ . (2.38)
In the same spirit, turning on the couplings J j in the tunneling term V will have as a result that the bare operators γ † will evolve into dressed operators Γ † . We therefore make the following ansatz for the dressed operators

Γ † lkσ = γ † lkσ + ud † σ + σvd -σ + ∑ jk ′ U(jk ′ ; lk)c † jk ′ σ + σV(jk ′ ; lk)c j-k ′ -σ . (2.39)
The dressed operator Γ † lkσ describes a quasiparticle being injected from the reservoir l with a momentum k and spin σ. The assumption is that it is still a BdG eigenstate of the total Hamiltonian H = H r + V, so that

H, Γ † lkσ = E Γ † lkσ . (2.40)
Calculating the commutators involved and matching the terms on the left-and on the right-hand side, one obtains a set of coupled equations for the amplitudes on the dot and on the reservoirs. Using

H, c † jkσ = ϵ k c † jkσ + ∆ * j c j-k-σ + J j d † σ H, c j-k-σ = -ϵ k c j-k-σ + ∆ j c † jkσ -J j d -σ V, d † σ = ∑ jk J j c † jkσ [V, d -σ ] = -∑ jk J j c j-k-σ (2.41)
one obtains:

E u σv = J l u lk -σv lk + ∑ jk ′ J j U(jk ′ ; lk) -σV(jk ′ ; lk) (2.42)
and

U(jk ′ ; lk) -σV(jk ′ ; lk) = J j E 2 -ϵ 2 k ′ -∆ j 2 E + ϵ k ′ -∆ j -∆ * j E -ϵ k ′ u σv . (2.43)
Then, the inhomogeneous LS equation for the amplitudes on the dot is:

E u v -∑ j J 2 j σ z g j (E)σ z u v = J l σ z u lk v lk (2.44)
where g j (ω) is the retarded Green's function of the superconducting reservoir j, obtained by integrating out the reservoirs, and σ z = 1 0 0 -1 is the third Pauli matrix. Using

∑ k ′ 1 ω 2 -ϵ 2 k ′ -∆ 2 + iη → +∞ -∞ dϵ ρ 0 ω 2 -ϵ 2 k ′ -∆ 2 + iη (2.45)
where ρ 0 is the density of states in the normal state of the superconductors one obtains

g j (ω) =            -πρ 0 √ ∆ 2 -ω 2 ω ∆ j ∆ * j ω |ω| < ∆ -iπρ 0 sign(ω) √ ω 2 -∆ 2 ω ∆ j ∆ * j ω |ω| > ∆ (2.46)
Defining Γ j = πρ 0 J 2 j , and g j ≡ σ z g j σ z we can write compactly,

J 2 j g j (ω) = Γ j iv F q(ω) ω -∆ j -∆ * j ω , (2.47)
where we have defined the quantity

v F q(ω) ≡ i ∆ 2 -ω 2 θ(∆ -|ω|) + sign(ω) ω 2 -∆ 2 θ(|ω| -∆). (2.48)
∆ is the magnitude of the gap which we take to be of equal value on both reservoirs ∆ = |∆ L | = |∆ R |. The amplitudes on the dot are then given by

E -∑ j Γ j iv F q(E) E -∆ j -∆ * j E u v = J l σ z u lk v lk . ( 2.49)
The right-hand side can be seen as a source term. The amplitudes on the dot can then be simply obtained by calculating the inverse of the operator on the left-hand side. This inverse operator can be identified as a Green's function of the dot or as the resolvent operator R(ω)

= (ω -H s -σ(ω) + iη) -1 .
Here the Hamiltonian of the dot is H s = ∑ σ ϵ d d † σ d σ and the self energy σ(ω) is obtained when integrating out the contribution from the reservoirs

σ(ω) ≡ ∑ j J 2 j σ z g j (ω)σ z . (2.50)
So far we have considered that the energy of the dot is ϵ d = 0, but the term can be easily restored if needed, in which case one finds

R(E) = E - ε d 0 0 -ε d -∑ j Γ j iv F q(E) E -∆ j -∆ * j E -1
.

(2.51)

We see that the resolvent has anomalous (off-diagonal) terms which mix the electron and hole amplitudes. This is a result of the proximity effect which induces superconducting correlations on the dot. Moreover, there are solutions for energies inside the superconducting gap |E| < ∆ : these are the Andreev bound states.

Andreev bound states

The poles of the resolvent give us the spectrum. For energies inside the gap |E| < ∆, they are the solutions of the equation

E 2 1 + ∑ j Γ j √ ∆ 2 -E 2 2 -ε 2 d - ∆ 2 ∆ 2 -E 2 Γ 2 L + Γ 2 R + 2Γ L Γ R cos(ϕ R -ϕ L ) = 0 (2.52) We can assume symmetric couplings Γ L = Γ R = Γ/2.
In fact, we can do this without any loss of generality, as has been shown in [START_REF] Kadlecová | Quantum dot attached to superconducting leads: relation between symmetric and asymmetric coupling[END_REF]. Denoting the phase difference ϕ Rϕ L = 2ϕ simplifies the above expression to

E 2 1 + Γ √ ∆ 2 -E 2 2 -ε 2 d - Γ 2 ∆ 2 cos 2 ϕ ∆ 2 -E 2 = 0. (2.53) 0 π 2 π 0 Δ -Δ ϕ R -ϕ L E Figure 2.1:
Andreev bound states as a function of the phase difference between the reservoirs, for Γ = ∆ and ϵ d = 0, ϵ d = ∆, ϵ d = 2 ∆ (black to light grey lines). Due to the symmetry of the BdG equation, states come in pairs mirrored around zero. Unlike other setups, the ABS of the S-QD-S system remain detached from the continua for all phases, even for large tunnel couplings. This will be important when we consider their dynamics.

The solutions can easily be found numerically. The ABS spectrum is plotted in Fig. 2.1 as a function of the phase difference ϕ, for different positions of the dot level. For ε d = 0, the above equation is simplified to

E ± ∆ cos ϕ + E √ ∆ 2 -E 2 Γ = 0. (2.54)
Therefore, the solutions tend to the energies of a perfectly transmitting junction E ≃ ±∆ cos ϕ in the limit that Γ ≫ ∆. As Γ is reduced, the ABS become detached from the superconducting continua for all phase differences. In the opposite limit of ∆ ≫ Γ, the energies are of the order of the tunnel couplings E ≃ ±Γ cos ϕ. The effect of a finite energy ε d is to open a gap between the positive and negative ABS, since it breaks the symmetric role of the electron-like and hole-like components in Eq. (2.51). Since for small tunnel couplings the solutions are detached from the gap with E ≪ ∆, we can approximate the square root by

√ ∆ 2 -E 2 ≈ ∆ -E 2 2 ∆
. Renormalizing all quantities of Eq. (2.53) by ∆ one finds

1 + Γ ∆ E 4 -ϵ 2 d + (Γ + ∆) 2 E 2 + ∆ 2 ϵ 2 d + Γ 2 ∆ 2 cos 2 (ϕ) = 0. (2.55)
Two of the four solutions obey the condition which gives us an analytic solution that works in the limit where ∆ is the largest energy in the system.

E 2 < ∆ 2 , E ABS ∆ 2 = ϵ 2 d + (Γ + ∆) 2 2∆(Γ + ∆) - ϵ 2 d + (Γ + ∆) 2 2∆(Γ + ∆) 2 - ϵ 2 d + Γ 2 cos 2 (ϕ) ∆(Γ + ∆) , (2.56)

The Andreev molecule

We now consider the three terminal S-QD-S-QD-S junction, which we could also call a bijunction, consisting of two quantum dots connected to three superconducting reservoirs, as shown in Fig. 2.2. We assume that the middle superconductor is large on the mesoscopic scale, so that it can be considered a reservoir, but that the distance between the dots remains finite. We set the phase in the middle reservoir to zero ϕ c = 0, which means that the phase difference through the first (second) junction is controlled by the phase of the left (right) superconductor ϕ a(b) . In the limit where the middle reservoir is short compared to the superconducting coherence length R ≲ ξ 0 the system is a realization of the Andreev molecule. To the best of my knowledge, the spectrum of this setup has not been presented anywhere in the literature before. Previous works have focused on the current-phase relation [START_REF] Kocsis | Strong nonlocal tuning of the current-phase relation of a quantum dot based Andreev molecule[END_REF][START_REF] Scherübl | Transport signatures of an Andreev molecule in a quantum dot-superconductor-quantum dot setup[END_REF], and on the effect of the Coulomb energy on the dots (which we do not consider here). Nevertheless, the results are similar with the spectra of Andreev molecules based on weak links [START_REF] Pillet | Nonlocal Josephson effect in Andreev molecules[END_REF][START_REF] Scherübl | Transport signatures of an Andreev molecule in a quantum dot-superconductor-quantum dot setup[END_REF][START_REF] Pillet | Josephson diode effect in Andreev molecules[END_REF], with the key difference being that, as we shall see, for the S-QD-S-QD-S system there are no leaky Andreev states.

For the bijunction we can repeat the procedure used for the single junction case presented in the previous section. By writing the Lippmann-Schwinger equations and integrating out the reservoirs, we derive an eigenvalue equation for the amplitudes on the two dots. We consider a setup as depicted in Fig. 2.2 with three superconducting reservoirs j = {a, b, c} connected to two dots i = {1, 2}. The total Hamiltonian in the case of the bijunction is

H = H r + V H r = ∑ jkσ ϵ k c † jkσ c jkσ + ∑ jk ∆ j c † jk↑ c † j-k↓ + ∆ * j c j-k↓ c jk↑ V = ∑ i∈dots ∑ jkσ J ji d † iσ (x i )c jkσ + J ji c † jkσ d iσ (x i ) (2.57)
with real tunnel couplings J ji = J ij . The creation operator on dot i can be written as d † iσ (x i ) = e ikx i d † iσ , where x i is the position of the dot. We note the distance between the two dots as x 2x 1 = R. For convenience we choose the positions of the dots to be x 1 = 0 and x 2 = R. As in the S-QD-S case, we define a dressed operator which must now include operators for the two dots,

Γ † lkσ = γ † lkσ + ∑ i∈dots u(i)d † iσ + σv(i)d i-σ + ∑ jk ′ U(jk ′ ; lk)c † jk ′ σ + σV(jk ′ ; lk)c j-k ′ -σ .
(2.58) Calculating the commutator

H, Γ † lkσ = E γ † lkσ (2.59)
is slightly more tedious, but gives

E u(i) v(i) = J li e ikx i u lk -v lk + ∑ jk J ji e ikx i U jk -V jk (2.60) and U jk -V jk = ∑ i ′ J ji ′ e -ikx i ′ σ z g j (E, k)σ z u(i ′ ) v(i ′ ) (2.61)
where

g j (ω, k) = 1 ω 2 -ϵ 2 k -∆ 2 +iη ω + ϵ k ∆ j ∆ * j ω -ϵ k
is the retarded Green's function of the j-th superconducting reservoir before integration over all the mo- menta k. We end up with a set of coupled inhomogeneous equations for the amplitudes on the two dots

E u(1) v(1) = ∑ jk J 2 j1 σ z g j (E, k)σ z u(1) v(1) + J j1 J j2 e -ikR σ z g j (E, k)σ z u(2) v (2)
+ J l1 σ z u lk v lk E u(2) v(2) = ∑ jk J j2 J j1 e ikR σ z g j (E, k)σ z u(1) v(1) + J 2 j2 σ z g j (E, k)σ z u(2) v (2)
+ J l2 e ikR σ z u lk v lk

(2.62)

Summing over the momenta gives local Green's functions, as well as non-local Green's functions that depend on the distance R between the dots. Written in the basis of (u(1), v(1), u(2), v(2)), the resolvent of the system is5

R(E) = E 1 4 -     σ 1 (E) Γ c g c (E, R) Γ c g c (E, R) σ 2 (E)     -1 = E 1 4 -      ∑ j=a,c Γ j1 g j (E) Γ c g c (E, R) Γ c g c (E, R) ∑ j=b,c Γ j2 g j (E)      -1
.

(2.63)

where

g j (ω) = 1 iv F q(ω) ω -∆ j -∆ * j ω ,
and

g j (ω, R) = e iq(ω)R cos(k F R)g j (ω) + sin(k F R)σ z (2.64)
The spectrum of the Andreev molecule is therefore given by the solutions of the equation

det R(E) = 0 ⇒ det(E 1 2 -σ 1 ) × det E 1 2 -σ 2 -Γ 2 c g c (E, R) 1 E -σ 1 g c (E, R) = 0 (2.65)
where each determinant on the second line involves 2 × 2 matrices, each one associated with one of the two dots. In practice, it is much easier to numerically calculate the spectral functions A 1 = -2 π Im R 11 (E) and A 2 = -2 π Im R 33 (E) of the first and second dot, respectively. A justification for these formulas will be given in the next chapter, see Eq. (3.47). The maxima of the spectral functions correspond to the solutions of Eq. (2.65).

We are now in a position to comment on the coupling between the dots: the term responsible for the coupling is the non-diagonal block of the resolvent. This term is the non-local Green's function of the reservoir S c and represents propagation of the quasiparticles in 1d. For propagation inside the gap |E| < ∆ its magnitude is controlled by a factor

e iq(E)R = e - √ ∆ 2 -E 2 R/v F = e - √ 1-( E ∆ ) 2 R/ξ 0 . (2.66)
The coupling therefore decreases exponentially as the interdot distance R becomes larger than the superconducting coherence length ξ 0 . At R ≫ ξ 0 one retrieves two decoupled junctions, with ABS given by the solutions of

det(E 1 2 -σ 1 ) × det(E 1 2 -σ 2 ) = 0.
The effect on the spectrum when decreasing the distance is shown in Fig. 2.3.

When the distance is large, R ≫ ξ 0 , the two junctions are decoupled and the spectrum of the second junction at x = R (red lines in Fig. 2.3a) is unaffected by a change of the phase difference at x = 0. In the molecular regime R < ξ 0 , the ABS on the two junctions are hybridized. This produces avoided crossings in the spectrum, as shown in Fig. 2.3b. In contrast with Andreev molecules based on short weak links [START_REF] Pillet | Nonlocal Josephson effect in Andreev molecules[END_REF], all states in the S-QD-S-QD-S case are within the gap. The phase difference through the second junction, ϕ b can be used to tune the hybridization at a certain interdot distance. As shown in Fig. 2.4, at ϕ b = 0 the junctions are almost decoupled, while increasing it at ϕ b = π/2 or ϕ b = π produces avoided crossings. For identical junctions (with equal tunnel couplings), the avoided crossings are always expected to be found at ϕ a = ±ϕ b [START_REF] Pillet | Nonlocal Josephson effect in Andreev molecules[END_REF].

Part II

F L O Q U E T -A N D R E E V R E S O N A N C E S F L O Q U E T -A N D R E E V R E S O N A N C E S I N M U LT I T E R M I N A L Q U A N T U M D O T S

introduction

The aim of this chapter is to present a formalism that can be used to study the Floquet spectrum of multiterminal Josephson junctions biased with commensurate voltages. In particular, we will first study a quantum dot coupled to two superconducting leads (S-QD-S system) before turning our attention to the Andreev molecule geometry (S-QD-S-QD-S system). voltage-biased s-qd-s junction. As we have seen in the previous chapter, the spectrum of a quantum dot coupled to two superconducting reservoirs contains a pair of discrete subgap states, the Andreev bound states (ABS). When a voltage difference V is applied to such a junction, the discrete states evolve into resonances with a finite width [START_REF] Mélin | Simple Floquet-Wannier-Stark-Andreev viewpoint and emergence of low-energy scales in a voltage-biased three-terminal Josephson junction[END_REF][START_REF] Mélin | Engineering the Floquet spectrum of superconducting multiterminal quantum dots[END_REF][START_REF] Douçot | Berry phase in superconducting multiterminal quantum dots[END_REF]. The finite lifetime of the associated quasiparticles is a result of the multiple Andreev reflection (MAR) processes that allow quasiparticles to eventually gain the energy needed to overcome the superconducting gap 2 ∆. Since MAR are higher-order processes, the width of the resonances is exponentially small in ∆/eV [START_REF] Mélin | Engineering the Floquet spectrum of superconducting multiterminal quantum dots[END_REF]. Note that it is important that we consider the specific case of quantum dots, because their equilibrium ABS are detached from the superconducting continua (for example, see Fig. 2.1 in the previous chapter). This means that the states on the dot are less hybridized with the continua of the reservoirs, and therefore the non-equilibrium states on the dot have a longer lifetime.

josephson junctions as floquet systems. A Josephson junction is often simply described as a perfect voltage-to-frequency converter. Due to the Josephson relation dϕ dt = 2e h V(t), an applied dc voltage means that the superconducting phase difference between the reservoirs evolves linearly with time. The Hamiltonian that describes the S-QD-S junction is then time-periodic H(t) = H(t + 2π ω 0 ), where ω 0 is the frequency associated with the drive. Such periodically driven systems can be studied using Floquet theory [1]. The immediate result of the existence of a discrete time symmetry is that the energy is only conserved up to a multiple of the drive frequency [START_REF] Shirley | Solution of the Schrödinger equation with a Hamiltonian periodic in time[END_REF][START_REF] Sambe | Steady states and quasienergies of a quantum-mechanical system in an oscillating field[END_REF]. One therefore talks of "quasi-energies" in direct analogy to the "quasi-momenta" of Bloch theory for systems with periodicity in space. For the spectrum of the S-QD-S system, this means that instead of just a pair of resonances, one expects to find sidebands which are separated by integer multiples of the drive.

the wannier-stark problem. In fact, it is instructive to draw the analogy to band theory further. Using a semiclassical picture, Zener [START_REF] Zener | A theory of the electrical breakdown of solid dielectrics[END_REF] showed that a particle in a periodic potential will oscillate under application of a static electric field, with a frequency of oscillations

ω B = dF/h
that is proportional to the period of the potential d and the electric force F. The existence of oscillations, known as Bloch oscillations, was verified many years later [START_REF] Waschke | Coherent submillimeter-wave emission from Bloch oscillations in a semiconductor superlattice[END_REF][START_REF] Ben Dahan | Bloch oscillations of atoms in an optical potential[END_REF]. The difficulty arose because their period is typically much larger than the scattering time of electrons on defects, making it difficult to observe them 1 . With a full quantum-mechanical treatment, Wannier showed [START_REF] Wannier | Wave functions and effective Hamiltonian for Bloch electrons in an electric field[END_REF] that the Bloch oscillations in the time-domain give rise to the Wannier-Stark ladder in the energy domain. Namely, if the Hamiltonian of the Bloch particle in a static external field has an eigenstate E α , then a whole ladder of eigenstates, separated by the Bloch frequency, exist in the spectrum

E α,m = E α + mhω B , m ∈ Z.
In reality, one should consider that the states are coupled to an energy continuum through Zener tunneling [START_REF] Zener | A theory of the electrical breakdown of solid dielectrics[END_REF]. The solutions are then complex [START_REF] Avron | Instability of the continuous spectrum: The N-band Stark ladder[END_REF] E α,m = E α + mhω B -iΓ α , giving rise to a Wannier-Stark ladder of resonances. See reviews [START_REF] Mendez | Wannier-Stark ladders and Bloch oscillations in superlattices[END_REF][START_REF] Glück | Wannier-Stark resonances in optical and semiconductor superlattices[END_REF][START_REF] Kolovsky | Bloch oscillations of cold atoms in optical lattices[END_REF] and references within. Coming back to the problem of a voltage-biased S-QD-S junction, the dc voltage gives rise to two sets of ladders, the Floquet-Wannier-Stark ladders (FWS) [START_REF] Mélin | Engineering the Floquet spectrum of superconducting multiterminal quantum dots[END_REF] which can be seen as a dynamical version of the Wannier-Stark ladder. The fact that they come in pairs of two is a consequence of the doubling that exists in superconducting systems where eigenvalues come in pairs that are symmetric around zero. We argue that the FWS ladders could be observed by probing the spectral function of the dot. This can be done by performing local tunneling spectroscopy on the dot.

multi-terminal junctions. The existence of FWS ladders of resonances holds more generally for multiterminal configurations [START_REF] Mélin | Simple Floquet-Wannier-Stark-Andreev viewpoint and emergence of low-energy scales in a voltage-biased three-terminal Josephson junction[END_REF]. However, in this case application of commensurate voltages is required in order to obtain a single basic frequency throughout the system. The simplest non-trivial case involves a three-terminal junction biased in the so-called quartet configuration, where two leads are biased at opposite voltages V a = -V b , and the third one is, without loss of generality, grounded, V c = 0. A particular realization of a three-terminal junction is the Andreev molecule. We will refer to the Andreev molecule biased in the quartet configuration as the driven Andreev molecule. Note that, in principle, the theory could be extended to the application of incommensurate voltages, where each additional incommensurate frequency will correspond to an extra dimension in the phase space of the problem. This has been used, for example, to study Anderson localization in three (synthetic) dimensions [START_REF] Casati | Anderson transition in a one-dimensional system with three incommensurate frequencies[END_REF].

In this chapter we will present a method for treating multi-terminal superconducting junctions with quantum dots biased with commensurate voltages. The method builds up on previous work [START_REF] Mélin | Simple Floquet-Wannier-Stark-Andreev viewpoint and emergence of low-energy scales in a voltage-biased three-terminal Josephson junction[END_REF][START_REF] Mélin | Engineering the Floquet spectrum of superconducting multiterminal quantum dots[END_REF]. We will show that the problem can be mapped to a tight-binding chain in Floquet space, where the time dimension is traded for a fictitious Floquet dimension. We take advantage of this tight-binding structure to derive the Green's function on the dot(s) with a continued fraction method. We find that the equilibrium discrete states on each dot are dressed by MAR processes and evolve into ladders of resonances. We will focus on two observables, namely the spectral function and the current. We show how to calculate these observables in the steady state. We will argue that the time-averaged spectral function can be probed with tunneling spectroscopy, which therefore can reveal the Floquet spectrum and the dynamical Wannier-Stark ladders. The method is first applied to the single S-QD-S junction and then to the driven Andreev molecule. The spectrum of the driven molecule exhibits level splitting when the separation between the dots is comparable to the superconducting coherence length, R ≲ ξ 0 . These Floquet-Andreev resonances leave their trace in the dc current. We therefore calculate the steady state current passing through one junction and see that the proximity of the second junction modifies the usual MAR steps, which accordingly exhibit splitting into substeps. Most of the results of this chapter have been published in [START_REF] Keliri | Driven Andreev molecule[END_REF].

method overview

We are interested in Josephson junctions made of two or more superconducting reservoirs. The reservoirs are modeled as BCS superconductors, described by the Hamiltonian:

H r = ∑ jkσ ϵ k c † jkσ c jkσ + ∑ jk ∆ j c † jk↑ c † j-k↓ + ∆ * j c j-k↓ c jk↑ . (3.1)
The c † jkσ and c jkσ operators create and annihilate an electron in the j reservoir with momentum k and spin σ, correspondingly. For simplicity, we consider that all superconductors have a gap of equal magnitude, ∆ j ≡ ∆e iϕ j .

Depending on the particular configuration, the reservoirs are then coupled to one or more small regions. We will specify the particular geometries in later sections. The small regions form quantum dots with a Hamiltonian

H s = ε i d † iσ (x i )d iσ (x i )
, where the operator d † iσ creates an electron with spin σ on the dot labeled i which is at some position x i . A gauge transformation can be used to get rid of the space dependence in H s so that d † (x i ) = e ikx i d † . For simplicity, we will model a quantum dot as a single discrete level at zero energy, i.e. ε i = 0.

When a dc voltage V j is applied to the superconductor j, its Hamiltonian acquires a time-dependence according to the Josephson relation ϕ j (t) = ϕ j + 2eV j t, where ϕ j is the superconducting phase of the j superconductor. The timedependence can be eliminated from H r with the gauge transformation c † j → c † j e -ieV j t c j → c j e ieV j t .

(3.

2)

The total Hamiltonian can therefore be written as a time-independent part

H 0 = H r + H s , plus a time-dependent part V (t), H(t) = H 0 + V (t). (3.3)
The time-dependent part V (t) describes the tunneling between a dot at position x i and the j reservoir

V (t) = ∑ i∈dots ∑ jkσ J ji (x i )e is j ω 0 t d † iσ c jkσ + J * ji (x i )e -is j ω 0 t c † jkσ d iσ (3.4)
The gauge transformation (3.2) has therefore transferred the time-dependence to the hopping amplitudes. The tunnel couplings are J ji (x i ) = J ji e ikx i , with real amplitudes J ji = J ij = J * ji . We will moreover use the notation V j = s j V, with s j ∈ {0, ±1} depending on the particular configuration. We choose commensurate voltages because this choice leads to a single basic frequency ω 0 = eV/h, and to a time-periodic Hamiltonian.

Using the basic idea of the Floquet method [START_REF] Oka | Floquet engineering of quantum materials[END_REF][START_REF] Shirley | Solution of the Schrödinger equation with a Hamiltonian periodic in time[END_REF][START_REF] Sambe | Steady states and quasienergies of a quantum-mechanical system in an oscillating field[END_REF][START_REF] Hänggi | Driven quantum systems[END_REF], quantities can be expanded into Fourier modes,

H(t) = ∑ m∈Z e -imω 0 t H m
where the harmonic index m can be thought of as a position on a fictional Floquet direction. Physically, it corresponds to the number of quanta absorbed or emitted by the system. One then obtains a time-independent tight-binding model in an extended Hilbert space [START_REF] Shirley | Solution of the Schrödinger equation with a Hamiltonian periodic in time[END_REF][START_REF] Grempel | Quantum dynamics of a nonintegrable system[END_REF]. A common procedure is to "project out" the contribution of sites n ̸ = m up to some large Floquet index n = N, and arrive to an effective Floquet Hamiltonian for the site m [START_REF] Oka | Floquet engineering of quantum materials[END_REF]. The dimensions of the obtained tight-binding model depend on the number of incommensurate drive frequencies [START_REF] Martin | Topological frequency conversion in strongly driven quantum systems[END_REF]. Here, we will have one basic frequency across the system, so we will obtain an effective 1D tight-binding model.

Dressed quasiparticle basis

In the absence of tunneling V (t) = 0, an eigenstate of the superconducting reservoir labeled l is

γ † lkσ = e -iE lk t u lk c † lkσ + σv lk c l-k-σ , (3.5)
which evolves in time according to Eq. (2.28b), and where u lk = x lk e iϕ l /2 , v lk = y lk e -iϕ l /2 and the coefficients x, y are the usual coefficients obtained by diagonalizing the BCS Hamiltonian

x 2 lk = E lk + ϵ k 2E lk and y 2 lk = E lk -ϵ k 2E lk . ( 3.6)
The energy E lk = ϵ 2 k + ∆ 2 is the excitation energy needed for adding an electron or a hole to the BCS ground state of the l superconductor. The Bogoliubovde Gennes (BdG) equation then has the solutions,

i d dt γ † lkσ = H r , γ † lkσ = E lk γ † lkσ . (3.7)
When the tunneling is turned on, the discrete ABS states on the dots become resonances due to the multiple Andreev reflection processes (MAR) that connect them to the superconducting continua. The main idea is that in contrast to a closed Floquet system, an open Floquet system does not, in general, thermalize, but may reach a non-equilibrium steady state [START_REF] Mori | Floquet states in open quantum systems[END_REF][START_REF] Tsuji | Floquet states[END_REF]. Moreover, subgap ABS states of quantum dots are detached from the superconducting continua; this protects the resonances from dissipation, resulting in long lifetimes. It is reasonable, then, to expect well-defined quasiparticles in this type of system. Therefore, we start by constructing dressed quasi-particle operators Γ † lkσ describing a quasi-particle being injected from a source reservoir l with momentum k and spin σ, to any of the reservoirs j, and quantum dot(s) i

Γ † lkσ =γ † lkσ + ∑ i∈dots u(i; lk)d † iσ + σv(i; lk)d i-σ + ∑ jk ′ U(jk ′ ; lk)c † jk ′ σ + σV(jk ′ ; lk)c j-k ′ -σ . (3.8)
In the absence of tunnel couplings, the Γ † operators reduce to the bare quasiparticle operators γ † . The amplitudes u(i; lk), v(i; lk) have respectively the meaning of an electron-or hole-like amplitude on the dot i, while capital letters U, V correspond to amplitudes in the reservoirs. It is important to keep track of the reservoir which serves as a source, and therefore we are obliged to add it as an index to all amplitudes. The amplitudes on the dots and on the reservoirs are time-periodic and can therefore be Fourier-expanded. As a result, the creation operator Γ † can be written as a Fourier series

Γ † lkσ (t) =γ † lkσ (t) + ∑ m∈Z ∑ i∈dots e -i(E lk +mω 0 )t u m (i; lk)d † iσ + σv m (i; lk)d i-σ + ∑ m∈Z ∑ jk ′ e -i(E lk +mω 0 )t U m (jk ′ ; lk)c † jk ′ σ + σV m (jk ′ ; lk)c j-k ′ -σ .
(3.9)

The dressed operators are Floquet solutions of the BdG equations:

i d dt Γ † lkσ (t) = H(t), Γ † lkσ , (3.10)
and obey the Floquet theorem,

Γ † lkσ (t + T) = e -iE lk t Γ † lkσ (t). (3.11)
We make the assumption that the dressed operators form an orthonormal and complete basis. A discussion on the completeness of scattering states can be found in [START_REF] Goldberger | Collision theory[END_REF]. Let us add that the completeness relation expresses the assumption that there are no bound states in the system, i.e. it expresses the physical intuition that all ABSs acquire a finite lifetime when the voltage is finite. The completeness relation allows to express all other operators (on the dot and on the reservoirs) as linear combinations of the dressed operators. In Heisenberg representation (see previous chapter, Eq. (2.30)), the "bare" operators can be written as:

d † iσ (t) = ∑ lk ∑ m e i(E lk +mω 0 )t u * m (i; lk)Γ † lkσ -σe -i(E lk +mω 0 )t v m (i; lk)Γ lk-σ d iσ (t) = ∑ lk ∑ m e -i(E lk +mω 0 )t u m (i; lk)Γ lkσ -σe i(E lk +mω 0 )t v * m (i; lk)Γ † lk-σ .
(3.12)

Similarly, for the reservoir operators

c † jk ′ σ (t) = ∑ lk ∑ m e i(E lk +mω 0 )t δ m0 δ jl δ kk ′ u * jk ′ + U * m (jk ′ ; lk) Γ † lkσ -σe -i(E lk +mω 0 )t δ m0 δ jl δ kk ′ v jk ′ + V m (jk ′ ; lk) Γ lk-σ c jk ′ σ (t) = ∑ lk ∑ m e -i(E lk +mω 0 )t δ m0 δ jl δ kk ′ u jk ′ + U m (jk ′ ; lk) Γ lkσ -σe i(E lk +mω 0 )t δ m0 δ jl δ kk ′ v * jk ′ + V * m (jk ′ ; lk) Γ † lk-σ (3.13)
Finally, it is natural to assume that the steady state is the state annihilated by any destruction operator Γ lkσ , that is

Γ lkσ |S⟩ = 0, (3.14)
for any set of {lkσ}.

Floquet-Lippmann-Schwinger equations

We insert Eq. (3.9) into the BdG equation (3.10). Calculation of the various commutators 2 leads to two sets of coupled equations for the amplitudes on the dots and the amplitudes on the reservoirs. Explicitly, we find:

(E lk + mω 0 + iη) u m (i; lk) v m (i; lk) = J li e ikx i δ m,-s l u lk -δ m,s l v lk + ∑ jk ′ J ji e ik ′ x i U m+s j (jk ′ ; lk) -V m-s j (jk ′ ; lk) (3.15)
and

U m (jk ′ ; lk) -V m (jk ′ ; lk) = ∑ i J ji e -ik ′ x i σ z g j (E lk + mω 0 , k ′ )σ z u m-s j (i; lk) v m+s j (i; lk) . ( 3.16)
We can therefore integrate out the reservoirs by performing the summation over the momenta

∑ k ′ J ji e ik ′ x i σ z U m (jk ′ ; lk) V m (jk ′ ; lk) = = ∑ i ′ ∑ k ′ J ji J ji ′ e ik ′ (x i -x i ′ ) σ z g j (E lk + mω 0 , k ′ )σ z u m-s j (i ′ ; lk) v m+s j (i ′ ; lk) = ∑ i ′ Γ j,ii ′ σ z g j,ii ′ (E lk + mω 0 )σ z u m-s j (i ′ ; lk) v m+s j (i ′ ; lk) (3.17)
2 Useful commutator relations:

V (t), c † jkσ = ∑ i J ji e ikx i e is j ω 0 t d † iσ V (t), d † iσ = ∑ jk J ji e -ikx i e -is j ω 0 t c † jkσ V (t), c j-k-σ = -∑ i J ji e ikx i e -is j ω 0 t d i-σ [V (t), d i-σ ] = -∑ jk J ji e -ikx i e is j ω 0 t c j-k-σ Moreover, H r , c † jkσ = ϵ k c † jkσ + ∆ * j c j-k-σ , and H r , c j-k-σ = -ϵ k c j-k-σ + ∆ j c † jkσ .
where we have defined Γ j,ii ′ = πρ 0 J ji J ji ′ , with ρ 0 being the density of states in the normal state of the superconductors. After integrating out the reservoirs, one is therefore led to a set of inhomogeneous Floquet-Lippmann-Schwinger (FLS) equations for the amplitudes on the dots:

(E lk +mω 0 + iη)u m (i; lk) = δ m,-s l J li (x i )u lk + ∑ ji ′ Γ j,ii ′ g 11 j,ii ′ (m + s j )u m (i ′ ; lk) + g 12 j,ii ′ (m + s j )v m+2s j (i ′ ; lk) (E lk +mω 0 + iη)v m (i; lk) = -δ m,s l J li (x i )v lk + ∑ ji ′ Γ j,ii ′ g 21 j,ii ′ (m -s j )u m-2s j (i ′ ; lk) + g 22 j,ii ′ (m -s j )v m (i ′ ; lk) (3.18)
In the above equation, g j,ii ′ δ ii ′ ≡ g j is the dimensionless Green's function for the j superconducting reservoir, defined as:

g j (ω) = 1 iv F q(ω) ω -∆ j -∆ * j ω , v F q(ω) ≡ i ∆ 2 -ω 2 θ(∆ -|ω|) + sign(ω) ω 2 -∆ 2 θ(|ω| -∆). (3.19)
Note the unusual notation g j = σ z g j σ z . Moreover, since the quasi-energy appears in the combination ω + mω 0 , it is convenient to use the shorthand f (m) instead of f (ω + mω 0 ), for any function f . In the case of more than one dot, there are also non-local Green functions that depend on the distance separating the dots. We will explicitly give the expression in the case of the driven Andreev molecule, see Section 3.5.

We can rewrite Eq. (3.18) in the basis of the Nambu spinor

Ψ m ≡ (u m (1), v m (1), u m (2), v m (2), ...) T
which collects the amplitudes on the dot(s), by defining a linear operator L that acts on the states Ψ m :

(LΨ) m ≡ M 0 m Ψ m -M + m+1 Ψ m+2 -M - m-1 Ψ m-2 = S m . (3.20)
Equation (3.20) defines a Floquet operator L. Written in a matrix representation, it is a tridiagonal block-matrix of dimension dot ⊗ Nambu ⊗ Floquet. In the tight-binding analogy, the matrix M 0 m describes an on-site energy at position m of the chain, while matrices M ± m±1 describe hopping to neighboring sites through local Andreev reflections. The term S m is the source term appearing on the RHS of Eq. (3.18), S li,m ≡ J li e ikx i δ m,-s l u lk , -δ m,s l v lk T . The recursive character of Eq. (3.20) makes it possible to write the Floquet chain operator in a continued fraction form [START_REF] Hänggi | Driven quantum systems[END_REF][START_REF] Dy | Exact solution for the resolvent matrix of a generalized tridiagonal Hamiltonian[END_REF].

Iterative construction of the resolvent operator

If the resolvent operator R is the inverse of the operator L, then knowledge of R allows straightforward calculation of the amplitudes on the dots:

u m (i; lk) = ∑ i ′ R e i e i ′ m,-s l J li ′ (x i ′ )x lk e iϕ l /2 -R e i h i ′ m,s l J li ′ (x i ′ )y lk e -iϕ l /2 , v m (i; lk) = ∑ i ′ R h i e i ′ m,-s l J li ′ (x i ′ )x lk e iϕ l /2 -R h i h i ′ m,s l J li ′ (x i ′ )y lk e -iϕ l /2 .
(3.21)

Finding the poles of the resolvent operator corresponds to finding the spectrum of the operator L. The resolvent R is an operator which lives in the extended dot ⊗ Nambu ⊗ Floquet space. Upper indices correspond to the combined dot ⊗ Nambu space, and lower indices correspond to the infinite Floquet space.

The in-homogeneous FLS equations for the resolvent elements are:

M 0 m R mn -M + m+1 R m+2,n -M - m-1 R m-2,n = δ mn 1. (3.22)
For a tridiagonal block-matrix Hamiltonian, such as the one we are dealing with, it follows that its resolvent can be generally written in continued fraction form [START_REF] Dy | Exact solution for the resolvent matrix of a generalized tridiagonal Hamiltonian[END_REF]123]. The continued fraction representation is equivalent to the usual Dyson equation [START_REF] Swain | Continued-fraction methods in atomic physics[END_REF].

Starting from Eq. (3.22) it is straightforward to construct the resolvent elements by iteration, assuming a source at some index n and a cutoff at some large Floquet index ±N, with |N| ≥ ∆ ω 0 . The latter is equivalent to assuming that the wavefunction on the dot decays exponentially at energies above the gap |ω + Nω 0 | >> ∆. Physically, the first values of m will correspond to multiple quasi-particle reflection processes, by which the quasi-particle gains energy equal to mω 0 . When m is large enough so that ω + mω 0 > ∆, the quasi-particle enters the superconducting continuum, thus macroscopically resulting into a dissipative quasi-particle flow with normal conductance values. The smaller the voltage value, the more Floquet harmonics we need to take into account. We have the following system consisting of 2N + 1 equations: We solve this system of equations by iteration, and find that the diagonal elements are resummed into a geometric series

M 0 -N R -N,n -M + -N+1 R -N+2,n - (

( ( ( ( ( ( ( ( h h h h h h h h h

M - -N-1 R -N-2,n = 0 . . . M 0 n-2 R n-2,n -M + n-1 R nn -M - n-3 R n-4,n = 0 M 0 n R nn -M + n+1 R n+2,n -M - n-1 R n-2,n = 1 M 0 n+2 R n+2,n -M + n+3 R n+4,n -M - n+1 R nn = 0 . . . M 0 N R Nn - X X X X X X X X M + N+1 R N+2,n -M - N-1 R N-2,n = 0 (3.23) R 0 m M + m+1 M - m+1 R 0 m R 0 m+2 R 0 m+4 M - m+3 M + m+3 M + m+1 M - m+1 R 0 m+2 R 0 m + + Σ + m = + • • •
R mm = M 0 m -M + m+1 Σ + m+2 M - m+1 -M - m-1 Σ - m-2 M + m-1 -1 (3.24)
with forward and backward self-energy matrices, Σ ± , that can be calculated recursively once boundary conditions are imposed, that is once Σ ± ±N is set to zero at some large number ±N. We find that

Σ + m = 1 M 0 m -M + m+1 Σ + m+2 M - m+1 , Σ - m = 1 M 0 m -M - m-1 Σ - m-2 M + m-1 . (3.25)
The non-diagonal resolvent elements R mn can be expressed using the selfenergy matrices

R mn = Σ + m M - m-1 . . . Σ + n+2 M - n+1 R nn if m > n, R mn = Σ - m M + m+1 . . . Σ - n-2 M + n-1 R nn if m < n. (3.26)
Note that the "Floquet self-energies" are diagonal in the Floquet index, while the resolvent is not. This is because the self-energies are associated to the leads whose Hamiltonian is static after gauge-transformation, while the resolvent describes the states on the dot(s), which are time-dependent. This means that the Green functions of the leads depend only on time differences (and therefore their Fourier transform depends only on one Floquet index), while the Green functions of the dot states need to be doubly Fourier transformed.

This description of the resolvent admits a simple diagrammatic representation. For example, by expanding the forward self-energy term Σ + into a series, we see that it regroups all paths that start from a point m and only return to it after having visited all sites m ′ > m, up to the boundary site N.

Σ + m = R 0 m 1 -R 0 m M + m+1 Σ + m+2 M - m+1 = R 0 m + R 0 m M + m+1 Σ + m+2 M - m+1 R 0 m + . . . = R 0 m + R 0 m M + m+1 R 0 m+2 M - m+1 R 0 m + . . . (3.27)
The expansion of the forward self-energy in terms of diagrams is illustrated in Fig. 3.1.

On the other hand, the backward self-energy term Σ -resums loops which pass through sites m ′ < m. Equation (3.26) then describes the shortest path connecting a site n to site m. The role of the self-energy terms Σ ± m is to renormalize the unperturbed diagonal elements of the resolvent R 0 m = 1/M 0 m by introducing a finite imaginary part, corresponding to virtual excursions to the superconducting reservoirs. This imaginary part is introduced in practice by the reservoir's Green's functions, contained in the self-energy, which become imaginary at energies larger than the gap |ω| > ∆. Physically, this corresponds to coupling the initial discrete levels (the ABSs) on the dot(s) to the superconducting continua through MAR. Then, Σ + corresponds to MAR processes which raise the energy of a quasiparticle above the gap ω > ∆, while Σ -corresponds to MAR processes which lower the energy below the gap ω < -∆. The truncation of the continued fractions at some cutoff index |N| > ∆ ω 0 is equivalent to considering that the self-energies become small Σ ± (±N) → 0 at large energies |ω ± Nω 0 | ≫ ∆ as MAR processes become less probable. Therefore, at voltages which are a significant fraction of the gap, one can greatly simplify the expressions of Σ ± , while at small voltages an increasingly greater number of Floquet harmonics need to be taken into account. This expansion is a locator-type expansion of the resolvent, commonly used in disordered systems [START_REF] Anderson | Absence of diffusion in certain random lattices[END_REF][START_REF] Ziman | Models of disorder: the theoretical physics of homogeneously disordered systems[END_REF], in the sense that the unperturbed part of the resolvent R 0 m locates a quasi-particle on site m, in contrast to the more usual "propagator" which describes the propagation of a particle. The resolvent R mn then represents the probability that a quasi-particle is localized on site m, given that it was originally on site n. The major difference between a propagator and a locator expansion is the restriction on repeated indices which is necessary in the later.

There is an alternative and equivalent formulation of the above iterative equations, depending on how one chooses to define the self-energy matrices. I will use both formulations interchangeably, depending on which is more convenient. The second way to write Eqs. (3.24-3.26) can be conveniently summarized as follows:

R mm = M 0 m -Σ + m -Σ - m -1 , R mn = M 0 m -Σ ± m -1 M ∓ m∓1 R m∓2,n if m ≷ n, Σ ± m = M ± m±1 M 0 m±2 -Σ ± m±2 -1 M ∓ m±1 .
(3.28)

As we will see in the following section, the diagonal terms of the resolvent give access to a spectral function, while non-diagonal elements are needed for the calculation of transport properties. The fact that we can express the nondiagonal elements in terms of diagonal elements (due to Eq. (3.26)) means that the current will also, in the end, be closely related to the spectral function.

observables

We will focus on two observables. The first one is the spectral function. We will show that it can be probed by tunneling spectroscopy with a normal tip. This should therefore give experimental access to the Floquet spectrum. The second observable we will examine is the dc component of the current (the MAR current) because it is sensitive to the structure of the Floquet ladder of resonances. Since the system is time-periodic, all observables in the steady state will possess the same periodicity. We therefore expect that the steady state current can be written as a Fourier series I(t) = ∑ p e ipω 0 t I p . The dc current will then correspond to the zeroeth harmonic I 0 .

Probing the Floquet spectrum with tunneling spectroscopy

In order to observe the Floquet spectrum, we imagine that we approach a normal probe to a dot which is coupled to superconducting leads. The probe consists of a metal tip which is coupled to the dot by a tunneling term. A voltage bias is used to tune the chemical potential of the probe, so that a tunneling current can flow between the dot and the probe. We will show that measuring the differential conductance gives direct access to the spectral function of the dot. A similar discussion (as well as a discussion on other kinds of spectroscopy) can be found in [START_REF] Pillet | Tunneling spectroscopy of the Andreev bound states in a carbon nanotube[END_REF]. The technique has been used for the observation of the (equilibrium) ABS [START_REF] Pillet | Andreev bound states in supercurrent-carrying carbon nanotubes revealed[END_REF]. A recent experiment reported observation of Floquet-Andreev states realized by continuous microwave irradiation of a graphene Josephson junction [START_REF] Park | Steady Floquet-Andreev states in graphene Josephson junctions[END_REF], although the validity of the interpretation of the replicas as Floquet states has been challenged [START_REF] Haxell | Microwave-induced conductance replicas in hybrid Josephson junctions without Floquet-Andreev states[END_REF]. I follow a standard derivation, using linear response theory: one expects that switching on a perturbation modifies the expectation values of operators according to the Kubo formula [START_REF] Kubo | Statistical-mechanical theory of irreversible processes. I. General theory and simple applications to magnetic and conduction problems[END_REF] which connects the linear response to a retarded Green's function. The appearance of the retarded Green's function is an expression of causality. In the case of Floquet-Green's functions, it has been proven that a non-negative spectral function exists and that it can be interpreted as a density of states [START_REF] Uhrig | Positivity of the spectral densities of retarded Floquet Green functions[END_REF]. Indeed, I will show that the spectral function can be defined through the diagonal (in Floquet) elements of the Green's function of the dot.

The total Hamiltonian of the system is given by the Hamiltonian of the S-QD-S system, the Hamiltonian of the probe, and the tunneling term between dot and probe: H = H(t) + H probe + H tun . The probe is modeled as a simple free fermion system

H probe = dω ωψ † (ω)ψ(ω) (3.29)
with populations defined by the usual relations

⟨ψ † (ω)ψ(ω ′ )⟩ = δ(ω -ω ′ )n(ω) ⟨ψ(ω)ψ † (ω ′ )⟩ = δ(ω -ω ′ )(1 -n(ω)) (3.30)
where n(ω) is the Fermi-Dirac distribution. If the probe is biased with a voltage V tun , then n(ω) = f FD (ω -eV tun ) = (1 + e β(ω-eV tun ) ) -1 . At zero temperature, this simply reduces to a step function n(ω) = θ(eV tunω).

The tunneling term between the probe and the dot can be written as

H tun = dω J tun (ω) d † ψ(ω) + ψ † (ω)d , (3.31)
and the tunneling current from the tip of the probe to the quantum dot is

I tun = -ie dω J tun (ω) d † ψ(ω) -ψ † (ω)d . (3.32)
Using linear response 3 we obtain the expectation value of the tunneling current

⟨I tun (t)⟩ = -e t -∞
ds dω dω ′ J tun (ω)J tun (ω ′ )

× d † (t)ψ(ω, t) -ψ † (ω, t)d(t), d † (s)ψ(ω ′ , s) + ψ † (ω ′ , s)d(s) . (3.33)
The expectation values appearing on the RHS can be calculated using Wick's theorem. One obtains

⟨I tun (t)⟩ = e t -∞ ds dω J 2 tun (ω) n(ω) ⟨d(s)d † (t)⟩ e -iω(t-s) + ⟨d(t)d † (s)⟩ e iω(t-s) + (n(ω) -1) ⟨d † (t)d(s)⟩ e -iω(t-s) + ⟨d † (s)d(t)⟩ e iω(t-s) .
(3.34)

3 At time t = -∞ the perturbation H ext is turned on. The system is initially in thermal equilibrium. In the perturbed system the expectation value of an operator A at some time t is given by linear response theory

⟨A(t)⟩ = ⟨A⟩ -i t -∞ ds ⟨[A(t), H ext (s)]⟩
where ⟨A⟩ is the expectation value of the operator in the absence of the perturbation.

A much more simple formula can be obtained if we consider the differential conductance. Since only the populations n(ω) depend on the tunneling voltage, we get

∂ ⟨I tun (t)⟩ ∂V tun = e ds dω J 2 tun (ω) ∂n(ω) ∂V tun × × ⟨ d † (t), d(s) ⟩ θ(t -s)e -iω(t-s) + ⟨ d(t), d † (s) ⟩ θ(t -s)e iω(t-s) .
(3.35)

In the above expression, one can recognize the components of the retarded Green's function in Nambu space:

G r,11 (t, s) = -i ⟨ d(t), d † (s) ⟩ θ(t -s) G r,22 (t, s) = -i ⟨ d † (t), d(s) ⟩ θ(t -s).
(3.36)

Time-independent Hamiltonian

If the Hamiltonian of the dot is static, the Green's function depends only on time differences. We can define its Fourier transform as

G(t -s) = dω 2π e -iω(t-s) G(ω), with G(ω) = dt e iωt G(t) (3.37)
One then gets 4

∂ ⟨I tun ⟩ ∂V tun = -2e dω J 2 tun (ω) ∂n(ω) ∂V tun Im G r,11 (ω)
= -e dω J 2 tun (ω) 11 (ω) + G r,22 (-ω) .

∂n(ω) ∂V tun Im G r,
(3.38)

The last equality in the above equation was obtained by considering that, first, the two Nambu components of the retarded Green's function are clearly related by complex conjugation G 11 (t) * = -G 22 (t) which then implies that, by considering their Fourier transform

G 11 (ω) * = dt e -iωt G 11 (t) * = -dt e -iωt G 22 (t) = -G 22 (-ω). (3.39)
We therefore have

G 11 (ω) -G 11 (ω) * = -G 22 (-ω) * + G 22 (-ω) which gives us the identity Im G 11 (ω) = Im G 22 (-ω). (3.40)
4 The Fourier transform of the delta function is useful here

δ(ω -ω ′ ) = 1 2π dt e -i(ω-ω ′ )t
From Eq. (3.38), one is lead to define the spectral function on the dot as

A(ω) = - 1 π Im G r,11 (ω) + G r,22 (-ω) = - 2 π Im G r,11 (ω) . (3.41)
At zero temperature, the derivative of the Fermi-Dirac distribution is a delta function

∂n(ω)
∂V tun = -eδ(eV tunω) so the differential conductance becomes particularly simple: 11 (eV tun ), temperature T=0.

∂ ⟨I tun ⟩ ∂V tun = 2e 2 J 2 tun (eV tun ) Im G r,
(3.42)

Time-periodic Hamiltonian

In the case of a voltage-biased junction, one can use a double Fourier transform for the Green's function on the dot

G(t, s) = ∑ m,n dω 2π e -i(ω+mω 0 )t e i(ω+nω 0 )s G mn (ω), (3.43)
with Fourier coefficients

G mn (ω) = 1 T ds +T/2 -T/2 dt e i(ω+mω 0 )t e -i(ω+nω 0 )s G(t, s). ( 3
A(ω) = - 1 π Im G r,11 00 (ω) + G r,22 00 (-ω) = - 2 π Im G r,11 00 (ω) (3.47)
For the S-QD-S system, G r can be simply identified with the resolvent operator, while for multidot systems, the Green's function on dot i can be identified with the part of the resolvent corresponding to the subspace of dot i. The quantity A(ω) = -1 π Im R 00 (ω) can be seen as a time-average of the spectral function over one period of the drive.

Current harmonics and MAR current

The total current passing through the dot i is related to the change of the electron number on the dot Ni (t

) = ∑ σ d † iσ (t)d iσ (t). In particular, ∑ j I j→i (t) = d dt Ni (t), (3.48)
where I j→i is the current flowing from the reservoir j to the dot i. Without loss of generality, we can take the dot i to be at position x i = 0. In the steady state, the average current is given by the commutator of the number operator and the tunneling term,

∑ j I j→i (t) = d dt ⟨ ∑ σ e iHt d † iσ d iσ e -iHt ⟩ = i ∑ σ ⟨[H(t), d † iσ d iσ ] ⟩ = i ∑ σ ⟨[V (t), d † iσ d iσ ] ⟩ = -2 ∑ jk ′ σ Im J ji e -is j ω 0 t c † jk ′ σ d iσ . (3.49)
We can isolate the contribution from the j reservoir, and write

I j→i (t) = -2 ∑ k ′ σ Im J ji e -is j ω 0 t c † jk ′ σ d iσ . (3.50)
Given that the dressed quasi-particle operators Γ † , Γ form a complete basis, we can express all other operators in this basis, see Eq. (3.12) and (3.13). The advantage of such decompositions is that one can then very easily derive expressions for expectation values in the stationary state |S⟩, which is simply defined as the state which is annihilated by the application of the Γ operator, meaning that Γ lkσ |S⟩ = 0 for any {lkσ}. Therefore, any average taken in the steady state will contain only contributions from ⟨S|ΓΓ † |S⟩ terms. We then find that

∑ σ e -is j ω 0 t c † jk ′ σ d iσ = 2 ∑ lk ∑ mn e -i(m-n+s j )ω 0 t δ m0 δ jl δ k ′ k v jk ′ + V m (jk ′ ; lk) v * n (i; lk) = 2 ∑ lk ∑ mn e -i(m-n)ω 0 t δ m,s j δ jl δ k ′ k v jk ′ + V m-s j (jk ′ ; lk) v * n (i; lk) (3.51)
and therefore the current is given by

I j→i (t) = -4 ∑ k ′ ,lk ∑ mn Im J ji e -i(m-n)ω 0 t δ m,s j δ jl δ k ′ k v jk ′ + V m-s j (jk ′ ; lk) v * n (i; lk) .
(3.52) Using the FLS equations (3.18), we can re-express the source term:

J ji δ m,s j v jk ′ = -(E jk ′ + mω 0 )v m (i; jk ′ ) -∑ lk J li V m-s l (lk; jk ′ ). (3.53)
This gives

I j→i (t) = 4 Im ∑ k ′ ∑ mn e -i(m-n)ω 0 t (E jk ′ + mω 0 )v m (i; jk ′ )v * n (i; jk ′ ) -4 Im ∑ k ′ ,lk ∑ mn e -i(m-n)ω 0 t J ji V m-s j (jk ′ ; lk)v * n (i; lk) -J li V m-s l (lk; jk ′ )v * n (i; jk ′ ) . (3.54)
Since 2i Im z = zz * , the first term on the RHS can be re-written in a more meaningful way:

4 Im ∑ mn e -i(m-n)ω 0 t (E + mω 0 )v m v * n = = -2i ∑ mn e -i(m-n)ω 0 t (E + mω 0 )v m v * n -e -i(m-n)ω 0 t (E + nω 0 )v * n v m = -2i ∑ mp e -ipω 0 t pω 0 v m v * m-p . (3.55)
This means that the first term goes to zero for the dc component of the current, i.e when the Fourier component p = 0. In fact, now the current can be decomposed into a series of harmonics:

I j→i (t) = -2i ∑ k ′ ∑ mp e -ipω 0 t pω 0 v m (i; jk ′ )v * m-p (i; jk ′ ) -4 Im ∑ k ′ ,lk ∑ mp e -ipω 0 t J ji V m-s j (jk ′ ; lk)v * m-p (i; lk)
term where {jk ′ } ⇌ {lk} ≡ ∑ p e -ipω 0 t I j→i,p .

(3.56)

The MAR current corresponds to the dc component and, therefore, to the zeroth term I j→i dc = I j→i,0 . The last step is to re-express the amplitudes on the RHS. The current will then contain only quadratic terms in the resolvent. Using Eq. (3.17) the reservoir amplitudes are integrated out, and only the amplitudes on the dot(s) are left. The dot amplitudes can in turn be expressed as a function of the resolvent operator and the source terms due to Eq. (3.18). The reservoir amplitudes give

∑ k ′ J ji V m (jk ′ ; lk) = -∑ i ′ Γ j,ii ′ g 21 j,ii ′ (m)u m-s j (i ′ ; lk) + g 22 j,ii ′ (m)v m+s j (i ′ ; lk)) . (3.57)
We can then write the p harmonic of the current as

I j→i,p = -2ipω 0 ∑ k ∑ m v m (i; jk)v * m-p (i; jk) +4 Im ∑ k ∑ i ′ ,l,m ⟨h| Γ j,ii ′ g j,ii ′ (m -s j ) u m-2s j (i ′ ; lk) v m (i ′ ; lk) v * m-p (i; lk) -Γ l,ii ′ g l,ii ′ (m -s l ) u m-2s l (i ′ ; jk) v m (i ′ ; jk) v * m-p (i; jk) . (3.58)
At this point we have introduced a Dirac bra-ket notation. The resolvent is an operator acting in the combined space dot ⊗ Nambu ⊗ Floquet. In order to have more compact expressions, we introduce the matrix σ j,ii ′ (ω) ≡ Γ j,ii ′ g j,ii ′ (ω), so that σ j is acting in the dot ⊗ Nambu space. Finally, we use the notation ⟨i, e, m|R|i ′ , h, n⟩ = R e i h i ′ mn . The amplitudes can be directly calculated if we know the resolvent elements: With the above notation we can re-write Eq. (3.59):

u m (i; lk) v m (i; lk) = R ii mn J li (x i ) δ n,-s l u lk -δ n,s l v lk + ∑ i ′ ̸ =i R ii ′ mn J li ′ (x i ′ ) δ n,-s l u lk -δ n,s l v lk (3.
u m (i; lk) v m (i; lk) = ∑ i ′ ⟨i, m| RC -s l i ′ , 0 S li ′ , (3.61)
where the source term

S li ′ ≡ J li ′ (x i ′ ) u lk , -v lk T .
We can then write the term:

u m-s j (i; lk) v m+s j (i; lk) = ∑ i ′ i, m C s j RC -s l i ′ , 0 S li ′ . (3.62)
With this notation, the p harmonic of the current becomes:

I j→i,p = -2ipω 0 ∑ m ∞ 0 dω 2π ⟨i, h| RC -s j m0 Q j (ω) RC -s j † m-p,0 |i, h⟩ +4 Im ∑ l,m ∞ 0 dω 2π ⟨i, h|σ j (m -s j ) C s j RC -s l m-s j ,0 Q l (ω)(RC -s l ) † m-p,0 |i, h⟩ -4 Im ∑ l,m ∞ 0 dω 2π ⟨i, h|σ l (m -s l ) C s l RC -s j m-s l ,0 Q l (ω) RC -s j † m-p,0
|i, h⟩ .

(3.63)

The source terms gives rise to a matrix which corresponds to the populations of the reservoir of injection

S li S † li ′ = Q l,ii ′ , with Q l,ii ′ (ω) = i(σ r l,ii ′ (ω) -σ a l,ii ′ ) θ(ω),
so that it is zero unless the excitation energy is ω > ∆. The dc component of the current is then given by

I dc j→i = 4 Im ∑ l,m ∞ 0 dω 2π σ j (m -s j ) C s j RC -s l m-s j ,0 Q l (RC -s l ) † m0 h i -term where j ⇌ l (3.64)
Equation (3.64) is at first sight not easy to calculate since it involves an integral over all quasi-particle excitation energies, as well as a summation over all the Floquet harmonics. Fortunately, the resolvent elements decay exponentially at large energies outside the gap, so that the integration can be drastically truncated. Moreover, at large enough voltages we observe a localization (analogous to the Wannier-Stark localization) which gives a rapidly convergent summation over the Floquet harmonics. However, Eq. (3.64) even though compact, is rather opaque physically speaking. We will try to rectify this in the next section when talking concretely about the MAR current in the simple case of an S-QD-S system. We will show that, by using the locator expansion of the resolvent, we can bring the MAR current in a form that is a function of the spectral function of the dot and a string of anomalous Green's functions (representing the MARs).

voltage-biased s-qd-s

We will start by applying the formalism we presented in the previous sections to the case of a single quantum dot coupled to two superconducting leads (S-QD-S). This is a well-known system, at least when it comes to the calculation of the MAR current, and for a review the interested reader can look at the relevant section in [START_REF] Martín-Rodero | Josephson and Andreev transport through quantum dots[END_REF].

The symmetric configuration of a resonant dot with energy ϵ d = 0 coupled to two reservoirs which are voltage biased at opposite voltages (V L = -V R ) is well-understood: the subgap structure of the current-voltage curve shows steps at voltage differences which are integer subdivisions of the gap, 2eV = 2 ∆ m , with m a positive integer. These steps in the subgap current appear whenever a new "MAR trajectory" becomes possible: a quasi-particle is Andreev reflected m times, changing its energy by e(V R -V L ) = 2eV with each reflection, until it has enough energy (equivalent to the size of the gap 2 ∆) to reach the superconducting continuum of states and give a contribution to the current. The presence of the resonant level gives a "resonant condition": only trajectories which pass through the resonant level contribute. This restricts m to being an odd integer, and suppresses even MAR processes [START_REF] Martín-Rodero | Josephson and Andreev transport through quantum dots[END_REF][START_REF] Levy Yeyati | Resonant tunneling through a small quantum dot coupled to superconducting leads[END_REF][START_REF] Jonckheere | Nonequilibrium supercurrent through a quantum dot: Current harmonics and proximity effect due to a normal-metal lead[END_REF][START_REF] Johansson | Resonant multiple Andreev reflections in mesoscopic superconducting junctions[END_REF]. Displacing the resonant level ε d ̸ = 0 produces additional structure [START_REF] Levy Yeyati | Resonant tunneling through a small quantum dot coupled to superconducting leads[END_REF]. In fact, we will argue that the structure in the vicinity of the MAR steps is connected to the positions of the Floquet ladders. Any asymmetry in the system that produces shifts in the Floquet spectrum then produces a modification of the MAR steps. The asymmetry can be produced by somehow breaking the mirror symmetry of the system. For example by moving the resonant level (or in general by hav-

ing |ε d -V L | ̸ = |ε d -V R |),
or by asymmetric couplings to the reservoirs. We will first consider the special case of a symmetrically biased junction V L = -V R as well as the asymmetrically biased junction V L = -V, V R = 0. Another reason for being interested in the (-V, 0) configuration is that it will be relevant in the study of the Andreev molecule. In the regime of a large interdot distance, the two junctions will be uncoupled. For an Andreev molecule biased in the quartet configuration (-V, 0, +V), the current through the first dot will then correspond to the current of the junction biased at (-V, 0).

FLS equations

Below we write explicitly the equations needed to study the S-QD-S case. The starting point is Eq. (3.18) which becomes

(E lk + mω 0 + iη)u m (lk) = +δ m,-s l J l u lk + ∑ j Γ j g 11 j (m + s j )u m (lk) + g 12 j (m + s j )v m+2s j (lk) (E lk + mω 0 + iη)v m (lk) = -δ m,s l J l v lk + ∑ j Γ j g 21 j (m -s j )u m-2s j (lk) + g 22 j (m -s j )v m (lk) (3.65)
where the choice of voltages is s j = ±1 in the symmetric case or s j = {0, -1} in the asymmetric case. In the basis of a two-component Nambu spinor

Ψ m = u m , v m T we write (LΨ) m ≡ M 0 m Ψ m -M + m+1 Ψ m+2 -M - m-1 Ψ m-2 = J l δ m,-s l u lk -δ m,s l v lk , (3.66)
where the matrices M 0 , M ± depend on the particular configuration.

floquet chain in the asymmetric case.

Taking s L = -1, s R = 0 we get 5 M 0 m = (E + mω 0 )1 2 - Γ L g 11 L (m-1)+Γ R g 11 R (m) Γ R g 12 R (m) Γ R g 21 R (m) Γ L g 22 L (m+1)+Γ R g 22 R (m) , M + m = 0 0 Γ L g 21 L (m) 0 , and
M - m = 0 Γ L g 12 L (m) 0 0 . (3.67)
By expressing the amplitudes on the dot as a function of the resolvent,

u m (lk) v m (lk) = R mn J l δ n,-s l u lk -δ n,s l v lk , (3.68)
we can make a few remarks: As can be seen from Eq. (3.26), once we know the source index s l , we can express all non-diagonal elements R m,±s l as a function of the diagonal elements R ±s l ,±s l . From the same equation, we see that for R m,±s l to be non-zero, m has to have the same parity as the source index ±s l . Therefore m can change by steps of 2. This doubling is related to the fact that Ψ m is coupled to Ψ m±2 by second-order Andreev processes.

floquet chain in the symmetric case. Taking s L = -1, s R = +1 we get the matrices

M 0 m = (E + mω 0 )1 2 - Γ L g 11 L (m-1)+Γ R g 11 R (m+1) 0 0 Γ L g 22 L (m+1)+Γ R g 22 R (m-1) , M + m = 0 Γ R g 12 R (m) Γ L g 21 L (m) 0 , and
M - m = 0 Γ L g 12 L (m) Γ R g 21 R (m) 0 . (3.69)
Here, because we have taken the voltage difference across the junction to be 2eV, we expect to have a doubling with respect to the previous case. Moreover, due to the high symmetry of this case, we find that the resolvent has a particular structure, both in Nambu and in Floquet space. First of all, we remark that the self-energy matrices are diagonal in Nambu space. 6 So are all elements R mm .

The rest of the elements R mn with m ̸ = n are either diagonal in Nambu space

(if m-n 2 = 2p, with p ∈ Z) or anti-diagonal (if m-n 2 = 2p + 1)
. In particular, this implies the existence of "selection rules" for the amplitudes u m (l), v m (l). Indeed, if we express these amplitudes as a function of the resolvent and the source term, a careful analysis will give us a selection rule for the index m : As in the previous case m must have the same parity as the source index

Ψ m (l) ̸ = 0 if m -s l = 2p, p ∈ Z.
(3.70) Moreover, the "contents" of the wavefunction alternate between electron-like and hole-like:

u m (l) ̸ = 0 if m + s l = 4p v m (l) ̸ = 0 if m -s l = 4p (3.71)
so that there is a doubling of the periodicity with respect to the previous case.

Floquet spectrum

Using the definition of the matrices M 0 , M ± and the iterative formulas of Eq.

(3.28) we can calculate the diagonal part of the resolvent R 11 00 (ω) whose imaginary part gives us the spectral function of the dot, see Eq. (3.47). All calculations are performed in Mathematica. In practice, it is sufficient to choose a cutoff N ∼ ∆ ω 0 for the calculation of the self-energies.

Single junction biased at (-V,+V).

As we have noted, in this case the resolvent elements R mm are always diagonal in Nambu space. It is also worth noting that they are independent of the phases of the superconducting leads. In general, due to the presence of particlehole symmetry, we have the relation R 11 00 (E, ω 0 ) = R 22 00 (E, -ω 0 ). If, moreover, the coupling strengths to the left and right leads are equal, Γ L = Γ R , there is a mirror symmetry in the system that results in R ii 00 (E, ω 0 ) = R ii 00 (E, -ω 0 ). The combination of the two symmetries then gives R 11 00 (E, ω 0 ) = R 22 00 (E, ω 0 ). Since Im R 11 00 (E) = Im R 22 00 (-E) always holds (Eq. (3.40)), then the two relations together imply that the spectrum always has a resonance at E = 0. By the translation property of the resolvent, it follows that the spectrum consists of resonances at all multiples of the basic frequency 2ω 0 . The result is that the Floquet spectrum consists of uncoupled Floquet-Wannier-Stark ladders which are situated at even multiples of the frequency

E = 2pω 0 .
This is illustrated in Fig. 3.2a where we plot the spectral function using the frequency ω 0 as a scaling parameter [START_REF] Shirley | Solution of the Schrödinger equation with a Hamiltonian periodic in time[END_REF]. The inverse voltage scaling makes apparent the periodicity of the system, T = 2π ω 0 , and gives us something akin to a "Floquet bandstructure". We therefore plot the spectrum in reduced units of E/ω 0 as a function of the inverse voltage ∆/ω 0 . The density plots used here essentially show the maxima of the spectral function, corresponding to the resonances of the resolvent element R 11 00 (E). The spectral function (without scaling) is shown at various voltages in Fig. 3.2b, illustrating its resonant structure, with peaks becoming sharper at small voltages, and increasingly smeared as the voltage frequency approaches the value of the superconducting gap.

The degeneracy of the Floquet ladders is, however, lifted as soon as some asymmetry is introduced in the system. For example, if the couplings are unequal, Γ L ̸ = Γ R , the bandstructure becomes dispersive and asymmetric around zero. This is in stark difference with the equilibrium case, where the symmetric coupling case turns out to be the most general, and all properties of the asymmetric system can be derived from the symmetric one [START_REF] Kadlecová | Quantum dot attached to superconducting leads: relation between symmetric and asymmetric coupling[END_REF]. As shown in Fig. 3.2c, when Γ L ̸ = Γ R , the Floquet ladders are no longer decoupled and avoided crossings appear, signaling coupling between them via Landau-Zener-Stückelberg-Majorana transitions [START_REF] Ivakhnenko | Nonadiabatic Landau-Zener-Stückelberg-Majorana transitions, dynamics, and interference[END_REF][START_REF] Giacomo | The Majorana formula and the Landau-Zener-Stückelberg treatment of the avoided crossing problem[END_REF]. Note that the coupling is between different ladders, and not between the electron-hole sectors (since the resolvent remains diagonal in Nambu space). The resonances are now slightly shifted away from even multiples of the frequency, i.e. resonances are found at

E p,± = 2pω 0 ± ϵ,
with a basic period of 4ω 0 . The parameter ϵ can be seen as an average of the equilibrium ABS energy over the phase variable [START_REF] Mélin | Simple Floquet-Wannier-Stark-Andreev viewpoint and emergence of low-energy scales in a voltage-biased three-terminal Josephson junction[END_REF]. This becomes exact away from avoided crossings, when the interladder coupling can be neglected. Similar expressions were found long ago by Wannier [START_REF] Wannier | Wave functions and effective Hamiltonian for Bloch electrons in an electric field[END_REF]. This situation is illustrated in Fig. 3.2c.

Figure 3.2d shows the spectral function A(E) = -2 π Im R 11 (E) at various frequencies (solid lines), together with the "hole" part A(-E) = -2 π Im R 22 (E) (dashed lines). We observe (not shown here) that the distance between the peaks increases when the difference between the couplings increases. The "electron" peak moves to the right if Γ L < Γ R and towards the left if Γ L > Γ R . The opposite is true for the "hole" peak. Clearly, in this case there is an asymmetry of the spectral function around zero since A(E) ̸ = A(-E). This prediction should be experimentally verifiable by measuring the differential conductance at opposite tunneling voltages.

Single junction biased at (-V,0).

In this case the spectrum consists of ladders of resonances at quasi-energies E p,± = 2pω 0 ± ϵ, with a basic period of 2ω 0 . This is illustrated in Fig. 3.3a which shows the spectral function of the system, plotted with inverse scaling. The Floquet ladders show avoided crossings which is a sign of coupling between them. Even though this looks similar to the previous case, here the resolvent elements R mm are no longer diagonal in Nambu space. This results in coupling between the electron and hole sectors. As a result, even though both electronic and hole parts of the resolvent have peaks at the same energies, their weight away from avoided crossings differs when varying the voltage. To illustrate this, in Fig. 3.3b we plot the difference between the electron and hole parts of the spectral function

A(E) -A(-E) = - 2 π Im R 11 00 (E) -R 22 00 (E) .
Near avoided crossings, a rapid change of the Floquet states is expected to happen. Accordingly, we see that the sign of the aforementioned quantity changes signs at avoided crossings, signaling the change in character between electronlike and hole-like states. Our results suggest that, in order to observe this asymmetry it would be necessary to break the mirror symmetry of the system. This could be achieved either by having asymmetric tunnel couplings in a symmetric bias situation V L = -V R , or by an asymmetric bias situation V L ̸ = -V R , with any configuration of tunnel couplings. One could then probe the asymmetry away from avoided crossings by measuring the conductance at opposite tunneling voltage values.

floquet engineering. The avoided crossings present in Fig. 3.2c or 3.3 could be used to find dynamical sweet spots of the system. These so-called sweet spots are optimal working points corresponding to the extrema in quasienergy differences and have been proposed as a way to protect qubits from noise. Contrary to the static case where few sweet spots are present, the extra dimension of time in periodically driven systems allows to find a manifold of dynamical sweet spots [START_REF] Huang | Engineering dynamical sweet spots to protect qubits from 1/ f noise[END_REF]. An added advantage is that one can tune the system to an avoided crossing by changing the drive in situ. Realization of a Floquet qubit has been proposed along this line [START_REF] Gandon | Engineering, control, and longitudinal readout of Floquet qubits[END_REF], where a periodic driving can be used to tune the system near one of the avoided crossings, and a second drive can be used to control transitions between the Floquet states. Typically, the smaller the quasi-energy dispersion, the more insensitive would the qubit states be to fluctuations. Moreover, tuning a fluxonium qubit to a dynamical sweet spot away from its half-flux bias static spot has been shown to increase coherence times [START_REF] Mundada | Floquetengineered enhancement of coherence times in a driven fluxonium qubit[END_REF], demonstrating the relevance of Floquet engineering to the qubit community.

Δ 9ω 0 7ω 0 ω 0 /Δ=0.16 Δ 5ω 0 ω 0 /Δ=0.24 |R -1-1 | |R -3-1 | |R -5-1 | |R -7-1 | Δ 3ω 0 ω 0 /Δ=0.4

Resonant structure of the resolvent

Transport properties such as the current require calculation of the non-diagonal elements of the resolvent R mn with m ̸ = n, for energies above the gap E > ∆.

We remind that the indices m, n represent positions on the Floquet chain, and an element R mn represents the probability to be at a position m on the Floquet chain, having started at some source point n, see Sec. 3.2.3. In the context of the current calculation, the source point n takes values that depend on the voltage biasing of the leads, n = V j /V. For example, in the symmetric case (-V, +V) the source is at n = ±1, and therefore the current requires calculation of the resolvent elements which connect the "source sites" on the Floquet chain to sites on positions m, R m,±1 . Here, we need to make several points:

1. Outside the gap the spectrum shows signs of dissipation because of strong hybridization with the reservoirs, i.e. the resolvent elements R 11,22 00 (E) decay exponentially for large energies E ≫ ∆. Depending on the configuration, Floquet resonances "survive" at, or around, frequencies ω 0 > ∆/m, with m the smallest even number such that the condition holds. The peaks are sharper at small voltages, and become smeared when increasing ω 0 .

2. In the previous section we have seen that the resolvent R 00 has resonances at (or around) even multiples of ω 0 . Then, the symmetry by translation of the resolvent R m,n (ω + pω 0 ) = R m+p,n+p (ω) means that elements R ±1,±1 (ω) = R 00 (ω ± ω 0 ) have resonance peaks at (or around) odd multiples of ω 0 .

3. Since the non-diagonal elements of the resolvent R m,±1 are obtained from the diagonal elements R ±1,±1 through Eq. (3.26) we can conclude that they are resonant only when the corresponding diagonal elements are resonant, i.e. their resonances are at (or around) odd multiples of the frequency.

4. Also as a consequence of Eq. (3.26), there is a hierarchy of peaks in m which depends on the voltage. This means that an element R m,±1 becomes dominant when mω 0 is the dominant peak above the gap, i.e. when m is the minimal odd integer for which mω 0 > ∆. For example, the element R -3,-1 is dominant when the second order MAR trajectory is dominant, ∆ 3 < ω 0 < ∆, the element R -5,-1 is dominant when the third order MAR trajectory is dominant, ∆ 5 < ω 0 < ∆ 3 , and so on.

These observations are illustrated in Figure 3.4 for the non-diagonal elements R m,-1 .

MAR current

We will focus on the MAR current of the S-QD-S system in the symmetric case (-V, +V). Using Eq. (3.64) we can write explicitly the dc component of the current, which simplifies to

I dc = 4 Im ∑ m ∞ ∆ dω 2π × × Γ L g 21 L (m + 1), g 22 L (m + 1) R 11 m+2,-1 R 12 m+2,1 R 21 m,-1 R 22 m,1 Q R R 21 m,-1 R 22 m,1 * -Γ R g 21 R (m -1), g 22 R (m -1) R 11 m-2,1 R 12 m-2,-1 R 21 m,1 R 22 m,-1 Q L R 21 m,1 R 22 m,-1 * . ( 3.72)
The matrix Q l can be seen as a spectral function of the reservoir of injection l We integrate Eq. (3.72) numerically in order to produce the current-voltage characteristics. The results of our calculation for different values of tunnel couplings Γ = Γ L = Γ R are depicted in Fig. 3.6a, and we verify that our method gives the expected results for the I-V curves at zero temperature, with jumps at odd subdivisions of the gap

Q l (ω) =2iΓ l (g r l (ω) -g a l (ω)) θ(ω) = 2Γ l √ ω 2 -∆ 2 ω -∆ l -∆ * l ω , ω > ∆.
Γ L =0.3Δ=Γ R Γ L =0.3Δ, Γ R =0.4Δ Γ L =0.3Δ, Γ R =0.5Δ
V R -V L = 2eV = 2 ∆ 2p + 1 ⇒ eV = ω 0 = ∆ 2p + 1 .
At each step, a new MAR trajectory becomes possible, resulting in an increase in the current. The increase is sudden due to the BCS density of states which diverges at the gap edges. The steps become wider with increasing coupling. It has perhaps been less commented on in the literature that the subgap current steps appear at exactly ω 0 = ∆ 2p+1 only when there is a mirror symmetry which happens when the tunnel couplings to the reservoirs are equal, Γ L = Γ R and the voltage-biasing is symmetric around the dot energy V L = -V R . In this mirror-symmetric case, electron-like and hole-like MAR trajectories are equally favorable. We have seen that the corresponding spectrum consists of decoupled ladders situated exactly at even multiples of ω 0 , see Fig. 3.2a. When this symmetry is broken, we found that the electron-like part -Im R 11 00 (ω) and the hole-like part -Im R 22 00 (ω) of the spectrum have peaks at different energies E p,± = 2pω 0 ± ϵ. We find that the current carries a trace of this characteristic of the spectrum: the MAR steps break into two sub-steps, positioned around the original ω 0 = ∆±ϵ 2p+1 frequencies. The exact shape of the steps (cusps or peaks) depends on the choice of the couplings [START_REF] Houzet | Multiple Andreev reflections in hybrid multiterminal junctions[END_REF], and the distance depends on the difference between the couplings |Γ L -Γ R |, as shown in Fig. 3.5b. We will clarify the connection between the current and the spectral function in the next section. Indeed, Eq. (3.72) is hardly transparent. We will now try to find an expression that gives some physical intuition. To do this, we will go back to an expression of the dc current with respect to the amplitudes on the dot (essentially rewriting Eq. (3.58)):

I dc j→dot = 4 Im ∑ l,k ∑ m Γ j g 22 j (m) v m+s j (l, k) 2 -Γ l g 22 l (m)|v m+s l (j, k)| 2 + Γ j g 21 j (m)u m-s j (l, k)v * m+s j (l, k) -Γ l g 21 l (m)u m-s l (j, k)v * m+s l (j, k) .
(3.74)

In the case we examine we can set s j = s L = -1, s l = s R = +1, so that

I dc L→dot = 4 Im ∑ k ∑ m Γ L g 22 L (m)|v m-1 (R, k)| 2 -Γ R g 22 R (m)|v m+1 (L, k)| 2 +Γ L g 21 L (m)u m+1 (R, k)v * m-1 (R, k) -Γ R g 21 R (m)u m-1 (L, k)v * m+1 (L, k) . (3.75)
The first two terms on the RHS containing the normal part of the Green's function are clearly non-zero only when Im g 22 L,R (m) ̸ = 0. This means that the energy must satisfy the condition |ω + mω 0 | > ∆. Since the energy is ω > ∆, this is always satisfied for positive m. However, such terms will not have a significant contribution to the current for mω 0 > ∆, because at large energies the resolvent decays exponentially. Rather, this term mostly contributes when mω 0 < -2 ∆, which implies that m is mostly negative. Moreover, recall that Eq. (3.71) says that the amplitudes v m (l) are non-zero when m has a distance 4p from the source voltage s l . Therefore both |v m-1 (R)| 2 and |v m+1 (L)| 2 are non-zero for m = 4p + 2 = . . . , -6, -2, +2, . . . . We see that the sum over m will reduce to a sum over only a few harmonics. This term then gives the usual MAR steps at odd subdivisions of the gap as expected, ω 0 = ∆ 2p+1 . The last two terms containing the anomalous Green's function, proportional to

u m+1 (R)v * m-1 (R) and u m-1 (L)v * m+1 (
L) also follow the same rules, although it is not as evident at first glance. As a result, the sum over m in Eq. (3.75) can then be restricted by the condition that

m 2 = 2p + 1, p ∈ Z.
contribution of floquet harmonics. The number of harmonics one needs to consider depends on the voltage. For strong driving, only a few harmonics need to be taken into account, and the sum in Eq. (3.75) can be drastically truncated. We therefore have a "localization" on the Floquet chain, analogous to the Wannier-Stark localization of electrons in solids at strong electric

Γ L = Γ R = 0.3 ∆.
fields. This localization is illustrated in Fig. 3.6. At large voltages ω 0 > ∆, the drive is strong enough to promote quasi-particles directly above the gap without any MAR processes, and we only need to sum over two harmonics m = ±2.

As we lower the voltage, we progressively need to add more harmonics, in correspondence to the MAR processes which are dominant. In the region of the first allowed MAR process, ∆ 3 < ω 0 < ∆, the current is well approximated by summing over three harmonics m = ±2, -6, in the next region of ∆ 5 < ω 0 < ∆ 3 , we need to add one more m = ±2, -6, -10, and so on.

In order to better understand the structure of the current, one needs to analyze the integrand of the current I dc L→dot = dω 2π ∑ m I (m) :

I (m) = 4 Im Γ L g 22 L (m)|v m-1 (R)| 2 -Γ R g 22 R (m)|v m+1 (L)| 2 + 4 Im Γ L g 21 L (m)u m+1 (R)v * m-1 (R) -Γ R g 21 R (m)u m-1 (L)v * m+1 (L) (3.76)
The sum over the positive m is nonzero, and it can be verified numerically that it gives an almost constant contribution to the current, without any structure near the MAR steps. We can then concentrate on the sum over the negative m. For a given m, we can approximate the amplitudes by keeping only the dominant contributions with respect to the resolvent elements. Note that this is not done in the numerical calculation of the current, where we kept all terms. However, it is useful for gaining some intuition. In order to do this, we rewrite the amplitudes appearing in Eq. (3.76) as

|v m (l)| 2 = R 21 m,-s l 2 Q 11 l + R 22 m,s l 2 Q 22 l + 2 Re R 21 m,-s l (R 22 m,s l ) * Q 21 l u m (l)v * n (l) = R 11 m,-s l Q 11 l (R 21 n,-s l ) * + R 11 m,-s l Q 12 l (R 22 n,s l ) * + R 12 m,s l Q 21 l (R 21 n,-s l ) * + R 12 m,s l Q 22 l (R 22 n,s l ) * (3.77)
The non-diagonal elements of the resolvent represent paths which connect the "source sites" at n = V L,R /V = ±1 on the Floquet chain to "measure sites" at positions m. These elements therefore correspond to processes where the system changes its energy by an amount |m -n|ω 0 , equivalent to absorbing or emitting an |m -n| number of "photons". We then expect that the dominant term in each of the above expressions is the one involving the shortest possible path. Explicitly, for m < 0 we approximate the amplitudes by

|v m-1 (R)| 2 ≃ R 21 m-1,-1 2 Q 11 R |v m+1 (L)| 2 ≃ R 22 m+1,-1 2 Q 22 L (3.78)
as well as,

u m+1 (R)v * m-1 (R) ≃ R 11 m+1,-1 Q 11 R (R 21 m-1,-1 ) * u m-1 (L)v * m+1 (L) ≃ R 12 m-1,-1 Q 22 L (R 22 m+1,-1 ) * . (3.79)
Taking into account Eq. (3.28), we can write

R m-1,-1 = 1 M 0 m-1 -Σ - m-1 0 Γ R g 12 R (m) Γ L g 21 L (m) 0 R m+1,-1 .
The change in the integrand near a MAR step can therefore be approximated by

δI (m) ω 0 → 2 ∆ + m ≃ -8Γ L Γ R ω √ ω 2 -∆ 2 Im g 22 R (m) R 22 m+1,-1 2 θ(ω -∆) + 8Γ 3 L Γ R ω √ ω 2 -∆ 2 Im g 22 L (m) g 21 L (m) (M 0 m-1 -Σ - m-1 ) 22 2 R 11 m+1,-1 2 θ(ω -∆) -8Γ 2 L Γ R ω √ ω 2 -∆ 2 g 21 L (m) 2 Im 1 M 0 m-1 -Σ - m-1 22 R 11 m+1,-1 2 θ(ω -∆) + 8Γ L Γ 2 R ω √ ω 2 -∆ 2 g 21 R (m) 2 Im 1 M 0 m-1 -Σ - m-1 11 R 22 m+1,-1 2 θ(ω -∆). (3.80)
It is clear that the first two terms in the above expression contribute to the current when Im g(m) ̸ = 0. As explained previously, this implies the condition

ω 0 > 2 ∆ m .
For the last two terms, in the vicinity of a MAR step ω 0 → 2 ∆ + m we can make the approximation

Im 1 M 0 m-1 -Σ - m-1 11,22 ≃ Im 1 M 0 m-1 + iη 11,22 ≃ Im 1 ω -(m -1)ω 0 -Γ R,L g 11,22 R,L (m) + iη . (3.81)
Although it was not evident from the start, the above expression means that the last two terms of (3.80) are also nonzero when the condition Im g(m) ̸ = 0 holds. The final step is to re-express the quantity R 11,22 m+1,-1 2 appearing in (3.80) through the locator expansion. We find:

R 22 m+1,-1 2 = (Γ L Γ R ) |m| 2 -1 R 22 -1,-1 2 × g 12 R (-2)g 21 L (-4) . . . g 12 R (m + 2)g 21 L (m) (M 0 -3 -Σ - -3 ) 11 . . . (M 0 m+3 -Σ - m+3 ) 11 (M 0 m+1 -Σ - m+1 ) 22 2 .
(3.82)

An analogous expression can be found for the R 11 Finally, in order to make the connection between the current and the spectral function, we can show that the absolute value of the diagonal elements appearing in (3.82) is connected to the spectral function. The proof goes as follows:

The imaginary part of the resolvent can be written as the difference between the resolvent and its hermitian conjugate, 2i

Im [R(ω)] = R(ω) -R † (ω). If we write the resolvent as R(ω) = [M -Σ] -1 then it follows that (R) -1 -(R † ) -1 = 2i Im(M -Σ). (3.83)
Multiplying the above expression with R from the left and with R † from the right we get the general relation

R † -R = 2iR Im(M -Σ)R † ⇒ -Im R = R Im(M -Σ)R † . ( 3

.84)

We define the matrix γ(ω) ≡ 2 Im(M 0 0 -Σ + 0 -Σ - 0 ). Now we can use the fact that the resolvent and the self-energy are diagonal in Nambu space to get the relation between the spectral function defined in Eq. (3.47) and the absolute value of the resolvent elements

A(ω) = - 2 π Im R 11 00 (ω) = γ 11 (ω) π R 11 00 (ω) 2 .
(3.85) Figure 3.7: Three terminal junction biased in the quartet configuration. When the distance of the dots is comparable to the superconducting coherence length, R ∼ ξ 0 , the system behaves like a molecule.

Equivalently, we can write

A(-ω) = - 2 π Im R 22 00 (ω) = γ 22 (ω) π R 22 00 (ω) 2 . (3.86)
Moreover, the translation property of the resolvent means that

R 11 ±1±1 (ω) 2 = π γ 11 (ω ± ω 0 ) A(ω ± ω 0 ), R 22 ±1±1 (ω) 2 = π γ 22 (ω ± ω 0 ) A(-ω ± ω 0 ). (3.87)
To summarize, we can replace absolute values R 11,22 m+1,-1 2 appearing in (3.80) by a spectral function times a function containing a string of Andreev reflections, as expressed in (3.82). The result is similar to Eq. 29 (for an S-S contact) in [START_REF] Cuevas | Hamiltonian approach to the transport properties of superconducting quantum point contacts[END_REF] or Eq. 3 (for a resonant level) in [START_REF] Levy Yeyati | Resonant tunneling through a small quantum dot coupled to superconducting leads[END_REF], with the major difference that, in order to avoid divergences we had to use a renormalized perturbation theory instead of a naive one. However, our derivation is not completely general, as it relied on the simpler structure of the resolvent in the V L = -V R case.

driven andreev molecule

We will now turn our attention to the Andreev molecule, which we have studied in the absence of voltage bias in the previous chapter, see Sec. 2.4.2. We consider a three-terminal junction setup as sketched in Fig. 3.7, where two quantum dots are each connected to a superconducting lead S a , S b , and both are coupled to the superconducting lead S c . The leads are voltage-biased in the quartet configuration with V a = -V, V b = +V and the middle lead is grounded, V c = 0. We take the first dot position to be x 1 = 0, and the second to be at some distance x 2 = R. Similar to the equilibrium case, we will show that the system behaves like a molecule when the interdot distance is comparable to the superconducting coherence length, R ≲ ξ 0 . This will cause a splitting of the Floquet-Andreev resonances and a modification of the dc current in the vicinity of the MAR steps.

FLS equations

The homogeneous part of Eq. (3.18) is 7 (suppressing the source index l):

(E + mω 0 ) u m (1) v m (1) = Γ c1 g c (m) + Γ a g 11 a (m-1) 0 0 g 22 a (m+1) u m (1) v m (1) + Γ a 0 g 12 a (m-1) g 21 a (m+1) 0 u m+2 (1) v m-2 (1) + Γ c g c (m, R) u m (2) v m (2) (E + mω 0 ) u m (2) v m (2) = Γ c2 g c (m) + Γ b g 11 b (m+1) 0 0 g 22 b (m-1) u m (2) v m (2) + Γ b 0 g 12 b (m+1) g 21 b (m-1) 0 u m-2 (2) v m+2 (2) + Γ c g c (m, R) u m (1) v m (1) . (3.88)
The new element with respect to the single dot case is a non-local Green's function which couples the the amplitudes on different dots g j,ii ′ (m) ≡ g j (m, R).

It depends explicitly on the distance between the two dots:

g j (m, R) = e iq(m)R cos(k F R)g j (m) + sin(k F R)σ z , (3.89)
where k F is the Fermi wave-vector. The reservoirs have been chosen as 1D for simplicity, but this should not influence the qualitative description. In equilibrium, a more careful analysis of the effect of dimensionality [START_REF] Kornich | Overlapping Andreev states in semiconducting nanowires: competition of one-dimensional and three-dimensional propagation[END_REF] and comparison with [START_REF] Pillet | Nonlocal Josephson effect in Andreev molecules[END_REF] shows that if the propagation in the middle reservoir is 3D instead of 1D, then the resulting hybridization of the Andreev molecule will be smaller at a given distance by an order of magnitude [START_REF] Benzoni | Hybridization of Andreev bound states in closely spaced Josephson junctions[END_REF]. Moreover, it is important that the middle superconductor S c is considered large on the mesoscopic scale, in the sense that its spectrum is continuous and its chemical potential is well-defined, but that the distance R between the two dots remains finite.

The non-local Green's function mediates the coupling between the junctions and oscillates on two very different length scales. For energies smaller than the gap |E| < ∆, the factor e iq(E)R = e - √ 1-(E/∆) 2 R/ξ 0 decays exponentially over distances larger than the superconducting coherence length ξ 0 ≡ v F /∆, while for energies above the gap, e iq(E)R oscillates without decay as long as there is no mechanism of decoherence in S c . A finite quasiparticle lifetime [START_REF] Kaplan | Quasiparticle and phonon lifetimes in superconductors[END_REF] will eventually produce decay of the quasiparticle propagation in S c over a mesoscopic coherence length that should be between two to three orders of magnitude larger than ξ 0 [START_REF] Mélin | Ultralong-distance quantum correlations in three-terminal Josephson junctions[END_REF]. These non-vanishing oscillations physically 7 The couplings to the superconducting leads have been defined as follows:

πρ 0 J 2 a1 ≡ Γ a , πρ 0 J 2 c1 ≡ Γ c1 , πρ 0 J 2 b2 ≡ Γ b , πρ 0 J 2 c2 ≡ Γ c2 ,
and

πρ 0 J c1 J c2 = √ Γ c1 Γ c2 ≡ Γ c
represent quasi-particle propagation in the continuum of the reservoirs, which is therefore not bound by the superconducting coherence length. On the other hand, the phase k F R oscillates rapidly at the scale of the Fermi wavelength λ F = 2π/k F which is typically much smaller than the superconducting coherence length, ξ 0 ≃ 10 3 λ F . The length-scale of the quasiparticle propagation depends on the quasi-particle energy E, while the phase k F R will produce geometric oscillations due to the trigonometric functions of Eq. (3.89). Since the two scales are very different, and we want to focus on new physics related to the energy-dependence rather than to any geometric effects, we will assume that the phase k F R is fixed. This choice should correspond to performing an average over the angle ψ = k F R, as is done in [START_REF] Mélin | Ultralong-distance quantum correlations in three-terminal Josephson junctions[END_REF].

In the basis of the Nambu spinor

Ψ m = (u m (1), v m (1), u m (2), v m (2)) T
we can rewrite Eq. (3.88) as

(LΨ) m ≡ M 0 m Ψ m -M + m+1 Ψ m+2 -M - m-1 Ψ m-2 . (3.90)
Explicitly, the 4 × 4 matrix

M 0 m = M (1) m -Γ c g c (m, R) -Γ c g c (m, R) M (2) m , (3.91)
contains the non-local coupling of the two dots in its antidiagonal blocks, as well as local terms on each of the dots, collected in the block-diagonal in the M

(1,2) m matrices:

M (1) m = (E + mω 0 ) 1 2 -Γ c1 g c (m) -Γ a g 11 a (m -1) 0 0 g 22 a (m + 1) , M (2)
m = (E + mω 0 ) 1 2 -Γ c2 g c (m) -Γ b g 11 b (m + 1) 0 0 g 22 b (m -1) . (3.92)
Finally, the matrices M ± m contain only local Andreev reflection terms

M - m =             0 Γ a g 12 a (m) 0 0 0 0 0 0 0 0 0 0 0 0 Γ b g 21 b (m) 0             , M + m =             0 0 0 0 Γ a g 21 a (m) 0 0 0 0 0 0 Γ b g 12 b (m) 0 0 0 0             . ( 3

Floquet spectrum

The resolvent can be calculated using the iterative formulas of Eq. (3.28). For the spectrum, we are interested in the diagonal elements R 00 . For example, if we perform local tunneling spectroscopy on the first dot, we will probe the imaginary part of the elements R 11,22 00 of the resolvent. We therefore define the spectral function of the first dot as

A 1 = - 2 π R 11 00 . Since R 00 = M 0 0 -Σ + 0 -Σ - 0 -1
, the influence of the second dot on the spectrum of the first is contained in the off-diagonal blocks both of M 0 0 and of the self-energy matrices. For the subgap part of the spectrum, |E| < ∆, the dominant interdot coupling term is the off-diagonal block of the matrix M 0 0 , i.e. a term which depends on the non-local Green's function, Γ c g c (0, R). When the interdot distance is comparable to the superconducting length R ∼ ξ 0 , one therefore expects that the two dots are hybridized, forming a non-equilibrium molecular state. The hybridization produces a splitting of the quasienergies, with four FWS ladders instead of two. This is illustrated in Fig. 3.8 which shows the maxima of the spectral function A 1 using inverse voltage scaling. The distance is set to R = ξ 0 /2.

Since the coupling element g c (0, R) decays exponentially when the distance between the two dots becomes much larger than ξ 0 , we expect that the level splitting will also show an exponential decrease, and that at large distances we will find the quasienergies corresponding to a single junction biased at (-V, 0) (i.e. Fig. 3.3). The lift of degeneracy as a function of the distance is shown in Fig. 3.9. In equilibrium, the level splitting of the ABS at a fixed interdot distance can be controlled by adjusting the phase difference through one of the junctions. In the driven case, the voltage can play an analogous role.

In reality, there is another regime when the distance is larger than the superconducting length but smaller than a mesoscopic coherence length ξ 0 < R < l ϕ . In this regime the system behaves like an interferometer, and there are oscillations of the spectral function around the single-junction value. We discuss this long-range interference effect in the next chapter and show that the oscillations of the spectral functions lead to an oscillatory I-V curve.

non-resonant dots In Eq. (3.88) we have set the dot energies at resonance, with ε 1,2 = 0. The effect of detuning can easily be investigated by adding an additional term to the matrix M 0 m :

M 0 m → M 0 m - ε 1 σ z 0 0 ε 2 σ z .
In Fig. 3.10 we show the effect of detuning on the spectral function of the first dot. The energy of the first dot is set at resonance ε 1 = 0, while the energy of the second dot varies, ε 2 = δ. We observe a progressive decoupling of the dots as the detuning δ is increased. The results suggest that the hybridization is more robust under detuning at larger voltage values.

MAR current

In order to calculate the MAR current we use Eq. (3.64). There, we introduced the matrix Q l,ii ′ acting in the dot ⊗ Nambu space. In order to proceed, we keep only the "local" part of Q, i.e. elements Q l,ii = 2iΓ li g r l (ω)g a l (ω) θ(ω). The reason that this approximation is reasonable is that the resolvent elements already contain the coupling between the dots to an arbitrary degree, since they are calculated by resumming the contribution of different paths. The current from the first dot to the middle superconductor is then:

I dc 1→c =4 ∑ m ∞ 0 dω 2π Im Γ c1 g 21 c (m), g 22 c (m) R 11 m,1 R 12 m,-1 R 21 m,1 R 22 m,-1 Q a R 21 m,1 R 22 m,-1 * -Γ a g 21 a (m + 1), g 22 a (m + 1) R 11 m+2,0 R 12 m+2,0 R 21 m,0 R 22 m,0 Q c R 21 m,0 R 22 m,0 * . (3.94)
The main result for the bijunction current is presented in Fig. 3.11, where we plot the results of numerical calculations of Eq. (3.94). For simplicity, we consider all tunnel couplings to be equal Γ = 0.3 ∆, and we fix k F R = π/4, in order to avoid oscillations on the scale of the Fermi wavelength. We see that for large distances between the dots (grey line) the I-V curve approaches that of a single S-QD-S junction biased at (-V, 0) (dashed black line). The single junction biased at (-V, 0) has a voltage difference which is halved compared to the (-V, +V) configuration which we studied in 3.4.3. Its MAR steps appear when mω 0 2 > ∆ ± ϵ ⇒ ω 0 > 2(∆±ϵ) m . The behavior of the current at the regime R ≫ ξ 0 will be a subject for the next chapter.

For distances comparable to the superconducting length (red and orange line), we are clearly in a "molecular junction" regime [START_REF] Toroker | On the relation between steady-state currents and resonance states in molecular junctions[END_REF], and the MAR steps break into four sub-steps. These sub-steps correspond to the splitting of the energy levels in the Floquet spectrum in the Andreev molecule regime, as discussed earlier (see Fig. 3.8). The steps are visible when the resonances are not overlapping, that is when their widths are smaller than their separation. Given that the width of a resonance coupled to a continuum of states increases when the coupling to the continuum is increased, one expects that small values of voltage bias and tunnel couplings give sharper features. Indeed, the new features due to the proximity of a second junction are more clear at voltages equal to higher-order subdivisions of the gap. However, the modification should still be visible around the 2 ∆ 3 or the 2 ∆ 5 MAR steps. Moreover, the influence of the tunnel couplings on the I-V curves is shown in Fig. 3.12. As it is expected we observe that the subgap features are softened when increasing the tunnel coupling to the reservoirs. contribution of floquet harmonics. As discussed in the case of the S-QD-S junction, depending on the region of the I-V curve the sum over the Floquet modes in Eq. (3.94) of the current can be drastically truncated at large voltage drives. This localization is illustrated in Fig. 3.13. At large voltages ω 0 > 2 ∆, the drive is strong enough to promote quasi-particles directly above the gap without any MAR processes, and we only need to sum over two harmonics m = ±1. As we lower the voltage, we progressively need to add more harmonics, in correspondence to the MAR processes which are dominant. In the region of the first allowed MAR process, 2 ∆ 3 < ω 0 < 2 ∆, the current is well approximated by summing over three harmonics m = ±1, -3, in the next region of 2 ∆ 5 < ω 0 < 2 ∆ 3 , we need to add one more m = ±1, -3, -5, and so on.

conclusions

In conclusion, we have studied the S-QD-S system and the Andreev molecule subjected to a dc voltage drive, which brings the systems out of equilibrium. The superconductivity of the leads allows to explore Floquet physics by simple voltage biasing. As is often done in Floquet systems [START_REF] Oka | Floquet engineering of quantum materials[END_REF], one can expand quantities, such as the electron/hole amplitudes on the dot(s), into Fourier modes and map the initial time-dependent BdG equation into a tight-binding chain with sites labeled by the Fourier harmonics. The sites are coupled to their nearest neighbors through absorption or emission of virtual photons which here correspond to Andreev reflections. The tight-binding analogy can be exploited to get an iterative solution that corresponds to a Dyson equation for the Green's function on the dot(s). We find that the initial discrete levels are renormalized by a self-energy which corresponds to the sum of two independent processes, the Σ ± of Eq. (3.25). These matrices describe processes that connect the initial discrete levels to the superconducting continua either above the gap, ω > ∆, or below the gap, ω < -∆, and therefore are a source of dissipation. An effective Hamiltonian for the amplitudes on the dot(s) will therefore be non-Hermitian, since the S-QD-S system in the presence of voltage is an open quantum system. In the case of the Andreev molecule, the resonances on the dots are then coupled through the middle reservoir by the nonlocal Green's function of Eq. (3.89). This coupling produces splitting of the energy levels and avoided crossings in the spectrum at distances between the dots comparable to the superconduct-ing coherence length R ∼ ξ 0 . It also modifies the subgap structure, produc- ing splitting of the MAR steps. This direct coupling decays exponentially with R/ξ 0 . Instead, at R ≫ ξ 0 the two dots are coupled through higher-order processes contained in the self-energy which become the dominant mechanism for coupling at long distances. This long-range mechanism involves local MAR processes on each dot and subsequent quasiparticle propagation through the middle superconductor at energies |ω| > ∆. The system in this regime behaves like an interferometer, in the sense that the current becomes an oscillatory function of the voltage. We will explore this effect in the following chapter.

F L O Q U E T -T O M A S C H E F F E C T

introduction

In this chapter, we will study the S-QD-S-QD-S bijunction in the limit where the middle superconductor is long, making the interdot distance R larger than the superconducting coherence length. This means that the molecular effects we saw in the previous chapter are no longer relevant, as at large interdot distances there is no hybridization between the resonances on the dots. However, at a mesoscopic regime of distances we will show that the system behaves like an interferometer, and that a long-range coupling is induced between the dots.

In particular, we will show that the subgap current oscillates as a function of the voltage and the interdot distance. This interference is the "Floquet version" of the geometrical interference effect first discovered by Tomasch in thick superconducting films [START_REF] Tomasch | Geometrical resonance in the tunneling characteristics of superconducting Pb[END_REF][START_REF] Tomasch | Geometrical resonance and boundary effects in tunneling from superconducting In[END_REF][START_REF] Tomasch | Energy spacing of geometrical resonance structure in very thick films of superconducting In[END_REF]. The Tomasch effect can be explained as an interference between electron-like and hole-like excitations which are degenerate in energy, but differ in their wavenumbers [START_REF] Wolfram | Tomasch oscillations in the density of states of superconducting films[END_REF]

k ± ≃ k F ± E 2 -∆ 2 /hv F .
McMillan and Anderson [START_REF] Mcmillan | Theory of geometrical resonances in the tunneling characteristics of thick films of superconductors[END_REF] pointed out that such a situation can arise when there is a variation of the gap parameter, such as at an N-S interface. In particular, Andreev scattering on the interface mixes incident excitations with k ± and reflected excitations with k ∓ which results to a quasiparticle interference effect. Consequently, electron tunneling measurements on some part of the superconductor at some distance d from the N-S interface and at energies larger than the gap show that the density of states and the tunneling current are a periodic function of

(k + -k -)d = 2d hv F (eV) 2 -∆ 2 ,
where V is the applied voltage due to the tunnel probe. A typical distance d in the Tomasch experiments was a few tens of micrometers, corresponding to a distance two orders of magnitude larger than a typical superconducting coherence length. The Tomasch effect is therefore a nonlocal effect where the local density of states measured on one end of the superconductor depends on the properties of the interface at the other end.

In the case of the voltage-biased S-QD-S-QD-S system, there is a Floquet analogue of the Tomasch effect due to processes that involve local MARs on each dot, followed by quasiparticle propagation in the middle superconductor at energies above the gap, |E| > ∆. The existence of this phenomenon was proposed recently in [START_REF] Mélin | Ultralong-distance quantum correlations in three-terminal Josephson junctions[END_REF], showing the effects on the quartet current, which is the dissipationless current carried by pairs of Cooper pairs [START_REF] Huang | Evidence for 4e charge of Cooper quartets in a biased multiterminal graphene-based Josephson junction[END_REF]. This chapter builds up on Ref. [START_REF] Mélin | Ultralong-distance quantum correlations in three-terminal Josephson junctions[END_REF]. In particular, we focus on the consequences of the Floquet-Tomasch effect firstly on the dc current-voltage characteristics, and secondly on the spectrum of the bijunction in the subgap region, |E| < ∆. In both cases the Floquet-Tomasch effect results in the appearance of oscillations of the curves with respect to the single junction curves. The corresponding pole structure of the resolvent is also drastically modified with respect to the resolvent of the single junction, and the number of poles found increases with the dot separation.

In order to understand the physics related to coupling two dots via a continuum, we study a simplified model of two discrete levels coupled through a normal wire with a linearized dispersion relation. We show the connection between the analytic properties of the resolvent and the temporal dynamics, in particular that the appearance of a great number of poles in the complex energy plane is translated into oscillations in the time domain. Finally, a relation between the resolvent and the scattering matrix is derived, which makes the connection between scattering and decay properties. This relation is known in the literature as the Mahaux-Weidenmüller formula [START_REF] Dittes | The decay of quantum systems with a small number of open channels[END_REF].

Inspired by the study of this toy model, we show that the modification of the bijunction's resolvent, as well as the resulting oscillations in the spectrum, can be described by an effective non-Hermitian two-level model of resonances coupled through a continuum.

As in the previous chapter, we study the bijunction biased on the "quartet line" [START_REF] Freyn | Production of nonlocal quartets and phase-sensitive entanglement in a superconducting beam splitter[END_REF], that is we consider the voltage configuration (V a , V c , V b ) = (-V, 0, +V). This means that the phases on each reservoir evolve with time as per the Josephson relation ϕ a (t) = ϕ a -2eVt, and ϕ b (t) = ϕ b + 2eVt (having chosen ϕ c ≡ 0 as a reference phase). As a result, the quartet phase, defined as the combination ϕ q ≡ ϕ a (t) + ϕ b (t) -2ϕ c = ϕ a + ϕ b , is static. We find that the Floquet-Tomasch oscillations of the spectral function depend on the quartet phase. Specifically, we find that the amplitude of oscillations could be potentially controlled non-locally by changing the phase drop on one of the dots. This chapter is mostly based on [START_REF] Keliri | Long-range coupling between superconducting quantum dots induced by periodic driving[END_REF].

MAR current in the interferometer regime

We start by considering the numerical results for the MAR current. We find that at interdot distances larger than the superconducting coherence length, the dc current oscillates as a function of the applied voltage. The oscillations are more pronounced in the vicinity of the MAR steps. The effect of an increasing distance between the dots is shown in Fig. 4.1, where we plot the current from the first dot to the middle reservoir according to Eq. (3.94). The plot is focused in the region of the MAR step which occurs around ω 0 = 2∆ 5 . At R > ξ 0 the curves approach the I-V curve of a single junction (dashed line). However, we observe that a large distance (grey lines at R = 5ξ 0 , 10ξ 0 , 25ξ 0 ) produces oscillations of the I-V curve itself around the single junction curve. We also see that the frequency of oscillations increases with the distance.

Density of states above the gap

Figure 4.2 shows the spectral function of the first dot at energies above the gap, calculated at a fixed voltage value ω 0 = 0.45 ∆. As in the case of the MAR current, we observe that at large distances the spectral function oscillates around the spectral function of the single junction. We will show that the spectral function is a periodic function of a Tomasch phase factor 2R

v F √ E 2 -∆ 2 .
It has been argued [START_REF] Mélin | Ultralong-distance quantum correlations in three-terminal Josephson junctions[END_REF] that the Floquet-Tomasch effect could be used to create correlations of Cooper pairs over long distances, which are orders of magnitude larger than the superconducting coherence length (in the Tomasch experiment the thickness was some tens of micrometers 10 -30µm ∼ 100ξ 0 .). It is a non-local effect over distances which are not achievable in the absence of the voltage drive since it is mediated by quasi-particles which, by MAR processes, reach the continuum of states of the middle superconductor S c where they can propagate over a long distance without being bound by the superconducting coherence length.

Using the formalism developed in the previous chapter, we will show that the spectral function of the first dot A(E) = -2 π Im R e 1 e 1 00 (E) oscillates as a function of the energy E and the interdot distance R at energies above the gap, as shown in Fig. 4.2. We will slightly change notation with respect to the previous chapter: we use indices i = 1, 2 to refer to the first and second dot, respectively. Within the subspace of the dot i we use indices e i , h i to refer to the Nambu indices. The resolvent can be calculated using the expression

R mm = M 0 m -Σ m -1 ,
where M 0 m is given by Eq. (3.91) and Σ m = Σ + m + Σ - m is the total self-energy matrix. The self-energy is a 4 × 4 matrix which can be written in block form as

Σ m ≡ Σ (1) m Σ 12 m Σ 21 m Σ (2) m . (4.1)
The resolvent of the bijunction is also written in block form

R mm = M (1) m -Σ (1) m -Γ c g c (m, R) -Σ 12 m -Γ c g c (m, R) -Σ 21 m M (2) m -Σ (2) m -1 . (4.2)
As seen in Fig. 4.2, the oscillations of Im R e 1 e 1 00 are superimposed on the singlejunction curve. We therefore define the resolvents of the single junctions (where there is no coupling between the dots, or one dot is considered at a distance R → ∞ from the other) as

R (1,2) mm = M (1,2) m -Σ (1,2) m -1 (4.3)
where all matrices are now calculated in the 2 × 2 Nambu subspace of dot 1 or dot 2. The self-energies Σ

(1,2) m represent local MAR processes which connect the states of dot 1 or 2 above the gap to those below the gap.

For energies above the gap E > ∆, the element which gives access to the spectral function, R 00 , can be approximated by

R 00 ≈ M (1) 0 -Σ (1) 0 -Γ c g c (0, R) -Γ c g c (0, R) M (2) 0 -Σ (2) 0 -1 . (4.4)
The self-energy block matrices Σ 12 , Σ 21 in the off-diagonal can be neglected, because they are of higher-order in the tunnel couplings compared to the direct coupling term Γ c g c (0, R). Equation (4.4) describes resonances on each dot which are formed due to local MAR processes and which are then coupled by propagation in the middle reservoir, represented by the non-local Green's function g c (0, R) in the off-diagonal. Since the resolvent is a block matrix, we can invert it blockwise. For the resolvent of Eq. (4.4), the upper-left block (corresponding to the subspace of the first dot) can be approximated by

1 [R 00 ] dot 1 ≃ R (1) 00 + Γ 2 c R (1) 00 g c (0, R)R (2) 00 g c (0, R)R (1) 00 + O Γ 4 c . (4.5)
This is a Fabry-Pérot-like process, with the first term in the series corresponding to a single junction resonance and each successive term representing an even number of reflections between the two dots. The first correction to the resolvent of the two coupled dots with respect to the resolvent corresponding to an uncoupled dot then is

[R 00 ] dot 1 -R (1) 00 ≈ Γ 2 c R (1) 00 g c (0, R)R (2) 00 g c (0, R)R (1) 00 . (4.6)
A term like the above has a clear physical interpretation: the resonances on each dot represented by R (1,2) are coupled via propagating quasi-particles in the continuum of states of the middle superconductor. The amplitude of the effect will therefore depend on the specific geometry and the amount of disorder of the middle reservoir (here we have considered a clean, one-dimensional superconducting wire). We see that a Tomasch phase factor will appear due to the non-local Green's function g c (0, R) which is proportional to the phase

e i √ E 2 -∆ 2 R/v F . A total phase e 2i √ E 2 -∆ 2 R/v F
is then accumulated by quasi-particles which travel from dot 1 to dot 2 and back. As a result, at a fixed distance R, the resolvent of the coupled system oscillates as a function of the energy around the single-dot resolvent R (1) , as shown in Fig. 4.2.

1 If M is the inverse of a block matrix M =     A B C D     -1 =     M 11 M 12 M 21 M 22    
then its upper-left block is equal to

M 11 = (A -BD -1 C) -1 = A -1 + A -1 B 1 D -CA -1 B CA -1 .

floquet-tomasch effect on the spectrum

In the following section we will show that the periodic driving affects the subgap region of the spectrum, and induces a long-range interdot coupling at distances R ≫ ξ 0 . The main difference with respect to the previous discussion is that for subgap energies the approximation performed in Eq. (4.4) is not valid. The reason is that for |E| < ∆ the non-diagonal self-energy blocks Σ 12 , Σ 21 cannot be ignored. In fact they are the dominant terms responsible for the interdot coupling. At equilibrium (in the absence of voltage bias), there are two competing mechanisms for the coupling of the two dots in the molecular regime: a) crossed Andreev reflection (CAR) processes, involving the Andreev reflection of two electrons, one from each dot, which then form a Cooper pair in the middle superconductor, and b) elastic cotunneling (EC) processes, involving normal transmission of quasiparticles through the middle superconductor [START_REF] Deutscher | Crossed Andreev reflections[END_REF][START_REF] Russo | Experimental observation of bias-dependent nonlocal Andreev reflection[END_REF]. In terms of the superconducting Green's functions, CAR corresponds to the anomalous propagators, while EC corresponds to the normal components.

An efficient way to tune the rate between these two processes has recently been proposed and demonstrated [START_REF] Liu | Tunable superconducting coupling of quantum dots via Andreev bound states in semiconductorsuperconductor nanowires[END_REF][START_REF] Dvir | Realization of a minimal Kitaev chain in coupled quantum dots[END_REF]. At equilibrium, separating the two dots at distances larger than ξ 0 results in trivially recovering the spectrum of two single dots, as both CAR and EC will be exponentially suppressed. Driving the system, however, allows for a new mechanism of interdot coupling, which we will now analyze.

We will examine the Floquet spectrum at large interdot distances, and compare it to the single junction spectrum. As explained in the previous chapter, once we know the "Floquet chain operator"

L mm = M 0 m -Σ + m -Σ - m (4.7)
the spectral function is obtained by the imaginary part of the inverse operator R = L -1 . Equivalently, the solutions of the characteristic polynomial det L 00 = 0 are the (complex) quasienergies. We can therefore plot the function det L 00 in the complex plane, and look for its zeroes. If a zero of det L 00 appears at an energy E 0iγ, then the spectral function will have a peak around the energy E 0 , with a width given by the imaginary part γ.

single junction. If the two dots are decoupled, then there is a single junction with V a = -V, V c = 0. We denote its Floquet chain operator by

L (1) mm = M (1) m -Σ +,(1) m -Σ -, (
Γ a = Γ c1 = ∆/2.
and its inverse by the resolvent operator R (1) = L (1) -1

, where

M

(1)

m = E + mω 0 -Γ c1 g c (m) -Γ a g ee a (m -1) 0 0 g hh a (m + 1) , ( 4.8)
and the self-energy matrices are:

Σ +,(1) m = Γ 2 a    0 0 0 g he a (m + 1) 1 M (1) m+2 -Σ +,(1) m+2 ee g eh a (m + 1)    , (4.9a) Σ -,(1) m = Γ 2 a     g eh a (m -1) 1 M (1) m-2 -Σ -,(1) m-2 hh g he a (m -1) 0 0 0     . (4.9b)
The self-energy adds dissipation to the diagonal elements of M (1) due to MAR processes. To simplify things, if we consider a voltage such that ∆/2 < ω 0 < ∆, we see that if an initial ABS level is close to an energy |E| ≪ ∆, then two Andreev reflections would be enough to connect the level to the continuum of states above or below the gap by absorption or emission of virtual photons of energy 2ω 0 . In this case the self-energies will have an imaginary part proportional to g eh a (±1)g bijunction. The numerical results for the spectral function on dot 1 of the bijunction are presented in Fig. 4.4. At large interdot distances, the real part of the resonances is not shifted with respect to the single-junction peaks, but oscillations appear, superimposed on the single-junction peaks due to coupling with the second dot. In the complex plane, the resulting behavior is a proliferation of the solutions of the characteristic polynomial det L 00 = 0. The frequency of oscillations of the resolvent and, correspondingly, the number of zeroes in the complex plane increase with the distance. This behavior of the zeroes of det L 00 is shown in Fig.

ω 0 = ∆/2.
distance.

The remaining of this chapter is an effort to understand the pole structure shown in Fig. 4.4 and Fig. 4.5. Making certain approximations, we can study the system analytically, and derive an effective non-Hermitian Hamiltonian of a two-level system.

Large voltage bias, large separation approximation

The Floquet chain operator involves continued fraction expressions of the selfenergy. This corresponds to resumming an infinite number of (certain kinds) of diagrams. Technically, one truncates the continued fractions at some cutoff index |N| > ∆/ω 0 by considering that the self-energies become small Σ ± ±N → 0 at large energies |E ± Nω 0 | ≫ ∆. Therefore, at voltages which are a significant fraction of the gap, one can greatly simplify the expressions of Σ ± , while at small voltages an increasingly greater number of Floquet harmonics need to be taken into account. We will here concentrate on the former regime, because it facilitates the analytical calculation of the L operator, while giving some insight on the involved mechanism of coupling. However, the Floquet-Tomasch mechanism of coupling that will be described in the next section occurs at smaller voltage values as well, albeit at higher MAR order and therefore at a higher order in the tunnel couplings.

We will therefore study the subgap spectrum of the bijunction in a regime of voltages which are a significant fraction of the gap ∆/2 < ω 0 < ∆. The assumption is that two reflections connect any subgap state to the continua |E ± 2ω 0 | > ∆, but one reflection is not necessarily sufficient if the energy is close to zero |E ± ω 0 | < ∆. Moreover, we will assume that the interdot distance is large, meaning that any molecular effects can be ignored.

The starting point for understanding the results of Fig. 4.4 and 4.5 is the Floquet chain operator at m = 0:

L 00 = M 0 0 -Σ + 0 -Σ - 0 = M (1) 0 -Γ c g c (0, R) -Γ c g c (0, R) M (2) 0 -Σ + 0 -Σ - 0 . (4.10)
The operator L is written in the basis of the four-component Nambu spinor

Ψ m = u m (1), v m (1), u m (2), v m (2)
T and is therefore a 4 × 4 matrix acting in dot ⊗ Nambu space. The diagonal blocks of L correspond to intradot processes, while the off-diagonal blocks correspond to interdot processes. Specifically, the dots 1 and 2 are each coupled by local reflections to their closest reservoirs. This information is contained in the block matrices:

M (1,2) m = (E+mω 0 )1 2 -Γ c,(1,2) g c (m) -Γ (a,b) g ee (a,b) (m + s (a,b) ) 0 0 g hh (a,b) (m -s (a,b) )
.

(4.11) The off-diagonal blocks of L 00 couple the two dots through processes involving nonlocal Andreev reflections. The off-diagonal coupling term in M 0 0 is the nonlocal Green's function of the middle reservoir g c (0, R), given by Eq. 3.89. The nonlocal Green's function is the dominant source of coupling at small distances R ≲ ξ 0 , but becomes exponentially small at large distances for processes inside the gap (i.e. for energies |E| < ∆). Therefore, in the regime of interest R ≫ ξ 0 , the coupling of the two dots will be contained entirely in the self-energy matrices Σ ± 0 :

L 00 ≃ M (1) 0 0 0 M (2) 0 -Σ + 0 -Σ - 0 , R ≫ ξ 0 (4.12)
The self-energy elements do not go to zero as e -R/ξ 0 , but are limited by a mesoscopic coherence length instead [START_REF] Mélin | Ultralong-distance quantum correlations in three-terminal Josephson junctions[END_REF]. For a large voltage bias ∆/2 < ω 0 < ∆, we can truncate the expressions for the self-energies such that only two Floquet harmonics are kept, Σ ± |m|≥2 → 0. The self-energies can be calculated Figure 4.6: Self-energy terms (a) Σ +,h 1 e 2 and (b) Σ -,e 1 h 2 couple a hole (electron) on dot 1 to an electron (hole) on dot 2 through propagation in the middle reservoir. An overall phase e ∓iϕ q e iq(±2)R is accumulated.

using the iterative definition of Eq. (3.28) together with Eq. (3.93). Explicitly, the self-energy matrices have the form:

Σ + 0 =        0 0 0 0 0 Γ 2 a g he a (1) 1 M 0 2 e 1 e 1 g eh a (1) Γ a Γ b g he a (1) 1 M 0 2 e 1 h 2 g he b (1) 0 0 Γ a Γ b g eh b (1) 1 M 0 2 h 2 e 1 g eh a (1) Γ 2 b g eh b (1) 1 M 0 2 h 2 h 2 g he b (1) 0 0 0 0 0        , ( 4.13a)
Σ - 0 =          Γ 2 a g eh a (-1) 1 M 0 -2 h 1 h 1 g he a (-1) 0 0 Γ a Γ b g eh a (-1) 1 M 0 -2 h 1 e 2 g eh b (-1) 0 0 0 0 0 0 0 0 Γ a Γ b g he a (-1) 1 M 0 -2 e 2 h 1 g he b (-1) 0 0 Γ 2 b g he b (-1) 1 M 0 -2 e 2 e 2
g eh b (-1)

         . (4.13b)
By inverting the 4 × 4 matrix M 0 ±2 , one can express the self-energies (and therefore the coupling between the dots) as a function of local and nonlocal Green's functions of the reservoirs. The inversion of the matrix M 0 m can be performed blockwise. If the matrix has the form:

M 0 m = M (1) m -Γ c g c (m, R) -Γ c g c (m, R) M (2) m , ( 4.14)
then we can decompose its inverse as (suppressing the indices m, R for brevity)

1 M 0 ≈
where we have made a perturbative expansion in the tunnel couplings and kept only terms up to O Γ 2 c . The coupling between the dots is contained entirely in the non-diagonal terms of the self-energy. These can be written as

Σ -,e 1 h 2 ≈Γ a Γ b Γ c g eh a (-1) 1 M (1) -2 hh g he c (-2, R) 1 M (2) -2 ee g eh b (-1), (4.16a) Σ +,h 1 e 2 ≈Γ a Γ b Γ c g he a (1) 1 M (1) 2 ee g eh c (2, R) 1 M (2) 2 hh g he b (1), (4.16b) Σ +,e 2 h 1 ≈Γ a Γ b Γ c g eh b (1) 1 M (2) 2 hh g he c (2, R) 1 M (1) 2 ee g eh a (1) (4.16c) Σ -,h 2 e 1 ≈Γ a Γ b Γ c g he b (-1) 1 M (2) -2 ee g eh c (-2, R) 1 M (1) -2 hh g he a (-1) (4.16d)
The above formulas can be interpreted as specific physical processes that couple the two dots through local and nonlocal Andreev reflections. All processes couple an electron (hole) at initial energy |E| ≪ ∆ on dot 1 to another hole (electron) at same energy E on dot 2. Initially the quasiparticle on dot 2 is Andreev reflected locally on reservoir S b whereby its energy is changed by E ± ω 0 . This is then followed by a nonlocal Andreev reflection through the middle superconductor S c at energies which are above the gap |E ± 2ω 0 | > ∆, so that the propagation is not limited by the superconducting coherence length. Finally, a local Andreev reflection on reservoir S a returns the quasiparticle to the initial energy E on dot 1. A graphical representation of Eq. (4.16) is sketched in Fig. 4.6. The coupling due to processes like the above involves three Andreev reflections, meaning it is of order O(Γ a Γ c Γ b ) in the tunnel couplings. We can also see that the three Andreev reflections will contribute a quartet phase factor e ±iϕ q , where ϕ q = ϕ a + ϕ b -2ϕ c with ϕ c = 0. Finally, an energy-dependent phase factor, which we could call the "Floquet-Tomasch phase factor",

e iq(±2)R = e ±i √ (E±2ω 0 ) 2 -∆ 2 R/v F , |E ± 2ω 0 | > ∆, (4.17)
is also accumulated due to the propagation in the middle superconductor.

interlude: coupling discrete levels through a continuum

The preceding analysis suggests that at the regime of interest the relevant physics is that of two resonances coupled through a continuum. A crucial element is that the self-energy processes (4.16) are energy-dependent. Even more importantly, the existence of oscillations as a function of the energy means that we cannot ignore the energy dependence of the self-energies.

Often when considering the coupling between a discrete level and a continuum, one can consider their interaction as constant. This results in a resonance of Lorentzian form and any state initially on the discrete level will eventually decay into the continuum as per the Fermi Golden Rule. It is useful to study this kind of problems with the use of the resolvent operator of the Hamiltonian R(z) ≡ 1 z -H . From R the time evolution operator is obtained through contour integration, giving access to the temporal dynamics of the system. Indeed, for positive times t > 0, the evolution operator is given by

U(t) = - 1 2πi ∞ -∞ dω e -iωt R(ω + iη), t > 0.
The time evolution therefore depends on the analytic properties of the resolvent.

In the case of a single discrete level coupled to a continuum of states, the part of the resolvent on the dot has the form

R(z) = 1 z -ε d -Σ(z)
.

The self-energy can be written as Σ(ω) = Σ ′ (ω) + iΣ ′′ (ω). The real part Σ ′ usually adds a slight shift to the energy of the unperturbed level. The most important qualitative change comes from the imaginary part Σ ′′ . The assumption leading to Fermi's Golden Rule, and which amounts to replacing the selfenergy by a complex valued constant, is that Σ(ω) is a slowly varying function, particularly in the vicinity of ϵ d . This assumption leads to the approximation

R(ω + iη) ≃ 1 ω -ε d + iΓ ,
in which case the resolvent has a simple pole at ω = ε d -iΓ. Taking the appropriate contour as in Fig. 4.7a, the time evolution operator is U(t) = e -i ε d t e -Γt , so that the probability to find a particle on the discrete level decays exponentially

P(t) = |U(t)| 2 = e -2Γt .
However, the exponential decay is, in many cases, modified. For example, in the previous discussion we have glossed over the fact that the continuum adds a branch cut in the complex plane. Then if the continuum has a finite width [START_REF] Kogan | Decay of discrete state resonantly coupled to a continuum of finite width[END_REF], a bandgap [START_REF] Lambropoulos | Fundamental quantum optics in structured reservoirs[END_REF], or a singular density of states [START_REF] Basko | Landau-Zener-Stueckelberg physics with a singular continuum of states[END_REF] (this is a non-exhaustive list) the exponential decay can be significantly modified, or even altogether suppressed [START_REF] Basko | Landau-Zener-Stueckelberg physics with a singular continuum of states[END_REF]. In this section we will take the wide-band limit which means we ignore any such effect. A more careful discussion can be found in [START_REF] Cohen-Tannoudji | Atom-photon interactions: basic processes and applications[END_REF] or [START_REF] Goldberger | Collision theory[END_REF].

We will instead focus on the interference effects due to the existence of two discrete levels, which are coupled indirectly through propagation in a continuum of states. We will first study the spectrum of the system and the time-evolution of the populations of the levels. We will show that in a strongcoupling regime (whose meaning will be given later) the resolvent has a great number of poles in the complex plane. In the time-domain the existence of the poles is translated into an oscillation of the populations: a state initially prepared on one level will travel through the wire and be reflected back a few times before gradually evaporating. We will then study the out of equilibrium situation where the wire is connected to reservoirs with different chemical potentials on either end. We derive a relation between the scattering matrix and the resolvent. In this case, the existence of multiple poles in the complex plane is translated into Fabry-Pérot oscillations of the transmission coefficient.

Spectrum and time-evolution

We will study a toy model of two discrete states which are indirectly coupled through a one-dimensional wire. In first quantization, we associate discrete states |d i ⟩ with quantum dots at positions x i , separated by a distance x 2x 1 ≡ L, and continuous states |k⟩ with the reservoirs. The discrete states and the continuous states live in different subspaces and are therefore orthogonal. The continuous states act as channels of decay for the bound states initially on the dots, and upon coupling dots and reservoirs, the bound states generally turn into resonances with a finite lifetime and eventually decay into the environment. The total Hamiltonian of this system is the Fano-Anderson Hamiltonian 18) and the wavefunction of the system can be written as a linear combination of states on the dots and on the wire

H = ∑ k ω k |k⟩⟨k| + ∑ i∈dots ϵ i |d i ⟩⟨d i | + ∑ k,i J k e -ikx i |k⟩⟨d i | + J * k e ikx i |d i ⟩⟨k| , (4.
|Ψ(t)⟩ = ∑ i∈dots g i (t) |d i ⟩ + ∑ k f (k, t) |k⟩ . (4.19)
We assume that initially, t = 0, there is only one electron in the system, localized on the first dot, so that we have the initial conditions: g 1 (0) = 1, g 2 (0) = 0, f (k, 0) = 0. Then, the Schrödinger equation i d dt |Ψ(t)⟩ = H |Ψ(t)⟩ gives a system of coupled equations for the amplitudes

i dg i (t) dt = ϵ i g i (t) + ∑ k J * k e ikx i f (k, t) i d f (k, t) dt = ω k f (k, t) + ∑ i J k e -ikx i g i (t). (4.20)
In order to resolve this system, it is convenient to use the Fourier transform

g(ω) = +∞ 0 dt e iωt g(t). (4.21)
For positive times t ≥ 0, the integral is convergent if Im ω > 0. Using

i +∞ 0 e iωt dg(t) dt dt = ωg(ω) -ig(0), (4.22)
Eq. (4.20) is transformed into

ωg 1 (ω) -i = ϵ 1 g 1 (ω) + ∑ k J * k e ikx 1 f (k, ω) ωg 2 (ω) = ϵ 2 g 2 (ω) + ∑ k J * k e ikx 2 f (k, ω) ω f (k, ω) = ω k f (k, ω) + J k e -ikx 1 g 1 (ω) + J k e -ikx 2 g 2 (ω) (4.23)
By eliminating the amplitudes of the reservoirs, we obtain the linear system:

ω -ϵ 1 -∑ k |J k | 2 ω -ω k g 1 (ω) -∑ k |J k | 2 ω -ω k e -ikL g 2 (ω) = i ω -ϵ 2 -∑ k |J k | 2 ω -ω k g 2 (ω) -∑ k |J k | 2 ω -ω k e ikL g 1 (ω) = 0. (4.24)
Formally, it is straightforward to compute the amplitudes by defining the resolvent operator as the inverse of the linear operator appearing in Eq. (4.24)

g 1 (ω) g 2 (ω) = R(ω) i 0 . ( 4.25)
The behaviour of the amplitudes as a function of time can be then calculated using the inverse Fourier transform of the amplitudes

g i (t) = 1 2π g i (ω)e -iωt dω , (4.26)
and the survival probability to find a particle on the dot i at some time t is given by the square of its amplitude

P i (t) = |g i (t)| 2 . ( 4.27)
Since t > 0, the integral of Eq. (4.26) is convergent if the contour is chosen to be a semi-circle of infinite radius in the lower half of the complex plane, Im ω < 0. However, the initial definition of g(ω) was for Im ω > 0. In practice this means one has to analytically continue the functions g i (ω), initially defined on the upper half complex plane, to the lower half complex plane. The amplitudes are then calculated by contour integration of the resolvent elements. The presence of multiple poles will modify the Fermi golden rule behavior of the survival probability. At long times, the main contribution to the integral will come from any poles near the real axis.

Linear dispersion

The solution (4.26) essentially reposes on knowledge of the analytical properties of the self-energy:

Σ(ω) ≡ ∑ k |J k | 2 ω -ω k 1 e -ikL e ikL 1 . (4.28)
In order to simplify things, the wire connecting the dots will be considered one-dimensional with a linear dispersion ω k = vk. For a normal metal the linearization (Tomonaga-Luttinger model) is justified for the description of the low lying excitations around the Fermi surface, which in 1d consists of two Fermi points at k = ±k F . The only complication is that there are two branches: a right-moving branch with energy ω k = vk and a left-moving branch with energy ω k = -vk. In the wide-band approximation, the coupling J k = J can be considered constant. Appropriately taking the contours then gives

Σ(ω + iη) = dk 2π |J k | 2 ω -vk + iη 1 e -ikL e ikL 1 = -iγ 1 e iωL/v e iωL/v 1 , (4.29)
where γ ≡ J 2 /v. Equation (4.26) becomes

ω -ϵ 1 + iγ iγe iωL/v iγe iωL/v ω -ϵ 2 + iγ g 1 (ω) g 2 (ω) = i 0 . ( 4.30)
The resolvent has two series of poles in the lower half of the complex plane. Assuming identical dots, ϵ 1 = ϵ 2 = ϵ 0 ≡ vk 0 , these are given by the solutions of the transcendental equation:

(ω -ϵ 0 + iγ) 2 + γ 2 e 2iωL/v = 0 ⇒ ω ± = ϵ 0 -iγ(1 ± e iωL/v ). (4.31)
One of the solutions ω + is found numerically and plotted in Fig. 4.8 for increasing interdot distances. The other solution ω -behaves in a similar way. For small distances, Eq. (4.31) has two solutions ω ± slightly shifted below the real axis near the dot energy ϵ 0 . As the length of the wire is increased, the number of solutions increases and the poles of the resolvent start approaching the real axis. When the coupling γ is increased, the proliferation of the poles is

-2 -1 0 1 2 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 Re(ω + ) Im(ω + ) Lϵ 0 /π v=1 -2 -1 0 1 2 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 Re(ω + ) Im(ω + ) Lϵ 0 /π v=1.5 -2 -1 0 1 2 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 Re(ω + ) Im(ω + ) Lϵ 0 /π v=2 -2 -1 0 1 2 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 Re(ω + ) Im(ω + )
Lϵ 0 /π v=2.5 The red dots at their intersection are the solutions ω + . We take ϵ 0 = 0.5, γ = 0.01. As the distance between the dots grows, poles proliferate and approach the real axis.

faster (in the sense that the same amount of poles will appear at smaller values of the distance). We will show that the behavior of the poles depends on the competition between the lifetime of the resonances, associated with the inverse of the coupling strength τ = γ -1 , and the time of propagation in the wire t p = L/v. The far left plot in Fig. 4.8 will turn out to correspond to a regime of weak-coupling γL/v ≪ 1. The next three plots, from left to right, show the transition from the weak-coupling regime to a strong-coupling regime.

Weak-coupling limit

To a first approximation, the coupling term e iωL/v ≃ e iϵ 0 L/v = e ik 0 L , and there are two solutions to Eq. (4.31),

ω ± = ϵ 0 -iγ(1 ± e ik 0 L ). (4.32)
It is easy to see that when the condition k 0 L = nπ, with n an integer, is satisfied, then one of the poles lies on the real axis at ω 1 = ϵ 0 = vk 0 . The mode is therefore "trapped" and does not decay into the continuum. It was von Neumann and Wigner that first proposed the existence of such solutions [START_REF] Neumann | Uber merkwürdige diskrete eigenwerte[END_REF], which are now called "bound states in the continuum" (BICs) [START_REF] Hsu | Bound states in the continuum[END_REF]. At the same time, the second solution acquires a maximum dissipation ω 2 = ϵ 0 -2iγ. This phenomenon is repeated whenever the phase shift k 0 L is an integer multiple of π, meaning that the round trip phase shift is a multiple of 2π (creation of a standing wave). When this standing wave is created, destructive interference causes one of the modes to essentially become uncoupled from the environment, recovering its original energy ϵ 0 , while the second state is more strongly coupled to the environment and becomes short-lived. Such Fabry-Pérot BICs are known to exist in a variety of systems, both classical and quantum-mechanical, whenever two resonances are coupled through the same continuum, see [START_REF] Hsu | Bound states in the continuum[END_REF], and references within. The appearance of a BIC is accompanied by a non-avoided crossing of the real parts. On the contrary, when k 0 L = nπ/2, with n an odd integer, the solutions have equal imaginary parts and a maximum separation of the real parts, ω 1,2 = ϵ 0 ± γiγ. This is illustrated in Fig. 4.9. To make this plot, we use Eq. (4.31) without making an approximation on the energy, and search for the two solutions with the smallest imaginary parts. For γL/v ≪ 1, the approximation that there are only two solutions given by Eq. (4.32) holds, but this starts breaking down as the distance increases.

Re(ω 1 -ω 2 )

Im(ω 1 )
Im(ω 2 ) survival probability at a non-avoided crossing and at an anti-crossing. At a non-avoided crossing, one root is on the real axis ω 1 = ϵ 0 , while the second ω 2 = ϵ 0 -2iγ decays with twice the original decay rate. It is then easy to calculate the survival probabilities P i (t) = |g i (t)| 2 using contour integration

0 1 2 3
g 1 (t) = 1 2π dω i(ω -ϵ 0 + iγ)e -iωt (ω + iη -ϵ 0 )(ω + iη -ϵ 0 + 2iγ) = 1 2 e -iϵ 0 t (1 + e -2γt ), g 2 (t) = 1 2π dω ±γe -iωt (ω + iη -ϵ 0 )(ω + iη -ϵ 0 + 2iγ) = ∓ 1 2 e -iϵ 0 t (1 -e -2γt
).

(4.33)

In the long time limit the probabilities therefore stay finite, meaning that a "part" of the particle stays trapped on the dots

lim t→∞ P 1,2 (t) = 1 ± e -2γt 2
L = nπ 2 ,
with n an odd integer, the survival probabilities (red and dashed red line) eventually go to zero, and P 1 decays faster than the Fermi Golden Rule (grey line). Same parameters as Fig. 4.9.

If, however, we place ourselves at an anti-crossing (that is whenever k 0 L = nπ/2, n an odd integer), then we find that the two poles are ω 1,2 = ϵ 0 ± γiγ and the time-dependent amplitudes are

g 1 (t) = 1 2π dω i(ω -ϵ 0 + iγ)e -iωt (ω -ϵ 0 -γ + iγ)(ω -ϵ 0 + γ + iγ)
= e -iϵ 0 t e -γt cos γt

g 2 (t) = 1 2π dω ±iγe -iωt (ω -ϵ 0 -γ + iγ)(ω -ϵ 0 + γ + iγ)
= ∓e -iϵ 0 t e -γt sin γt

(4.34)
Therefore, the survival probability of both states at long times goes to zero. Since |cos γt| 2 < 1, a particle initially on the first dot dissipates faster than Fermi's golden rule. 2 The conclusions of this paragraph are summarized on Fig. 4.10, with black lines representing the behaviour when the levels cross and red lines for anti-crossings.

2 If there is no coupling between the two dots, the amplitude on the first dot is simply

g 1 (t) = 1 2π dω ie -iωt ω -ϵ 0 + iγ = e -iϵ 0 t e -γt ,
and the survival probability is given by Fermi's Golden Rule,

P 1 (t) = |g 1 (t)| 2 = e -2γt .

General solution

However, the above was only a weak-coupling approximation to the solution of Eq. (4.31). As seen from the numerical solution, Fig. 4.8, this approximation only captures the behavior of the solutions which are near the real axis, essentially ignoring the series of poles below. Starting from Eq. (4.30), we can rewrite the amplitudes on the dots as

g 1 (ω) = i 2 1 ω -ϵ 0 + iγ(1 -e iωL/v ) + 1 ω -ϵ 0 + iγ(1 + e iωL/v ) g 2 (ω) = i 2 1 ω -ϵ 0 + iγ(1 -e iωL/v ) - 1 ω -ϵ 0 + iγ(1 + e iωL/v ) (4.35)
Using the series expansion

1 ω -ϵ 0 + iγ(1 ± e iωL/v ) = +∞ ∑ n=0 (∓iγ) n e inωL/v (ω -ϵ 0 + iγ) n+1 , (4.36)
we get

g 1 (ω) = i ∑ n even (iγ) n e inωL/v (ω -ϵ 0 + iγ) n+1 , g 2 (ω) = i ∑ n odd (iγ) n e inωL/v (ω -ϵ 0 + iγ) n+1 . (4.37)
The inverse Fourier transform for the amplitude on the first dot follows by contour integration. 3

g 1 (t) = ∑ n even γ n n! (t -nL/v) n e -iϵ 0 (t-nL/v) e -γ(t-nL/v) θ(t -nL/v) = ∑ n even (γL/v) n n! (t/t p -n) n e -ik 0 L(t/t p -n) e -γL/v(t/t p -n) θ(t/t p -n). (4.38)
The amplitude on the second dot g 2 (t) can be found in the same way, with the only difference being that the sum is over odd integers.

3 For the pole of order n + 1 the inverse Fourier transform will require calculation of a term

i 2π dω (iγ) n e -iω(t-nL/v) (ω -ϵ 0 + iγ) n+1 .
Since the pole is in the lower-half of the complex plane, the integral is As seen from Eq. (4.38), there are two competing time lengths in this problem: the first one is the time of propagation in the wire, t p = L/v, while the second is the lifetime of the resonances on the dots, given as the inverse of their width, τ = 1/γ. The coupling strength γL/v = t p /τ therefore expresses the competition between the two timescales.

zero if t < nL/v. If t > nL/v, we get i 2π dω (iγ) n e -iω(t-nL/v) (ω -ϵ 0 + iγ) n+1 = (iγ) n n! (-i(t -nL/v)) n e -iϵ 0 (t-nL/v) e -γ(t-nL/v) θ(t -nL/v) P 1 (t) P 2 (t) 0
• At the limit where γL/v < 1, we are at a weak-coupling limit, and only the first few terms of the sum in Eq. (4.38) contribute. We will then recover the results of the previous discussion. In this case only the two poles of the resolvent near the real axis contribute to the time evolution.

• If the lifetime of the resonances is much shorter than the time of propagation, γL/v > 1, the image is that of a series of pulses: a wavepacket starting from dot 1 at time t = 0 arrives at dot 2 at time t p and returns on dot 1 at time 2t p . The wavepacket therefore oscillates between the two dots, and progressively dissipates. This is shown in Fig. 4.14. We see that the great number of poles in Fig. 4.8 results in the appearance of oscillations in the time domain.

Connection with transport

There is a connection between the S-matrix and the resolvent of the system (given by what is called the Mahaux-Weidenmüller formula in the literature) which permits to make a connection between scattering and decay problems [START_REF] Dittes | The decay of quantum systems with a small number of open channels[END_REF]. We imagine that the wire is connected, on either end, to a reservoir with a well-defined chemical potential. If a voltage bias V is applied between the two reservoirs, we can define µ L = ε F + eV and µ R = ε F . The imbalance between the reservoirs creates scattering states which propagate in the wire. If we assume a linearized dispersion ω = vk, there are two branches describing motion in the leads, a left-moving branch with ω L = -vk and a right moving one with ω R = +vk. The corresponding parts of the wavefunction are written as ψ L,R (x).

Then the Schrödinger's equation is

(ω -ϵ 1 )g 1 = J(ψ R (x 1 ) + ψ L (x 1 )) (ω -ϵ 2 )g 2 = J(ψ R (x 2 ) + ψ L (x 2 )) (4.39)
and

ωψ R = -iv ∂ψ R ∂x + Jδ(x -x 1 )g 1 + Jδ(x -x 2 )g 2 ωψ L = +iv ∂ψ L ∂x + Jδ(x -x 1 )g 1 + Jδ(x -x 2 )g 2 (4.40)
In order to construct the scattering matrix we need to build diffusive states.

With the ansatz that ψ L,R are plane waves away from the positions of the scatterers at x 1 , x 2 , we write

ψ R (x) =      α 1 e ikx , x < x 1 α 2 e ikx , x 1 < x < x 2 α 3 e ikx , x 2 < x , ψ L (x) =      β 1 e -ikx , x < x 1 β 2 e -ikx , x 1 < x < x 2 β 3 e -ikx , x 2 < x (4.41)
The situation is summarized in Fig. 4.12. By regularizing the wavefunction at the positions of the dots

Ψ(x 1 ) = 1 2 Ψ(x - 1 ) + Ψ(x + 1 ) Ψ(x 2 ) = 1 2 Ψ(x - 2 ) + Ψ(x + 2 ) (4.42)
Eq. (4.41) and (4.39) give

(ω -ϵ 1 )g 1 = J 2 (α 1 + α 2 )e ikx 1 + J 2 (β 1 + β 2 )e -ikx 1 (ω -ϵ 2 )g 2 = J 2 (α 2 + α 3 )e ikx 2 + J 2 (β 2 + β 3 )e -ikx 2 . ( 4.43)
Moreover, from Eq. (4.40) we get an expression for the discontinuities at

x 1 , x 2 α 2 -α 1 = -i J v e -ikx 1 g 1 α 3 -α 2 = -i J v e -ikx 2 g 2 β 2 -β 1 = i J v e ikx 1 g 1 β 3 -β 2 = i J v e ikx 2 g 2 .
(4.44)

x = x 1 x = x 2 α 2 β 2 α 1 β 1 α 3 β 3
Figure 4.12: Schema for the scattering states.

Substituting Eq. (4.44) in Eq. (4.43), we get an expression of the amplitudes on the dots as a function on the ingoing amplitudes

ω -ϵ 1 + i J 2 v i J 2 v e ikL i J 2 v e ikL ω -ϵ 2 + i J 2 v g 1 g 2 = J e ikx 1 e -ikx 1 e ikx 2 e -ikx 2 α 1 β 3 (4.45)
The inverse of the matrix on the left hand side is the resolvent R(ω), so that

g 1 g 2 = JR(ω)W(k) α 1 β 3 , with W(k) ≡ e ikx 1 e -ikx 1 e ikx 2 e -ikx 2 . (4.46)

Scattering matrix

The scattering matrix gives an expression of the outgoing amplitudes β 1 , α 3 , as a function of the in-going amplitudes α 1 , β 3 :

α 3 β 1 = S α 1 β 3 . (4.47)
By inspection, the scattering matrix can be defined as

S = t r ′ r t ′ (4.48)
where t and t ′ are transmission amplitudes and r, r ′ are reflection amplitudes. Using Eq. (4.44) the outgoing amplitudes can be written as a function of the amplitudes on the dots

α 3 β 1 = α 1 β 3 - iJ v e -ikx 1 e -ikx 2 e ikx 1 e ikx 2 g 1 g 2 (4.49)
which, together with Eq. (4.45), leads to

α 3 β 1 = α 1 β 3 -iγW † (k)R(ω)W(k) α 1 β 3 . ( 4

.50)

The scattering matrix is therefore connected to the resolvent by the relation

S(k) = 1 -iγW † (k)R(k)W(k). (4.51)

Transfer matrix

The transfer matrix expresses states on the right of a scatterer as a function of the states on its left. It provides a convenient expression of the total transmission amplitude as a sum of terms involving an increasing number of reflections in the region between the two dots. The transfer matrix at x 1 can be defined as

α 2 β 2 = M 1 α 1 β 1 = 1 ω -ϵ 1 ω -ϵ 1 -iγ -iγe -2ikx 1 iγe 2ikx 1 ω -ϵ 1 + iγ α 1 β 1 . (4.52)
At x 1 we can also define the scattering matrix

S 1 ≡ t 1 r ′ 1 r 1 t ′ 1 , through the relation α 2 β 1 = S 1 α 1 β 2 .
Rewriting the system in Eq. ( 4.52), one can get a relationship between the elements of the transfer matrix and the scattering matrix at x 1 , namely

t 1 r ′ 1 r 1 t ′ 1 = 1 M 22 1 1 M 12 1 -M 21 1 1 = 1 ω -ϵ 1 + iγ ω -ϵ 1 -iγe -2ikx 1 -iγe 2ikx 1 ω -ϵ 1 .
(4.53)

Similarly, we define the transfer matrix at x 2

α 3 β 3 = M 2 α 2 β 2 = 1 ω -ϵ 2 ω -ϵ 2 -iγ -iγe -2ikx 2 iγe 2ikx 2 ω -ϵ 2 + iγ α 2 β 2 , ( 4.54)
as well as the scattering matrix

t 2 r ′ 2 r 2 t ′ 2 = 1 ω -ϵ 2 + iγ ω -ϵ 2 -iγe -2ikx 2 -iγe 2ikx 2 ω -ϵ 2 . ( 4.55)
The transfer matrix construction is convenient, because it is straightforward to combine the individual transfer matrices in order to get the total transfer matrix: Since

α 3 β 3 = M 2 α 2 β 2 = M 2 M 1 α 1 β 1 , the total transfer matrix is equal to M = M 2 M 1 .

Transmission coefficient

One quantity of interest is the transmission coefficient, T(ω) = |t| 2 . In the Landau-Bütiker formalism, the average current under a voltage bias eV = µ L -µ R is given by integrating the transmission coefficient over the available states [START_REF] Datta | Electronic transport in mesoscopic systems[END_REF]

I = - 2e h T(ω)[ f L (ω) -f R (ω)] dω = e h ε F +eV ε F T(ω) dω , (at temperature T = 0). (4.57)
Then, the differential conductance at zero temperature has a rather simple form:

G(V) = dI dV = 2e 2 h T(ε F + eV). (4.58)
From Eq. (4.56), the transmission amplitude can be written as a geometric series

t = 1 M 22 = 1 M 2 M 1 22 = t 1 t 2 1 -r ′ 1 r 2 = t 1 t 2 +∞ ∑ n=0 (r ′ 1 r 2 ) n = t 1 t 2 + t 1 r ′ 1 r 2 t 2 + t 1 r ′ 1 r 2 r ′ 1 r 2 t 2 + • • • (4.59)
The physical interpretation is straightforward. The first term of the geometric series corresponds to the amplitude for transmission through the wire without any reflections. The second term corresponds to transmission after two reflections, the third after four reflections, and so on. The system is therefore a realization of an electronic Fabry-Pérot. We find:

t = (ω -ϵ 1 )(ω -ϵ 2 ) (ω -ϵ 1 + iγ)(ω -ϵ 2 + iγ) + γ 2 e 2iωL/v (4.60)
Assuming a detuning between the dots ϵ 1,2 = ϵ 0 ± δ, and setting ϵ 0 = 0, the transmission amplitude is

t = ω 2 -δ 2 (ω + δ + iγ)(ω -δ + iγ) +∞ ∑ n=0 (-1) n γ 2n e 2inωL/v (ω + δ + iγ) n (ω -δ + iγ) n . (4.61)
This expression can be scaled with respect to the natural energy scale related to a wire of length L, ω L ≡ v/L. Then, the transmission amplitude can be written in reduced energy units of ω = ω/ω L :

t( ω) = ω 2 -(δL/v) 2 ( ω + δL v + i γL v )( ω -δL v + i γL v ) +∞ ∑ n=0 (-1) n (γL/v) 2n e 2in ω ( ω + δL v + i γL v ) n ( ω -δL v + i γL v ) n . ( 4
.62) As we saw in the previous section, the strength of the coupling is given by the competition between the lifetime of the resonances τ and the propagation time t p :

t p /τ = γL/v.
This parameter explicitly controls the strength of the reflection processes in Eq. (4.62). weak coupling limit. In the weak coupling limit where γL/v < 1, the higher order terms in the sum can be ignored, and the transmission is simply given by a direct transmission through the system T = |t 1 t 2 | 2 . The physics is then described by the overlap between two Fano resonances. If the dots are identical with no detuning δ = 0, then the total transmission amplitude is

t ≃ ω ω + 2iγ .
As we saw in the previous section, in the weak coupling regime in the absence of detuning and for ε 0 = 0, one resonance becomes a BIC with zero width while the other becomes twice as dissipative. Comparing t with the transmission amplitude through a single dot t 1 = ω ω+iγ , one can conclude that the overall effect of adding a second dot is a doubled width 2γ. This is shown in Fig. 4.13a. As shown in Fig. 4.13b, introducing a small detuning between the dots leads to the system becoming transparent in a very narrow region [START_REF] Mingaleev | Coupledresonator-induced reflection in photonic-crystal waveguide structures[END_REF], while at strong detuning δ ≫ γ, one recovers two isolated resonances at ±δ. These phenomena have been experimentally observed, mostly in photonic systems, see for example the review [START_REF] Miroshnichenko | Fano resonances in nanoscale structures[END_REF] and the references within. I note that in the literature of the Fano resonances the limit of overlapping resonances is usually called a "strong coupling" limit. This is because the scales that are usually involved are the width of the resonances with respect to their detuning. Here, however, the situation is different since we have to compare three different energy scales. strong coupling limit. In the strong coupling limit γL/v > 1, the transmission is a sum of Feynman paths involving an increasing number of reflections between the two dots. The interference caused is expressed by the phase e 2inωL/v , accumulated after an even number 2n of reflections. The interference produces oscillations of the transmission coefficient as a function of the energy, as shown in Fig. 4.14. This will result in steps in the current. The position of the transport peaks corresponds to the position of the resolvent poles in the complex plane. This is more evident in Fig. 4.15 where a "phase diagram" of the system is presented, showing the transition from the weak to the strong coupling regime. For each selection of detuning and coupling strength, we plot the transmission coefficient together with the poles of the resolvent in the complex plane, as well as the resulting current.

discussion. We presented a simplified model that captures the transition from a regime of isolated resonances to an interferometer regime. This model shares some features with the driven S-QD-S-QD-S system since both models involve the coupling of two quantum dots through a continuum of states. In particular, they both show a proliferation of the resolvent poles at large distances and concomitant oscillations in the observables. However, we must note that there are important differences:

• We considered a normal wire instead of a superconducting one. Any effect related to the superconducting gap, as well as the singularity of the BCS density of states at the gap edges was not taken into account.

• Superconductivity mixes electrons and holes. The Floquet-Tomasch effect involves two different processes, represented by the self-energy matrices Σ ± that couple the electrons on one dot to the holes on the other. These processes have a more complicated dependence on the energy and on the applied voltage. In particular, they can become resonant for a Floquet-Andreev state if by MAR processes they connect it to the edge of the superconducting gap.

reduction to a two level system

With an appropriate transformation, we can show that the basic physics of the system at large interdot distances is that of two resonances coupled through a continuum. The resulting effective Hamiltonian is non-Hermitian, which is a result of the fact that we have focused on the Hamiltonian of a subsystem.

The linear operator L 00 can be transformed into the basis where the matrices M

(1,2) 0 of the uncoupled dots are diagonal. We will assume identical dots for simplicity, so that the matrices M

(1,2) 0 have the same pair of eigenvalues ±E 0 . We will take into account that due to the particle-hole symmetry of the BdG states the roots of the characteristic polynomial det L 00 = 0 come in pairs (if E is an eigenvalue, so is -E * ). We can then focus only on the positive sector of energies, assuming that the coupling between positive and negative energy states is small.

Effective Hamiltonian

For large voltage bias and small tunnel couplings we can approximate the matrices

M (1,2) 0 , assuming |E| ≪ ω 0 , ∆, by M (1,2) 0 ≈   E ∓ Γ (a,b) ω 0 √ ∆ 2 -ω 2 0 -Γ c,(1,2) -Γ c,(1,2) E ± Γ (a,b) ω 0 √ ∆ 2 -ω 2 0   . ( 4.63)
The roots of the characteristic polynomial det M

(1,2) 0 = 0 are

E ± (1,2) = ± Γ 2 c,(1,2) +   Γ (a,b) ω 0 ∆ 2 -ω 2 0   2 . ( 4.64)
For simplicity, we will assume that the dots are identical, meaning that we take the couplings Γ a = Γ b and Γ c1 = Γ c2 = Γ c . Then, the two matrices M

(1,2) 0 have the same pair of eigenvalues

±E 0 ≡ ± Γ 2 c +   Γ a ω 0 ∆ 2 -ω 2 0   2 ,
but different eigenvectors. Indeed, the structure of the matrices M

(1,2) 0 means that the corresponding eigenvectors satisfying M (j) 0 ψ ± j = 0 can be parametrized as

ψ + 1 = cos θ sin θ , ψ - 1 = -sin θ cos θ , (4.65a) ψ + 2 = sin θ cos θ , ψ - 2 = -cos θ sin θ . (4.65b)
We see that the "electron" and "hole" components of the eigenvectors are, in fact, reversed. This is because the two dots are coupled to reservoirs S a and S b which are biased with opposite voltages. Moreover, we can derive a simple relation for the angle θ involving the tunnel couplings and the voltage frequency ω 0 :

θ = 1 2 arctan   Γ c Γ a ∆ 2 -ω 2 0 ω 0   . (4.66)
Within our approximation that ∆ 2 < ω 0 < ∆, we can deduce from the above relation that the principal value of the angle is θ ∈ [0, π/4[. This angle therefore controls the electron/hole content of the eigenvectors.

Having found the eigenvectors, we can construct change of basis matrices P(θ) and Q(θ) that diagonalize M

0 = P(θ)DP(θ) -1 = cos θ -sin θ sin θ cos θ E -E 0 0 0 E + E 0 cos θ sin θ -sin θ cos θ , M (2)
0 = Q(θ)DQ(θ) -1 = sin θ -cos θ cos θ sin θ E -E 0 0 0 E + E 0 sin θ cos θ -cos θ sin θ . (4.67)
Using these, we can transform the initial Floquet operator L 00 in the basis of the eigenvectors

ψ + 1 , ψ - 1 , ψ + 2 , ψ - 2 : L = P(θ) -1 0 0 Q(θ) -1 L P(θ) 0 0 Q(θ) = P -1 M
(1)

0 P 0 0 Q -1 M (2) 0 Q - P -1 0 0 Q -1 Σ + 0 + Σ - 0 P 0 0 Q . (4.68)
By permutation of the basis vectors ψ - 1 ⇌ ψ + 2 , we can rewrite L in order to make apparent the two blocks which correspond to positive and negative eigenvalues. To lowest order of perturbation in the tunnel couplings, we can neglect the non-diagonal blocks in L. This amounts to neglecting the coupling between positive and negative eigenvalue sectors. For the upper-left block of L we then obtain We see that the parameter θ (defined in Eq. 4.66) controls the relative strength of the self-energy processes Σ ± , where Σ -connects the dots through S c at energies below the gap E -2ω 0 < -∆, and Σ + connects the dots through quasiparticle propagation in the middle superconductor at energies above the gap E + 2ω 0 > ∆. The parameter θ itself can be controlled by the voltage and the couplings which change the relative weights of the electronlike and holelike components of the eigenvectors.

L ++ = E -E 0 0 0 E -E 0 -cos 2 θ Σ -,e 1 e 1 0 Σ -,e 1 h 2 0 Σ -,h 2 e 1 0 Σ -,h 2 h 2 0 -sin 2 θ Σ +,h 1 h 1 0 Σ +,h 1 e 2 0 Σ +,e 2 h 1 0 Σ +,e 2
The resulting effective operator is of the form

L eff = E -E 0 + iγ 1 iγ 12 iγ 21 E -E 0 + iγ 2 . ( 4.70)
The above relation describes the coupling of two discrete levels initially at E 0 , which are coupled through a continuum of states. The overall action of the continuum is, as expected, to add a small shift to E 0 equal to the real part of the diagonal self-energy elements and a width equal to their imaginary part. Moreover, the two resonances are then coupled through the non-diagonal elements of the self-energies. However, the particularity of this system is that the self-energy processes that couple the dots are energy-dependent, and the imaginary parts γ, defined below, are themselves functions of the energy E, the voltage bias ω 0 and the distance R between the resonances. Explicitly, the diagonal components of the self-energy matrices will add a finite lifetime to the discrete levels at E 0 , given by:

iγ 1 = -Γ 2 a cos 2 θ • g eh a (-1) 1 M 0 -2 h 1 h 1 g he a (-1) -Γ 2 a sin 2 θ • g he a (1) 1 M 0 2 e 1 e 1 g eh a (1), (4.71a
)

iγ 2 = -Γ 2 b cos 2 θ • g he b (-1) 1 M 0 -2 e 2 e 2 g eh b (-1) -Γ 2 b sin 2 θ • g eh b (1) 1 M 0 2 h 2 h 2 g he b (1), (4.71b)
and γ 1 = γ 2 ≡ γ for identical dots. The non-diagonal components that couple the two resonances can be written in a form that makes apparent the dependence on the quartet phase ϕ q = ϕ a + ϕ b -2ϕ c , and the Floquet-Tomasch phase factors (defined in 4.17) γ 12 = αe iϕ q e iq(-2)Rβe -iϕ q e iq(2)R , (4.72a) γ 21 = αe -iϕ q e iq(-2)Rβe iϕ q e iq(2)R , (4.72b) where

α = Γ a Γ b Γ c ∆ 3 cos(k F R) cos 2 θ [∆ 2 -(ω 0 -E) 2 ] (2ω 0 -E) 2 -∆ 2 1 M (1) -2 hh 1 M (2) -2 ee , (4.73a) β = Γ a Γ b Γ c ∆ 3 cos(k F R) sin 2 θ [∆ 2 -(ω 0 + E) 2 ] (2ω 0 + E) 2 -∆ 2 1 M (1) 2 ee 1 M (2) 2 hh . ( 4.73b)
Finding the resulting eigenvalues due to the coupling between the two resonances requires finding the zeroes of the characteristic polynomial of L eff . The characteristic polynomial will be a transcendental equation, generally requiring a numerical solution. For identical dots, (E -E 0 + iγ) 2 + α 2 e 2iq(-2)R + β 2 e 2iq(2)R -2αβ cos 2ϕ q e iq(2)R e iq(-2)R = 0 (4.74)

The solutions of the above equation are found numerically and plotted on the complex plane in Fig. 4.16. At small distances, we find two solutions around the initial level E 0 , slightly shifted in the complex plane. As the distance between the dots grows, however, there are more solutions which appear around the two initial ones. The number of the solutions increases with the distance since the factors e iq(±2)R become more rapidly oscillating. Figure 4. [START_REF] Mahmood | Selective scattering between Floquet-Bloch and Volkov states in a topological insulator[END_REF] shows that the effective model roughly captures the expected behavior, i.e., the number of poles increases with increasing interdot distance, in agreement with Fig. 4.4 and Fig. 4.5 that were produced by numerically calculating the full operator L 00 .

Oscillations of the spectral function

From Eq. (4.70) we can calculate the corresponding effective resolvent operator R eff = L -1 eff :

R 11 eff = E -E 0 + iγ (E -E 0 + iγ) 2 + γ 12 γ 21 = ∞ ∑ n=0 (-1) n (γ 12 γ 21 ) n (E -E 0 + iγ) 2n+1 . (4.75)
This is reminiscent of the expressions we found for the toy model in the previous section. Here, the strength of the Floquet-Tomasch effect on the resolvent is controlled by the product γ 12 γ 21 = α 2 e 2iq(-2)R + β 2 e 2iq(2)R -2αβ cos 2ϕ q e iq( 2)R e iq(-2)R .

The system therefore resembles a Fabry-Pérot interferometer, but there are two types of propagating waves involved, corresponding to the two self-energy processes Σ ± . We can associate the parameter α with the backward self-energy Σand with propagation at E -2ω 0 < -∆. Meanwhile, the parameter β comes from the forward self-energy Σ + and is associated with propagation at E + 2ω 0 > ∆. The amount of interference between the two processes is controlled by the quartet phase.

In the weak coupling limit the interference effects can be ignored. Keeping only the first term of the geometric series gives a single resonance at E 0 . What we call a weak-coupling limit is more evident if we scale the resolvent with respect to the energy scale related to propagation at distance R, ω R = v F /R. In reduced energy units

R 11 eff = 1 ω R ∞ ∑ n=0 (-1) n (γ 12 γ 21 R 2 /v 2 F ) n ( E -E 0 R v F + i γR v F ) 2n+1 . (4.76)
The coupling strength is therefore controlled by the rescaled quantities (γ 12 R/v F ) and (γ 21 R/v F ). The amount of terms we need to keep in the series therefore increases with the interdot distance, resulting in the proliferation of poles as shown in Fig. 4.16.

For illustrative purposes, we can consider a voltage value ∆ -E 0 < 2ω 0 < ∆ + E 0 , where only the forward self-energy Σ + contributes around E 0 > 0. Then the expression for the resolvent at real energies close to E 0 simplifies to:

R 11 eff = ∞ ∑ n=0 (-β 2 ) n e 2inq(2)R (E -E 0 + iγ) 2n+1 ≈ 1 E -E 0 + iγ - β 2 e 2iq(2)R (E -E 0 + iγ) 3 + • • • (4.77)
The effect on the spectral function -2 π Im R 11 eff (E) will then be a Breit-Wignerlike resonance around E 0 , coming from the first term (E -E 0 + iγ) -1 , and smaller oscillations superimposed on the resonance due to the second term on the right-hand side. The spectral function will therefore oscillate as a periodic function of a Floquet-Tomasch factor:

2q(2)R = 2R v F (E + 2ω 0 ) 2 -∆ 2 .
Since the interdot coupling term β is proportional to Γ a Γ b Γ c , we expect that the Floquet-Tomasch oscillations are larger in amplitude when increasing the couplings. At the same time, the width of the resonances given by γ is also proportional to the tunnel couplings. Then one expects that the resonances are smeared out with increasing Γ.

Quartet phase

In Eq. (4.74), the quartet phase ϕ q appears in the last term as cos 2ϕ q . This can be related to "octet" processes, as discussed in [START_REF] Mélin | Ultralong-distance quantum correlations in three-terminal Josephson junctions[END_REF]. A sketch of an octet process is shown in Fig. 4.17b. Here, Eq. (4.74) gives us bounds for the appearance of the octets. At large distances, the Floquet-Tomasch phase factors e iq(±2)R , and therefore the corresponding self-energy processes Σ ± 0 , are non-zero at the current order in the tunnel couplings only if the condition |E ± 2ω 0 | > ∆ is satisfied. Then, around a resonance E 0 > 0 the e iq(2)R term contributes when the voltage is 2ω 0 > ∆ -E 0 , while the e iq(-2)R term contributes when 2ω 0 > ∆ + E 0 . As a result, the octet term can only contribute when both processes are present, that is, when 2ω 0 > ∆ + E 0 . There will therefore be a regime of voltages ∆ -E 0 < 2ω 0 < ∆ + E 0 , where only the forward self-energy Σ + contributes around E 0 > 0, while, making an analogous argument, only the backward self-energy Σ -will contribute around -E 0 < 0. When the voltage is increased above 2ω 0 > ∆ + E 0 both processes contribute, but with different weights, since the process with the larger absolute value of energy |E 0 ± 2ω 0 | will start to ex- ponentially decay at energies much larger than the gap. In the regime that the octet term is relevant, the quartet phase can nonlocally control the interdot coupling, since ϕ q can be tuned by changing the phase ϕ b across the second junction, while measuring the spectrum on the first. Equation (4.74) suggests that the amplitude of oscillations is enhanced at ϕ q = 0, and minimized at ϕ q = π/2. Moreover, the quartet phase does not affect the frequency of oscillations, which is rather a function of the energy, the voltage, and the interdot distance. These observations are verified numerically in Fig.

4.17a which shows the variation

Im[

R e 1 e 1 00 -R R (1) . Different colors represent a different quartet phase, ϕ q ∈ [0, π 2 ], as explained in the inset. Other parameters are set to ω 0 = 0.7 ∆, Γ j = 0.5 ∆, R = 25 ξ 0 . (b) Sketch of an octet process Σ -,e 1 h 2 0 × Σ +,e 2 h 1 0 , which couples a hole at energy E 0 on dot 1 to an electron at the same energy on dot 2 and accumulates a phase e 2iϕ q .

without making any approximations on the distance or the voltage, i.e., by calculating the operator L 00 . It is worth noting that Fig. 4.17a implies that the Floquet-Tomasch oscillations will not be smeared out if an average is taken over the quartet phase. The oscillations should therefore be observable even if the quartet phase should drift with time.

conclusions

We have studied two driven superconducting quantum dots connected to a common superconductor. In the limit where the superconductor is long and subgap transport is governed by MARs, we showed that a long-range coupling develops between the dots. We showed that the system can be described by an effective non-Hermitian model of two resonances coupled through higherorder processes that involve local MARs on each dot, followed by a nonlocal Andreev reflection through the common superconductor at energies above the gap. The induced interdot coupling modifies the Floquet spectrum, producing oscillations in the spectral function. The amplitude of these oscillations can be controlled nonlocally by changing parameters like the phase drop across one of the dots. This amounts to tuning the oscillations with the quartet phase, and we have found bounds for which the quartet phase is involved.

It remains to be seen if control of the quartet phase is feasible experimentally at finite voltage bias. This is the topic of recent work, which proposes an interferometric setup sensitive to quartet processes [START_REF] Mélin | Quantum interferometer for quartets in superconducting three-terminal Josephson junctions[END_REF]. It is therefore an open question whether the quartet phase can be used to control the amplitude of the Floquet-Tomasch oscillations. Regardless, our results suggest that the oscilla-tions as a function of the energy will not be smeared out, even if we consider a quartet phase that drifts with time.

Our approach is relevant for well-defined Floquet resonances on the dots, that is, for tunnel couplings to the reservoirs that are not very large Γ < ∆, and for subgap voltage values. We assumed a large subgap voltage bias ∆/2 < ω 0 < ∆ that allows to simplify the analytical part of this work since at strong driving only a few Floquet harmonics need to be taken into account. However, the mechanism that results in a long-range interdot coupling should exist at smaller subgap voltages as well, although at higher order in the tunnel couplings.

Part III

M A S T E R E Q U AT I O N S M A S T E R E Q U AT I O N S F O R N O N -I N T E R A C T I N G F E R M I O N I C S Y S T E M S
This chapter contains preliminary work on the derivation of a master equation in the Floquet basis of states for the S-QD-S system.

introduction

In Chapter 3 we have seen that a finite dc voltage bias turns the equilibrium Andreev bound states into Floquet-Andreev resonances with a finite lifetime. In this chapter, we ask the following questions:"Is it possible to derive a master equation in the basis of the Floquet-Andreev states? What are the populations of these states?"

The standard derivation of a Markov master equation is done by taking a weak tunneling limit. In our case, however, this would imply that MAR processes are washed out. Instead we want the coupling strength to be small with respect to the superconducting gap, i.e., Γ < ∆, but not too small (not Γ ≪ ∆). If we followed the standard perturbative derivation of a master equation (see for example [START_REF] Manzano | A short introduction to the Lindblad master equation[END_REF] or the classic textbooks [START_REF] Breuer | The theory of open quantum systems[END_REF][START_REF] Gardiner | Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics[END_REF]), this would imply that higher-order terms in the coupling strength need to be taken into account. We will therefore take an alternative route, which consists in calculating the reduced density matrix ρ s , as well as correlation functions, with the help of the resolvent. In fact, the time evolution of ρ s (t) is determined by the pole structure of the resolvent. We will replace the demand of Γ ≪ ∆ by the demand that the coupling to the reservoirs is sufficiently small so that the Floquet-Andreev resonances are narrow and, moreover, not overlapping. In other words, we will make a quasiparticle approximation, where the quasiparticles associated to the Floquet-Andreev states will be supposed well-defined.We will first apply the method to the case of a multilevel dot coupled to non-interacting fermionic leads. We will then consider the case of a driven dot, where we will derive the master equation and populations in a Floquet basis of states. Finally, we will consider a dot coupled to superconducting leads which are voltage biased with commensurate voltages, and we will derive an expression for the populations of the Floquet-Andreev states.

Open quantum systems

Up until now the analysis is exact and we have done nothing but a rewriting of the Schrödinger equation for the case of a statistical mixture of states described by the density matrix ρ. We suppose, however, that we are studying an open system, i.e., a system interacting with an environment. We therefore consider a small quantum system described by some Hamiltonian H s and density matrix ρ s , a reservoir described by some Hamiltonian H r and density matrix ρ r and a coupling between them V. The dynamics of the total system is given by the von Neumann equation (5.5) as before. We are interested in the dynamics of the small quantum system, given by the reduced density matrix ρ s (t) = Tr r (ρ), obtained by tracing out the degrees of freedom of the reservoir. The time evolution of the reduced density matrix is no longer unitary, but is given by a propagator V(t, t 0 ) :

ρ s (t) = V(t, t 0 )ρ s (t 0 ) (5.8)
with the initial condition V(t 0 , t 0 ) = 1. The Markovian approximation consists in assuming that the evolution superoperator V can be written in exponential form V(t, t 0 ) = e L(t-t 0 ) . This results in a Markovian master equation for the reduced density matrix [START_REF] Manzano | A short introduction to the Lindblad master equation[END_REF][START_REF] Breuer | The theory of open quantum systems[END_REF][START_REF] Gardiner | Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics[END_REF]:

ρs (t) = L(t)ρ s (t) (5.9)
Note that this is not the quantum Liouville equation (5.7), since the density matrix involved is the density matrix of the subsystem and not of the total system. As a result, the superoperator L is not simply given by the commutator with a Hamiltonian, but will contain additional terms that describe the fact that an open system has dissipation. The Markovian master equation in Lindblad form is the most general form that preserves desired properties of the density matrix (its trace and positivity). The Lindblad equation (or Gorini-Kossakowski-Sudarshan-Lindblad equation) [START_REF] Lindblad | On the generators of quantum dynamical semigroups[END_REF][START_REF] Gorini | Completely positive dynamical semigroups of N-level systems[END_REF] can be written as

ρs (t) = -i[H eff , ρ s (t)] + D[ρ s (t)] = -i[H eff , ρ s (t)] + ∑ i γ i L i (t)ρ s (t)L † i (t) - 1 2 L † i (t)L i (t), ρ s (t) .
(5.10)

The commutator corresponds to the unitary part of the dynamics and usually contains a contribution from the reservoirs (a Lamb shift of the energy levels). H eff is therefore not necessarily equal to the original Hamiltonian H s of the reduced system. The dissipator part represents stochastic processes due to the reservoirs, encoded in the jump operators L i .

method summary

Knowledge of the reduced density matrix suffices to calculate expectation values of system operators, for example ⟨A(t)⟩ = Tr(A(t)ρ s (t)). However, in order to have complete knowledge of how the system behaves, we would like to be able to compute correlation functions between observables at different times, for example: .11) This is possible if we know the propagator V :

C AB (t, t ′ ) = A(t)B(t ′ ) . ( 5
C AB (t, t ′ ) = Tr AV(t, t ′ )BV(t ′ , 0)ρ s (0) , t > t ′ ≥ 0.
(5.12)

The convention is that superoperators act on all operators to their right, so the notation VBρ actually means V[Bρ].

The usual derivation of the master equation (5.10) is to assume weak-coupling between the system and the reservoirs and to apply a perturbative treatment (for a review, see [START_REF] Manzano | A short introduction to the Lindblad master equation[END_REF]). The reason we do not wish to follow this approach is that in the case of a voltage-biased superconducting dot, MAR processes would require calculation of high-order terms in perturbation theory. However, since we are dealing with non-interacting fermionic systems, Wick's theorem holds. This means that in principle we are able to calculate two-time correlation functions of the form C AB (t, t ′ ) = ⟨A(t)B(t ′ )⟩ using the time-evolution imposed by the poles of the resolvent of the system.

We assume that we have a small subsystem which consists of a finite number of discrete states described by a Hamiltonian H s and coupled to some noninteracting fermionic reservoirs. For example, the small system can be a multilevel quantum dot, a chain of dots, a driven quantum dot, or a superconducting dot.

In the absence of coupling to the reservoirs, the retarded Green's function on the dot can be written as

G 0 s (ε) = 1 ε -H s + iη = dim H s ∑ a=1 |ψ a ⟩⟨ψ a | ε -ε a + iη , ( 5.13)
where H s is the Hilbert space of the small system and ⟨ψ a |ψ b ⟩ = δ ab . When the coupling is switched on, the reservoirs modify the states on the dot. This is expressed through a non-zero self-energy σ(ε). If the coupling is not too strong, the eigenvalues ε a are slightly pushed below the real axis and we can replace them by ε a (ε)iγ a (ε), where both ε a (ε) and γ a (ε) are real quantities. The self-energy is a non-Hermitian operator, expressing the fact that the reservoirs induce dissipation on the states of the quantum dot, which therefore acquire a finite lifetime. For a non-Hermitian matrix, however, it is no longer necessarily true that an orthonormal set of eigenvectors can be constructed. This means that if

(H s + σ(ε)) |ψ a (ε)⟩ = (ε a (ε) -iγ a (ε)) |ψ a (ε)⟩ then, generally, ⟨ψ a |ψ b ⟩ ̸ = δ ab . Indeed, it is easy to check that ⟨ψ a |ψ b ⟩ = -2i ⟨ψ a |Im σ|ψ b ⟩ (ε a -ε b ) + i(γ a + γ b )
, where 2i Im σ = σσ † . In order to proceed, one introduces the right eigenstates of the Hermitian adjoint

(H s + σ(ε)) † |ϕ a (ε)⟩ = (ε a (ε) + iγ a (ε)) |ϕ a (ε)⟩ .
The above relation implies that ⟨ϕ a | is a left eigenvector of (H s + σ(ε)) with eigenvalue ε a (ε)iγ a (ε). One can then show that an orthogonality relation exists between the left and right eigenstates [START_REF] Brody | Biorthogonal quantum mechanics[END_REF][START_REF] Ashida | Non-hermitian physics[END_REF] The resolvent can then be written in this biorthogonal basis as

ϕ a (ε) ψ b (ε ′ ) = δ ab δ(ε -ε ′ ),
G s (ε) = ∑ a |ψ a (ε)⟩⟨ϕ a (ε)| ε -ε a (ε) + iγ a (ε) (5.14)

Narrow resonance approximation

We assume a weak coupling to the reservoirs, so that the real part of the poles ε a (ε) is only slightly modified with respect to the eigenvalues ε a of H s , so that the equation ε = ε a (ε) will have solutions E a close to the unperturbed eigenvalues ε a . Moreover, γ a (ε) is expected to be small so that the poles of the resolvent in the lower half of the complex plain remain near the real axis. In other words, we assume that the quasiparticles associated with the poles of the resolvent are well-defined. We can then make a narrow resonance approximation, which amounts to approximating the resolvent by

G s (ε) ≃ ∑ a Z a |ψ a (E a )⟩⟨ϕ a (E a )| ε -E a + iΓ a (5.15)
The resolvent is defined by three quantities:

• The real part of the poles, defined by a self-consistent equation

E a = ε a (E a )
Generally, this equation requires a numerical calculation.

• The damping given by the imaginary part, related to the inverse lifetime • A renormalization factor

Γ a = Z a γ a (E a )
Z a = 1 - dε a (ε) dε -1 ε=E a .
In real time, one therefore obtains a time-evolution according to the propagator .16) This allows calculating the time evolution of observables and correlation functions in the basis of the residues of the resolvent. It is important to note that the approach breaks down in the presence of avoided crossings when the separation between any two resonances becomes comparable to their width.

G s (t, t ′ ) = -i θ(t -t ′ ) ∑ a Z a e -iE a (t-t ′ ) e -Γ a (t-t ′ ) |ψ a (E a )⟩⟨ϕ a (E a )| . ( 5

multilevel dot coupled to normal leads

The first example we have in mind is to study the transport through a multilevel dot as illustrated in Fig. 5.1. The quantum dot is considered to be a small noninteracting system composed of a finite number of discrete levels, described by the Hamiltonian H s . The system is then coupled to some metallic reservoirs which are described by a Hamiltonian H r

H r = ∑ j=L,R dε ε |εj⟩⟨εj| , (5.17)
and the continuum states are normalized with ⟨εj|ε ′ j ′ ⟩ = δ jj ′ δ(εε ′ ).

We are interested in the out of equilibrium case where the reservoirs' chemical potentials are different, µ L > µ R and at least some discrete levels of the system have an energy that is between the conduction and the valence band. The coupling of the system to the reservoirs is

V = ∑ j dε v j (ε) εj + εj v j (ε) , (5.18)
where v j (ε) ∈ H s . In order to proceed, we construct dressed states

|Ψ(εi)⟩ = |εi⟩ + |Ψ s (εi)⟩ + |Ψ r (εi)⟩ , (5.19)
which are solutions of the Lippmann-Schwinger equation

|Ψ(εi)⟩ = |εi⟩ + G 0 (ε)V |Ψ(εi)⟩ (5.20)
where G 0 (ω) = (ω -H s -H r + iη) -1 is the retarded Green's function when the coupling is turned off V = 0. Plugging Eq. (5.19) into Eq. (5.20) one gets two coupled equations:

(ε -H r + iη) |Ψ r (εi)⟩ = ∑ j dε ′ ε ′ j v j (ε ′ ) Ψ s (εi) (ε -H s + iη) |Ψ s (εi)⟩ = ∑ j dε ′ v j (ε ′ ) δ ij δ(ε -ε ′ ) + ε ′ j Ψ r (εi) (5.21)
The first equation gives

|Ψ r (εi)⟩ = ∑ j dε ′ |ε ′ j⟩ v j (ε ′ ) Ψ s (εi) ε -ε ′ + iη . ( 5

.22)

Introducing the self-energy operator

σ(ε) = ∑ j dε ′ v j (ε ′ ) v j (ε ′ ) ε -ε ′ + iη , ( 5.23)
the states on the dot are obtained through

|Ψ s (εi)⟩ = G s (ε) |v i (ε)⟩ , (5.24)
where the resolvent is defined as

G s (ε) = (ε -H s -σ(ε) + iη) -1 .
(5.25)

Reduced density matrix

We assume that the basis of dressed states is complete:

1 = ∑ i dε |Ψ(εi)⟩⟨Ψ(εi)| .
Moreover, we assume that in the steady-state the populations of the reservoirs fix expectation values of the total system's operators. This implies that the total density matrix of the system can be written as

ρ = ∑ i dε n i (ε) |Ψ(εi)⟩⟨Ψ(εi)| . ( 5

.26)

The reduced density matrix, corresponding to the part on the small system can then be written as

ρ s = ∑ i dε n i (ε) |Ψ s (εi)⟩⟨Ψ s (εi)| . ( 5

.27)

The states on the dot will make appear the resolvent, which now involves the influence of the reservoirs through the self-energy. Using Eq. (5.24) we can write

ρ s = ∑ i dε n i (ε)G s (ε) |v i (ε)⟩⟨v i (ε)| G † s (ε).
(5.28)

Using the narrow resonance approximation for the resolvent and inserting into Eq. (5.28), we obtain

ρ s = ∑ i dε n i (ε) ∑ a,b Z a Z b |ψ a (E a )⟩ ⟨ϕ a (E a )|v i (ε)⟩ ε -E a + iΓ a ⟨v i (ε)|ϕ b (E b )⟩ ⟨ψ b (E b )| ε -E b -iΓ b
(5.29) In the weak coupling limit, the integrand will have sharp peaks at energies E a and E b with corresponding widths Γ a and Γ b . The dominant term in the sum will then correspond to a resonant condition a = b. The Lorentzian distribution which appears will then approach a Dirac delta function lim

Γ a →0 1 (ε -E a ) 2 + Γ 2 a = π Γ a δ(ε -E a ),
resulting in the expression

ρ s ≃ ∑ a ∑ i πn i (E a )Z 2 a Γ a |⟨ϕ a (E a )|v i (E a )⟩| 2 |ψ a (E a )⟩⟨ψ a (E a )| . ( 5.30)
We introduce the dressed basis

∑ i πZ 2 µ ϕ µ (E µ ) v i (E µ ) 2 Γ µ ψ µ (E µ ) ψ µ (E µ ) ≡ f µ f µ (5.31)
which in the narrow resonance approximation is complete ∑ dim H s µ=1 f µ f µ = 1 and orthonormal f µ f ν = δ µν . We can then rewrite the reduced density matrix in the dressed basis in the simpler form

ρ s = ∑ µ p µ f µ f µ (5.32)
where p µ are the populations defined as

p µ ≡ ∑ i n i (E µ ) ϕ µ (E µ ) v i (E µ ) 2 ∑ i ϕ µ (E µ ) v i (E µ ) 2 .
(5.33)

Correlation functions

For non-interacting fermionic systems Wick's theorem holds. This is true in the case of a non-interacting dot connected to metallic leads, but also holds for the case of superconducting leads, since in the case of superconductivity the fact that the mean-field Hamiltonian is quadratic in the fermionic creation and annihilation operators means that it is effectively non-interacting. We can therefore calculate any higher-point correlation function by reducing it to a sum of two-point correlators. For example, consider the correlation function

⟨c † a (t)c a (t)c † b (t ′ )c b (t ′ )⟩ = ⟨c † a (t)c a (t)⟩ ⟨c † b (t ′ )c b (t ′ )⟩ + ⟨c † a (t)c b (t ′ )⟩ ⟨c a (t)c † b (t ′ )⟩ -⟨c † a (t)c † b (t ′ )⟩ ⟨c a (t)c b (t ′ )⟩ (5.34)
where the c, c † are operators of the small system. The anomalous terms ⟨c † a (t)c † b (t ′ )⟩ and ⟨c a (t)c b (t ′ )⟩ are nonzero only in the superconducting case. In the normal case, we therefore have to calculate the first two terms on the right-hand side. For a two-time correlation function

C AB (t, t ′ ) ≡ A(t)B(t ′ ) , (5.35)
with single-particle operators

A = ∑ µ,ν A µν c † µ c ν and B = ∑ µ ′ ν ′ B µ ′ ν ′ c † µ ′ c ν ′ , Wick's theorem gives C AB (t, t ′ ) = ⟨A⟩ ⟨B⟩ + ∑ µν ∑ µ ′ ν ′ A µν B µ ′ ν ′ ⟨c † µ (t)c ν ′ (t ′ )⟩ ⟨c ν (t)c † µ ′ (t ′ )⟩ . ( 5.36)
If in second quantization, we associated creation operators c † µ to the states f µ , the expectation values can be calculated using

⟨c † µ (t)c ν ′ (t ′ )⟩ = e iE µ t-Γ µ t e -iE ν t ′ -Γ ν t ′ ⟨c † µ c ν ′ ⟩ = e iE µ (t-t ′ ) e -Γ µ (t-t ′ ) p µ δ ν ′ µ , ⟨c ν (t)c † µ ′ (t ′ )⟩ = e -iE ν (t-t ′ ) e -Γ ν (t-t ′ ) (1 -p ν )δ νµ ′ , (5.37)
resulting in

C AB (t, t ′ ) = ⟨A⟩ ⟨B⟩ + ∑ µν e i(E µ -E ν )(t-t ′ ) e -(Γ µ +Γ ν )(t-t ′ ) p µ (1 -p ν )A µν B νµ (5.38)
Under the Markovian assumption, the correlation functions can be written as

C AB (t, 0) = Tr(AV(t)[Bρ s ])
with V(t) = e Lt in the time-independent case. We will derive the master equation by demanding that the evolution of the density matrix due to the Lindbladian matches the result of Eq. (5.38). We will present a weak derivation, considering an infinitesimal evolution. It is, however, possible to derive the evolution superoperator directly. We will present the direct calculation of V(t) in the Appendix A.

We can use Wick's theorem to calculate the expectation value

c † µ c ν c † κ c λ = δ µν δ κλ p µ p κ + δ µλ δ νκ p µ (1 -p ν ).
(5.47)

The first term therefore cancels, and we get the expected result

i ⟨[H eff , A]B⟩ = i ∑ µν (E µ -E ν )p µ (1 -p ν )A µν B νµ .
(5.48)

Dissipator

For the dissipative part, we define the dissipator superoperator in the basis of the c, c † as follows

D[•] = ∑ λ α λ c λ • c † λ - 1 2 c † λ c λ , • + β λ c † λ • c λ - 1 2 c λ c † λ , • ≡ ∑ λ α λ L - λ + β λ L + λ [•].
(5.49)

where the jump operators L ± correspond to an incoherent creation and destruction of particles. With this notation, (5.52) which are solutions of the LS equation. In the absence of coupling between the dot and the reservoirs, the dressed state |Ψ(εi)⟩ reduces to a state with energy ε in reservoir i, that is |Ψ(εi)⟩ = |εi⟩ . The time-dependent Schrödinger equation gives two coupled equations

Tr(AD[Bρ s ]) = ∑ λ α λ Tr Ac λ Bρ s c † λ - 1 2 A c † λ c λ , Bρ s + ∑ λ β λ Tr Ac † λ Bρ s c λ - 1 2 A c λ c † λ , Bρ s = 1 2 ∑ λ α λ [c † λ , A] c λ B + c † λ [A, c λ ] B
ε + i ∂ ∂t -H r |Ψ r (εi, t)⟩ = ∑ j dε ′ ε ′ j v j (ε ′ ) Ψ s (εi, t) , ε + i ∂ ∂t -H s |Ψ s (εi, t)⟩ = ∑ j dε ′ v j (ε ′ ) ε ′ j εi + ε ′ j Ψ r (εi, t) .
(5.57)

Due to the time-periodicity we can write the dot Hamiltonian, as well as the states on the dot, using a Fourier decomposition Using the Fourier transform of Ψ s (εi, t), one finds the states on the reservoirs |Ψ r (εi, t)⟩ = ∑ j,m dε ′ ε ′ j v j (ε ′ ) χ m (εi) e -imω 0 t ε + mω 0ε ′ + iη , (5.59) which will give rise to a self-energy operator This can be interpreted as a matrix equation in an enlarged Hilbert space H s ⊗ T, where H s is the Hilbert space of the isolated dot and T is the vector space of the Fourier modes [START_REF] Sambe | Steady states and quasienergies of a quantum-mechanical system in an oscillating field[END_REF]. If an element of this Hilbert space can be written as The cost of extending the Hilbert space is that there is a redundancy in the solutions of Eq. (5.61), so that shifting the quasienergy by multiples of the frequency results in the same physical state. Within the narrow resonance approximation, we will see that this means that all injection energies in the reservoirs whose energies differ by mω 0 produce the same time-evolution of correlation functions (within a phase). This will be interpreted as multiple conducting channels giving rise to the same response of the small system. It is convenient to re-write Eq. (5.61) with the help of a linear operator acting in the extended Hilbert space H s ⊗ T : In practice, L F is calculated numerically, for example by truncation and diagonalization, or using a continued-fraction method if the drive is harmonic. 2 In the absence of coupling to the reservoirs, the self-energy is zero, and the Floquet spectrum consists of Wannier-Stark ladders of dim H s . If H s has a set of eigenenergies ε a , this means that the Floquet spectrum consists of the quasienergies ε a,m = ε amω 0 .

σ(ε) = ∑ j dε ′ v j (ε ′ ) v j (ε ′ ) ε -ε ′ + iη . ( 5
L F (ε) ≡ ∑ m (ε + mω 0 -σ(ε + mω 0 )) ⊗ P m -∑ n H s,n ⊗ T n (5.
The corresponding eigenvectors are denoted |ψ a,m ⟩ . In the presence of coupling to a continuum, the quasienergies evolve into resonances with a finite imaginary part γ a,m . The resolvent operator can be written similarly to the static case (5.14) using a biorthogonal basis Assuming the resonances are narrow, we expect that the equation ε = ε a,m (ε) will have solutions E a,m = E amω 0 which are close to the eigenenergies of the uncoupled system ε amω 0 . We can then approximate the resolvent by From (5.69) it also follows that if the vector of the enlarged Hilbert space |ψ a,0 (E a )⟩ corresponds to the vector |ψ a (t)⟩ in the physical space, then |ψ a,m (E a -mω 0 )⟩ = T m |ψ a,0 (E a )⟩ corresponds to another vector |ψ ′ a (t)⟩ = e -imω 0 t |ψ a (t)⟩ . These properties will be useful in the calculation of the reduced density matrix: In the second step we have used the resonance condition

G s (ε)
E a -E b = (m -n)ω 0 .
In analogy to the static case, we define the dressed basis of states

f µ (t) f µ (t) ≡ ∑ i,m πZ 2 µ Γ µ ϕ a,0 (E µ ) v i (E µ -mω 0 ) ⊗ |-m⟩ 2 ψ µ (t) ψ µ (t)
(5.72) so that ρ s (t) = ∑ µ p µ ψ µ (t) ψ µ (t) , (5.73) where the Floquet populations are p µ = ∑ i,m n i (E µ -mω 0 ) ϕ a,0 (E µ ) v i (E µ -mω 0 ) ⊗ |-m⟩ 2 ∑ i,m ϕ a,0 (E µ ) v i (E µ -mω 0 ) ⊗ |-m⟩ 2 .

(5.74)

The states in the Floquet case are therefore populated very differently with respect to the equilibrium case. Instead of a Fermi-Dirac distribution associated to each reservoir, one gets infinite copies of Fermi-Dirac distributions with a chemical potential which is shifted by multiples of the driving frequency. This suggests the point of view that each state associated to the same quasienergy E a is populated due to infinitely many channels in the reservoirs. As explained in both reference [START_REF] Seetharam | Controlled population of Floquet-Bloch states via coupling to Bose and Fermi baths[END_REF] and [START_REF] Iadecola | Occupation of topological Floquet bands in open systems[END_REF], one could take advantage of the structure of Eq.

(5.74) to selectively control the population of the Floquet states by engineering the reservoir density of states encoded in the function n i (E). This result agrees with similar expressions found in the literature [START_REF] Seetharam | Controlled population of Floquet-Bloch states via coupling to Bose and Fermi baths[END_REF][START_REF] Matsyshyn | Fermi-Dirac staircase occupation of Floquet bands and current rectification inside the optical gap of metals: an exact approach[END_REF][START_REF] Kohler | Driven quantum transport on the nanoscale[END_REF].

With the definition of the dressed basis, the two-times correlation function can be written as The correlation functions factorize into a part corresponding to micromotion within a driving period, governed by Floquet modes f µ (t) , and a slowly oscillating part (which therefore determines the long-time dynamics) that is governed by the lifetime of the quasiparticles. If the two scales are well-separated, one can derive a master equation similarly to the static case. The main difference will be that the master equation will have periodic coefficients since the operators c † µ corresponding to the states f µ (t) are time-periodic.

superconducting dot

Knowledge of the density matrix completely suffices to characterize the state of a system. In the presence of superconductivity, the density matrix must therefore include information about anomalous expectation values of the type ⟨d j d i ⟩ and ⟨d † j d † i ⟩ , where d, d † are operators on the subspace of the dot. A density ma- trix in Nambu space can be defined in analogy to the definition of the density matrix elements in Fock space ρ ij ≡ d † j d i , but generalized by replacing the bare operators d, d † by Nambu spinors Ψ ≡ d d † . We are interested in a singlelevel dot coupled to some superconducting leads i, biased with commensurate voltages s i = V i /V. We are therefore interested in calculating the reduced density matrix on the dot subspace (5.76) The creation and annihilation operators on the dot can be expressed as linear combinations of the dressed operators

d † σ (t) = ∑ α ∑ m∈Z e i(E α +mω 0 )t u * αm Γ † ασ -σe -i(E α +mω 0 )t v αm Γ α-σ d σ (t) = ∑ α ∑ m
e -i(E α +mω 0 )t u αm Γ ασσe i(E α +mω 0 )t v * αm Γ † α-σ .

(5.77)

Here the index α is shorthand for the reservoir i with momentum k and spin σ, α = {ikσ}. Since the dressed annihilation operator destroys the steady-state we only have to compute terms of the form ΓΓ † . For example, ⟨d σ (t)d † σ ′ (t ′ )⟩ = ∑ αβ ∑ mn e -i(E α +mω 0 )t e i(E β +nω 0 )t ′ u αm u * βn ⟨Γ ασ Γ † βσ ′ ⟩ = ∑ α ∑ mn e -iE α (t-t ′ ) e -i(mω 0 t-nω 0 t ′ ) u αm u * αn δ σσ ′ .

(5.78)

Similarly, we find:

⟨d † σ (t)d σ ′ (t ′ )⟩ = ∑ α
∑ mn e -iE α (t-t ′ ) e -i(mω 0 t-nω 0 t ′ ) v αm v * αn δ σσ ′ ⟨d σ (t)d σ ′ (t ′ )⟩ = σ ∑ α ∑ mn e -iE α (t-t ′ ) e -i(mω 0 t-nω 0 t ′ ) u * αm v αn δ σ,-σ ′ ⟨d † σ (t)d † σ ′ (t ′ )⟩ = -σ ∑ α ∑ mn e -iE α (t-t ′ ) e -i(mω 0 t-nω 0 t ′ ) v αm u * αn δ -σ,σ ′

(5.79)

For the calculation of the density matrix we set t = t ′ , so this reduces to

ρ σ (t) = ∑ α ∑ m,n
e -i(m-n)ω 0 t u αm σv αm ⊗ u αn σv αn *

(5.80)

The coefficients u, v can be calculated if the resolvent operator is known by multiplying the resolvent with the source term (see Ch. x α e iϕ α /2 -y α e -iϕ α /2

(5.81)

The 2 × 2 matrix resulting from the product u αm v αm ⊗ u α,n v α,n * can then be written as a function of the resolvent elements and the reservoir populations: e -i(m-n)ω 0 t (RC -s i ) m,0 Q i (ε)(RC -s i ) † n,0 (5.84) This is a generalization of Eq. (5.28) to the case of a multiterminal superconducting dot. Using the narrow resonance approximation, we expect that the structure of the obtained reduced density matrix will be similar with the previous section, since for the derivation of Eq. (5.71) we assumed a general Floquet Hamiltonian without specifying its structure. The new element due to superconductivity is the extra structure in Nambu space.

∑ i,k u m (i, k) v m (i, k) ⊗ u n (i, k) v n (i, k) * = = ∑ i,k J 2 i R

Floquet-Andreev populations

In the narrow resonance limit, the resolvent is sharply peaked near energies ε = E a,p = E a + pω 0 , where a = ± is the label of the two Wannier-Stark ladders that appear in the superconducting case within one period. Near those energies, we can make the approximation R(ε) mn ε→E a,p ≃ Ψ a,m (E a,p ) Φ a,n (E a,p )

εpω 0 -E a + iΓ a .

(5.85)

As in the normal case, Ψ is a right eigenvector of the R -1 operator, while Φ is the corresponding left eigenvector. However, since we are working in Nambu space, for any given index m, |Ψ a,m ⟩ is a two-component column vector, while ⟨Φ a,m | is a row vector, so that |Ψ a,m ⟩ ⟨Φ a,m | is a 2 × 2 matrix. We therefore find ρ(t) ≃ ∑ a=± ∑ i,p π Γ a Φ a,0 (E a,p ) C -s i Q i (E a + pω 0 )C s i Φ a,0 (E a,p ) |Ψ a (t)⟩⟨Ψ a (t)|

(5.86) Compared to the case of normal reservoirs, Q i plays the role of the reservoir distribution function n i . The expression is slightly more complicated, because the action of the operator C is conditioned on the Nambu index, so that if we write the two component vector Φ = u v , then

C s i Φ † a,0 (E a,p ) = u * a,-s i (E a,p ) v * a,s i (E a,p ) = 1 N ⊗ T -p u * a,-s i +p (E a ) v * a,s i +p (E a )
(5.87)

The expectation value in Eq. (5.86) can be rewritten as Φ a,0 (E a,p )C -s i Q i C s i Φ a,0 (E a,p ) † = u a,-s i +p (E a ) v a,s i +p (E a )

T Q i u * a,-s i +p (E a ) v * a,s i +p (E a )

. (5.88) This quantity can be identified with the occupation of the Floquet-Andreev state a imposed by the reservoir labeled i. Indeed, if we write the reduced density matrix in a dressed basis of states as before ρ(t) = ∑ a p a |F a (t)⟩⟨F a (t)| then the Floquet-Andreev populations are given by p a ≡ ∑ i,p Φ a,0 (E a,p ) C -s i Q i (E a + pω 0 )C s i Φ a,0 (E a,p ) ∑ i,p Φ a,0 (E a,p ) Φ a,0 (E a,p ) (5.89)

Comparing to Eq. (5.74), one sees that the addition of superconductivity results in the Fermi-Dirac distribution being replaced by the matrix Q, defined as

Q i (E) = 2πρ 0 J 2 i √ E 2 -∆ 2 E -∆e iϕ i -∆e -iϕ i E θ(E -∆), (5.90)
so that the population of the Floquet-Andreev state with energy E a depends on an ensemble of step functions θ(E a + pω 0 -∆), representing copies of the superconducting reservoirs with a gap value that appears shifted by multiples of the voltage bias.

conclusions

We have presented work towards a method for deriving Markovian master equations, where the long-time evolution of correlation functions is controlled by the energies and linewidths of the resonances of the system's resolvent. The method relies on the assumption that resonances are narrow and wellseparated, so that for any two resonances the condition |E 1 -E 2 | ≫ Γ 1,2 holds, and on the use of Wick's theorem for the calculation of correlation functions. For a dot coupled to superconducting reservoirs, we have derived a density matrix description in the basis of the Floquet-Andreev resonances, and given an expression for the calculation of the populations associated to the resonances. The derivation of the master equation in the superconducting case is left for future work.

The main theme of this thesis is to study the Floquet spectrum of voltagebiased superconducting junctions. We focused on two particular systems, the S-QD-S junction and the S-QD-S-QD-S bijunction. In Chapter 3, we presented a method that builds up on previous work [START_REF] Mélin | Simple Floquet-Wannier-Stark-Andreev viewpoint and emergence of low-energy scales in a voltage-biased three-terminal Josephson junction[END_REF][START_REF] Mélin | Engineering the Floquet spectrum of superconducting multiterminal quantum dots[END_REF][START_REF] Douçot | Berry phase in superconducting multiterminal quantum dots[END_REF] that permits to extract observables from the resolvent of the particular system. We presented an efficient algorithm for the calculation of the resolvent elements based on a continued fraction method. In particular, we argue that the Floquet spectrum could be probed by local tunneling spectroscopy. Such an experiment has yet to be performed, but it would allow the exploration of Floquet physics with simple application of constant voltage bias. Contrary to closed Floquet systems that often suffer from thermalization, the voltage-biased superconducting junction should reach a Floquet steady state because a) the system is open so it does not have to thermalize [START_REF] Mori | Floquet states in open quantum systems[END_REF], and b) the existence of the superconducting gap offers protection to the resonances. The use of a quantum dot is therefore critical, because its states are typically detached from the superconducting continua. To the best of my knowledge, there is only one experiment which has measured the Floquet spectrum in Josephson junctions, under microwave irradiation [START_REF] Park | Steady Floquet-Andreev states in graphene Josephson junctions[END_REF], although the interpretation of this experiment has been recently questioned [START_REF] Haxell | Microwave-induced conductance replicas in hybrid Josephson junctions without Floquet-Andreev states[END_REF]. We studied the bijunction in two regimes depending on the distance R between the QDs. The first regime R ≲ ξ 0 results in hybridization of the states on the dots. In equilibrium, this hybridization results in the formation of a macroscopic Andreev molecule and in a non-local Josephson effect [START_REF] Pillet | Nonlocal Josephson effect in Andreev molecules[END_REF]. We studied the Andreev molecule in the presence of commensurate voltage bias, and shown how the Floquet spectrum and the MAR current are modified due to the hybridization.

In Chapter 4 we studied the bijunction in the opposite regime R ≫ ξ 0 . This regime was first studied in Ref. [START_REF] Mélin | Ultralong-distance quantum correlations in three-terminal Josephson junctions[END_REF], proposing that a long-range coupling will develop between the dots. We showed that the system in this regime is an interferometer and, as a result, the MAR current oscillates as a function of the applied voltage. We focused on the influence of the Floquet-Tomasch effect on the spectrum of the bijunction at large R, and showed that, in this limit, the system is effectively described by two resonances coupled indirectly through a continuum of states. Moreover, we found that the amplitude of the Floquet-Tomasch oscillations of the spectral function depends on the quartet phase and could potentially be controlled non-locally by changing the phase drop on one of the dots. A different point of view is taken in Chapter 5, where we present (preliminary) work on the derivation of a master equation in the Floquet basis of states.

The method relies on the assumption that if the coupling to the reservoirs is weak, the Floquet-Andreev resonances will be sharp, so that associated quasiparticles are well-defined. The lifetime of these quasiparticles would then control the long-time dynamics. We apply the method first to a dot coupled to normal leads and then to a dot with a time-periodic Hamiltonian. In the case of a superconducting dot, we derived the populations of the Floquet-Andreev states, but the derivation of the master equation is left for future work.

An immediate extension of the work presented in the last chapter is to treat the case of degeneracies or near degeneracies in the Floquet spectrum. In particular, the method we presented breaks down near avoided crossings, so this case needs to be considered separately. An interesting direction for future work would be to consider whether a Floquet-Andreev qubit is possible. Since the Andreev level qubit has already been proposed and experimentally realized, it would be interesting to see if a periodically driven qubit can be realized and manipulated with the use of a second drive. In fact, it has been proposed that Floquet qubits present certain advantages with respect to their equilibrium counterparts: for example, a Floquet qubit would offer multiple dynamical "sweet spots" which can be chosen in situ by changing the drive frequency [START_REF] Huang | Engineering dynamical sweet spots to protect qubits from 1/ f noise[END_REF][START_REF] Gandon | Engineering, control, and longitudinal readout of Floquet qubits[END_REF][START_REF] Mundada | Floquetengineered enhancement of coherence times in a driven fluxonium qubit[END_REF].

A third point is to consider a more realistic description of the S-QD-S junction. Electrons in a quantum dot are confined in a small space, therefore making the effect of Coulomb interactions important. A starting point would be to work in the Fock space of the dot, where one can derive an effective four-level system in the absence of voltage bias [START_REF] Meng | Self-consistent description of Andreev bound states in Josephson quantum dot devices[END_REF]. Assuming that, upon turning on the voltage, the four levels evolve into four Floquet-Andreev resonances, a possible way to move forward would be to try to derive a master equation in the Floquet basis of this Fock space [START_REF] Avriller | Andreev bound-state dynamics in quantumdot Josephson junctions: a washing out of the 0-π transition[END_REF][START_REF] Damanet | Controlling quantum transport via dissipation engineering[END_REF].

Finally, incommensurate voltages could be considered. This could be done as a perturbation on the commensurate case, or by mapping the system to a Floquet tight-binding lattice of higher dimensions [START_REF] Martin | Topological frequency conversion in strongly driven quantum systems[END_REF]. It is, in fact, possible to derive the superoperators L and V(t, t ′ ) directly, and not for infinitesimal times as was done in 5.3.3. Given the resolvent, we can associate creation and annihilation operators c λ , c † λ to its poles at energies E λ + iΓ λ . We expect a time-independent Liouville operator L, resulting in an exponential form for the evolution superoperator V(t, t ′ ) = e (t-t ′ )L .

(A.1)

We write L = -iad H + D, with

H = ∑ λ E λ c † λ c λ , D[•] = ∑ λ α λ c λ • c † λ - 1 2 c † λ c λ , • + β λ c † λ • c λ - 1 2 c λ c † λ , • ≡ ∑ λ α λ L - λ + β λ L + λ [•].

5 )

 5 Using the Fourier transform, the time-dependent Schrödinger equation is mapped to an eigenvalue problem

Figure 1 . 1 :

 11 Figure 1.1: (a) A conventional Josephson junction. (b) A multiterminal Josephson junction.

Figure 1 . 2 :

 12 Figure 1.2: Andreev bound state formation. (a) Transfer of a Cooper pair through an S-N-S junction due to successive Andreev reflections. (b) The density of states in the normal region contains a pair of Andreev bound states within the superconducting gap.

Figure 1 . 3 :

 13 Figure 1.3: (a) Sketch of a multiple Andreev reflection of order n = 3. In this MAR ladder picture, the leads are placed at the same chemical potential and quasiparticles "climb" the ladder with each Andreev reflection. (b) Current-voltage characteristics in the S-QD-S case.

Figure 1 . 4 :

 14 Figure 1.4: Sketch of the Andreev molecule.

  It is therefore of interest to understand how two Josephson junctions might interact when they are closely spaced. In an analogy to the formation of a molecule, bringing two ABS carrying junctions close enough should result in a hybridization of the ABS wave-functions, thus creating an "Andreev molecule". The hybridization results in non-local effects in the Josephson current, whereby changing the phase at one junction would change the current of Cooper pairs flowing through the other[START_REF] Pillet | Nonlocal Josephson effect in Andreev molecules[END_REF]

Figure 2 . 2 :

 22 Figure 2.2: Sketch of the Andreev molecule.

2 Figure 2 . 3 :

 223 Figure 2.3: (a) Single junction regime (R > ξ 0 ). The second junction (red lines) is uncoupled and thus not affected by the phase difference of the first junction. (b) Andreev molecule regime (R ≲ ξ 0 ). The states on the first dot (black lines) and the states on the second dot (red lines) are hybridized. All couplings are set to Γ = ∆. ϕ b = 3π/5 and k F R = π/4.

Figure 2 . 4 :

 24 Figure 2.4: Andreev molecule spectra. Hybridization of the ABS of each junction produces avoided crossings in the spectrum. Parameters used for the calculation: R = ξ 0 /2, Γ = ∆, k F R = π/4.

Figure 3 . 1 :

 31 Figure 3.1: Diagrammatic representation of the forward scattering self-energy.Σ + m resums loops to the right of site m.

59 )

 59 It is convenient to introduce a translation operator C which acts on the resolvent by translating the Floquet indices by an integer number ±s, and is conditioned on the Nambu indices C s ≡ ∑ m |i, e, m + s⟩⟨i, e, m| + |i, h, m -s⟩⟨i, h, m| . (3.60)

Figure 3 . 2 :

 32 Figure 3.2: (a) Density plot of the spectral function of the dot in the S-QD-S system, biased with V L = -V R . The couplings to the reservoirs are equal Γ L = Γ R = 0.3 ∆. This configuration produces uncoupled FWS ladders. (b) Cuts at various voltages. For equal tunnel couplings the peaks are always at even multiples of ω 0 . (c) The more generic case of unequal tunnel couplings produces avoided crossings between ladders. (d) For Γ L ̸ = Γ R , the electron-like part of the spectral function (solid lines) and the hole-like part (dashed lines) have peaks at different energies 2mω 0 ± ϵ. The couplings are Γ L = 0.3 ∆, Γ R = 0.4 ∆.

Figure 3 . 3 :

 33 Figure 3.3: (a) Spectral function of the S-QD-S system biased with V L = -V, V R = 0. (b) Particle-hole asymmetry of the spectral function A(E) -A(-E). Red represents positive values of the difference (electron-like) and blue represents negative values (hole-like). The signs change at avoided crossings. The couplings are Γ L = Γ R = 0.3 ∆.

Figure 3 . 4 :

 34 Figure 3.4: Resonant structure of the non-diagonal resolvent elements R m,-1 at energies above the gap and for Γ L = Γ R = 0.3 ∆. Resonances appear at odd multiples of ω 0 . There is a hierarchy of peaks which depends on the voltage, and the resonance of an element R m,-1 dominates when m is the minimum integer that satisfies |m|ω 0 > ∆.

Figure 3 . 5 :

 35 Figure 3.5: Subgap current of a single resonant dot. (a) Subgap structure for the highly symmetric left-right configuration Γ = Γ L = Γ R and (V L = -V, V R = +V) with the known MAR steps at odd subdivisions of the gap. Logarithmic scaling has been used for better visibility of the features. (b) Comparison of subgap current for equal (dashed black line) and unequal (red and orange lines) tunnel couplings (Γ L ̸ = Γ R ), showing the modification of the MAR steps when there is asymmetry in the tunnel couplings.

Figure 3 . 6 :

 36 Figure 3.6: Localization at strong drives. The amount of harmonics needed increases with decreasing voltage. Γ L = Γ R = 0.3 ∆.

  above expression makes more explicit a diagrammatic view of the MAR current, where the major contribution is given by the shortest path between the source site at the left reservoir s L = -1 and the measure site at m + 1 through |m + 2|/2 electron-hole conversions. The resonances of R 11,22 m+1,-1 2 now explain the relationship between the Floquet ladders and the current. For a symmetric junction and Γ L = Γ R all resonances are at some odd multiple of the frequency ω = |2p + 1|ω 0 , as shown on Fig. 3.4. All the terms in (3.80) consequently cre- ate MAR steps at odd subdivisions of the gap. For Γ L ̸ = Γ R , however, the terms R 11 m+1,-1 and R 22 m+1,-1 give peaks that are a distance 2ε apart. The MAR steps will then be split and the distance of the steps ought to increase with an increase in the difference |Γ L -Γ R |. Breaking of the mirror symmetry therefore produces the modification of the MAR steps that are shown in Fig. (3.5b).

. 93 )Figure 3 . 8 :

 9338 Figure 3.8: (a) Density plot of the spectral function of the first dot. The distance between the dots is taken to be R = ξ 0 /2 and all couplings are Γ = 0.3 ∆.

Figure 3 . 9 :

 39 Figure 3.9: Density plots of A 1 (E) -A 1 (-E) as a function of the distance between the two dots, for different values of voltage bias. At R ≫ ξ 0 the quasienergies approach those of a single junction. At R ∼ ξ 0 there is a lift of degeneracy. Red represents positive values of the difference (electron-like) and blue represents negative values (holelike). We set k F R = π/4 and all couplings Γ = 0.3 ∆.

Figure 3 . 10 :

 310 Figure 3.10: Effect of detuning and robustness at strong driving. The distance between the dots is taken to be R = ξ 0 /2 and all couplings are Γ = 0.3 ∆.

Figure 3 . 11 :Figure 3 . 12 :

 311312 Figure 3.11: Andreev molecule I-V curve calculated at various distances between the two dots. Parameters used: Γ = 0.3 ∆, k F R = π/4.

Figure 3 . 13 :

 313 Figure 3.13: "Localization" on the Floquet chain means only a few of the harmonics need to be taken into account when calculating the current. The amount of harmonics required increases with decreasing voltage. Parameters: R = ξ 0 , Γ = 0.3 ∆, k F R = π/4.

Figure 4 . 1 :

 41 Figure 4.1: Long-range Floquet-Tomasch effect on the subgap structure. Zoom near the second MAR step, around ω 0 = 2∆ 5 , compared to the single junction case (blue line). At large distances oscillations of the I-V curves appear around the single junction current. Coupling strength Γ = ∆/2.

Figure 4 . 2 :

 42 Figure 4.2: Spectral function at energies above the superconducting gap, at ω 0 = 0.45 ∆. Coupling strength Γ = ∆/2.

Figure 4 . 3 :

 43 Figure 4.3: Single junction spectrum. Multiple Andreev reflections turn the initial ABS on the dot into resonances with a finite width. The lower panel is a contour plot of log det L (1) 00 showing the zeroes of det L (1) 00 in the complex plane. The upper panel shows the corresponding spectral function. Drive frequency is ω 0 = ∆/2 and Γ a = Γ c1 = ∆/2.

Figure 4 . 4 :

 44 Figure 4.4: Floquet-Tomasch effect on the bijunction spectrum. The interdot distance is R = 50ξ 0 . Contour plot of log(|det L 00 |) in the com- plex plane showing the zeroes of det L 00 (lower panel) and the corresponding spectral function of dot 1, -2 π Im R e 1 e 1 (upper panel). Floquet-Tomasch processes produce oscillations in the DOS. All couplings are set to Γ = ∆/2 and the frequency of the drive is ω 0 = ∆/2.

4 . 5 5 Figure 4 . 5 :

 45545 Figure 4.5: Evolution of the zeroes of det L 00 for different couplings to the reservoirs when the interdot distance is increased. The number of poles increases with the distance. The frequency of the drive is set to ω 0 = ∆/2.

Figure 4 . 7 :

 47 Figure 4.7: Contours of integration for the evolution operator. (a) A simple pole leads to exponential decay. (b) Interference of multiple poles can modify the Fermi Golden Rule.

Figure 4 . 8 :

 48 Figure 4.8: Solutions ω + of Eq. (4.31) in the complex plane. Dashed blue lines are the solutions of the imaginary part of Eq. (4.31) and black lines are the solutions of the real part.The red dots at their intersection are the solutions ω + . We take ϵ 0 = 0.5, γ = 0.01. As the distance between the dots grows, poles proliferate and approach the real axis.

4 -Figure 4 . 9 :

 449 Figure 4.9: Plot of the solutions of Eq. (4.32) as a function of the distance Lbetween the two dots. At periodic intervals L = nπv ϵ 0 there are nonavoided crossings; one of the roots becomes purely real (bound state with no dissipation), while the second state has a maximum dissipation Im ω = 2γ. Parameters γ = 0.01, ϵ 0 = 0.5.

Figure 4 . 11 :

 411 Figure 4.11: In the strong coupling regime γL/v > 1 a particle oscillates between the two dots before dissipating.

Figure 4 . 13 :

 413 Figure 4.13: In the weak coupling limit γL/v ≪ 1, the physics is that of overlapping Fano resonances. (a) Identical dots with no detuning (black solid line). Transmission for a single dot (dashed line) (b) Small detuning δ = 0.1 (black line) and strong detuning δ = 1 (grey line).

Figure 4 . 14 :

 414 Figure 4.14: Strong coupling limit γL/v > 1. Transmission coefficient for γL/v = 2 (grey lines), and γL/v = 5 (black lines). (a) In the absence of detuning. (b) Detuning δL/v = 2. As the coupling increases, oscillations become more prominent. In both plots the functions are calculated in reduced units of ω/ω L .

1 Figure 4 . 15 :

 1415 Figure 4.15: Phase diagram.

Figure 4 . 16 :

 416 Figure 4.16: Contour plots of the solutions of Eq. (4.74) in the complex plane.The solutions are found at the intersections of Re(det L eff ) = 0 (black lines) and Im(det L eff ) = 0 (blue dashed lines). The roots are marked with red dots. Choice of parameters: Γ j = 0.2 ∆, ω 0 = 0.5 ∆.

]Figure 4 . 17 :

 417 Figure4.17: (a) Variation of the spectral function Im R-R(1)

  and no orthogonality relation exists in general between the left and right eigenvectors themselves ⟨ψ a |ψ b ⟩ ̸ = δ ab , ⟨ϕ a |ϕ b ⟩ ̸ = δ ab . One can, however, impose a normalization condition ⟨ψ a |ψ a ⟩ = 1, or ⟨ϕ a |ϕ a ⟩ = 1, but not both at the same time. The completeness of the biorthonormal basis is expressed by 1 = ∑ a |ψ a ⟩⟨ϕ a | .

Figure 5 . 1 :

 51 Figure 5.1: Multilevel dot.

1 2 1 2 ∑

 112 c λ , A] c † λ B + c λ [A, c † λ ] B(5.50)Using the commutatorsc † λ , c † µ c ν = -c † µ δ λν , and c λ , c † µ c ν = c ν δ λµ ,(5.51)we getTr(AD[Bρ s ]) = -1 2 ∑ µν (α µ + α ν )A µν c † µ c ν B + 1 2 ∑ µν (β µ + β ν )A µν c ν c † µ B + α ν )A µν δ µν p µ ⟨B⟩ + p µ (1p ν )B νµ + ∑ µν (β µ + β ν )A µν δ µν (1p µ ) ⟨B⟩p µ (1p ν )B νµ = ∑ µ β µ -(α µ + β µ )p µ A µµ ⟨B⟩ µν α µ + β µ + α ν + β ν p µ (1p ν )A µν B νµ .

H

  s (t) = ∑ m e -imω 0 t H s,m |Ψ s (εi, t)⟩ = ∑ m e -imω 0 t |χ m (εi)⟩ .(5.58)

. 60 )

 60 The states on the dot are then determined through an infinite set of equations(ε + mω 0σ(ε + mω 0 )) |χ m ⟩ -∑ n H s,n |χ m-n ⟩ = |v i (ε)⟩ δ m0 .(5.61)

  |ψ a ⟩ = ∑ m |χ a,m ⟩ ⊗ |m⟩ ,(5.62)then the corresponding vector in the physical space is|ψ a (t)⟩ = ∑ m e -imω 0 t |χ a,m ⟩(5.63)

64 )

 64 where the translation operator T acts on the Floquet states by translating them in Fourier space T n |m⟩ = |m + n⟩ . Moreover, P m |m ′ ⟩ = δ mm ′ |m⟩ is the projec- tor on the m Fourier mode. With this definition, inverting L F gives a generalization of Eq. (5.24) to the time-periodic case|Ψ s (εi)⟩ = L -1 F (ε) |v i (ε)⟩ ⊗ |m = 0⟩ .(5.65)

  m (ε)⟩⟨ϕ a,m (ε)| εε a,m (ε) + iγ a,m (ε) . (5.66)

ε∼E a, m ≃

 m Z a |ψ a,m (E a,m )⟩⟨ϕ a,m (E a,m )| ε -E a,m + iΓ a . (5.67)The linear operator L F has the covariance propertyL F (ε)T n = T n L F (ε + nω 0 ). (5.68)Due to this translation property we can show thatT n |ψ a,m (ε)⟩ = |ψ a,m+n (εnω 0 )⟩ ⟨ϕ a,m (ε)| T -n = ⟨ϕ a,m+n (εnω 0 )| (5.69)

  ρ s (t) = ∑ i dε n i (ε) |Ψ s (εi, t)⟩⟨Ψ s (εi, t)| .(5.70)Within the narrow resonance approximation the reduced density matrix isρ s (t) ≃ ∑ i dε n i (ε) ∑ |ψ a,m (E a,m )⟩⟨ϕ a,m (E a,m )| (|v i (ε)⟩ ⊗ |0⟩) ε -E a + mω 0 + iΓ a × (⟨0| ⊗ ⟨v i (ε)|) |ϕ b,n (E b,n )⟩⟨ψ b,n (E b,n )| ε -E b + nω 0 -iΓ b πn i (E a -mω 0 )Z 2 a Γ a ⟨ϕ a,0 (E a )| (|v i (E a -mω 0 )⟩ ⊗ |-m⟩) 2|ψ a (t)⟩⟨ψ a (t)|(5.71)

  E C T D E R I VAT I O N O F T H E T I M E -E V O L U T I O N S U P E R O P E R AT O R

(A. 2 )

 2 In order to have a compact notation, we have introduced the superoperator adH (•) = [H, •].We want to calculate V(t)[Bρ] = e Lt [Bρ]. It is convenient to introduce a super- operator M B which acting on the operator ρ multiplies it on the left by the operator B, i.e. M B [ρ] ≡ Bρ. In the stationary state, we expect that Lρ = 0. We can take advantage of this fact to equate L(Bρ) = LM B (ρ) = [L, M B ](ρ). We are interested in the case of the operator B being any single-particle operator acting on the subsystem statesB = ∑ µν B νµ c † ν c µ . (A.3)For the non-dissipative part of the Liouvillian it is straightforward to calculate-i[ad H , M B ]ρ = -i[H, Bρ] + iB[H, ρ] = -i[H, B]ρ. (A.4) Explicit calculation of the commutator c † λ c λ , c † ν c µ = c † ν c µ (δ λνδ λµ ), gives -i[H, B]ρ = ∑ µν i(E µ -E ν )B νµ c † ν c µ ρ. (A.5)

  The many-body wavefunction should describe a coherent state created by adding pairs of electrons with opposite momenta and spins. A Cooper pair is therefore created by the pair operator c † -k↓ c † k↑ , where the "bare" operators c † kσ create an electron with momentum k and spin σ (σ can be either spin-up ↑ or spin-down ↓), and their hermitian conjugates c kσ annihilate the vaccum state c kσ |0⟩ = 0. A fitting wavefunction for the superconducting ground state is then 2

  This relation reduces to the usual Heisenberg equation of motion when the Hamiltonian H and the operator A are time-independent

	i	d dt	A H (t) = A H , H .	(2.22)

  Figure 4.10: Survival probabilities P 1,2 at a non-avoided crossing (black and dashed black line), showing that a particle initially on the first dot does not completely dissipate into the continuum. At anticrossings k 0

	1.0	P 1 (t) at crossing
	0.8	P 2 (t) at crossing
	0.6	P 1 (t) at anti-crossing
		P 2 (t) at anti-crossing
	0.4	Fermi's Golden Rule
	0.2
	0.0	1.0	2.0	3.0
		t/τ
				2	=	1 4	.

  C AB (t, t ′ ) = ⟨A(t)⟩ ⟨B(t ′ )⟩ i(E µ -E ν )(t-t ′ ) e -(Γ µ +Γ ν ))(t-t ′ ) p µ (1p ν )× × f µ (t) A f ν (t) f ν (t ′ ) B f µ (t ′ )

	+ ∑ µ,ν	(5.75)

e

  We can write this more compactly with the use of a translation operator in Fourier space, conditioned by the Nambu components(RC -s i ) mn = 1 N ⊗ ⟨m| R |ns i ⟩ ⟨n| 0 0 |n + s i ⟩ ⟨n| |n⟩ ⊗ 1 N .(5.83)With this definition, the density matrix takes the form ρ(t) = ∑

		11 m,-s i R 12 m,s i R 21 m,-s i R 22 m,s i	x 2 ik -x ik y ik e -iϕ i	-x ik y ik e iϕ i y 2 ik	R 11 n,-s i R 12 n,s i R 21 n,-s i R 22 n,s i	†	(5.82)
	= ∑ i	dε	R 11 m,-s i R 12 m,s i R 21 m,-s i R 22 m,s i	Q i (ε)	R 11 R 21 n,-s i R 22 n,s i n,s i n,-s i R 12	†

i dε ∑ m,n

Floquet-Andreev resonances in multiterminal quantum dots 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Helium has a boiling point of 4.2K and Onnes managed to get to a record temperature of about 1.5K

There is freedom in the definition since one can first add the k ↑ electron and then its pair -k ↓, or vice versa. Indeed, in the original BCS paper the pair creation operator is defined using the ordering c † k↑ c † -k↓ , while we have used the opposite order c † -k↓ c † k↑ . One simply has to substitute v k → -v k in order to retrieve the wavefunction in[START_REF] Bardeen | Theory of superconductivity[END_REF].

Another way to derive Eq. (2.27) is to start from the condition that the quasiparticle annihilation operator destroys the ground stateγ(0) |ψ(0)⟩ = 0 ⇒ U(t)γ(0) |ψ(0)⟩ = 0 ⇒ U(t)γ(0)U † (t) |ψ(t)⟩ = 0 ⇒ γ(t) = U(t)γ(0)U † (t).

Γ ji ≡ πρ 0 J 2 ji , Γ c ≡ √ Γ c1 Γ c2 = πρ 0 J c1 J c2

It was the development in semiconducting superlattices that made the observation possible. Not much later, it also became possible to observe the Bloch oscillations in cold atom systems. There, the crystal lattice is replaced by an optical lattice, making it possible to increase relaxation time significantly.

Note that in all cases we have set the dot energy to zero, ε d = 0. This can be rectified simply by adding a (-ε d σ z ) term to the M 0 m matrices.

The total self-energy matrix Σ = Σ + + Σ -is actually diagonal in the asymmetric case as well.

In order for the eigenvalues of the density matrix to be interpreted as probabilities, ρ must fulfill the following two conditions

If the drive is harmonic, for example H s (t) = 2ε d cos(ω 0 t) |d⟩⟨d| , there are only two modes contributing to the Fourier decomposition: H s,n = ε d (δ n,+1 + δ n,-1 ), so that m couples only to m ± 1. The operator L F then becomes tridiagonal and a continued-fraction method gives the resolvent elements, as described in Ch.3.

master equations for non-interacting fermionic systems

Density matrix theory and master equations

When dealing with a closed system, quantum mechanics tells us that the state of the system is fully described by a vector |ψ(t)⟩ belonging to a Hilbert space. The dynamics of the system are then described by the Schrödinger equation. Equivalently, we can represent this evolution by an operator that propagates the state vector from an initial time t ′ to a final time t : |ψ(t)⟩ = U(t, t ′ ) ψ(t ′ ) .

(5.1)

Inserting into the Schrödinger equation we get an equivalent equation for the evolution operator

The evolution operator is unitary

Moreover, if the Hamiltonian is time-independent we have the well-known expression U(t, t ′ ) = e -iH(t-t ′ ) .

(5.3)

In most situations, the system is not in a pure state, but in a mixed state, meaning that there is a set of states {|ψ i ⟩} that the system can be in with some associated probability p i . We can then describe the system by a density matrix ρ ≡ ∑ i p i |ψ i ⟩⟨ψ i | . 1 Once the density matrix ρ is defined, the expectation value of any observable is given by the relation ⟨A⟩ = Tr(Aρ).

(5.4)

Since the system still evolves according to the Schrödinger equation, it follows that the equation of motion for the density matrix is given by the von Neumann equation ρ(t) = -i[H(t), ρ(t)] (5.5) which can equivalently be written as

In analogy to the classical Liouville equation, we can rewrite the von Neumann equation as

where we have defined the Liouville superoperator (an operator that acts on operators) as

master equations for non-interacting fermionic systems

(Weak) derivation of the master equation

Starting from Eq. (5.38), we take the time derivative:

(5.39)

We then expect that the above expression is equal to

and the dissipative part

(5.43)

Non-dissipative part

We can check that the non-dissipative part is reproduced by a Hamiltonian which is diagonal in the basis of the c † , c operators:

(5.44) Indeed, using the cyclic property of the trace,

and using the commutator

(5.46)

master equations for non-interacting fermionic systems

Setting α µ + β µ = 2Γ µ , the second term gives the expected dissipative part. The first part vanishes by imposing that populations satisfy the relation

We then find

(5.54)

We therefore have a master equation in the basis of the creation and annihilation operators corresponding to the poles of the resolvent:

(5.55)

driven dot

We now consider a dot which is driven so that its Hamiltonian H s is a periodic function of time with the frequency ω 0 = 2π T , i.e., H s (t) = H s (t + T). In order to simplify the discussion, we will consider that the Hamiltonian of the reservoirs H r and the coupling between dot and reservoirs V are time-independent and given by Eq. ( 5.17) and (5.18), respectively. We will use Floquet theory to show that, at least within the narrow resonance approximation, it is possible to derive a master equation in the basis of the Floquet states on the dot.

We start by constructing the dressed states

For the dissipative part, for any superoperator of the form

We therefore have

for α, β given by Eq. (5.54), we find

and the evolution superoperator follows with an exponentiation

It remains to study the special case of µ = ν for which

We can postulate V(t)(c † µ c µ ρ) ≡ g(t)ρ + f (t)c † µ c µ ρ. Using the property d dt V(t) = LV(t) we get ġ(t)ρ + ḟ (t)c † µ c µ ρ = L(g(t)ρ + f (t)c † µ c µ ρ) = f (t)L(c † µ c µ ρ) (A.12) from which we get ḟ (t) = -2Γ µ f (t) ⇒ f (t) = c 1 e -2Γ µ t ġ(t) = 2Γ µ p µ f (t) ⇒ g(t) = c 2p µ e -2Γ µ t (A.13) Since V(0) = 1 we get the conditions f (0) = 1 and g(0) = 0, which in turn give c 1 = 1, c 2 = p µ . We therefore have