Superconducting nanoscopic systems represent one of the most interesting classes of systems   . This includes specifically the possibility to study and manipulate entangled quantum systems    as well as the possibility to generate even higher order correlations (i.e. multiple Andreev reflections)    . Andreev reflection refers to the retroreflection of a hole leaving a Cooper pair behind. Such processes may happen multiple times if the two leads contacted are both superconductors.
The properties of superconductor-quantum dot-superconductor (S-QD-S) systems are particularly relevant and rich due to the manifold energy scales being present in the system (superconductor gap , phonon frequency , Coulomb interaction U, tunneling rate as already analyzed in previous works   which requires to take simplications to get situations close to experiment without overcomplifications  ).
The highly nonlinear transport characteristics of the S-QD-S setup we wish to analyze have been studied in the non-interacting setup  and noise properties have been analyzed in non-equilibrium in less complicated systems  , but a full analysis of the system in question including all aspects of the current flow is missing so far.
2. Description of the System
The quantity of our primary interest is the full counting statistics (FCS) in terms of the cumulant generating function (CGF) χ. It represents a very convenient tool for the calculation of a variety of transport properties. It is directly related to the probability distribution function P(Q) to transfer Q elementary charges during a fixed very long measurement time τ. By simple derivation with respect to some parameters (counting fields) χ gives all cumulants (irreducible moments) of P(Q). This allows to measure and understand in depth the different transport properties contributing to the CGF. The Hamiltonian for the system under consideration is given by
The superconducting electrodes on the left and right can be described by BCS Hamiltonians with the gap of the superconducting terminal. We assume both gaps to be identical in order to simplify the analysis as no large-scale new features are to be expected.
In equilibrium the respective Hamiltonians can be written as ( )
using . The voltage is applied symmetrically so that . The corresponding Green’s functions are given in  and are abbreviated as respectively. The quasi-particle density of states is strongly energy dependent . The tunnel Hamiltonian has to take into account the time-dependent phase of the SCs.
For the resonant level model applied here  it is convenient to work with dressed electronic states by applying a polaron transformation.
which leads to a Hamiltonian where the electron-phonon interaction characterized by the coupling and resonance frequency is completely absorbed in the tunnel part of the Hamiltonian. We are left with
replacing and . We absorb the polaron shift for the dot energy by a redefinition of .
The general solution of the system described by Equation (1) seems out of question and approximations to proceed are generally necessary.
We want to first consider (without electron phonon interaction). We want to consider the case of strong onsite interaction and an odd number of electrons on the quantum dot. In this case the most prominent effect is the emergence of a spin 1/2 Kondo resonance around the Fermi edge for temperatures below the Kondo temperature . In normal conducting systems the Kondo temperature is directly related to the onsite interaction via 
where refers to the tunnel rate between the normal conductor and the quantum dot. Due to the additional energy scale in the problem, two scenarios may occur: for large the Kondo resonance couples to the quasiparticles in the superconductor leading to a behavior similar to the one for normal conducting systems  whereas for small the Kondo resonance is weakly coupled to the superconductors due to the absence of mobile electrons at the Fermi edge  . In this case we can greatly simplify the problem as an effective model in this case is the resonant level model with small hybridisations of the superconductor with the dot (which we assume to be equal for simplicity).
We evaluate the cumulant generating function using the generalized Keldysh technique  . is the counting field and is the phase difference of the two superconductors. in the case of charge transport is given by
where is the probability of charge q to be transferred during a given (long) measurement time .
Partial derivations of with respect to give direct access to the cumulants (irreducible moments). For our system the connection to the Hamiltonian is given by
where denotes the expectation value of with the substitution . C means the Keldysh contour and means time ordering on it. The counting field has to be both time and contour dependent. It changes sign on the different branches of the Keldysh contour to account for the charge transfer. Additionally, is nonzero only during the measurement interval . We use the standard expression
to find the CGF as the counting field derivative.
Compared to the case of contacts between normal metals and superconductors the counting field derivative has normal but also anomalous contributions leading to
We proceed by defining the exact-in-tunneling λ-dependent dot Green’s functions and the free electrode Green’s functions in Keldysh space:
Equation (10) can now be integrated with respect to λ to access the CGF. The normal contributions give rise to MARs and quasi-particle tunneling whereas the anomalous part gives rise to Josephson tunneling. We discuss both parts separately.
3. Josephson Tunneling
We first evaluate the second part of the expression in Equation (10). Following  and assuming the limit of the Kondo resonance being weakly coupled the result for the CGF is
The inelastic tunnel processes associated with the emission and absorption of phonons are described by the function being the Fourier transform of the phonon-phonon correlation function. For this correlation function we assume the phonons are thermally distributed, which may be due to coupling to a thermal environment given by the substrate of backgate  . The effect of coupling to an external bath can be characterized by an additional coupling constant and in the following we assume the bath to be fully Ohmic. In this case the phonon spectral density has Lorentzian shape
The phonon correlation function can now be analytically to be:
where . The functions are given by
In the case of small treated here we find the first harmonic for :
which corresponds to the dc and ac Josephson current depending on whether a voltage is applied or not.
Typically we would therefore expect a reduced Josephson current given the low transparency of the contact. We would not expect pronounced features of the phonons given that only the phase would depend on voltage.
4. Multiple Andreev Reflection
We go over to the evaluation of Equation (10). The calculation can be carried forward as in  . In the case of low tunneling transparency the electronic transport can be described by a sequential tunneling picture with a product of tunneling rates for each MAR
The result for the CGF for is then
where is the Fermi function.
We observe that multiple Andreev reflections are strongly suppressed as . The suppression can only be overcome in case conductance is boosted by a resonant phonon.
For voltages above the gap we recover the previously discussed situation for normal conducting systems as in  .
To conclude, we have derived the cumulant generating function for a superconducting quantum dot involving both Josephson tunneling as well as multiple Andreev reflections. Whereas we discovered that the Josephson effect is typically strongly suppressed in these setups and only slightly affected by resonant phonons, multiple Andreev reflections are strongly affected and show clear signs when resonant phonons are involved.
We expect that these results pave the way towards future usage of these junctions as transistors as they show a clearly nonlinear behavior. We would also believe that the corresponding current and noise characteristics can be observed in experimental setups soon. The corresponding analysis of the correspondence between the model developed here and the experimental results would clearly mark the next step also given the approximations in the Kondo limit taken here.
The author would like to thank Ingo Deppner for numerous discussions.