Understanding resonant charge transport through weakly coupled single-molecule junctions

Off-resonant charge transport through molecular junctions has been extensively studied since the advent of single-molecule electronics and is now well understood within the framework of the non-interacting Landauer approach. Conversely, gaining a qualitative and quantitative understanding of the resonant transport regime has proven more elusive. Here, we study resonant charge transport through graphene-based zinc-porphyrin junctions. We experimentally demonstrate an inadequacy of non-interacting Landauer theory as well as the conventional single-mode Franck–Condon model. Instead, we model overall charge transport as a sequence of non-adiabatic electron transfers, with rates depending on both outer and inner-sphere vibrational interactions. We show that the transport properties of our molecular junctions are determined by a combination of electron–electron and electron-vibrational coupling, and are sensitive to interactions with the wider local environment. Furthermore, we assess the importance of nuclear tunnelling and examine the suitability of semi-classical Marcus theory as a description of charge transport in molecular devices.


Supplementary Note 1 Chemical structures
The chemical structures of the five molecules deposited from solution to form single-molecule transistors are given in Supplementary Figure 2. The porphyrin central unit (P) is functionalised with different groups designed to act as π-anchors, interfacing the molecule with the graphene electrodes. The synthesis of all five compounds has been reported previously. 4 (TDP) 2 P is molecule M, and is the molecule in devices A-D in the main text, and device E in Section Supplementary Note 6. Devices F-N are displayed in Section Supplementary Note 7 and contain either one of (1-pyrene) 2 P, (2-pyrene) 2 P, (TBF) 2 P or (HBC) 2 P.
As stated in the main body of this work, during the oxidation of these compounds the charge density is predominantly withdrawn from both the anchor groups and the porphyrin ring, as opposed to the aryl side groups. This is shown in Supplementary Figure 3 which shows the single-particle molecular orbitals in various charge states most likely involved in the investigated resonant charge transport (where N corresponds to the neutral molecular species).

Supplementary Note 2 Theory
In this section we outline the theoretical approach used to describe the transport properties of the molecular junctions. We describe the molecular structure within the junction as a single 'site' that is coupled to two (source S, and drain D) electrodes and is additionally interacting with a collection of (molecular and broader environmental) vibrational modes, as schematically pictured in Supplementary Figure 4. The molecular system is modelled as a single energy level with the energy µ (which depends on the gate voltage V g ). As described in the main body of this work, we assume that the vibrational modes interacting with the molecular energy level can be found in their thermal equilibrium state at all times.

Reduction and oxidation rates
As discussed in the main body of this work, we model the overall charge transport as a sequence of electron hoppings (oxidations and reductions) at the source and drain interfaces. We denote the rates of oxidation and reduction processes as γ l ox and γ l red , respectively, the molecular electrode couplings are given by Γ l , where l = S, D. The reduction and oxidation rates at the source and drain interfaces are given by: respectively, as given in the main body of this work. In the above, f l (ε) denotes the Fermi-Dirac distribution: where µ l is the chemical potential of the lead l. The parameter Ω depends on the degeneracy of the molecular electronic level and takes the value of 0 in the case of even N and Ω = 1 otherwise. To understand the physical meaning of the (2 − Ω) and (1 + Ω) factors let us consider a system which (for a given V g ) is found in an even-N charge state at zero bias. As V b is increased, the resonant transport (charging and de-charging of the molecule) becomes possible. Since the N + 1 charge state is doubly degenerate (it possesses an odd number of electrons), the reduction rate is twice what would have been in the absence of this degeneracy -both spin-up and spin-down electrons can charge the molecule. The efficiency of an electron hopping off the molecule, on the other hand, is not affected since the N charge state is singly-degenerate. The opposite is true for an odd N.
The physical meaning of Supplementary Equation 1 and Supplementary Equation 2 is as follows. They describe the effective hopping rate between a donor and an acceptor system, one of which is represented by a single (molecular) energy level with the other one comprising a continuum of energy levels. To obtain an effective hopping rate, one must therefore integrate over this continuum of energy levels (each of which can act as a donor/acceptor) weighted by the Fermi distribution [or 1 − f l (ε)] which describes the population of the above donor [acceptor] levels.

6/32
In Supplementary Equation 1, Supplementary Equation 2, k red (ε) and k ox (ε) are the molecular density of states (DOS), and are given by The above rates were derived using a second-order quantum master equation in the polaron-transformed frame. In order to go beyond the conventional Born approximation (which does not capture the lifetime broadening) the free system evolution was replaced with an effective one obtained through the equations of motion, see Ref. 6 for details. Therein, Γ l = 2π|V l | 2 ρ l where V l is the molecule-lead coupling strength (with the electrode l), ρ l is the constant density of states in the lead l (wide-band approximation), and τ −1 = (Γ S + Γ D )/2h is the lifetime of the molecular energy level. Overall, the damping term e −t/τ induces lifetime broadening in our description (the origins of which can be traced back to the uncertainty principle which relates the uncertainty in energy of a quantum state to its lifetime). Finally, B(t) is the phononic correlation function: where J (ω) is the phononic spectral density which describes the distribution of the vibrational modes weighted by the strength of the electron-vibrational coupling. It is formally defined as: where g q is the strength of the coupling between the molecular electronic energy level and the mode with frequency ω q . The phononic correlation function B(t) describes the nuclear (vibrational) dynamics accompanying the nonadiabatic electron transfer. Semi-classically, it can be understood as a time-dependent Franck-Condon factor, see Ref. 7 for a detailed discussion.

Rate Equation model
In the present case, the overall quantum master equation describing the dynamics of the system pictured in Supplementary Figure 4 can be reduced to a simpler rate equation model for the population of the empty P N and charged molecular state P N+1 : 6 We solve the above in the steady-state limit: dP N dt = dP N+1 dt = 0. This yields:

(Supplementary Equation 10)
The overall current is given by: which gives the expression given in the main body of this work:

Origin of Ω factors
To additionally justify the presence of the Ω prefactors in Eqs. (2) and (3) in the main body of this work let us consider the N/N + 1 transition in detail (where N corresponds to a non-degenerate neutral charge state with an overall spin of zero). The set of rate equations for this scenario should account for the populations of three possible states: the non-degenerate N charge state (for simplicity referred to as P 0 ), and two degenerate N + 1 charge state with an overall spin of (+1/2) and (−1/2), denoted by P ↑ and P ↓ , respectively. Then, this set of rate equations can be compactly written as: where γ is the overall reduction rate (for the spin-up or spin-down electrons): andγ is the overall oxidation rate (again for the spin-up or spin-down electrons): We note that P 0 + P ↑ + P ↓ = 1 and that in the absence of magnetic field, P ↑ = P ↓ . The steady-state populations can now be easily found: The electric current can be found at either molecule-lead interface. Considering, for instance, the source interface, the electric current is given by: which trivially leads to an expression equivalent to the one in the main body of this work (for Ω = 0).

High-temperature limit: Marcus theory
A number of simplifications to the above expressions can be made in the high-temperature limit. 6 Firstly, we note that the damping term e −t/τ in Supplementary Equation 3, Supplementary Equation 4 induces only broadening of the IV characteristics. The broadening of the Fermi distributions in the leads will have the same effect on the transport properties of the junction, and we can therefore ignore this damping factor for τ −1 k B T . Secondly, we can take the high-temperature limits in the phononic correlation function in (Supplementary Equation 5): coth (ω/(2k B T )) ≈ 2k B T /ω; sin(ωt) ≈ ωt, and cos(ωt) − 1 ≈ −ω 2 t 2 /2. This yields is the reorganisation energy. The above assumptions give the well-known Marcus expressions 8-10 for the energydependent hopping rates: as discussed in the main body of this work, see Ref. 6 for details.

Supplementary Note 3 DFT calculations
Gaussian 09 5 was used to carry out geometry optimisations and frequency calculations for different charge states of the molecules studied; B3LYP/6-31G(d) functional/basis set combination was used. From the results it is possible to calculate ab initio values for electron-vibration couplings and generate a spectral density that accounts for the inner sphere reorganisation of the molecules we have studied in nanogaps. For each vibrational mode an electron-vibration coupling constant, Λ q , can be obtained from where ω q is the vibrational frequency of mode q, and K is the Duschinsky shift vector. 11 The shift vector is obtained from (Supplementary Equation 23) where x 0 and x 0 are optimised coordinates in final and initial states, respectively. L is the (3N × 3N − 6) normal mode matrix of the initial state outputted by Gaussian09, m is a (3N × 3N) matrix that contains the masses of the atoms along the diagonal three times, and V is a (3N − 6 × 3N − 6) matrix that contains the reduced masses of each vibration along the diagonal. The purpose of m and V is to mass-weight the normal mode matrix; B is the axis-switching matrix. 12 The total inner sphere reorganisation energy, λ i , obtained by this method is given by the sum over the reorganisation energy for each mode: where S q = Λ 2 q is the Huang-Rhys parameter.

Molecular structures
For the DFT work presented here, the molecular structures have been simplified in two ways. Firstly, the 3,5bis(trihexylsilyl)phenyl groups (used to enhance the solubility of our porphyrins) have been replaced by hydrogens. Secondly, the long alkoxy chains on the anchors of (TDP) 2 P and (HBC) 2 P have also been replaced by -OH groups. The 3,5-bis(trihexylsilyl)phenyl groups are twisted to an angle of 70 • with respect to the porphyrin ring, and therefore only weak conjugation will be present between the π-systems of these groups. The electrons of the alkyl chains are localised in σ -bonds. Therefore the substitutions will have a minimal influence on the electronic properties of the porphyrins, whilst greatly reducing the computational cost of the ab initio calculations presented in the following sections. One additional step was added to the analysis -it was assumed that the anchor groups of the molecule are bound to graphene and this 'clamping' suppresses vibrations that exhibit large out-of-plane motions of the anchor groups. Therefore the vibrations that contain large out-of-plane motions of the anchor atoms (z > 0.02 Å amu 1/2 ) are filtered from the vibrational analysis. 13 The sliding energy (the change in energy when moving the molecule parallel to the substrate) is significantly lower than the change in energy associated with moving the molecule in the perpendicular direction (DFT calculations suggest the sliding energy is around 1% of the binding energy in the perpendicular direction) and therefore vibrations in which there are in-plane motions (e.g. x & y) are retained. 13,14 This filtering did not make a difference to the (pyrene-2) 2 P, (HBC) 2 P, (TBF) 2 P coupling constants, however filters some vibrations from the analysis of (pyrene-1) 2 P and (TDP) 2 P as there are changes in the angle between the π-anchors and central porphyrin unit depending upon the charge state of the molecule.

Ab initio electron-vibration coupling constants
In Supplementary Figures 5 to 11 we plot the electron-vibration coupling constants of the above DFT calculations for the molecular compounds pictured in Supplementary Figure 2. Supplementary Table 1 shows the resulting (inner-sphere) reorganisation energies.

10/32
Anchor )  TDP  67  55  Pyrene-1  88  124  Pyrene-2  72  56  TBF  27  50  HBC  132  61   Supplementary Table 1 In the (TBF) 2 P structure there is an sp 3 -hybridised carbon between the acetylene and the π-anchor, and the anchor rotates with respect to the porphyrin upon geometry optimisations of different charge states. This leads to unphysically large values of Λ q and λ i . Supplementary Figure 10 shows the HOMO/LUMO are localised upon the porphyrin moiety, therefore we assume geometric changes will also be localised to this portion of the molecule, and to calculate the electron-phonon couplings of the charge transitions of (TBF) 2 P we can replace the TBF anchors with hydrogen atoms.

Examples of vibrational modes
We finish this section by schematically plotting in Supplementary Figure 12 two relatively strongly coupled modes calculated for the N − 1/N transition of the molecule M [(TDP) 2 P]. The low-frequency mode (at 6 meV) corresponds to a collective twist of the anchor groups with respect to the porphyrin core. The higher-frequency mode is a combination of various bond stretches, atom rockings and wagings. We note, however, that the exact nature of the considered vibrational modes is inconsequential to the conclusions of this work.

The effect of bulky side-groups
In the above, to make our calculations tractable, we have ignored the bulky side groups (shown in Section Supplementary Note 1). This can be justified by the fact that, as shown in Section Supplementary Note 1, the frontier orbitals have very little electron density on the aforementioned side groups. Consequently, the charging of the molecule is not expected to alter the equilibrium nuclear positions of atoms in these bulky side-groups.

Supplementary Note 4 Substrate reorganisation energy
As discussed in the main body of this work, besides the molecular vibrational modes, the molecular electronic energy level can be also coupled to the phononic modes of the substrate. Here, we limit ourselves to a phenomenological description of these interactions, and assume that they can be accounted for the following spectral density: where λ o is the reorganisation energy for this (outer-sphere) interaction, and ω c is the cut-off frequency. The average frequency is given by ω = 4ω c . One should anticipate that the cut-off frequency an intrinsic property of the substrate while the reorganisation energy (which quantifies the molecule-substrate vibrational coupling strength) depends on the orientation and position of the zinc-porphyrin on the substrate and should be expected to vary between different devices.
The aforementioned (surface-induced) reorganisation energy can be estimated (for a given molecule-surface distance d) as: 15,16 In the above, ε r0 is the static dielectric constant of the insulator (SiO 2 ) film: ε r0 = 3.9; ε r∞ is the optical dielectric constant: ε r∞ = 2.3.Q F,M are the effective charges (fractions of the electric charge e) defined in (B8) of Ref. 16 which depend on the porphyrin-substrate separation, d. We approximate the porphyrin as a 7 × 7 Å square with a uniformly distributed charge (assuming the additional charge density is not localised on the anchor groups). Furthermore, a tot is the total width of the insulator film, a tot = 300 nm. Consequently, the second term on the right-hand side of (Supplementary Equation 26) can be largely ignored. We estimate that the porphyrin-substrate distance is roughly d ∼ 0.6 nm (as inferred from the AFM measurements of the porphyrin compounds deposited on highly oriented pyrolytic graphite, see SI of Ref. 4 ). Using formula (Supplementary Equation 26) we therefore estimate the outer-sphere reorganisation energy as λ o ∼ 200 meV.
When modelling the 77 K data, we set ω c to 25 meV which results in an excellent match of the model used here to the experimental data for all considered devices (vide infra), as well as yielding λ o values which correspond to porphyrin-substrate distances of between roughly 0.5 and 1 nm (in agreement with the AFM studies), see Section Supplementary Note 7.
We note however that the effectiveness of our model (quality of the fit of the IV characteristics) is largely unaffected by the exact choice of ω c (within the range of roughly 10 to 50 meV). In Supplementary Figure 13 we compare the experimental and fitted IV characteristics for different values of ω c . A reasonable qualitative agreement can be reached regardless of the exact value of ω c . Quantitatively, the model used here performs best for ω c ∼ 25 meV. This is further shown in Supplementary Figure 14  Furthermore, Supplementary Figure 16 shows the experimental, as well as calculated, zero-bias conductance of device A at T = 5 K. The two theoretical results were obtained with and without lifetime broadening respectively, 6,18 and using the parameter values extracted in the main body of this work. As can be seen in Supplementary Figure 16, the zero-bias conductance is overestimated by ∼ 30% in the absence of lifetime broadening.  In the above, the distance of the porphyrin from the surface (SiO 2 ) was approximated in accordance with Refs. 15

Supplementary Note 8 Alternative intermediate-temperature models
As discussed in the main body of this work, Marcus theory (MT) fails to correctly describe the transport behaviour of considered molecular junctions at the liquid-nitrogen temperature, T = 77 K. Here, we attempt to account for the empirical data with two low-temperature extensions of the conventional Marcus theory.
Firstly, we consider the Low-Temperature Corrected MT (LTC-MT) derived in Ref. 6 . It is derived from the phononic correlation function in (Supplementary Equation 5) by retaining higher-order terms while taking the high-temperature limit. The molecular DOS are given by where the low-temperature correction is introduced through the presence of additional parameter which accounts for the coupling to high-frequency vibrational modes for which the high-temperature assumption of the conventional Marcus theory is not justified. Secondly, we shall make use of Marcus-Jortner theory. 19 It separates the low-frequency environmental coupling (which is treated as in the conventional Marcus theory with the reorganisation energy λ o ) and the high-frequency molecular modes. The latter are approximated by a single effective mode with frequency ω 0 which couples to the electronic degrees of freedom with the Huang-Rhys parameter S. The molecular DOS are then given by: 19

30/32
Supplementary Figure 23 shows the experimental IV traces on resonance as well as the theoretical fits to the quantum approach (from the main body of this work), LTC-MT and Marcus-Jortner theory. LTC-MT yields comparable fits to the (more sophisticated) quantum approach. On the other hand, the Marcus-Jortner approach tends to, in general, perform worse than the above methods as it gives rise to artefacts akin to those of the usual Marcus treatment. The artefacts are visible most clearly in 5 out of 12 devices -D, F, G, I, and L. Despite this, it still performs better than the conventional Marcus approach although at the cost of two additional fitting parameters. We then use the parameters from these fits to reconstruct the full stability diagrams and compare these to experimental data, this is shown in Supplementary Figure 24. Despite reproducing the IV traces well, LTC-MT overestimates the degree of vibrationally-induced broadening of the IV characteristics off resonance, in agreement with earlier predictions. 6 . This can be seen most evidently in the stability diagrams of devices such as B, C, F, G, and H. The performance of Marcus-Jortner approach is comparable to that of the conventional Marcus theory as it again predicts early plateaus in the IV characteristics, c.f. Supplementary Figure 21. Overall, Marcus-Jortner and LTC-MT approaches both have drawbacks when fitting certain devices to IV traces and stability diagrams respectively, as described above. Therefore, by considering both Supplementary Figure 23 and Supplementary Figure 24, we conclude our quantum approach provides the most robust description of the device B-M dataset.