Important issues facing model-based approaches to tunneling transport in molecular junctions

Extensive studies on thin films indicated a generic cubic current-voltage $I-V$ dependence as a salient feature of charge transport by tunneling. A quick glance at $I-V$ data for molecular junctions suggests a qualitatively similar behavior. This would render model-based studies almost irrelevant, since, whatever the model, its parameters can always be adjusted to fit symmetric (asymmetric) $I-V$ curves characterized by two (three) expansion coefficients. Here, we systematically examine popular models based on tunneling barrier or tight-binding pictures and demonstrate that, for a quantitative description at biases of interest ($V$ slightly higher than the transition voltage $V_t$), cubic expansions do not suffice. A detailed collection of analytical formulae as well as their conditions of applicability are presented to facilitate experimentalists colleagues to process and interpret their experimental data by obtained by measuring currents in molecular junctions. We discuss in detail the limits of applicability of the various models and emphasize that uncritically adjusting model parameters to experiment may be unjustified because the values deduced in this way may fall in ranges rendering a specific model invalid or incompatible to ab initio estimates. We exemplify with the benchmark case of oligophenylene-based junctions, for which results of ab initio quantum chemical calculations are also reported. As a specific issue, we address the impact of the spatial potential profile and show that it is not notable up to biases V somewhat larger than V_t, unlike at higher biases, where it may be responsible for negative differential resistance effects.


Introduction and background
The roughly parabolic shape of the conductance G(V ) ≡ ∂ I/∂V , or the related cubic dependence of the current (I) on bias (V ), was considered a prominent characteristic of transport via tunneling. This conclusion emerged from extensive studies on a variety of macroscopic thin film junctions of oxides, insulators, superconductors up to relatively large biases. [1][2][3][4][5][6] A quick glance at I − V measurements in a variety of molecular junctions may convey the impression that this cubic dependence is satisfactory also for such systems. 7 As illustration, we have chosen in Fig. 1a a measured I −V curve, 8,72 for which fitting with a cubic polynomial looks particularly accurate. The view based of such third-order Taylor expansions (cf. eqn (1)) was able to qualitatively describe a series of interesting aspects rea Theoretische Chemie, Universität Heidelberg, Im Neuenheimer Feld 229, D-69120 Heidelberg, Germany. † Electronic supplementary information (ESI) available.
See DOI: 10.1039/C5CP02595H ‡ E-mail: ioan.baldea@pci.uni-heidelberg.de. Also at National Institute for Lasers, Plasmas, and Radiation Physics, Institute of Space Sciences, RO 077125, Bucharest-Mȃgurele, Romania lated to charge transport in molecular junctions 9 Above, G ≡ lim V →0 I/V is the low bias conductance. 10 However, if the approximation of eqn (1) were quantitatively adequate, it would be deceptive for a model-based description of transport in molecular junctions. Resorting to simple phenomenological models enabling expedient data processing and interpretation of charge transport through nanoscale devices represents common practice in molecular electronics. The validity of a certain model/mechanism for transport in a given system/device is often assessed by considering its ability to fit the measured current-voltage I −V curves. Adopting this "pragmatic" standpoint, the problem encountered is that an approximate generic parabolic shape of the conductance could be inferred at not too high biases regardless the tunneling model. 7,9 Whether the electron wave function tunnels across a structureless average medium modeled as a tunneling barrier or through tails of densities of states of one (frontier) or a few offresonant molecular orbital levels, one could "appropriately" adjust a few parameters, and the third-order expansion of the current formula I = I(V ), often available in closed analytical forms from literature, 1,3,4,6,11,12 could relatively easy be "made" to fit the usually featureless I − V curves measured. Then the quality of the fitting 1-14 | 1 alone cannot be invoked in favor of one tunneling mechanism out of other possible mechanisms underlying the various phenomenological models.
The inspection of Fig. 1a may indeed convey the impression of an overall very good agreement between experiment and eqn (1). However, a more careful analysis reveals that the range of the highest biases that can be sampled experimentally is quantitatively not so satisfactorily described. Typically, these V -values are only slightly larger than the transition voltage V t . 13 The difference in V t -values extracted from the experimental and fitting curves of Fig. 1b exceeds typical V t -experimental errors.
The range V ∼ V t is of special interest. As discussed in Section 2, V t represents an intrinsic property of a junction out of equilibrium; it characterizes the (high) bias range exhibiting significant nonlinearity. Therefore, even more than an overall fitting of I − V transport data, it is important for theory to correctly understand/reproduce V t . And, as illustrated in Fig. 1b and c, to this aim, the theoretical description should go beyond the cubic expansion framework provided by eqn (1). Emphasizing this aspect in the analysis of the various transport models considered below, which are among the most popular in the molecular electronics community, [1][2][3]5,6,11,[15][16][17] represents one main aim of the present study.
Emerging from idealized description of reality, the models utilized in transport studies, like those used elsewhere, inherently have a limited validity. The various model parameter values are subject to specific restrictions, which represent intrinsic limitations for the model in question; merely succeeding to provide good quality data fitting is meaningless if the the parameter values deduced in this way fall outside the range of model's validity. In some cases, these conditions of validity simply reflect restrictions on parameter values justifying certain mathematical approximations made to express a certain result in closed analytical form. Checking whether the various parameters deduced from data fitting do indeed satisfy the corresponding (mathematical) restrictions or are consistent with ab initio estimated values is absolutely necessary. Nevertheless, not rarely current data analysis misses this important consistency step, which may be related to an insufficient discussion in the literature of the physical background of the various phenomenological models and the limited validity of the pertaining analytical formulas. Exposing the limits of applicability of such models frequently used in molecular electronics represents another aim of the present paper. In doing that, this paper is intended as an effort to make the community more aware of these limitations. It should by no means be understood as challenging the utility of model-based studies in gaining valuable conceptual physical insight into charge transport.
Throughout, we will consider the case of zero temperature, as appropriate for off-resonant tunneling and symmetric coupling to electrodes (equal width parameters Γ L = Γ R = Γ due to "left" and "right" electrodes). In agreement with the fact that conductances of typical molecular junctions are much smaller than the quantum conductance G ≪ G 0 = 2e 2 /h = 77.48 µS, we will assume throughout Γ's much smaller than relevant (molecular orbital) energy offsets relative electrodes' Fermi energy ε B ≡ E B − E F (E F ≡ 0). We will further assume ε B > 0, because, for tunneling barriers, the corresponding formulas are simpler to write for electrons (n-type conduction); in all formulas presented below, ε B should be understood as |ε B | in cases of holes (p-type conduction) tunneling across negative barriers, as easily obtained by performing the charge conjugation transformation e → −e, m → −m, see, e.g., ref. 18). Likewise, since the tight-binding Hamiltonians of eqn (S1) and (S3) are invariant under particle-hole conjugation (e.g. ref. 19 and citations therein), the results for +ε B and −ε B coincide.
While discussing in general the aspects delineated above, we will examine the benchmark case of conducting probe atomic force microscope (CP-AFM) molecular junctions based on oligophenylene dithiols (OPDs) 8 as a specific example.

Physical properties envisaged
In this paper, we will mainly consider two physical properties (zerobias conductance G and transition voltage V t 13 ) and focus our attention on how these properties vary with the molecular length d or size N across homologous series of molecules described within schematic models. Typical for nonresonant tunneling transport is a dependence where β (β ) characterizes the exponential decay of the zero-bias conductance G with the molecular size (length). "Historically", 13 the transition voltage V t was introduced as the bias at the minimum of the Fowler-Nordheim quantity log I/V 2 ; the initial claim was that of a mechanistic transition from direct tunneling across a trapezoidal energy barrier (V < V t ) to Fowler-Nordheim (field-emission) tunneling across a triangular barrier (V > V t ). 13 As already noted in previous studies, 17,20,21 such a Fowler-Nordheim transition does not occur in molecular junctions. Physically, V t has nothing to do with the Fowler-Nordheim tunneling theory. [22][23][24] In particular, in molecular junctions the quantity log I/V 2 does not linearly decrease with 1/V at higher biases; this was a result deduced for extraction of electrons from cold metals by intense electric fields. 23 Still, transport measurements in molecular junctions do yield curves for log(I/V 2 ) exhibiting a minimum. Mathematically, the bias V = V t at the minimum of log I/V 2 coincides with the bias where the differential conductance is two times large than the nominal conductance 25 ∂ Eqn (3) can be taken as alternative mathematical definition of V t revealing at the same time its physical meaning. Using the minimum of log I/V 2 plotted a function of V (or 1/V ) 13 merely represents a mathematical trick to extract the V t -value of eqn (3). To eliminate the confusion that continues to exist in the literature on a mechanistic (Fowler-Nordheim) transition occurring at V = V t , instead of Fowler-Nordheim diagrams log I/V 2 vs. V , it might be more appropriate to use diagrams V 2 /I vs. V , which have maxima exactly at the same V = V t given by eqn (3); 14 so, instead of "transition voltage spectra" ("TVS", Fig. 1b), one can speak of a "peak voltage spectra" ("PVS", Fig. 1c).
As expressed by eqn (3), V t represents a genuine nonequilibrium property characterizing the charge transport through a molecular junction out of equilibrium. So, V t is qualitatively different from  Fig. 1 (a) Fitting the red I −V curve measured for a CP-AFM Ag/OPD1/Ag molecular junction 8 with a cubic polynomial, eqn (1), depicted by the blue line may appear to be quite satisfactory. (b) Nevertheless, for the same curves, the difference in the positions of the minimum of the Fowler-Nordheim (FN) quantity log I/V 2 defining the transition voltage V t is substantially larger than experimental uncertainties (δV t < 0.1 V 8 ). In the case illustrated here, the transition voltage obtained via eqn (1) is even beyond the bias range accessed experimentally. (c) As alternative to the FN plot, one can inspect a "peak voltage spectrum", namely the quantity V 2 /I plotted vs. V , 14 which exhibits a maximum located at V = V t . (We show only the range V > 0 because the measured I −V curve is symmetric to a very good approximation.) the zero-bias conductance G; in fact, the latter (as well as other properties commonly targeted in measurements, e.g., thermopower Seebeck coefficient) can be expressed by properties of a device (let it be molecular junction or else) at equilibrium via fluctuation-dissipation theorem. 26 Interestingly, as revealed by recent studies [27][28][29] and supported by measurements comprising thousands of molecular junctions 30,31 , it is even more justified to consider V t (rather than G) as a junction's characteristic property; in a given class of molecular junctions, V t is much less affected by inherent stochastic fluctuations than individual I −V traces or low bias conductances G. 30 The aforementioned should be taken as motivation why V t represents a quantity on which the present transport study focuses.

Models and relevant analytical formulas for symmetric I − V curves
In this section we will examine a series of models widely employed in studies on charge transport through molecular devices. For some of these models, the dependence I = I(V ) or other useful formulas can be given in closed analytical form. Worthy to be remembered, describing idealized situations, these analytical results hold only if certain conditions are satisfied. Whenever the case, the conditions of validity will be also given below. Other models examined in this paper rely on tight-binding Hamiltonians, for which the exact dependence I = I(V ) can be obtained without substantial numerical effort.
For reasons presented later, the expansion coefficients c 2 and c 4 of the current up to the fifth order, eqn (4), will also be provided for each case They can be directly inserted in eqn (5) and (6) to compute the transition voltages V t3 and V t5 within the third and fifth order expansion, Notice that the expansion coefficients can be expressed as where V 2 and V 4 have dimensions of voltages and allows one to express the expansion in terms of dimensionless voltages Next we briefly describe the models considered here and, whenever possible, give the relevant analytical formulas for the low bias conductance G (G 0 ≡ 2e 2 /h = 77.48 µS is the conductance quantum), the expansion coefficients c 2 and c 4 entering eqn (4) as well as the transition voltage V t along with the conditions for their validity.
(i) For biases e|V | < ε B , within the Simmons WKB-based approximation for electron tunneling across a tunneling barrier of effective height ε B without lateral constriction (Simmons's model), the current I is given by eqn (7). 1,5 The Taylor expansion of the RHS of eqn (7) allows one to deduce formulas for the zero-bias conductance G and the coefficients c 2 and c 4 entering eqn (4). They are expressed by 1-14 | 3 eqn (8), (10), and (11), respectively Validity condition: z ≡β d > ∼ 4 (12) Above, A stands for the junction's transverse area and α 0 ≡ 2 √ 2m 0 /h, m 0 and m * being the Sommerfeld's constant 32 , free and effective electron mass, respectively. Eqn (12) (like eqn (18) below) defines a dimensionless junction width z. The transition voltage for the Simmons and Simmons-based models has been investigated in a series of recent works. 18,20,25,[33][34][35] Eqn (7) and (8) represent mathematical approximate results; in addition to the restriction on the bias range specified, they only hold for barriers satisfying eqn (12), requiring tunneling barriers sufficiently high and wide. 15 In fact, the applicability of Simmons' approximation is more restrictive than required by eqn (12), as revealed by a more detailed analysis. 5,36,37 (ii) For electron tunneling across a tunneling barrier with lateral constriction at biases e|V | < ∼ ε B , the current I can be expressed by eqn (14), 18 which applies provided that eqn (18) is satisfied. In addition to the quantities G, c 2 and c 4 obtained by expanding the RHS of eqn (14), the transition voltage V t can also be expressed in closed analytical form; see eqn (15), (16), and (17) Validity condition: z ≡β d > ∼ 8 (18) Noteworthy, the restriction imposed by eqn (18) to the results of eqn (14) and (15) is more severe than that expressed by eqn (12), which refers to a situation less appropriate for molecular electronic devices. Like eqn (12), eqn (18) expresses the fact that the above results for model (ii) hold for physical situations where tunneling barriers are sufficiently high and wide.
(iii) The highly off-resonant sequential tunneling 11 (superexchange limit 38 ) across a molecular bridge consisting of a wire with N sites and one energy level (ε B ) per site represents a limiting case of the situation depicted in Fig. 2a. Provided that eqn (22) is satisfied, the current is expressed by eqn (19); eqn (20) and (21) follow via straightforward series expansion.
Here, t h is the hopping integral between nearest neighboring sites. Eqn (19) and (20) emerge from the Landauer formula, eqn (25), using the transmission given in the RHS of eqn (S2) for arbitrary N. Notice that the power dependent transmission expressed by the latter yields a strict proportionality 73 where the dimensionless quantity u t is the solution of the parameterfree algebraic eqn (24), which results by inserting eqn (19) into eqn (3) Eqn (22) expresses the physical fact that the above results for model (iii) hold for situations corresponding to strongly off-resonant tunneling.
(iv) A molecular chain consisting of N monomers, each monomer being characterized by a single orbital is schematically depicted in Fig. 2a. Model (iii) presented above represents a limiting case (cf. eqn (22)) of this model. To illustrate this, we give in eqn (S2) of the ESI † the transmission function computed exactly and within the sequential tunneling approximation for the case N = 2 Adapted to oligophenylene chains, ε B should be taken as a model for the HOMO energy of an isolated phenylene unit, and t h as the coupling between HOMO levels of adjacent rings. (This is the motivation for choosing the subscript h.) Under applied bias V , two limiting cases can be considered. One limit is that of an applied field completely screened out by delocalized electrons ("metallic" molecule, Fig. 2c); the potential is constant across the molecule (V r = 0), and the site energies ε r B are the same ε r B ≡ ε B . Another limit (no screening, "insulating" molecule) is that of a potential V r varying linearly from site to site between the values +V /2 and −V /2 at the two electrodes ( Fig. 2d), yielding r-dependent on-site energies ε r B = ε B − eV r .
(v) Molecular wires consisting of six-site rings with a single (πelectron) energy level per site (CH-unit) depicted in Fig. 2b represent a tight-binding description of oligophenylene molecules. The corresponding second-quantized Hamiltonian is given in the ESI †. As illustrated by Fig. 2b, t i and t stand for intra-and inter-ring hopping integrals, respectively.  Although for models (iv) and (v) the dependence I = I(V ) cannot be expressed in closed analytical form, it can easily be obtained numerically within the Landauer framework 39 Analytical formulas for the transmission function T (E) are lengthy even for Γ much smaller than ε B and t h 's, 40,41 so they are of little interest in the present context. The numerical effort to compute T (E) via the Landauer trace formula is insignificant: it consists of a numerical integration and a matrix inversion, which is needed to compute the retarded Green's function. 39 (vi) We will also examine the case of a generalized exponential transmission Results for this case, including an approximate analytical formula for the current valid for biases not considerably exceeding V t and ε B ≫ ∆ that can be deduced resorting to a Stratton-like approximation 25,42 for transmission are presented below. The conductance G of eqn (29) and the coefficients c 2,4 of eqn (31) and (32) can be straightforwardly obtained from the fifth-order expansion of the RHS of eqn (27) The particular case δ = 2, which corresponds to a Gaussian transmission, has occasionally been studied in the literature 9,17 ; results for this case are presented in the ESI †. The relation given below holds in general (not only for the Gaussian transmission δ ≡ 2) Here, V approx t, 3 and V approx t,5 represent estimates obtained from eqn (5) and (6)  has been given previously. 9 Notice that above the superscript approx refers to physical situations of strongly off-resonant tunneling, mathematically expressed by the inequality ∆ ≪ ε B .

Results for symmetric I − V curves using generic model parameter values
Detailed results for current (I) and transition voltages computed exactly and within the third-and fifth-order expansions 3 , and V t,5 , respectively) are collected in Fig. 3, S1, 4, 5, S2, S3, S4, 6, 7, S5 and 8. Below we will emphasize some main aspects of the results obtained for the corresponding numerical simulations done with generic parameter values.
The transition voltage V t,exact computed by using the "exact" current for the Simmons model [termed model (i) in the main text], eqn (7) (panel (a)), along with the approximate estimates V t, 3 and V t,5 obtained by inserting the expansion coefficients c 2 and c 4 of eqn (10) and (11) into eqn (5) and (6). 5 agrees with V t,exact within a few percents (z ≡β d is a dimensionless barrier width).
The transition voltage V t,exact computed from the current of eqn (19) corresponding to the superexchange mechanism across a molecular wire modeled as chain of sites having a single orbital of energy ε B [termed model (iii) above], along with the approximate estimates V t, 3 and V t,5 obtained by inserting the expansion coefficients c 2 and c 4 of eqn (21) into eqn (5) and (6) (panel (a)). Panel (b) shows that V t, 3 deviates from V t,exact by up to 73% (much larger than typical experimental uncertainties in V t of ∼ 10% 8,31 ), while V t,5 agrees with V t,exact within at most 14%. Notice that, within the sequential tunneling approximation, V t,exact , V t, 3 , and V t, 5 do not depend on the inter-site hopping integral t h .

The need to go beyond the parabolic conductance approximation
As visible in the aforementioned figures, the estimate V t,5 , based on the fifth-order expansion of eqn (4), represents a good approximation for the transition voltage V t,exact ≡ V t computed exactly for each model considered. For reasonably broad ranges of the model parameter values, the relative deviations V t,5 −V t /V t usually fall within typical experimental errors ( < ∼ 10%). 8,31 This does not apply to the (over)estimate V t,3 , whose deviations from the exact values are considerably larger that experimental errors. As visible in Fig. 4 and S8, V t,3 could be almost two times larger that V t,exact . This demonstrates that the cubic approximation for current (or parabolic approximation for conductance), eqn (1), can only be used for qualitative purposes.

The spatial potential profile as possible source of negative differential resistance
The impact of the spatial potential profile across a molecular junction is another issue worth to mention. Presumably, realistic potential profiles fall between the two situations depicted in panels (c) and (d) of Fig. 2 (a constant potential and a linearly varying potential, respectively); so examining the differences in the results obtained in these two cases may be taken to assess how important is to exactly know the actual potential profile in a real case. The results presented in Fig. S2, S3 and S4 indicate a weak effect at biases not much larger than V t . Obviously, this is a further plus of transition voltage spectroscopy.
The manner in which the potential varies across a molecular junction is important only at biases substantially higher than V t . This is shown in Fig. S2, S3, S4 and especially in Fig. 6, which, to the best of our knowledge, demonstrates a qualitatively new effect. As visible there, instead of a current plateau occurring for a potential flat across the molecule beyond bias values at which molecular orbital  Fig. 5 Transition voltage V t computed exactly for a molecular wire whose monomers are modeled as single sites characterized by a single orbital of energy ε B . For arbitrary nearest-neighbor hopping integral values t h , this is termed model (iv) in the main text, which becomes model (iii) in the limit t h ≪ ε B . Notice the almost insignificant change V t → V t + δV t when the spatial potential profile across the chain changes from constant (Fig. 2c) to a linear drop (Fig. 2d). The differences between the dashed and solid lines in panel (a) represent the absolute changes δV t . The relative changes δV t /V t are shown in panel (b). energies (ε 1,2,... ) become resonant to the Fermi level of one electrode (eV res = 2ε 1,2,... ) 39 , the current decreases for a linearly varying potential. We are not aware of a previous indication on such a negative differential resistance (NDR) effect related to a specific potential variation across a molecular junction. Noteworthy, this is an NDR effect for uncorrelated transport computed within the limit of wideband electrodes. Previous work on uncorrelated transport indicated finite electrode bandwidths and energy-dependent electrode density of states as possible sources of NDR. 43

Gaussian transmission versus Lorentzian transmission
The comparison between the cases of Lorentzian (obtained as a particular case of model (iv) for N = 1, discussed in detail previously, e.g., ref. 12) and Gaussian transmission (a particular case of eqn (26)) reveals that, even systems wherein the transport is dominated by a single energy level, may have qualitatively different physical properties.
As visible in Fig. 7 and S5 and eqn (S7), V t strongly depends on the width parameter ∆ and is inversely proportional to ε B . By contrast, in the case of Lorentzian transmission, V t linearly increases with ε B and is nearly independent of the width parameter Γ as long as 12,21,29 In fact, both the proportionality to ε B and the insensitivity to Γ-variations for Γ ≪ ε B is a common feature shared by the above models (iii-v), while a V t nearly proportional to ε 1/2 B at large ε B for models (i) and (ii) is the consequence of the fact that the corresponding transmissions are similar to that of eqn (26) with δ = 1/2 (cf. eqn (17) and (30)). Common features for these two types of transmissions are the inflection point of the I − V curves ("current steps", Fig. 8a) located at biases V res = 2ε B /e where the level becomes resonant to electrodes' Fermi energy, and the fact that the associated maximum in the differential conductance has a width proportional to the transmission widths; δV res ∝ ∆ (Fig. 8b), similar to δV res ∝ Γ for Lorentzian transmission.
The opposite dependence of V t on ε B for Lorentzian (V t ∝ ε B ) and Gaussian (V t ∝ 1/ε B ) transmission is relevant for the correctly inter-preting gating effects [44][45][46] or the impact of electrodes' work function. 27 Concerning the approximate estimates for V t,exact , Fig. 7b and S5b indicate a behavior similar to that found in the above cases: V t,5 represents an accurate estimate for the exact V t,exact , while V t,3 is quantitatively unsatisfactory.

Asymmetric I − V curves described within the Newns-Anderson model
Let us briefly examine the case where, unlike those assumed throughout above, the I −V curves are not symmetric under bias reversal, i.e., I(−V ) = −I(V ). It is clear that the asymmetry with respect to bias polarity reversal cannot be satisfied within the framework provided by eqn (1) and (4); even powers in V should also be added By including only the lowest even power (i.e. c 1 = 0, c 3 = 0) 9 , it is possible to account for an asymmetry I(V ) = −I(−V ). However, doing this, i.e. using eqn (35), has an important drawback. By inserting eqn (35) in eqn (3) one easily gets i.e., transition voltages of equal magnitudes for positive and negative bias polarities, and this holds no matter how large is  . 6 Tunneling current computed for molecular wires [model (iv), eqn (S1)] for a constant (Fig. 2c) and linearly varying (Fig. 2d) potential. The linear potential drop has an impact at higher biases more pronounced than that visible in Fig. S2, S3 and S4; displacements of the current step positions (which are located at on-resonance situations, eV = 2ε j , in the case of a constant potential) and negative differential effects. ε j 's ( j ≥ 1) represent the absolute values of the highest occupied orbital energies measured from the electrodes' Fermi energy. For a two-site chain, these values are ε 1,2 = ε B ∓ t h (cf. inset of panel (a)). Notice the logarithmic scale on the y-axis. relevant; this clearly demonstrates that the "generic parabolic" dependence of the conductance of eqn (35) is an unsatisfactory approximation.
To find the counterpart of eqn (6) for the asymmetric case one has to solve an algebraic fourth order equation obtained by inserting eqn (36) into eqn (3) This yields asymmetric transition voltages V t,5+ = −V t,5− . Their analytical expressions are too long and will be omitted here. Starting from the models discussed above, asymmetric I − V curves and V t,+ = −V t,− can be obtained by allowing a bias-induced energy shift ε B → ε B (V ) = ε B + γeV (39) where −1/2 < γ < 1/2 is a voltage division factor (see, e.g., ref. 12 and citations therein).
Although calculations are straightforward, being too lengthy, the counterpart of the formulas given for symmetric I − V curves in the preceding section will not be given here. Instead, we will restrict ourselves to the Newns-Anderson (NA) model 48-52 within the wide-band limit, which assumes a single level of energy ε B (V ) and Lorentzian transmission. For (off-resonant) situations of practical interest (Γ ≪ ε B , biases up to ∼ 1.5ε B /e) 29,31,46,53-57 , exact formulae for the current and V t± have been deduced 12 ; see eqn (S10) and (S11) in the ESI †. The fifth-order expansion in eqn (S10) straightforwardly yields the expansion coefficients entering eqn (36) wherec j ≡ c j (ε B /e) j ( j = 1 to 4). By inserting the above expression in eqn (37) and (38) (c) Fig. 9 Results for a single level and Lorentzian transmission. Whether computed exactly, using the third-or fifth-order expansions eqn (S10), (35), and (36), respectively) I −V curves are asymmetric with respect to the origin for a nonvanishing voltage division factor γ (panel (a)). The transition voltages (minima of the Fowler-Nordheim quantity in panel (b)) computed exactly from eqn (S11) for opposite polarities are of different magnitudes (V t+ = |V t− |). The fifth-order expansion correctly describes this inequality (V t,5+ = |V t,5− |), while the third-order expansion incorrectly predicts V t,3+ = |V t,3− |.
bias polarities can be obtained. The comparison with the exact V t± obtained from eqn (S11) is depicted in Fig. S6 and 10.
To end this subsection, lower order expansions (eqn (35) and (36)) for the asymmetric I −V curves considered above appear to be more problematic than for symmetric cases. Fig. 10 shows that even the fifth-order expansion provides a satisfactory description for both bias polarities only for weak asymmetries (|γ| < ∼ 0.1).

Case of molecular junctions based on oligophenylene dithiols
Analytical expressions for current and conductance like those given above are very often utilized by experimentalists to fit their measurements. However, and this is one important point on which the present work wants to emphasize, these formulas are only valid under well defined conditions/restrictions, and utilization beyond these conditions makes no sense. To exemplify, in this section we will analyze recent transport data obtained for molecular junctions based on oligophenylenes using silver electrodes Ag/OPDs/Ag containing up to four phenyl rings (1 ≤ N ≤ 4). 8 Experimental studies of OPD-based molecular junctions with silver electrodes containing 1 ≤ N ≤ 4 phenyl rings yielded a value of the conductance attenuation factor β ≃ 1.56 per ring; 8 for a phenyl ring size d 1 = 4.3Å, this corresponds toβ ≃ 0.36Å −1 8 . This βestimate can be invoked to immediately rule out the tunneling barrier picture underlying models (i) and (ii) discussed above as valid framework to analyze the charge transport in these junctions. Out of the phenylene dithiol species studied in ref. 8, it is only the fourmember species, d → d 4 = 4d 1 , that satisfies condition of eqn (12), β d = 6.24 > 4. In fact, as discussed in detail recently 18,25 , the Simmons model in its original formulation 1 does not account for the fact that, in molecular junctions, electron motion is laterally confined. When lateral constriction is incorporated into theory, the condition becomes more restrictive. Instead of eqn (12), one should apply eqn (18), which is invalidated even for N = 4.
The presently employed parameter ε B represents an effective Fig. 10 Results for a single level and Lorentzian transmission. Transition voltages for positive and negative biases computed exactly V t± from eqn (S11), and using the the third-or fifth-order expansions [V t,3∓ , eqn (35), and V t,5± , eqn (36), respectively]. Notice the incorrect prediction V t,3+ = |V t,3− |. The fifth-order expansion provides good estimates V t,5± ≃ V t± only for small γ; at larger γ, only one polarity is satisfactorily described (e.g., V t,5− ≈ V t− for γ > 0).
barrier height, which also embodies image effects. According to eqn (13) and (9), to "explain" a certain β -value adjusted to fit experimental data, both ε B and the the effective mass m * can empirically be "adjusted" to make theory to "agree" with experiment. However, what cannot be "manipulated" in this way is the β -value itself. The restriction expressed by eqn (18) represents the mathematical condition for a valid description within the tunneling barrier model, and it cannot be modified whatever "appropriate" is the choice of ε B and m * . The argument based on the (too small) β -value also demonstrates that a description based on the superexchange limit of tunneling, which underlines model (iii), 11 is impossible. For β = 1.56, eqn (20) yields t h /ε B = 0.458, which is not much smaller than unity, as required by eqn (22). Again, whatever "appropriate" the adjustment of the parameters ε B and t h , they cannot be chosen to satisfy the condition required by theory for a valid superexchange mechanism, because the experimental data "fix" the value of β and thence via eqn (20) the ratio t h /ε B . So, one should go beyond the superexchange limit.
The most straightforward generalization is to use model (iv) for a chain with having one level per monomer, for which theory does not impose restrictions to the ratio between ε B and t h . Within a microscopic description based on model (iv), t h represents the hopping integral coupling the two HOMOs of two adjacent benzene rings. It can be deduced from the difference between the ionization energies I 1 and I 2 of benzene and biphenyl: 2t h = I 1 − I 2 . Quantum chemical calculations based on refined ab initio methods, EOM-CCSD (equation-of-motion coupled clusters singles and doubles), as implemented in CFOUR 58 and OVGF (outer valence Green's functions) 59-61 done with GAUSSIAN 09, 62 allowed us to reliably determine I 1 and I 2 . The values are I 1 = 9.219 eV and 9.197 eV, and I 2 = 8.217 eV and 8.204 eV, respectively. They allow to compute the hopping integral needed for calculations based on model (iv): t h = 0.501 eV and 0.497 eV, respectively. The agreement between these two t h -values demonstrates the reliability of the ab initio estimates. The direct ab initio determination of ε B (relative alignment between HOMO and electrodes' Fermi level) is an issue too difficult to be addressed here. It can be determined from the experimental value V t | N=1 ≃ 1.15 V for Ag/OPD1/Ag. 8,27 The equa-tion ε B = eV t √ 3/2 12 yields ε B = 0.996 eV. Exact currents for the Hamiltonian of eqn (S1) can be computed numerically, and this allows to obtain V t for various N's. Results based on model (iv) using this ab initio estimated t h = 0.5 eV and adjusting ε B to the value ε B = 0.996 eV fixed by the experimental value V t | N=1 ≃ 1.15 V 8 are completely unacceptable; not even an exponential decay of the conductance with increasing size N can be obtained.
In view of this disagreement, we have attempted to keep only one of the above parameters (either t h = 0.5 eV or ε B = 0.996 eV) and to determine the other by fitting the experimental value β = 1.56. 8 Fixing ε B = 0.996 eV yields t h = 0.387 eV, while fixing t h = 0.5 eV requires a value ε B = 1.288 eV. The results obtained in this way are depicted in Fig. S6 and 11, respectively, and they show that none of these empirical changes can make model (iv) to provide a satisfactory description. Differences between the cases of a flat potential and a linearly varying potential are insignificant ( Fig. S6c and 11c); so, the lack of information of the real spatial potential profile across junctions cannot be advocated as possible source of this disagreement.
A description based on model (v), which considers the phenyl rings explicitly, represents the highest reasonable refinement of a tight-binding approach to charge transport in OPDs junctions. To determine, t i , which within this model is a characteristic of a ring, we use the lowest singlet excitation energy E exc ; it can be expressed as E exc = 2t i . Using aug-cc-pVDZ basis sets, we have determined E exc from EOM-CCSD 58 and SAC-CI (symmetry adapted cluster/configuration interaction) 62 calculations; the values are E exc = 5.154 eV and E exc = 5.051 eV, respectively. They are in good agreement with the experimental values: 4.9 eV 63,64 and 5.0 eV. 65 So, we use the value t i = 2.5 eV, which also agrees with the overall description of the excitation spectrum of benzene as well as with polyacetylene data 66 . To determine t, we use again the difference between the lowest ionization energies of benzene and biphenyl. (Notice that the parameter t of model (v), which is the hopping integral between neighboring C sites belonging to adjacent phenyl rings (Fig. 2b), is different from the parameter t h of model (iv), which is the hopping integral between HOMOs of adjacent phenyl rings (Fig. 2a).) To reproduce I 1 − I 2 ≃ 1 eV (see above), a value t = 3.677 eV is needed. The very fact that t = 3.677 eV is larger than t i = 2.5 eV is unphysical; the hopping integral t associated with a single C-C bond cannot  Fig. 11 Results for model (iv) using the ab initio value t h = 0.5 eV and the value ε B = 1.288 eV obtained by fitting the experimental conductance tunneling attenuation coefficient β = 1.56 (panel (a)) 8 The agreement between the calculated V t -values (panel (b)) and experiment is poor. The spatial potential profile across these junctions do not notably affect V t , as illustrated by the results obtained for constant and linearly varying potentials (panel (c)). The points (panels (b) and (c)) and error bars (panel (b)) represent experimental results for CP-AFM Ag/OPDs/Ag junctions. 8 be larger than the hopping integral t i associated to neighboring carbon atoms in an aromatic ring characterized by C-C distances shorter than of a single C-C bond (Fig. 2b). It is a clear expression that, parameters of tight-binding models (however refined they are) cannot be adjusted to satisfactorily oligophenylene molecules. This finding is in line with recent studies demonstrating limitations of tight-binding approaches for molecules of interest for molecular electronics. 67 Using the ab initio values t = 3.677 eV and t i = 2.5 eV and adjusting ε B to fit V t | N=1 ≃ 1.15 V yields ε B = 3.433 eV. This parameter set is not even able to qualitatively describe the exponential decay with N of the conductance. Therefore, we have also considered two alternative parameter sets: ε B = 6.74 eV, t i = 2.5 eV, and t = 3.677 eV (Fig. 12) and ε B = 3.433 eV, t i = 2.5 eV, and t = 0.964 eV (Fig. S7). To get the first set, instead of fitting V t | N=1 ≃ 1.15 V (as done above), we have adjusted ε B by imposing β = 1.56. To obtain the second set, we have adjusted ε B and t to fit the experimental values β = 1.56 and V t | N=1 ≃ 1.15 V. Notice that the value t = 0.964 eV is much smaller not only than the ab initio estimate t = 3.677 eV and also smaller than t i = 2.5 eV, which cannot be explained by the fact that a C-C single bond is slightly shorter than an aromatic one. 68 As visible in Fig. 12 and S7, none of these two sets provides a satisfactory description.
The (in)accuracy of V t -estimates based on the third-and fifth-order expansions of eqn (1) and (4) using parameter values specific for OPDs junctions, which is illustrated in Fig. S8, is comparable to that encountered in the previous situations (Fig. 3, S1, 4, 7 and S5).
To conclude this part, neither representing OPDs junctions as chains of phenyl rings merely considering coupled HOMO's of neighboring rings within model (iv) nor a full treatment of the coupled rings at tight binding level is able to provide an acceptable description of the transport experiments. 8 Limitations of tight binding descriptions of molecules of interest for molecular electronics have been recently pointed out. 67 The data collected in 1 demonstrate why a description based on a single level and Gaussian transmission should also be ruled out for OPDs junctions. The difference between the values ε B deduced from transport data 8 and from ultraviolet photoelectron spectroscopy (UPS) 27 for Ag/OPD1/Ag is enormous in the only case (N = 1) where the latter is available. In addition, ε B , which represents the HOMO energy offset relative to the Fermi level, is found to increase with increasing N, which is completely unphysical. 8

Additional remarks
To avoid possible misunderstandings related to the present paper, a few remarks are in order, however. The above considerations refer to existing models, which are currently used for interpret experimental data for the charge transport via tunneling in molecular junctions. Our main aim was to present a list of formulae that can be used for data processing along with the pertaining applicability conditions. This paper is intended as a working instrument enabling to check whether certain experimental transport data are compatible or not with one of the existing models. By presenting the benchmark case of OPD-based junctions, we mainly aimed at illustrating how experimental data could/should be utilized to (in)validate a certain theoretical model, not merely checking whether I − V measured curves can be fitted without checking whether the values of the fitting parameters are acceptable. An exhaustive analysis of transport data available in the literature for various molecular junctions within the presently considered models is beyond the present aim; this may make the object of a (review) paper of interest on its own.
The presently considered models disregard important physical effects (e.g., details of interface molecule-electrode couplings or metalinduced gap states). It may be entirely possible that none of these models is entirely satisfactory just because such effects are signif-  Fig. 12 Results for the conductance (panel (a)) and transition voltage (panel (b)) computed within model (v) using the values t i = 2.5 eV and t = 3.677 eV deduced ab initio and adjusting ε B = 6.74 to fit the tunneling attenuation factor β = 1.56. In panel (b), the points and error bars represent experimental results for CP-AFM Ag/OPDs/Ag junctions. 8 icant. A corresponding extension of the theoretical models to provide experimentalists with simple formulae enabling data processing/interpretation, obviating demanding microscopic transport calculations remains a desirable task for the future.
In this paper we did not discuss the (off-resonant limit (Γ ≪ ε B ) of the) Newns-Anderson model in much detail. Relevant formulae (eqn (S10) and (S11)) are given in the ESI †. As shown in recent studies 8,69-71 , cases where this description applies and is justified microscopically exist. There, the transition voltage can be used to directly extract the orbital energy offset using a formula (eV t = 2ε B / √ 3 = 1.155 eV t for γ = 0 12 ) close to that initially claimed (eV t = ε B ) 13 . However, the results presented for the various models analyzed above have demonstrated that V t cannot be uncritically used to straightforwardly deduce the energy alignment of the dominant molecular orbital. This does not diminish the importance of V t . As emphasized in Section 2, V t is an important property characterizing the nonlinear transport. By studying its behavior across homologous molecular classes (i.e., varying the molecular size d or n) or under gating (i.e., varying ε B ), the (in)applicability of a certain model can be concluded. Moreover, V t turned out to be a key quantity, as it allowed to reveal that charge transport across different experimental platforms and different molecular species exhibit a universal behavior, which can be even formulated as a law of corresponding states (LCS) free of any empirical parameters 71 .
In order to demonstrate that the cubic polynomial approximation is insufficiently accurate to describe V t , we have presented in Fig. 1a a raw I − V trace measured on a CP-AFM junction 8 . Making the point in connection with Fig. 1 is only possible by using neat, smooth curves; typical raw I −V curves measured for single-molecule junctions are blurry and therefore inadequate for this purpose. A comparison between junctions consisting of a single molecule and a bundle of molecules is of interest on its own; still, we note that, as long as the charge transfer occurs via off-resonant tunneling, the similarity of the transport properties of these two types of junctions is expressed by the aforementioned LCS 71 . Considering transport through CP-AFM junctions within the models discussed above amounts to neglect proximity effects due to other (identical) molecules in the bundle. This assumption may certainly fail in some cases. However, the fact that OPD-based CP-AFM junctions obey this LCS 71 can be taken as a strong indication that the failure of the various models in case of these molecular devices discussed in Section 6 is not related to the approximate description in terms of of a bundle of independent molecules.

Conclusion
With the manifest aim of providing experimental colleagues a comprehensive working framework enabling them to process and interpret measurements of transport by tunneling in molecular junctions, this paper have presented a detailed collection of analytical formulae, emphasizing on the fact (less discussed in the literature) that these formulae only hold if specific conditions of applicability are satisfied, which often impose severe restrictions on the model parameters. From a more general perspective, the theoretical results reported above have demonstrated that: (i) The often accepted idea of a generic parabolic V -dependence of the conductance as fingerprint of transport via tunneling emerged from studies based on high and wide energy barriers. This picture misses a microscopic foundation for tunneling across molecules characterized by discrete energy levels. At biases of experimental interest, within all the models examined here, the description based on a third-order expansion I = I(V ) of transport measurements in molecular junctions is insufficient for quantitative purposes. To illustrate, we have shown this by systematically analyzing the fifth-order expansions, which turned out to be reasonably accurate at least for biases up to the transition voltages.
(ii) Merely adjusting model parameters (e.g., within a tunneling barrier picture by claiming renormalization effects due to image charges or effective mass) does not suffice to "made" a model valid for describing a specific molecular electronic device. On one side, models are normally valid only under specific conditions, which these parameters should satisfy. On the other side, the parameter values deduced from fitting experimental data should be consistent with ab initio estimates. Parameters for tight-binding models can be easily estimated via reliable ab initio quantum chemical calculations, as illustrated in this study.
(iii) The experimentally measured values of the β -tunneling coefficient of the low bias resistance can be inferred for quickly assessing 12 | 1-14 the inapplicability of certain tunneling models. This turned out to be the case for OPDs junctions, whose (too small) value (β = 1.56 orβ = 0.36Å −1 ) deduced from experiments is incompatible with descriptions based on tunneling barrier or superexchange mechanism. Because molecular junctions based on other aromatic aromatic species have often even smaller β 's (β ≈ 0.2Å −1 27,28 ), descriptions based on those models should be excluded. In such cases, uncritical application of the mathematical formula of I vs. V is meaningless, even if the measured I −V curves can be satisfactorily fitted.
(iv) For biases not too much higher than the transition voltage, a realistic description of the spatial potential profile across a molecular junction appears considerably less important that usually claimed. This is no longer the case at high biases, where the spatial potential profile may be responsible for qualitatively new phenomena, e.g., negative differential resistance (cf. Fig. 6).
(v) The fact that the estimate V t,5 ≈ V t based on the fifth-order I(V )-expansion turned out to be acceptable in many cases can be of practical help; eqn (4) (or eqn (36)) can be employed to process noisy experimental I −V curves, for which straightforwardly redrawing measurements as log(I/V 2 ) vs. V (or V 2 /I vs. V ) may yield substantial uncertainties to determine the minimum (or maximum) position. Fitting I − V curves with fifth-order polynomials can be used to extract the coefficients c 2,3,4 and to estimate V t ≈ V t,5 by means of eqn (6) or (38).
(vi) Finally, we refer to situations where a successful description based on a single level and Lorentzian transmission (also known as the Newns-Anderson model) has been concluded 8,69,70 . The fact that the I −V curves could be very well fitted within the Newns-Anderson framework was not the only argument leading to that conclusion; equally important was that the energy offset ε B deducing from fitting the I −V data has been correlated with ab initio estimates 8,70 and/or independent experimental information. 69 Even if the transport would have been dominated by a single level, the conclusion of those studies would have not emerged in case of, e.g. a Gaussian transmission, because of the completely different dependence of V t on ε B .
With the (few) examples mentioned in the preceding paragraph, we want to end by reiterating an idea already presented in Introduction: while attempting to make the community more aware of the limits of applicability of the various models utilized, the present paper did by no means intend to challenge the overall usefulness of modelbased studies in gaining conceptual insight into the charge transport by tunneling at nanoscale.