Counterintuitive issues in the charge transport through molecular junctions

Whether at phenomenological or microscopic levels, most theoretical approaches to charge transport through molecular junctions postulate or attempt to justify microscopically the existence of a dominant molecular orbital (MO). Within such single level descriptions, experimental current-voltage I-V curves are sometimes/often analyzed by using analytical formulas expressing the current as a cubic expansion in terms of the applied voltage V, and relate possible V-driven shifts of the level energy offset relative to the metallic Fermi energy \varepsilon_{0} to an asymmetry of molecule-electrode couplings or to an asymmetric location of the"center of gravity"of the MO with respect to electrodes. In this paper, we present results demonstrating the failure of these intuitive expectations. For example, we show how typical data processing based on cubic expansions yields a value of \varepsilon_0 underestimated by a typical factor of about two. When compared to theoretical results of DFT approaches, which typically underestimate the HOMO-LUMO gap by a similar factor, this may create the false impression of"agreement"with experiments in situations where this is actually not the case. Further, such cubic expansions yield model parameter values dependent on the bias range width employed for fitting, which is unacceptable physically. Finally, we present an example demonstrating that, counter-intuitively, the bias-induced change in the energy of an MO located much closer to an electrode can occur in a direction that is opposite to the change in the Fermi energy of that electrode. This is contrary to what one expects based on a"lever rule"argument, according to which the MO"feels"the local value of the electric potential, which is assumed to vary linearly across the junction and is closer to the potential of the closer electrode.


Introduction
In spite of significant advances [1][2][3][4][5][6][7][8][9][10][11][12][13] , charge transport across molecular junctions continues to remain a nonequilibrium problem difficult to understand 14,15 . Resorting to a single (Newns-Anderson) model 12,[16][17][18][19][20][21][22][23][24][25][26][27] to describe the transport within a picture assuming the existence of a dominant molecular orbital is a common procedure in the field, also allowing to rationalize more sophisticated microscopic transport calculations 28 . In fact, this single-level picture turned out to excellently explain a series of transport measurements beyond the ohmic bias range [29][30][31] and to back the model parameters extracted from fitting experimental data with high-level quantum chemical calculations 25,27,32 . In the present work, we consider two issues related to the analysis of the transport data within this framework: (i) Typical current-voltage I − V characteristics measured a Theoretische Chemie, Universität Heidelberg, Im Neuenheimer Feld 229, D-69120 Heidelberg, Germany. ‡ E-mail: ioan.baldea@pci.uni-heidelberg.de. Also at National Institute for Lasers, Plasmas, and Radiation Physics, Institute of Space Sciences, Bucharest-Mȃgurele, Romania in molecular junctions are featureless curves. An example is depicted in Fig. 1 below. Due to their general appearance, although a general analytic formula I = I(V ) is available from literature 18,23,28,33,34 , using third-order expansions instead of the exact expression (see eqn (2) below) appears to be a convenient and reasonable simplification and was used in earlier studies 12,24 .
Such third-order expansions are inspired by studies on a variety of macroscopic and mesoscopic junctions up to biases of current experimental interest, wherein it was considered a prominent characteristic of transport via tunneling [35][36][37][38][39][40] .
(ii) Like those shown in Fig. 1, experimental I − V curves may exhibit a more or less pronounced asymmetry upon bias polarity reversal I(−V ) = −I(V ), which is particularly desirable for achieving current rectification using molecular devices 13 . The most common way to embody this asymmetry in analytic transport approaches is either to relate it to asymmetric molecule-electrode couplings 24,28 or to assume that the ("center of gravity" of the) dominant molecular orbital is located asymmetrically relative to the two (say, "substrate" s and "tip" t) electrodes 21,33 .
The analysis presented below will demonstrate that, although the aforementioned assumptions seem to be justified intuitively, they are in fact of rather limited applicability. Examples will be presented showing cases where the opposite is true.

Results and discussion
To provide the reader with a convenient reference when reading the text that follows, the definition of the variables utilized below is given in Table 1. Table 1 List of the main variables utilized in the present paper.

Symbol
Meaning MO (dominant) molecular orbital ε 0 (V ) MO energy offset under applied bias (V = 0) ε 0 ≡ |ε 0 (V )| V =0 MO energy offset without applied bias or for junctions with symmetric I −V curves γ (dimensionless) Stark effect strength MO broadening functions due to coupling to electrodes s ("substrate") and t ("tip") Γ a arithmetic average of Γ s and Γ t Γ g geometric average of Γ s and Γ t Γ h harmonic average of Γ s and Γ t V t ≡ V p transition voltage (alias peak voltage, cf. ref. 41 and 15) for junctions with symmetric I −V curves transition voltage for positive/negative biases; transition voltage for positive or negative bias within the cubic approximation δ asymmetry of the MO-electrode couplings (0 < δ < 1) G low bias conductance G 0 conductance quantum µ s,t electrodes' Fermi energy G f it , ε f it 0 conductance and MO energy offset deduced by fitting using cubic expansions I vs. V V b bias range used for fitting ε f it MO energy offset deduced by fitting I −V curves obtained experimentally of via eqn (9) using cubic expansions in the bias range −V b < V < V b δ ε 0 bias-driven MO shift z molecular axis d molecular length z 0 MO "center of gravity" location E z electric field along the molecular axis δ ε 0 (z) bias-driven MO shift expected according to the "lever rule" (cf. ref. 20)

Basic working equations
By assuming electrode bandwidths much larger than all other characteristic energies (wide band limit), the current mediated by a single, possibly bias-dependent energy level having an energy offset ε 0 (V ) relative to the electrodes equilibrium Fermi energy can be written in a compact form The above formula results by recasting the more familiar expression 18,23,28,33,34 with the aid of the trigonometric identity Above, the biased ("substrate" s and "tip" t) electrodes are assumed to have Fermi energies µ s,t = ±eV /2, N is the effective number of molecules in junction, and G 0 = 2e 2 /h is the conductance quantum. Γ h,a,g stand for the harmonic, arithmetic, and geometric averages of the level broadening functions Γ s,t arising from the couplings between molecule and electrodes. In the zero-bias limit (V → 0), the conductance G has the form By assuming a bias-independent level energy offset, the third-order expansion O (V ) 3 in terms of V of eqn (2) or (1) reads In the off-resonance limit (Γ a ≪ |ε 0 |) which characterizes the vast majority of experimental situations, the above expression acquires the form An applied bias V can shift the energy of the dominant orbital, ε 0 → ε 0 (V ). By assuming a linear dependence ε 0 (V ) = ε 0 + γeV (6) a series of experiments could be successfully analyzed. In this case, the counterpart of the third-order expansions of eqn (4) and (5) read (7) and respectively. In off-resonance cases (Γ a ≪ |ε 0 ± eV /2|), an expression for the current not limited to low-order expansions in V can be deduced from eqn (2) 29 To end this section, we briefly refer to a quantity useful for the subsequent analysis, namely the transition voltage V t , defined as the bias at the minimum of the Fowler-Nordheim quantity log |I|/V 2 , or the equivalent peak voltage V p (≡ V t ), defined as the bias at the maximum of V 2 /|I|. The latter has been recently introduced 15,41 to emphasize that no mechanistic transition (e.g., from direct tunneling to field-emission tunneling, as initially claimed 42 In off-resonant situations described by eqn (9) and (5) or (8), simple expressions for the transition (peak) voltages for both polarities (V 1 t V 2 t < 0) can be derived analytically 15,29 Eqn (9) ⇒ eV 1,2 We checked that the off-resonance limit applies in all the cases presented below. Therefore, using the simplified eqn (9) instead of eqn (2) or (1), and eqn (5) or (8) instead of eqn (4) or (7) is legitimate.

Exact Newns-Anderson description versus cubic expansion
By fitting experimental I −V -curves using eqn (9) and (8), the values of the fitting parameters ε 0 , γ, and G entering these equations can be deduced. The analysis presented in this subsection will reveal surprising differences between the values estimated with the aid of these equations. The black symbols of Fig. 1 depict a typical, moderately asymmetric (I(−V ) = −I(V )) curve measured in molecular junctions 43 . Fitting the experimental data shown in Fig. 1 with the aid of eqn (9) yields a curve (red line in Fig. 1) in virtually perfect agreement with experiment. Although not so "perfect" as the red line, the green curve, obtained by fitting the experimental data using eqn (8), is in fact quite satisfactory. Still, much more importantly than the quality of the two fits, the two fitting procedures yield substantially different parameter values (see caption of Fig. 1). Particularly noteworthy is the fact that the cubic approximation drastically underestimates the HOMO-energy offset ε 0 , which represents ∼ 60% from the exact estimate. (Let us mention that in the junctions considered conduction is mediated by the highest occupied molecular orbital (HOMO) 43 ).
What is wrong with the green curve of Fig. 1 is the fact that the cubic expansion leading to eqn (8) is a legitimate approximation of eqn (2) or (1), or eqn (9) (because in the present off-resonant limit (Γ a ≪ |ε 0 |, see caption of Fig. 1) eqn (2) or (1) reduce to eqn (9)) only at sufficiently low biases. If this were the case, the differences between the I − V curves computed using eqn (8) (green curve in Fig. 1) and eqn (9) (blue curve in Fig. 1) at the same parameter values would be negligible. However, the inspection of Fig. 1 clearly reveals that this is not the case. The cubic expansion does not hold in the whole experimental V -range; it is legitimate only up to biases |V | ∼ 0.3 V, wherein the differences between the green and blue curves are small.
Using experimental transport data from ref. 43  As visible in the figure, the blue dashed line, obtained by using eqn (9) and the parameters corresponding to the green line, substantially deviates from the green curve, which would not be the case if the cubic approximation of eqn (8) were justified. The red curve could hardly be distinguished within the drawing accuracy from that obtained using eqn (1) or (2), because of the small values Γ s ≈ Γ t ∼ 10 −2 |ε 0 | deduced from eqn (3) with N ≈ 100 43 . The negative ε 0 -values reflect the HOMO-mediated conduction in the junctions considered 43 .
1-10 | 3 its inability to account for the experimental fact that the transition ( Fig. 2a) or, alternatively, the peak ( Fig. 2b) voltage spectra have minima or maxima located asymmetrically (V t,+ = −V t,− ). This result, based on experimental measurements, gives further support to a similar finding emerging from a theoretical simulation presented recently 15 .
Typical transport measurements on molecular junctions sample bias ranges slightly exceeding the transition voltage (V > ∼ V t ) 42,44-46 . In Fig. 3 we present results of a numerical simulation. There, the (red) curve has been computed using eqn (9) at γ = 0 (i.e., for a curve I(−V ) = −I(V ) symmetric about origin V = 0) for such a bias range (−1.25V t < V < 1.25V t ) along with the (green) curve obtained by fitting the red curve by means of eqn (5). Similar to Fig. 1, the quality of the fit is very good; nevertheless, the fitting parameter ε f it 0 (63% of the actual value ε 0 ) is drastically underestimated. For convenience, the results presented in Fig. 3 are presented in dimensionless variables obtained by using the "natural" bias and current units V t and I t ≡ I (V t ) introduced recently 41 .
Like those presented in Fig. 1, the results presented 3a reemphasize why attempting to fit transport data in bias ranges sampled in typical experiments by using the cubic approximation represents an inadequate procedure: this V -range is beyond the applicability of the cubic expansion. If the cubic expansion were legitimate, differences between curves computed via eqn (5) or (8) and the (practically) exact eqn (9) using the same parameters deduced via cubic fitting would be negligible. (In Fig. 1 the differences between the red curve and experimental data 43 are insignificant, so here we could refer to the results computed via eqn (9) as the "experimental" results.) The comparison between solid green lines (cubic fitting) and the blue dashed lines (exact equation + parameters from cubic fitting) of Fig. 1 and 3a, which depict the two aforementioned curves, shows that the opposite is true. These differences are small at low biases only; the contributions of the higher order terms neglected in eqn (5) and (8)) is witnessed as significant differences between the green and blue curves at higher biases.
In this vein, one can attempt to employ cubic expansions for data fitting in narrower bias ranges, where higher order terms are indeed negligible. Simulations of this kind are presented in Fig. 3b, 4a and 4b. For simplicity, in these figures we have chosen the case γ = 0 (I(−V ) = −I(V )), so fittings in a bias range −V b < V < V b and 0 < V < V b yield identical results. In Fig. 4a and 4b, we present results obtained from fitting by means of the cubic approximation of eqn (5) symmetric I − V curves (corresponding to γ = 0) computed using the "exact" eqn (9) (which mimic the "experimental" curves in these simulations) within bias ranges |V | < 0.6V t and |V | < 0.8V t , respectively. Differences (∼ 5% and 10%, respectively) between the exact values of MO energy offsets (ε 0 ) and those (ε f it 0 ) estimated in this way (ε f it 0 ≃ 0.95 ε 0 and ε f it 0 ≃ 0.90 ε 0 , respectively) are reasonable and comparable to experimental inaccuracies (e.g., ref. 27). The examples depicted in Fig. 4a and 4b pass the self-consistency test: the differences between the fitting curves (green curves, cubic expansions) and the (blue) curves computed via the exact eqn (9) with the parameter values deduced by fitting (parameters of the green curves) are reasonably small in the whole bias range |V | < 1.25V t shown (which mimics the bias range of experimental interest). should be sufficiently narrow (as is the case in Fig. 4). This may seem unexpected: in principle, a better fit may be expected when more data are sampled. In fact, this is surprising only at the first sight; the data to be fitted here are ideal data resulting from computations via eqn (9) that mimic (and actually very accurately reproduce) measurements (as visible in Fig. 1 and 2, or elsewhere 25,27,29,30,51 ), but are not affected by (statistical or measurement 41 Fig. 3b is related to the possibility to accurately "detect" slight deviations from linearity in the data to be fitted. This poses no problem in cases where "ideal" data not affected by errors are used (like those utilized to generate Fig. 3b). However accurate eqn (9) is, in contrast to the situation analyzed in Fig. 3b, real I −V measurements are affected by inherent experimental errors, and data in a sufficiently broad V b -range are needed for reliable fitting. The noise of experimental I −V -curves, which is the typical situation for STM setups 12,24,48 , acts detrimentally when too narrow bias ranges are employed for fitting.
The unacceptably strong dependence of the fitting parameters on the bias range encountered above for cases where γ = 0 becomes even more problematic in the case of asymmetric curves (I(−V ) = −I(V ), γ = 0). To illustrate this fact, in Fig. 5 we depict results showing the dependence on the bias range (V b ) of the γ-parameter obtained by fitting the experimental I −V curve 43 with the aid of the "exact" eqn (9) and the cubic approximation, eqn (8). They are shown as green and blue symbols, respectively. The difference between the two sets is obvious; while the V b -dependence of γ-values obtained via eqn (9) is insignificant, γ-values obtained via eqn (8) vary by a factor ∼ 3. It is worth noting in this context that accurate estimates of the parameter γ are needed to adequately describe the current asymmetry upon bias polarity reversal ("current rectification"). Fig. 1 presents a situation where the cubic approximation seems reasonable for one bias polarity (small differences between the green and blue curves for V < 0) but is totally unsatisfactory for the other bias polarity (V > 0). Cases  Fig. 1 via eqn (8) and (9) (blue and green symbols, respectively) within bias ranges is the variable on the abscissa. Notice the scattering of the blue symbols around average (much more pronounced than those of the green symbols, which are within experimental errors), which reflects the fact that the cubic approximation, eqn (8), represents an unsatisfactory description: model parameter values should not depend on how broad is the bias range employed for fitting. of (inadequate) methods able to describe one bias polarity while failing for the opposite bias polarity have been presented earlier; see Strong dependencies of the model parameters obtained by fitting using cubic expansions similar to eqn (7) have been found earlier in ref. 24 and ascribed to a limited applicability of the single-level model. Certainly, such limitations cannot be ruled out in some cases. However, the present investigation suggests a different possibility: the single-level description may apply (eqn (2) or (9)) while the related cubic approximations (eqn (7) or (8)) fail because are employed for too broad V -ranges where terms beyond the third order are important.

Bias-driven molecular orbital energy shift. I. State of the art
The most common view of current rectification is that the applied bias yields an energy shift of the dominant molecular orbital according to eqn (6).
To describe this bias-driven shift, a series of studies resort to a simplification 24,28,52 ; namely, they relate the asymmetric shift of the molecular orbital energy γ to the molecule-electrode couplings Γ s,t While this procedure reduces the number of fitting parameters, one should be aware that this is an ad hoc hypothesis without microscopic justification 52 . The fact that current rectification is not (necessarily 53,54 ) a result of the asymmetry of molecule-electrode couplings (Γ s = Γ t ) and, contrary to what eqn (12) claims, γ should be treated as a model parameter independent of Γ s,t has been emphasized earlier in a series of works 30,[53][54][55][56][57] . If eqn (12) applied, the parameter δ ≡ Γ t /Γ a (13) which quantifies the asymmetry of the molecule-electrode couplings, and the parameter γ would depend on each other 58 Values of the parameters γ and δ have been estimated by quantitatively analyzing various experimental data measured under STM platforms 53,56 ; the values found there do not satisfy eqn (14).
Another category of works ascribed the bias-driven shift of the energy level to an asymmetric location of the relevant molecular orbital in the space between electrodes. In this picture, the potential V (z) is assumed to drop linearly between electrodes. By assuming that the left contact located at z = −d/2 has the potential V (−d/2) = +V /2 and the right contact located at z = +d/2 has the potential V (+d/2) = −V /2 (Fig. 6b), the potential profile across the junction can be expressed as Therefore, to lowest order (δ ε 0 = O(V )) the energy correction for an MO having its "center of gravity" at z = z 0 is Formulated in words the "lever rule" 20 of eqn (16) expresses the fact that, upon applied bias, the MO energy changes according to the change in the local value of the electric potential, which is assumed to vary linearly across the junction and is closer to the potential of the closer electrode. For the case depicted in Fig. 6b, the MO center of gravity is closer to the right electrode (z 0 > 0); the MO is shifted upwards, following the upward change (+eV /2) in the Fermi energy of the right electrode, by an amount determined by the MO fractional position z 0 /d.

Bias-driven molecular orbital energy shift. II. A counterintuitive example backed by quantum chemical calculations
It is worth emphasizing that the "lever rule" (schematically depicted in Fig. 6b), which justifies the term of potential profile asymmetry or voltage division factor 21 used for the parameter γ, assumes a linear potential drop across the junction. In some cases, data could be quantitatively analyzed within this picture validated by inserting in a controlled way spacers in the molecules embedded in junctions 10 . However, by assuming a linear potential profile, screening effects are neglected. To demonstrate that a negligible screening is a fact that can by no means be taken for granted, we present below (Fig. 6) results for the HOMO energy of the alkanethiol molecule CH 3 (CH 2 ) 7 SH placed in an external field E z along the (z-)molecular axis.
These results have been obtained via genuine ab initio quantum chemical calculations (OVGF and CCSD, vide infra). Our aim is to bring surprising aspects to experimentalists' attention when they process molecular transport data. Therefore, to make this subsection accessible to a broader audience some relevant details will be given below in order to justify why such ab initio quantum chemical calculations beyond the widely employed density functional theory (DFT) are needed and how they are performed.
DFT calculations are very useful to obtain a variety of ground state properties. For geometry optimization, such calculations based on the B3LYP functional as implemented in GAUSSIAN 09 59 have also been performed in this paper. However, as well documented 60,61 , the Kohn-Sham (KS) "orbitals" utilized in DFT calculations are not physical molecular orbitals. Less problematic conceptually is the HOMO. If one knew the exact exchange-correlation functional (the key DFT quantity), the KS-HOMO energy would correspond to the lowest ionization energy 62 . However, for typical molecules used to fabricate molecular junctions this is not the case; even the HOMO energy is poorly described within the DFT 14 .
To avoid this issue, here we present results for the HOMO energies obtained via genuine ab initio quantum chemical calculations based on the outer valence Green's function (OVGF) method 63,64 . The OVGF method is a diagrammatic many-body approach 65 , wherein the self-energy entering the electronic Dyson equation includes full (i) second-and (ii) third-order terms of electron-electron interaction. Moreover, (iii) it is augmented by a geometrical approximation (physically associated to a screening factor) to also partially include fourth-and higher-order corrections 66 . The ionization energies are determined from the poles of the Green's function computed in this way. In Fig. 6c the labels 2P, 3P, and OVGF refer to the lowest ionization energy with reversed sign (HOMO energy) corresponding to the methods denoted above by (i), (ii), and (iii), respectively.
To obtain the dependence on the electric field E z applied along the molecule of the HOMO energy (Fig. 6c), calculations at the OVGF/6-311++g(d,p) level of theory have been performed. Such quantum chemical calculations are known to be accurate not only for medium size molecules like the presently considered alkanethiol molecule but also for larger molecules (like C 60 67 ). As expected 14 , differences between the OVGF HOMO energies and the Kohn-Sham "energies" (also shown in Fig. 6c) are very large. Differences between the OVGF energies and the Hartree-Fock (HF) values and those obtained within the second-order pole approximation (2P) are also significant, while the third-order (3P) pole approximation appears to be accurate in this case.
We applied the OVGF method because this is the most accurate approach to compute ionization (and electron attachment relevant for LUMO-mediated conduction) energies in the presence of an external electric field implemented in existing (or, at least, to our disposal) Kohn-Sham (KS) HOMO energies. Notice that although the HOMO density is concentrated near the thiol group SH, it is the more distant electrode that prevails in shifting the HOMO energy. quantum chemical packages. However, the example presented below suggests that there is no practical need to resort to even more elaborate many-body approaches to compute the lowest ionization energy (which would be the HOMO energy with reversed sign if the oneparticle description applied).
CCSD (coupled-cluster (CC) singles (S) and doubles (D)) 68,69 is such an approach; it represents the state-of-the-art of quantum chemistry to treat many-electron systems of medium size molecules. The CC technique constructs multi-electron wave functions by applying exponential cluster operators on the the Hartree-Fock (HF) (molecular orbital) wave function. In the specific case of CCSD, the cluster operator is truncated to single and double excitations (in short, "singles" and "doubles"). Within the CCSD framework, the lowest ionization energy can be computed either by applying the equation of motion method (EOM-IP-CCSD) 68 or by subtracting the total CCSD ground state energies of the cationic and neutral molecular species (∆-CCSD 14 ). As a further check of the OVGF approach, we mention that the OVGF-value (9.089 eV) agrees well agreement with the values thus obtained (8.996 eV using EOM-IP-CCSD and 8.945 eV using ∆-CCSD) without applied electric field.
CCSD calculations of this work have been performed with CFOUR 70 , a package also utilized to compute the spatial distribution of the HOMO (see ref. 14 for details). Again, to avoid issues related to KS orbitals, we have calculated the natural orbital expansion of the reduced density matrix at the EOM-IP-CCSD level, as the most reliable approach to characterize the spatial distribution of the extra hole (or electron in cases of LUMO-mediated conduction). By inspecting the natural orbital expansion, we found that the extra hole is almost entirely (∼ 97%) concentrated in a single natural orbital ("HOMO"). It is this spatial distribution that is shown in Fig. 6a and b.
The most important finding of this subsection emerges from the comparison of panels b and c of Fig. 6: the OVGF method (as well as the other methods related to it discussed above) predicts HOMO energies clearly exhibiting a trend opposite to the "lever rule" expectation. In view of this fact, namely, that cases exist, where the "lever rule" may fail, rather than voltage division factor or potential profile asymmetry, Stark effect strength 71 may be a possible, more appropriate term when referring to the parameter γ.
Although a linear dependence of ε 0 on the applied bias is often assumed in transport studies (also in electrochemical context 53,72 ), we are not aware of any quantum chemical study reporting such a result for molecules often used to fabricate molecular junctions. Therefore, we believe that the strict linearity of the dependence of the HOMO energy on the applied field/bias represents an important result of the present paper. Noteworthy, the results of the OVGF-computations (depicted by points in Fig. 6c), which perfectly lie on a straight line, correspond to electric field values up to 2 V/nm. These values safely cover the typical experimental range for molecular devices, which is in most cases up to about 1 V/nm, since beyond this value field ionization may become significant.

Conclusion
An important finding of the present paper is the demonstration that current transport data processing based on cubic expansions of the current as a function of voltage is inappropriate. First, this typically underestimates the energy offset of the dominant molecular orbital by a factor of about two. Because DFT calculations typically underestimate the HOMO-LUMO gap by a similar factor, in the light of the present finding, "agreements" between experiments and theories using Kohn-Sham orbital energies uncorrected by employing more accurate quantum chemical methods and/or image charge effects could/should be reconsidered. Second, the application of the cubic expansion for bias ranges of experimental interest (almost inherently) yields parameter values depending on how broad is the bias range employed for fitting; this may easily be interpreted as an unphysical result, creating the impression that the single level description is invalid. In reality, more plausible is that the cubic expansion rather than the single-level description is inadequate. A third drawback of the cubic expansion is its inability to quantitatively describe asymmetric I(V ) = −I(−V ) curves. This is revealed by the fact that, contrary to experiments 43 , it yields transition (peak) voltages of equal magnitude for both bias polarities (cf. Fig. 2); this is an aspect on which a theoretical simulation presented in ref. 15 already drew attention.
Another important finding reported in this paper is the fact that the bias-driven shifts of molecular orbital energies are necessarily determined neither by the asymmetry of the molecule-electrode couplings nor to the asymmetric location of the "center of gravity" of the molecular orbitals relative to the two electrodes. The latter aspect is important also because it emphasizes that even if a single orbital dominates the charge transport through a certain molecular junction, other molecular orbitals can have indirect contributions via subtle screening effects that may yield counterintuitive effects of the kind presented above. We chose to present a single (counter-)example, namely, the case of an isolated benchmark molecule (octanethiol) in external field. We could present more (counter-)examples (e.g., a fuller class of alkanethiols), but this would not add any further evidence; we do by no means claim that the "lever rule" fails in all cases. To avoid ambiguity, we considered an isolated molecule. If we presented a molecule linked to electrodes, never ending questions might arise, e.g. on the contacts' geometry (atop, bridge, hollow) or nature (chemisorption vs. physisorption). What is important for the present purpose is to show that the (upward or downward) MO shift due to an applied electric field is not necessarily directly related to the MO location.