Magic number theory of superconducting proximity effects and Wigner delay times in graphene-like molecules

When a single molecule is connected to external electrodes by linker groups, the connectivity of the linkers to the molecular core can be controlled to atomic precision by appropriate chemical synthesis. Recently, the connectivity dependence of the electrical conductance and Seebeck coefficient of single molecules has been investigated both theoretically and experimentally. Here we study the connectivity dependence of the Wigner delay time of single-molecule junctions and the connectivity dependence of superconducting proximity effects, which occur when the external electrodes are replaced by superconductors. Although absolute values of transport properties depend on complex and often uncontrolled details of the coupling between the molecule and electrodes, we demonstrate that ratios of transport properties can be predicted using tables of 'magic numbers,' which capture the connectivity dependence of superconducting proximity effects and Wigner delay times within molecules. These numbers are calculated easily, without the need for large-scale computations. For normal-molecule-superconducting junctions, we find that the electrical conductance is proportional to the fourth power of their magic numbers, whereas for superconducting-molecule-superconducting junctions, the critical current is proportional to the square of their magic numbers. For more conventional normal-molecule-normal junctions, we demonstrate that delay time ratios can be obtained from products of magic number tables.

For example, such junctions are often described using the Landauer formula = 0 ( ), where 0 =

METHODS
To illustrate how these phase-dependent phenomena can be predicted using magic number theory,  In a typical experiment using mechanically controlled break junctions or STM break junctions [13][14][15][16][17][18]21 , fluctuations and uncertainties in the coupling to normal-metallic electrodes are dealt with by measuring the conductance of such molecules many thousands of times and reporting the statistically-most-probable electrical conductance. If is the statistically-most-probable conductance of a molecule such as 1 (see Figure 1), with connectivity i-j and is the corresponding conductance of a molecule such as 2 (see From a conceptual viewpoint, magic ratio theory views the shaded regions in Fig. 1 as "compound electrodes", comprising both the anchor groups and gold electrodes, and focuses attention on the contribution from the core alone. As discussed in 28 , the validity of Eq. (1) rests on the key foundational concepts of weak coupling, locality, connectivity, mid-gap transport, phase coherence and connectivity-independent statistics. When these conditions apply, the complex and often uncontrolled contributions from electrodes and electrode-molecule coupling cancel in conductance ratios and therefore a theory of conductance ratios can be developed by focussing on the contribution from molecular cores alone.
The term "weak coupling" means that the central aromatic subunit such as anthanthrene should be weakly coupled to the anchor groups via spacers such as acetylene , as shown in Fig. 1. Weak coupling means that the level broadening Γ and the self energy Σ of the HOMO and LUMO should be small compared with the HOMO-LUMO gap . Any corrections will then be of order Γ/ or Σ/ , which means that such terms can be ignored, provided the Fermi energy lies within the gap. Clearly a central condition for the applicability of the Landauer formula and therefore magic-number theory is that the molecular junction is described by a time independent mean-field Hamiltonian. Coulomb interactions can be included in such a Hamiltonian, at the level of a self-consistent mean field description such as Hartree, Hartree-Fock or DFT. The concept of 'mid-gap transport' is recognition of the fact that unless a molecular junction is externally gated by an electrochemical environment or an ele ctrostatic gate, charge transfer between the electrodes and molecule ensures that the energy levels adjust such that the Fermi energy E F of the electrodes is located in the vicinity of the centre of the HOMO-LUMO gap and therefore transport takes place in the co-tunnelling regime. In other words, transport is usually 'off-resonance' and the energy of electrons passing through the core does not coincide with an energy level of the molecule. Taken together, these conditions ensure that when computing the Gree n's function of the core, the contribution of the electrodes can be ignored. The concept of 'phase coherence' recognises that in this co-tunnelling regime, the phase of electrons is usually preserved as they pass through a molecule and therefore transport is controlled by QI. 'Locality' means that when a current flows through an aromatic subunit, the points of entry and exit are localised in space. For example , in molecule 1 (see Figure 1), the current enters at a particular atom i and exits at a particular atom j. The concept of 'connectivity' recognises that through chemical design and synthesis, spacers can be attached to different parts of a central subunit with atomic accuracy and therefore it is of interest to examine how the flow of electricity depends on the choice of connectivity to the central subunit. The condition of "connectivity-independent statistics" means that the statistics of the coupling between the anchor groups and electrodes should be independent of the coupling to the aromatic core. To be more precise, we note that in an experimental measurement of single-molecule conductance using for example a mechanically-controlled break junction, many thousands of measurements are made and a histogram of logarithmic conductances is constructed. This statistical variation arises from variability in the electrode geometry and in the binding conformaton to the electrodes of terminal atoms such as the nitrogens in figure 1. The assumption of "connectivity-independent statistics" means that this variability is the same for the two different connectivities of figure 1. When each of these conditions applies, it can be shown 2,25,26 that in the presence of normal-metallic electrodes, the most probable electrical conductance corresponding to connectivity i-j is proportional to | ( F )| 2 where ( F ) is the Green's function of the isolated core alone, evaluated at the Fermi energy of the electrodes. In the absence of time -reversal symmetry breaking, ( F ) is a real number. Since only conductance ratios are of interest, we define magic numbers by where A is an arbitrary constant of proportionality, chosen to simplify magic number tables and which cancels in Eq. (1). Magic ratio theory represents an important step forward, because apart from the Fermi energy F , no information about the electrodes is required. The question we address below is how is the theory modified in the presence of superconducting electrodes and how can the theory be extended to describe Wigner delay times?
In the presence of normal-metallic electrodes, many papers discuss the conditions for destructive quantum interference (DQI), for which ≈ 0 9,18,27,29-33 . On the other hand, magic ratio theory aims to describe constructive quantum interference (CQI), for which may take a variety of non-zero values. If is the non-interacting Hamiltonian of the core, then since ( F ) = ( F − ) −1 , the magic number table is obtained from a matrix inversion, whose size and complexity reflects the level of detail contained in .
The quantities were termed "magic" 2,25,26 , because even a simple theory based on connectivity alone yields values, which are in remarkable agreement with experiment 25 . For example, for molecule 1 (see Figure 1), the prediction was = −1, whereas for molecule 2, = −9 and therefore the electrical conductance of molecule 2 was predicted to be 81 times higher than that of 1, which is close to the measured value of 79. This large ratio is a clear manifestation of quantum interference (QI), since such a change in connectivity to a classical resistive network would yield only a small change in conductance. To obtain the above values for and , the Hamiltonian was chosen to be where the connectivity matrix of anthanthrene is shown in Fig. 2. In other words, each element was chosen to be -1 if , are nearest neighbours or zero otherwise and since anthanthrene is represented by the bipartite lattice in which odd numbered sites are connected to even numbered sites only, is block off-diagonal. The corresponding core Green's function evaluated at the gap centre = 0 is therefore obtained from a simple matrix inversion (0) = − −1 . Since and therefore − −1 are block off-diagonal, this yields = ( 0 ̅ ̅ 0 ) ∝ (0), where is the magic number table of the polycyclic aromatic hydrocarbons (PAHs) core. The connectivity matrix and off-diagonal block of the magic number table for anthanthrene are shown in Figure 2b and c respectively. As noted above, for molecule 1 (see Figure   1), with connectivity 22-9, 22,9 = −1, whereas for molecule 2, with connectivity 12-3, 12,3 = −9. Magic number tables such as Figure 2c are extremely useful, since they facilitate the identification of molecules with desirable conductances for future synthesis. Conceptually, tables obtained from Hamiltonians such as (3) are also of interest, since they capture the contribution from intra-core connectivity alone (via the matrix , comprising -1's or zeros), while avoiding the complexities of chemistry.
The fact that magic number theory predicts experimentally-measured conductance ratios for a range of molecules 25 demonstrates that at least for PAHs measured to date, magic number theory is valid.
In what follows, the studied molecules are chosen, because their behaviour in the presence of normalmetal (ie gold ) contacts they exhibit a sizeable connectivity dependence in their electrical conductances, as demonstrated experimentally and theoretically in ref 25 . The tight binding parameters describing the molecular core are chosen following the philosophy in refs 2,25 , where the aim is to highlight the role of connectivity in determining the transport properties of these molecular cores. For this reason, the hopping integrals are set to unity and the site energies 0 are set to zero. In other words, the unit of energy is the hopping integral and the site energy is the energy origin. This means that the Hamiltonian is simply a connectivity matrix and therefore all predicted effects are a result of connectivity alone.
Remarkably, as demonstrated in 2,25 , this approach yields the experimentally-measured conductance ratios of a range of PAHs.
In the wide band limit, the only other parameter is the coupling between the terminal sites and the electrodes. Our aim is to compute ratios of transport properties corresponding to different connectivities to the electrodes. As shown in refs 2,25 , transport ratios do not depend on these couplings, provided they are sufficiently weak. Comparison with experiment in these refs shows that this weak -coupling criterion is satisfied by acetylene linker groups connecting the aromatic core to pyridyl anchor groups, which in turn bind to electrodes, as shown in Fig. 1.

ANALYTIC RESULTS
For convenience, we first state our main analytic results. Details of the derivations and discussion of the results will be given in later sections and in the Supplementary Information. The first outcome of the study is that for normal-molecule-superconducting (N-M-S) junctions Eq. (1) is replaced by where and label the atoms in contact with the normal electrode, while and label atoms in contact with the superconducting electrode. Equation (4) where − is the superconducting phase difference between the right and left superconducting electrodes. In equation (6), labels the atom in contact with the normal electrode, while ( ) label the atom in contact with the superconducting electrode S.
According to equation (6) the current through the normal lead is sensitive to the phase difference − . As we shall see in Sec. 4c, this phenomenon can be understood as an interference effect between two transport paths. Finally, for N-M-N junctions (see section 5), the ratio of Wigner delay times corresponding to connectivities , and , is

RESULTS FOR MOLECULAR JUNCTIONS WITH ONE OR MORE SUPERCONDUCTING ELECTRODES
In the presence of superconductivity, electron transport is controlled by both normal electron scattering and Andreev scattering. When the normal region between electrodes is a diffusive metal, one traditionally treats transport using quasi-classical equations 34 . Such an approach is not appropriate for phase-coherent transport through molecular cores, where the arrangement of atoms within the central region is deterministic. Instead, one should use a scattering approach which preserves such atomic-scale details 35 Calculating the conductance using the Landauer formula, one can indeed arrive s at equation (4). To demonstrate the validity of equation (4) we consider the normal-molecule-superconducting junctions with a pyrene core, shown in Figure 3.

Figure 3:
Normal-pyrene-superconducting junction with normal electrode coupled to site 1' and with superconducting contact at site 3. γ is the hopping amplitude between atoms, ∆ is the pairing potential in the superconducting contact and γ c is the coupling between the molecule and the superconducting contacts.
In our numerical calculation we used the tight binding model to describe the non-interacting Hamiltonian of the molecule with hopping amplitude γ. Following the reasoning of reference 42 , the conductance can be calculated using the normal and Andreev reflection amplitudes by the formula: where is the number of propagating electron-like channels in the normal lead, 0, is the normal and , is the Andreev reflection coefficient associated with a normal lead connected to site , and with a superconducting lead connected to site of the molecule. The second equality is valid for low biases ≪ |∆| , where ∆ is the superconducting pair potential. These reflection coefficients can be calculated via the Green's function theory of reference 43 .
Using the numbering convention of the sites given in Figure 3, one can easily show that the magic numbers are non-zero only between non-primed i and primed j' sites. The magic numbers between sites i and j' are shown in Table 1. We assumed that the coupling γ c between the molecule and the superconducting contacts is weak. Then it is reasonable to take the pairing potential ∆ to be non-zero only in the contact.  Figure 3 for the meaning of indexes i and j'. Magic numbers connecting two sites both labeled by a prime (or both labeled without a prime) are zero.
We calculated numerically the electrical conductances between different pairs of sites and '. The superconducting pairing potential in our calculations was |∆| = 1 meV, the hopping amplitude was γ = 2.4 eV and coupling amplitude was γ c = 0.45 γ. The onsite potential was zero at all sites. Our numerical results are summarized in Table 2. As one can see, our theoretical prediction in Eq. (4) is confirmed; namely the ratio of the calculated conductance between different sites and ′ agrees very well by the fourth power of magic numbers given in Table 1.  trodes, which do not affect conductance ratios. Furthermore conductance ratios are independent of the parameters used to define the electrodes. The above value of Δ was chosen, because superconducting gap of common superconductors such as niobium, tantalum and mercury is on the scale of a meV. The value of γ_c was chosen to ensure that the level broadening due to the contacts is small compared to the HOMO-LUMO gap, which is in the experimentally-relevant regime. Just as conductance ratios are independent of the parameters used to define the normal lead, they are also independent of the parameters used to define the superconducting electrode, including the size of the superconducting energy gap.
Magic number theory is a weak coupling theory and eventually, as the coupling to the electrodes increases, there will be a difference between magic number theory and a full tight binding calculation. To illustrate this point and at the same time to show that the deviati on is small, we performed the same calculation for γ c = 0.85 γ, which is larger than in the previous calculation. The results are given in Table 3.
Note that the ratios of the calculated conductance deviate from the fourth power law of magic numbers given in equation (4). However, for γ c < 0.85 γ the deviation from the weak coupling limit remains of order 10 %.

4b. S-M-S' junctions
In mesoscopic superconductivity it is well known that in phase biased Josephson junctions, at low temperatures, under rather general conditions, the critical current is inversely proportional to the normal state resistance of the junction 44 , i.e., ~1~. Therefore one may expect that in S-M-S' junctions the critical current flowing between the superconducting electrodes connected to sites i and j of the molecule is proportional to the transmission amplitude between these points leading to ~2. Equation  Using our numerical method presented in Ref. 45 we calculated the critical current between the different pairs of sites i and j ' . We used a superconducting pairing potential |∆| = 1 meV, a hopping amplitude γ = 2.4 eV and coupling amplitude γ c = 0.45 γ. The onsite potential was zero at all sites. The calculated critical current are summarized in Table 4. Clearly our theoretical prediction in equation (5) is confirmed, namely the ratio of the calculated critical currents between different sites i and j ' agree very well by the squares of the magic numbers given in Table 1.  We performed the same calculation for larger γ c . The molecular Josephson effect in the strong coupling limit was also studied in a recent work 46 . For stronger coupling, the critical current becomes much larger than in the weak coupling limit, which is consistent with the results of reference 46 . (The critical current increases by a factor of 10 4 compared to the weak coupling limit.) The results for the critical currents are given in Table 5. Note that the ratio of the calculated critical currents starts to deviate from the ratios of the corresponding magic numbers.

4c. N-M-SS' junctions
We now examine the conductance of an N-M-SS' Andreev interferometer as a function of the superconducting phase difference − (see Figure 5) between the left and right superconducting electrodes.
When electrons go through the path connecting the normal and superconducting leads, the electronic states acquire a phase that depends on the superconducting pair potential (due to the Andreev reflection at the superconducting surface). For N-M-S junctions we have seen that the transmission amplitude related to an electron incoming from the normal electrode and reflected back as a hole is proportional to  According to Equation (10), in general the conductance at the normal electrode is expected to show a periodic interference pattern as a function of the phase difference between the left and right superconducting leads. To verify our analytic expression, we compared the predictions of equation (10) to numerical tight binding simulations, see Figure 6. To calculate the conductance at the normal lead numerically, we make use of equation (10) generalized for three-terminal systems. Consequently, the Green's function, and also the normal and the Andreev reflection coefficients in equation (10)   One can see in Figure 6 that the interference pattern strongly depends on the position of the contacts. In agreement with the numerical results, we found that the interference pattern can only be observed if the connectivity between each pair of the normal and superconducting el ectrodes is finite [see Figure 6(a)-(c)], otherwise one of the interfering paths from formula (10) would be missing. Such a situation is shown in Figure 6(d), where the conductance is indeed independent of the phase difference φ R − φ L .

RESULTS FOR THE WIGNER DELAY TIME IN N-M-N JUNCTIONS
The Wigner delay time function was proposed by Wigner in 1955 for a single scattering channel derived from a Hermitian operator based on the scattering amplitude and then generalized by Smith in 1960 to the multichannel scattering matrices [47][48][49] . Delay times from different molecular orbitals can be measured experimentally, as described in 50 . More generally, within a molecular junction driven by an ac applied voltage, they are related to the off-diagonal elements of the admittance matrix, which describes the current response to such a time-varying voltage 51 .
Consider a scatterer, with one-dimensional leads connected to sites and , whose transmission amplitude is ( ) = | ( )| × ( ) . The corresponding Wigner delay time is define by If the scatterer is connected to single-channel current-carrying electrodes by couplings and , then it can shown that 1 where, if is the Hamiltonian describing the isolated molecular core, = ( − ) −1 and In deriving this expression, the electrodes are assumed to be one -dimensional tight-binding chains, with nearest neighbor hopping elements -, (where > 0) with a dispersion relation = −2 cos , which relates the energy of an electron travelling along the electrode to its wave vector , where 0 ≤ ≤ .
The group velocity of such electrons within the electrodes is therefore = = 2 sin .
In equation (13), , are the couplings between molecule and the left and right electrodes respectively and g ab is the , element of the core Green's function . Since we are interested in the contribution to the delay time from the molecular core, we shall consider the 'wide band limit', where is independent of energy in the energy range of interest, between the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) of the scattering region formed by the molecule. When H is real is real and therefore the delay time is obtained from the phase of the complex number ∆ = 1 + ∆ 1 + ∆ 2 .
As an example, consider the mathematically simple ballistic limit, where the scatterer is a linear chain On the other hand, we are interested in the opposite limit of a scatterer, which is weakly coupled to the leads, such that ≪ 1 and ≪ 1 and transport is off-resonance, such that the energy lies within the HOMO-LUMO gap. (The case of on-resonance transport is discussed in appendix 1.) In this case, β ≪ , and ≪ 1 so the delay time reduces to This equation shows that the total delay time is a sum of independent times due to each contact.
where we have defined an intrinsic core delay time to be: which is independent of the coupling to the leads. Since = ∑ , this yields Since the local density of states is given by . This demonstrates that is proportional to the local density of states at atom of the isolated molecule.
In the case where the couplings to the leads ( and ) are identical, then the ratio of delay times corresponding to connectivities , and , is This delay time ratio is a property of the core Green's function alone. It is interesting to note that as illustrated by all the above examples, in the weak coupling limit, the delay time is always positive.   a-e of Figure 7, we calculate the maximum and minimum delay times for each core as a function of the number of rings. For structures shown in Figure 7a-e, Figure 9 shows the maximum and minimum of the Wigner delay times, corresponding to the connectivities marked red and blue respectively. For example, in Figure 7b, for naphthalene, the maximum delay time is corresponds to atoms number 2, 4, 9, 7 and atoms 3 and 8 have the minimum value   Table 6: Maximum and minimum core delay times for the molecules of figure 7.
The above behavior is clearly reflected in the local density of states of the molecules, shown in Figure   10Figure 10.

DISCUSSION
To understand how superconducting proximity effects manifest themselves in molecular-scale junctions, one needs to reconcile the vastly different energy and length scales associated with superconductivity and molecular-scale transport. For the former, typical superconducting coherence lengths are on the scale of hundreds of nanometres and energy scales on the order of meV, whereas for the latter, molecular lengths are typically a few nanometres and energies on the scale of 1-5 eV. To investigate this question experimentally, there is also a need to identify signatures of the interplay between molecularscale transport and superconductivity, which are resilient to the sample-to-sample fluctuations arising from variability in the atomic-scale contacts between molecules and electrodes. In this article we have addressed this issue by investigating the connectivity dependence of superconducting proximity effects in various molecular structures connected to one or two superconducting contacts. We found that under certain conditions (weak coupling, locality, connectivity, mid-gap transport, phase coherence and connectivity-independent statistics) the electrical transport properties of the molecular junctions can be well described by a magic-number theory, which focusses on the connectivity between the individual sites of the molecules. For normal-molecule-superconducting junctions, for example, we find that the electrical conductance is proportional to the fourth power of their magic numbers, whereas for molecular Josephson junctions, the critical current is proportional to the square of their magic numbers. We also studied interference effects in three-terminal Andreev interferometers, where the interference pattern was driven by the superconducting phase difference. Our analytical predictions were in good agreement with the performed numerical simulations for all the studied systems.
We also investigated the connectivity dependence of Wigner delay times. At first sight, it seems unreasonable that the core Green's function and corresponding magic number table can yield information about delay times, because in the absence of a magnetic field, the core Hamiltonian and corresponding Green's function = ( − ) −1 are real, whereas delay times are associated with the phase of the complex transmission amplitude. Nevertheless we have demonstrated that delay time ratios can be obtained from the core Green's function or equivalently from the associated magic number tables. (4)

Calculation of the Wigner delay times for on-resonance transport
In the absence of external gating, electron transport through molecules under ambient conditions is usually off-resonance. On the other hand, if a molecule is gated such that the energy of electrons passing through the molecule is close to an energy level of the molecule, then in principle transport could be on resonance. To illustrate the properties of in this limit, we now examine a number of examples, under the condition that the scatterer is weakly coupled to the leads, such that ≪ 1 and ≪ 1.

Example 1
In the limit β ≪ , equation (14) of the main text reduces to Example 1a. As an example, consider the case where In this case, = 0 and In this expression, For ≪ Γ ab , this yields ≫ 1/Γ ab , which means that the electron spends a long time on the pendant orbital.
In this case, which reflects the fact that vanishes at the Fano resonance, where = .
Derivation of equations (4) and (6) of the main text.

S4
As discussed in [1], the low-bias, low-temperature electrical conductance of a molecule connected to normal-metal electrode and one or more superconducting electrodes is given by Here is the Andreev reflection coefficient describing the probability that an electron from the normal electrode traverses the molecule, Andreev reflects as a hole and then exits through the normal electrode. This process is accompanied by a Cooper pair entering the superconducting electrode and therefore carries a current. The electrons and holes of such a structure are described by the Bogoliubov de Gennes Hamiltonian In this equation, Δe iϕ is the superconducting order parameter, with uniform phase , which is non-zero in the superconducting electrode only and provides the coupling between electrons and holes. Conceptually, this doubling of the Hamiltonian, means that the pi system of the central core should be views as comprising both electronic and hole degrees of freedom, as shown in the figure below In the absence of the electrodes, the mid-gap Greens function of the electronic system is = −ℎ −1 , whereas the Green's function of the hole system is ℎ = ℎ −1 . When the superconductor is weakly attached, the Green's function matrix element connecting the electronic pi orbital at to the hole pi orbital at is, to lowest order in the coupling,

Derivation of equation (5) of the main text
According to Ref. [2], at zero temperature the critical current through a phase biased Josephson junction flowing through lead R (see Figure 4 in the main text) can be calculated by:   ( 3 )) .

S7
In the derivation of Eq. (16) we also took advantage of the block-diagonal structure of the Green's function of the central molecule, namely = ⅆ ( , ℎ ) = 3 ⊗ , where 3 is the third Pauli matrix.
Usually the main contribution to the critical current originates from the Andreev bound states with energy located within the superconducting gap (0 > > −Δ). In this energy regime the Green's function of the central molecule can be considered to be energy independent, and the above expression can be further simplified to: ( 3 )) .
Hence, we arrive at