# Secondary ”Smile”-gap in the density of states of a diffusive Josephson junction for a wide range of contact types

## Abstract

The superconducting proximity effect leads to strong modifications of the local density of states in diffusive or chaotic cavity Josephson junctions, which displays a phase-dependent energy gap around the Fermi energy. The so-called minigap of the order of the Thouless energy is related to the inverse dwell time in the diffusive region in the limit , where is the superconducting energy gap. In the opposite limit of a large Thouless energy , a small new feature has recently attracted attention, namely, the appearance of a further secondary gap, which is around two orders of magnitude smaller compared to the usual superconducting gap. It appears in a chaotic cavity just below the superconducting gap edge and vanishes for some value of the phase difference between the superconductors. We extend previous theory restricted to a normal cavity connected to two superconductors through ballistic contacts to a wider range of contact types. We show that the existence of the secondary gap is not limited to ballistic contacts, but is a more general property of such systems. Furthermore, we derive a criterion which directly relates the existence of a secondary gap to the presence of small transmission eigenvalues of the contacts. For generic continuous distributions of transmission eigenvalues of the contacts, no secondary gap exists, although we observe a singular behavior of the density of states at . Finally, we provide a simple one-dimensional scattering model which is able to explain the characteristic ”smile” shape of the secondary gap.

###### pacs:

75.76.+j, 74.50.+r, 75.50.Xx, 75.78.-n## I Introduction

One of the most striking impacts of a contact with a superconductor (S) onto a small piece of normal metal (N) is the modification of the local density of states (LDOS). This effect, known as the superconducting proximity effect, is related to the induction of superconducting correlations resulting in a finite value of the pair amplitude on the normal side deutscher:69 (). In the absence of phonon mediated attraction between electrons on the normal side, decoherence between electronlike and holelike amplitudes leads to an exponential decay of the pair amplitude with distance from the contact, with a characteristic length scale exceeding the superconducting coherence length.

Modification of the LDOS on both sides of the contact strongly depends on the scattering properties of the contacts (described in terms of transmission eigenvalues) and the properties of the normal region (geometry, size, and impurity concentration). In the case of diffusive systems, it was predicted theoretically that the LDOS can even be fully suppressed in a specific energy range around the Fermi energy which is known as the minigap mcmillan:68 (). The minigap width is of the order of the inverse dwell time in the normal structure, which is given by the Thouless energy , where denotes the mean level spacing of the normal region and is the total conductance of the structure which is assumed large compared to the conductance quantum . In the decades after its discovery it has in detail been studied theoretically golubov:88 (); beenakker:92 (); belzig:96 (). The development of more elaborate experimental techniques with high spatial resolution made variations of the LDOS in this energy range accessible to experiments gueron:96 (); scheer:01 (); moussy:07 (); mourik:12 (); churchill:13 (); garcia:13 (); cherkez:14 (), which was found to be in agreement with theoretical calculations to a high degree leseur:08 (); pillet:10 (); wolz:11 ().

Much interest was concentrated on systems built up of a finite normal region sandwiched between two superconductors: a Josephson junction josephson:62 (). In such systems, another parameter, i.e., the phase difference between the superconducting order parameters comes into play and leads to a phase-dependent minigap belzig:99 (); golubov:04 (). Classical ballistic systems lodder:98 () were investigated as well as diffusive systems belzig:96 () and the crossover between both pilgram:00 (). It turns out that not only diffusive systems exhibit a minigap, but also ballistic systems with a chaotic classical motion melsen:97 (); lodder:98 (); VavilovLarkin:03 (); beenakker:05 (); kuipers:10 (); kuipers:11 ().

At this point, one might think that such structures are sufficiently explored and all relevant properties are understood. However, recently Levchenko reported the finding of a dip in the LDOS close to the gap edge for short diffusive Josephson junctions with ideal contacts levchenko:08 (). Actually, this dip was already seen in former publications golubov:96 (); belzig:96 (); wilhelm:00 (); heikkilae:02 (); bezuglyi:05 (); hammer:07 (), however, no special attention was paid to it. In a previous work reutlinger:14 () we found the peculiar result that the suppression of the LDOS at is not limited to a dip, but a secondary gap of finite width appears for a diffusive system or chaotic cavity with the normal region connected through ballistic contacts to the superconductors. This secondary gap has a finite width as a function of the superconducting phase difference symmetrically around zero and closes with the characteristic shape of a ”smile”. It is situated directly below the superconducting gap edge for large . For decreasing the upper edge of the secondary gap detaches from and the gap vanishes completely below a critical value of . We have furthermore shown that the secondary gap is robust against asymmetries in the setup, comprising a difference in ballistic couplings or a weak spatial dependence.

In this work, we investigate a wide range of possible nonballistic contacts and show that the secondary ”smile”-gap is not only an exotic feature which appears for ballistic contacts, but is a more general property of short diffusive or chaotic Josephson systems. Using quasiclassical Green’s functions in the form of the quantum circuit theory, we begin by generalizing our ballistic calculations to contacts with constant transmission eigenvalues , for which we calculate the density of states as a function of and . We find a secondary gap which scales for large like in the ballistic case. In a specific example of contacts described by different constant transmission eigenvalues we show that one is not limited to a single secondary gap, but this gap can be split up into multiple subgaps. Numerical considerations of continuous transmission distributions (diffusive, dirty, double ballistic contacts) suggest that the secondary gap below the superconducting gap edge vanishes if the contacts include channels with close-to-zero transmission coefficients. We prove this conjecture by an analytical calculation. By considering asymmetric setups with a tunnel contact on one and a ballistic contact on the other side we show however, that in this case a secondary gap can exist at slightly smaller energies. By considering a 3-node system we show that although the LDOS varies at different nodes the secondary gap appears either in all nodes or in none of them. It should as well be observable in the integrated DOS of the normal part. Finally, we provide a simple one-dimensional (1D) model in order to describe transmission through the normal region, which is able to explain the ”smile” shape of the secondary gap.

## Ii Model

In order to calculate the LDOS in the normal region, we make use of the retarded Green’s function in the quasiclassical approximation. In the diffusive or dirty limit, the angle-averaged Green’s functions are described by the nonlinear diffusive Usadel equation which has the form of a continuity equation for coherence functions including the leakage of coherence due to the finite energy difference between electrons and holes. Since the spatial dependence of the Green’s function is not important for our needs (for more details see Sec. III. F), we can solve the problem by applying the so-called quantum circuit theory nazarov:94 (); qt (); nazarov:99 (). We can discretize the system and reduce the equations to an algebraic problem. A sketch of the investigated system is shown in Fig. 1. The superconductors have equal energy gaps , however in general the phases of the order parameters can be different. Since the global phase is of no significance, only the phase difference enters our calculation and we can assign the phase to the left and right superconductors, respectively. The Green’s function in the normal node is determined by the constraint of matrix current conservation, including the currents to the two superconductors ( being the index denoting left and right lead) as well as the leakage current related to the volume of the normal region through :

(1) |

The scattering properties of the contacts are contained in the expressions for the matrix currents nazarov:99 ()

(2) |

in terms of a set of transmission eigenvalues . In general the transmission eigenvalues can be different on both sides. For continuous transmission distributions the sums must be replaced by integrals over the particular distributions. The Green’s functions in the leads are those of a bulk superconductor, given by with the spectral functions and being given by for and by for , being the Pauli matrices in Nambu space of electrons and holes. In the normal node the Green’s function can be parametrized as . and are related via the normalization condition for quasi-classical Green’s functions , which is equivalent to . In the general case with different contacts on both sides, this corresponds to the solution of two equations for two complex variables. Expanding (III.6) in Pauli matrices and comparing the coefficients provides two independent equations

(3) | |||||

(4) |

Note that and as well as are complex valued in general. All information on the contacts is contained in the characteristic functions given by

(5) |

with and being the conductance of the particular side. Again, for a continuous transmission distribution, the sums must be replaced by integrals over the particular distributions . For a symmetric setup and . Equation (4) becomes trivial and only one equation in one complex variable remains. From Eq. (3) we find

(6) |

The density of states is finally obtained from through , being the density of states at the Fermi energy of the normal state.

## Iii results

In previous analysis reutlinger:14 () this setup was investigated for ballistic contacts with all . It turned out that the secondary ”smile”-gap which appears in the symmetric case is stable under asymmetries . For an asymmetric setup, two further gaps, complementary to the usual minigap and the ”smile”-gap, appear symmetrically around . In this work, we want to extend these calculations and consider a wider range of contact types, corresponding to a wider range of characteristic functions , either described by discrete transmission eigenvalues or by continuous distributions . The idea of this work is to investigate the stability of the secondary gap under deviation from the ballistic limit. Especially, we want to determine which contact properties define the existence of the secondary gap, since it is known that for tunnel contacts no secondary gap is found. For this reason, we consider symmetric as well as asymmetric setups in the intermediate regime between the tunnel and ballistic limits.

### iii.1 Constant transmission

A natural generalization of the ballistic contact is to stick to constant transmission eigenvalues, however, to allow for . As approaches 0, the secondary gap is expected to disappear and the tunnel result for the LDOS should be reproduced. We begin by considering symmetric contacts and thus solve Eq. (6) with the characteristic function

We find a secondary gap in the LDOS similar to the ballistic result, which survives even for small but finite . The numerical results for the critical phase , for which the gap closes, as well as for the upper and lower gap edges and at , are shown in Fig. 2. The colored regions denote the gap. Above a special value , which scales linearly with and is given by , the upper gap edge is fixed to and the lower edge approaches for increasing following a power law. The linear scaling of with follows from Eq. (6) for and . The dependence of the lower gap edge on for is derived in the following.

Below , the upper edge is detached from and approaches the lower edge until the gap disappears at some critical value of which as well seems to scale linearly with . The maximum of the critical phase at which the secondary gap closes [Fig. 2 (b)] does not depend on , however, it is shifted to smaller Thouless energies with decreasing . The dependence of the critical phase on seems to scale linearly with . In the limit , the gap disappears and the tunnel result without secondary gap is reproduced. However, for each finite value of the secondary gap exists, if the Thouless energy is made large enough. To get further insight to the analytic properties we linearize Eq. (6) in the energy range below and around in the limit . We find

(7) |

with being the dimensionless energy relative to .

This equation can be solved analytically, however, the expression for the general solution is quite long and will not be given here. In Fig. 3, it is plotted for various values of . Figure 3 (a) shows the width of the secondary gap approaching with decreasing ; Figure 3 (b) shows the structure of the density of states above the upper edge of the minigap. Note the different scales of the energy axes in the two plots. For decreasing , the LDOS approaches the tunnel limit without secondary gap but with the usual singularities golubov:89 () at and above the minigap.

Considering provides an analytical expression for the critical phase

Similarly for we find an analytical expression for the critical energy in the limit of large describing the width of the gap

being a lengthy expression related to the solution of a quartic equation. For it has the form . The position of the minigap edge (Fig. 3 (b)) can as well be calculated analytically. For it is given by

### iii.2 Combination of transmission eigenvalues

A generalization of the calculations from the previous section can be achieved by considering not only one constant transmission eigenvalue, but a whole set of different transmission eigenvalues, each weighted with a specific weight . We stick to a symmetric system with only one set of transmission eigenvalues and weights describing both sides. From the huge variety of possible sets, which could be analyzed, we pick only one in order to demonstrate that the secondary gaps structure in principle is not limited to only a single gap: An even finer subdivision of the LDOS below can be observed for certain contact types. We calculate the LDOS for one representative set given by

0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1 | |
---|---|---|---|---|---|---|---|---|---|---|

10 | 5 | 5 | 2 | 3 | 20 | 50 | 70 | 90 | 200 |

A plot of the numerical results in the energy range below is shown in Fig. 4. We find that for certain sets of transmission eigenvalues not only one secondary gap appears, but the DOS acquires an even finer structure with multiple subgaps. The number of subgaps depends strongly on the set under consideration and on . In the presented case we find three gaps at . Similar to the previously found secondary gaps they are symmetric around and vanish at some critical phase, which is not the same for different subgaps.

### iii.3 Continuous transmission distributions

So far, we investigated systems where scattering in the contacts is described by constant transmission eigenvalues. However, in systems experimentally accessible scattering is rather described by a continuous transmission distribution than by discrete transmission eigenvalues . It is thus of crucial interest whether the secondary gap appears as well if continuous transmission distributions are considered. Here, we investigate three generic contact types, each characterized by a distribution of transmission eigenvalues of the form with . The normalization constant is determined by the condition . The distributions for and correspond to a dirty and a diffusive connector schep:97 (); dorokhov:82 () , respectively. The distribution for is equivalent to two ballistic connectors with equal conductances in series qt ().

Again we stick to symmetric setups with equal contacts on both sides. The characteristic functions corresponding to the considered distributions have a relatively compact form belzig:00 (). For a dirty connector () it is given by . In case of a diffusive contact () the characteristic function is . And for the double ballistic contact qt () we have .

Figures 5(a)–(c) show the numerical results for the LDOS calculated from Eq. (6) for the three cases. Since in previous calculations the suppression of the LDOS around was strongest at , only this case is presented. Figure 5 (a) contains the numerical results for dirty contacts () for different values of . No signature for a suppression of the LDOS at the superconducting gap edge is found. The plots in Fig. 5 (b) are the results for diffusive contacts (). For , a weak suppression at can be seen. The inset with a higher resolution of the energy range of interest, however, shows, that this suppression is no gap. At first glance, this seems to disagree with Ref. levchenko:08 () but is possibly due to differences in the considered geometries. The lowest Fig. 5 (c) shows the results for the third contact type () with the highest weight at transmission eigenvalues around of all three distributions. We find a strong suppression at for all considered values of . The inset of this plot confirms that the LDOS is suppressed to 0 at .

In summary we find that with decreasing weight of at the suppression of energy levels at is reinforced. This observation supports the idea that the existence of Andreev bound states with energies directly below is related to the tunnel character of the boundaries, i.e., the transport channels with transmissions at .

In order to check this, a fourth type of transmission distribution is considered. Each contact in this system is built up of two ballistic contacts in series ( and ) having different conductances. Note that the total setup is still symmetric. The corresponding transmission distribution fundamentally differs from the previously considered distributions in the sense that it has no contribution at small transmissions, i.e., below a critical value given by . For it follows which corresponds to the previously investigated distribution with . Since we again stick to symmetric setups we drop the lead index in the following. The characteristic function is given by qt ()

The results are plotted in Fig. 6. Compared to the previous distributions, a gap of finite width appears directly below (inset of Fig. 6). As expected, the gap appears when there is no contribution of the transmission distribution around . This agrees with the results of Sec. III. A. In the following we derive a criterion that relates the existence of the secondary gap directly below to the weight of the transmission distribution around .

### iii.4 Analytical criterion

In this section, we show that the existence of the secondary gap below is indeed directly related to the distribution of the transmission eigenvalues in the vicinity of . We show that the existence of a minimal transmission results in a secondary gap below . To reach this conclusion, we consider Eq. (6) for a symmetric setup containing one general function and . It reads as

We linearize this equation in in the limit . The left side is expanded in which must be small in order to be valid. This has to be verified for the solution. Using the limiting form of and at small we get

In the leading order in , the general expression for given after Eq. (5) yields: . Introducing rescaled variables we get an equation without explicit dependence on . With the definitions and we have

(8) |

The LDOS is related to the real part of , so a purely imaginary solution for at small means having a gap at energies E close to . To show this we consider the sum in Eq. (5), which is an integral in the continuous case. For an arbitrary distribution , which can be normalized to satisfy , it reads

where the minimal transmission in was used to replace the lower boundary in the integral. Considering small and assuming with (later verified by the solution) we can neglect all other terms in the denominator compared to . We find

being a constant factor defined as . It is constrained to the interval , corresponding to a ballistic setup. Its value depends on the exact form of in the whole interval, particularly on . For a diffusive connector which is cut at we have , which becomes for and approaches as approaches . With this Eq. (8) reduces to a quadratic equation

with solutions . Both solutions are purely imaginary for signifying a gap in the LDOS. However only the solution for is consistent with the previously made assumption with for small . Since this solution is finite for small the second assumption, which assumed to be large, can always be fulfilled for by choosing sufficiently large. This is in agreement with our numerical results which predict no secondary gap below some critical value of . This critical depends on the value of and thus on the whole transmission distribution . The condition is related to the existence of a , since only with leads to . To conclude this section, we have shown that for an arbitrary transmission distribution without contribution in a finite interval above a secondary gap appears directly below , if is sufficiently large.

### iii.5 Asymmetric setup

In this section, we demonstrate that having the distribution function of transmission coefficients at below some is a sufficient but not necessary condition for having a secondary gap. We consider a device with two nonidentical junctions and show that the secondary gap may exist even if the latter condition on is violated. The secondary gap, however, does not appear directly below but is shifted to slightly smaller energies, thus it appears in a different regime than considered in Sec. III D.

The most interesting case to which we confine ourselves here is the one which combines the two extremal contact types: A tunnel contact with for all transport channels on one side and a ballistic contacts with for all channels on the other side. The conductances of both sides enter our calculation via the relation , where denotes the conductance of the tunnel contact and corresponds to the ballistic contact. For the role of the tunnel contact is negligible and the result of a symmetric ballistic system at showing a secondary gap below reutlinger:14 () should be reproduced. It is of particular interest how the transition from the ballistic limit to the tunnel limit occurs. Results of the previous section indicate that there should be no secondary gap just below : The condition on for the appearance of the secondary gap is violated by the presence of a tunnel junction even for .

Figure 7 shows the results of our numerical calculation. In Fig. 7 (a), the phase is fixed at and the vanishing of the secondary gap is shown for different values of as increases. At small , the upper gap edge is close to its ballistic value and decreases as increases. Similarly, the lower gap edge is close but slightly below its value of the symmetric ballistic case and decreases with increasing . At some critical value of which depends on the secondary gap disappears. This critical value is smaller for smaller . In Fig. 7(b), the Thouless energy is fixed at and the phase dependence of the secondary gap is plotted for different values of . We find that for the considered system the secondary gap has its maximum width not at as one might expect, but at . With decreasing , the result of the symmetric ballistic system is approached at all phases. Compared to previous findings for asymmetric ballistic contacts reutlinger:14 (), where a similar behavior with decreasing was found, no band of finite DOS separates the gap at from the gap at .

### iii.6 Spatial dependence

In order to achieve a spatial resolution of the local density of states (LDOS) we consider a symmetric model of three normal islands connected to two superconductors at . Due to symmetry the Green’s functions in the left and right normal nodes are equal. Both nodes are thus called in the following, the central node is called . The nodes are connected via a ballistic conductance to the superconductor and via to . In each normal node, electron-hole-decoherence is described through Thouless energies , respectively. The system setup is sketched in Fig. 8.

Matrix current conservation in one of the nodes and in the node determines the Green’s function and the LDOS in the particular node nazarov:94 (); nazarov:99 ():

The matrix currents and are defined according to the definition (2). In general, the system is described by three parameters , , and . We fix one parameter by considering only systems with equal mean level spacings in all three normal nodes . Furthermore, we fix . In this case for the LDOS in and are equal and correspond to the result of Figure 1(a) in reutlinger:14 (), and for we have the BCS-DOS in and again the result of Figure 1(a) of reutlinger:14 () in . Figure 9 shows the LDOS in the two nodes for [Fig. 9 (a)] and for [Fig. 9 (b)]. Taking this into account, and can be expressed in terms of . We find

In Fig. 10, we show the numerical results for the LDOS in both normal nodes in dependence of and energy in the secondary gap region below the superconducting gap edge . Figure 10 (a) shows the result in the outer nodes , and Fig. 10 (b) shows the result for the inner node . The white region in Fig. 10 (a) denotes . The main finding of our calculations is that the LDOS in the two nodes differs only where in both nodes. Whenever it is zero in the central node it is also zero in the outer nodes and vice versa. We thus find a behavior of the secondary gap similar to what is already known from the usual minigap pilgram:00 (). The width of this gap is not position dependent, only the LDOS above/below the particular gap edge varies with position.

We thus expect the secondary gap not only in the LDOS of a singular point, but as well in the integrated DOS of a finite region. Depending on the parameters, not for every system does a secondary gap appear. However, if it appears in one point, it exists also in every other point of the normal part. The previously used model with only a single normal node between the superconductors is thus sufficient if the main interest concerns the existence of the secondary gap and its properties. However, with this method we cannot calculate a position-resolved LDOS and thus cannot make statements about the integrated DOS in the energy interval between minigap and secondary gap.

## Iv 1D Scattering Model

The secondary gap we found for diffusive Josephson systems was calculated using Green’s function techniques in the quasiclassical approximation. Whereas this method is very powerful in calculating expectation values of physical observables, it does not provide a simple intuitive explanation for the absence of Andreev levels in the secondary gap region and the dependence of these levels on the phase difference between the superconductors. In this section, we investigate a simple 1D scattering model which is able to explain qualitatively the secondary ”smile”-gap. However, since we deal with diffusive or chaotic scattering systems with large conductance, we should not expect to reproduce the details of 3D solutions.

### iv.1 Single-trajectory Andreev level

We consider a semiclassical path between the left () and the right () superconductor [Fig. 11] and first recall the characteristics of Andreev bound states between superconductors on a ballistic trajectory. The bound state energies follow from the semiclassical quantization condition:

(9) |

Here, the first term is the phase difference acquired between electron and hole upon traversing the normal region. is essentially the inverse traversal time, which could also be due to ballistic motion for a trajectory of length . The second factor is twice the energy-dependent Andreev reflection phase and the third term the phase difference between the superconducting order parameters. All terms together have to add to an integer multiple of .

This equation reproduces two limiting cases. In long junction limit , we replace and find the usual spectrum of Andreev levels . In this case, levels move up and down in energy linearly with the phase difference . The lowest positive energy states have the energies . The levels split with and cross 0 at , corresponding to the closing of the minigap. In the opposite, short junction limit , we neglect the first term in Eq. (9) and find the Andreev levels .

The most interesting case is the ”not-so-short” junction limit . Assuming the energy is close to , we can replace by in the first term of Eq. (9), and taking we obtain

(10) |

Thus, we obtain two states shifted in phase by the (small) parameter . They touch the gap at the critical phase and the maximal distance to (at ) is . This is in quantitative agreement with the characteristics of the secondary gap found previously within the quasiclassical Green’s function theory. Note that in the present approximation, the two levels cross at . We can expect that finite backscattering will lead to an anticrossing and the phase dependence of the level resembles the ”smile” shape of the secondary gap.

### iv.2 Single-trajectory Andreev level with scattering

We investigate a simple model for the anticrossing and calculate the Andreev bound-state energies for a 1D-model with impurity scattering modeled by a scattering matrix. Although this model takes only backward scattering into the same trajectory into account and neglects the complex interference effects of three-dimensional impurity scattering which are covered by our original Green’s function calculations, the results provide an understanding of the phase-dependent Andreev level density of states. The bound-state energies are obtained from the scattering matrices in the normal region qt (). We consider the geometry shown in Fig. 11. The normal scattering matrix encompasses the back scattering at the impurity as well as the dynamical phases along the trajectory to the superconductor and is given by

where accounts for the position of the impurity along the path and is the transmission probability. The normal region scattering matrix for holes is related through .

The scattering matrices for electron-hole conversion at the interface to the superconductors are given by and , respectively. Note that the -space is not Nambu space. An electron arriving at either superconductor is reflected as a hole traveling towards the normal region from the same side, thus Andreev reflection is described by a diagonal matrix. The condition for a bound state reads as

(11) |

The bound-state energies in dependence of are plotted in Fig. 12 for different values of . Without backscattering in the normal region, at the two Andreev levels are degenerate [Fig. 12 (a)]. Taking into account impurity scattering in the normal part [Fig. 12 (b)] this degeneracy is lifted (the exact curve depends on the position where scattering occurs, i. e., on the parameter ). This results in the characteristic shape of the minigap and the secondary ”smile”-gap below . Figure 12(b) shows the -averaged results for Andreev bound states with one scattering event with (weak scattering). It is worth mentioning that only channels without scattering contribute to the zero-energy Andreev states at (not shown). For paths with one or more scattering events (more scattering matrices in the normal part), these levels are shifted to higher energies. Thus, we have shown that the secondary gap can be understood from the phase-dependence of the Andreev level when the junction length exceeds a length of the order of the superconducting coherence length, given by . The ”smile” shape can be traced back to the effect of backscattering.

## V Conclusion

To summarize, we have calculated the local density of states for diffusive Josephson systems for a wide range of contact types with attention to the energy range below , in which a secondary gap can appear. We have generalized previous calculations for ballistic contacts reutlinger:14 () and shown that the secondary ”smile”-gap is a robust feature in the proximity density of states for large Thouless energies. We thus suggest that this feature should be accessible to an experimental detection by means of high-resolution scanning tunneling spectroscopy and want to encourage research in this direction.

Acknowledgments. J. R. and W. B. were supported by the DFG through SFB 767 and BE 3803/5 and by the Carl Zeiss Foundation. Y. N. and L. G. thank the Aspen Center for Physics, supported in part by NSF Grant No. PHYS-1066293, for hospitality. Work at Yale is supported by NSF DMR Grant No. 1206612.

### References

- G. Deutscher and P. G. de Gennes, in Superconductivity, edited by R. D. Parks (Dekker, New York, 1969), Vol. 2, p. 1005.
- W. L. McMillan, Phys. Rev. 175, 537 (1968).
- A. A. Golubov and M. Yu. Kuprianov, J. Low Temp. Phys. 70, 83 (1988).
- C. W. J. Beenakker and H. van Houten, in: Single-Electron Tunneling and Mesoscopic Devices, ed. by H. Koch and H. Lübbig (Springer, Berlin, 1992)
- W. Belzig, C. Bruder, and G. Schön, Phys. Rev. B 54, 9443 (1996).
- S. Gueron, H. Pothier, N. O. Birge, D. Esteve, and M. H. Devoret, Phys. Rev. Lett. 77, 3025 (1996).
- E. Scheer, W. Belzig,Y. Naveh, M. H. Devoret, D. Esteve, and C. Urbina, Phys. Rev. Lett. 86, 284 (2001).
- N. Moussy, H. Courtois, and B. Pannetier, Europhys. Lett. 55, 861 (2007)
- V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Science 336, 1003 (2012).
- H. O. H. Churchill, V. Fatemi, K. Grove-Rasmussen, M. T. Deng, P. Caroff, H. Q. Xu, and C. M. Marcus, Phys. Rev. B 87, 241401 (2013).
- L. Serrier-Garcia, J. C. Cuevas, T. Cren, C. Brun, V. Cherkez, F. Debontridder, D. Fokin, F. S. Bergeret and D. Roditchev, Phys. Rev. Lett. 110, 157003 (2013)
- V. Cherkez, J. C. Cuevas, C. Brun, T. Cren, G. Menard, F. Debontridder, V. S. Stolyarov, and D. Roditchev, Phys. Rev. X 4, 011033 (2014)
- H. le Sueur, P. Joyez, H. Pothier, C. Urbina, and D. Esteve, Phys. Rev. Lett. 100, 197002 (2008).
- J.-D. Pillet, C. H. L. Quay, P. Morfin, C. Bena, A. Levy Yeyati, and P. Joyez, Nat. Phys. 6, 965 (2010).
- M. Wolz, C. Debuschewitz, W. Belzig, and E. Scheer, Phys. Rev. B 84, 104516 (2011).
- B. D. Josephson, Phys. Letters 1, 251 (1962).
- W. Belzig, F. K. Wilhelm, C. Bruder, G. Schön, and A. D. Zaikin, Superlatt. Microstruct. 25, 1251 (1999).
- A. A. Golubov, M. Yu. Kupriyanov, and E. Ilichev, Rev. Mod. Phys. 76, 411 (2004).
- A. Lodder and Yu. V. Nazarov, Phys. Rev. B 58, 5783 (1998).
- S. Pilgram, W. Belzig, and C. Bruder, Phys. Rev. B 62, 12462 (2000).
- M. G. Vavilov and A. I. Larkin, Phys. Rev. B 67, 115335 (2003).
- C. W.J . Beenakker, Lect. Notes Phys. 667, 131 (2005).
- J. A. Melsen, P. W. Brouwer, K. M. Frahm, and C. W. J. Beenakker, Physica Scripta T69, 233 (1997)
- J. Kuipers, D. Waltner, C. Petitjean, G. Berkolaiko, and K. Richter, Phys. Rev. Lett. 104, 027001 (2010).
- J. Kuipers, T. Engl, G. Berkolaiko, C. Petitjean, D. Waltner, and K. Richter, Phys. Rev. B 83 , 195316 (2011).
- A. Levchenko, Phys. Rev. B 77, 180503(R) (2008).
- J. C. Hammer, J. C. Cuevas, F. S. Bergeret and W. Belzig, Phys. Rev. B 76, 064514 (2007).
- T. T. Heikkilae, J. SÃ¤rkkÃ¤ and F. K. Wilhelm, Phys. Rev. B 66, 184513 (2002).
- F. K. Wilhelm and A. A. Golubov, Phys. Rev. B 62, 5353 (2000).
- A. Golubov and M. Kupriyanov, Physica C 259, 27 (1996).
- E. V. Bezuglyi, A. S. Vasenko, V. S. Shumeiko, G. Wendin, Phys. Rev. B 72, 014501 (2005)
- J. Reutlinger, L. Glazman, Yu. V. Nazarov and W. Belzig Phys. Rev. Lett. 112, 067001 (2014).
- Yu. V. Nazarov, Phys. Rev. Lett. 73, 134 (1994).
- Yu. V. Nazarov and Ya. M. Blanter, Quantum Transport (Cambridge University Press, Cambridge, 2009)
- Yu. V. Nazarov, Superlatt. Microstruct. 25, 1221 (1999).
- A. A. Golubov and M. Yu. Kupriyanov, ZhETF 96, 1420 (1989) [Sov. Phys. JETP 69, 805 (1989)].
- K. M. Schep and G. E. W. Bauer, Phys. Rev. Lett. 78, 3015 (1997); Phys. Rev. B. 56, 15860 (1997).
- O. N. Dorokhov, Solid State Commun. 51, 381 (1984).
- W. Belzig, A. Brataas, Yu. V. Nazarov, G. E. W. Bauer, Phys. Rev. B 62, 9726 (2000)