Analyzing Molecular Current-Voltage Characteristics with the

Feb 28, 2007 - Defect Scaling with Contact Area in EGaIn-Based Junctions: Impact on ..... Neil Bennett , Gengzhao Xu , Louisa J. Esdaile , Harry L. An...
0 downloads 0 Views 426KB Size
Download PDF
J. Phys. Chem. C 2007, 111, 4431-4444

4431

Analyzing Molecular Current-Voltage Characteristics with the Simmons Tunneling Model: Scaling and Linearization Ayelet Vilan† Ilse Katz Center for Nanotechnology, Ben Gurion UniVersity, Israel ReceiVed: October 18, 2006; In Final Form: December 31, 2006

Use of the Simmons model for analyzing tunneling transport across molecular junctions is reviewed, and its inherent limitations are examined, specifically for cases where there are no molecular length-dependent data (to extract a decay parameter), be it for experimental reasons or because of changes of the molecular energetics or packing with molecular length. The potential barrier across a molecular junction is shown to be strongly bias-dependent, much more so than is assumed in the commonly used version of the Simmons model. The means to distinguish true tunneling from conduction via pinholes (or hot spots) are also considered. Power expansion to the Simmons model shows that I/V vs V2 plots should be linear over that range, thus providing a simple and standardized parameter extraction. From such plots, we can extract values for the equilibrium conductance and for the “shape factor”, a complementary parameter that describes the shape of the I-V relations. The applicability of these two parameters for describing actual transport is illustrated by analyzing data for three different types of molecular junctions. The linearity of the I/V vs V2 plots can be used to evaluate if, and if so, at which bias (and if at all) direct tunneling occurs and under which conditions this is not the case. The extracted equilibrium conductance and the “shape factor” provide an empirical method for quantifying electronic transport across practical molecular junctions where the exact packing of the molecules is rather uncertain. As such, the analysis can be used to weigh and classify the effect of chemical modifications to molecular junctions or compare contacting methods, so as to allow a deeper understanding of transport via molecular junctions.

1. Introduction Molecular electronics is a fascinating scientific frontier, which already has demonstrated proof of concept for such phenomena as rectification,1 negative differential resistance,2 and others.3 Yet, basic understanding of electronic charge transport across molecular systems is still rudimentary. One basic question is that of structure-function relations,4 viz., which of the junction’s parameters (e.g., conjugation, position, density and symmetry of energy levels, bonding to the electrodes) affect the current and to what extent? An extreme example for a seemingly minor structural change corresponding to a huge functional effect was given by Tomfohr and Sankey5 for alkyl thiols on gold. On the basis of complex band structure calculations, they showed that varying the adsorption position from a hollow lattice site to an on-top Au atom site shifts the HOMO level of the adsorbed molecule by ∼1 eV, which is expected to exponentially affect the molecular conductance.5 Although this is only one example of what is nowadays a very significant endeavor,6-9 transport modeling and simulations generally rely on extremely well-defined atomic spatial positions, far beyond our experimental accuracy. This is important because, although it is accepted that charge transport across molecular junctions is dictated by both the spatial structure (thickness, uniformity) and the energetics (energy levels and effective mass or coupling), analyses of molecular currentvoltage (I-V) curves most commonly rely on schematic, † Present address: Department of Chemical Research Support, Weizmann Institute of Science, Rehovot, 76100 Israel.

nominal molecular structures, ignoring possible structural uncertainties (predominantly pinholes). Indeed, although the example, given above, suggests that chemical control over the exact binding site can have farreaching implications for molecular electronics, it is still an elusive goal. We would like to devise and perform sets of experiments to measure current-voltage (I-V) characteristics that respond in a controllable and reproducible manner to a controlled change in conditions. However, even if this is achieved, the experimentalist still faces the problem of how to quantify the observations and assign them to specific functional properties such as energy level positions (which dictate the barrier height), carrier delocalization (or coupling between energy levels), or, for molecular monolayers, defect density, something that may well affect not only the energy levels but also the contact area for transport. The limited knowledge on the exact molecular arrangement in actual molecular junctions impedes the use of accurate theories for quantification of experimental measurements. Furthermore, reducing detailed theoretical models into few averaged characteristic parameters (e.g., barrier height or effective mass) enables decomposing experimentally measured current-voltage (I-V) curves into a few basic parameters. Therefore, approximate, effective models are at present a necessary step for quantifying the increasing amount of molecular transport data. Within the largely accepted simplifications of nonresonant tunneling and that of the WKB (semiclassical) approximation to nonresonant tunneling, one of the most ubiquitous of such parameters is the decay parameter (β)10-12 (in Å-1, see equations

10.1021/jp066846s CCC: $37.00 © 2007 American Chemical Society Published on Web 02/28/2007

4432 J. Phys. Chem. C, Vol. 111, No. 11, 2007 a1 and a2 of the Appendix).13 This parameter describes the exponential decay of the current, I, with the length of the molecule (l): I ∼ exp(-βl). A major disadvantage of the decay concept is that it requires a set of data with systematically varying tunnel barrier width (molecular length), while keeping all other molecular parameters the same. Practically, such a strategy is limited mostly to simple alkyl chains and related molecules, because, if functional groups are introduced (e.g., conjugated moieties, polar groups, redox centers), some variation of energy levels or change in packing with molecular length is nearly inevitable. Alternatively, measuring current as a function of temperature can be very informative. It can tell if transport is thermally activated (pure tunneling is temperature-independent) and, if so, what is the activation energy.14,15 This is a factor that needs to be considered in the choice of experimental measurement setup, as several types (e.g., Hg drop and electrochemical junctions) allow only very limited temperature variation. Furthermore, monolayers can undergo a liquid-solid phase transition16 which can naturally influence their charge transport characteristics. Other molecules, such as proteins, are mostly incompatible with vacuum or dry conditions required for cooling, again limiting the range of temperature variation experimentally. The question that we will address here is the following: Can we extract from experimental current-voltage characteristics functional parameters that allow meaningful comparisons between experiments, in a way that is independent of tunnel barrier thickness, i.e., without the need for experimental current-thickness relations? In molecular electronics8,9,14,17-21 Simmons’ inelastic tunneling model22-24 (see section 2 for a summary) is used often, due to its simplicity, even though, in most cases it is likely oversimplified as it ignores important contributions such as super-exchange17,21 or hopping.9,25 Fitting the Simmons model to I-V data (instead of I-l, current-thickness data) was found to be useful for evaluating the contribution of chemical modifications at bilayer interfaces17 and was shown to produce parameters in good agreement with I-l extracted parameters.14 Wang et al. also reported on negligible temperature dependence of the measured current, indicating transport by tunneling rather than hopping14,15 (on tunneling vs hopping see also ref 25). Aswal et al. distinguished between three different bias ranges (low, medium, and very high), using simplifications for each of these ranges.20 The present work focuses on the intermediate bias range. However, as noted by Cui et al.18 and Wang et al.,14,15 fitting over different voltage intervals yields quite different results. A discussion on the problematic interdependence between the different parameters controlling the Simmons relations was given in ref 26. In addition, our own experience27 shows that a visually reasonable fit to the Simmons model can be reached fairly easily, but the resulting extracted parameters are extremely unstable and sensitive to the initial choice of starting values. Such parameters can hardly be regarded as “characteristic”. This is the reason for reconsidering the limitations inherent in, and the possibilities for using, the Simmons approach (2.1). We can understand the problem of multiple fits of the Simmons relations to a data set, within which no length-dependent data are available, by using a “scaled” analysis (2.2). The idea is to separate between the magnitude of the current (2.2.3) and the shape of the I-V curve (2.2.2). The magnitude is characterized by the equilibrium conductance, G0 (in Ω-1), rather than by an absolute current value (in Ampere, at a specific bias).

Vilan The shape of the I-V curve is characterized by a parameter which we will call the “shape factor” (F), the ratio between the width and the height of the tunneling barrier. These parameters reflect the product of the height and the width of the tunnel barrier, which dictates the current magnitude, G0, (see section 2.2.3, below) and the ratio between these two, (F), which dictates the shape of the I-V curve (see section 2.2.2). Note that G0 and F are functional parameters within which are contained several physical parameters, such as the junction’s barrier height, width, and contact area. In this manuscript, we will occasionally refer to them as lumped parameters so as to contrast them to the more familiar parameters that describe the junction. Although the combination of these two parameters enables us to unravel the known coupling between barrier height and effective mass,18 we show that this is possible only if the contact area is known or if data with systematic length variation are available (i.e., at least one spatial parameter is required). Another basic limitation of the Simmons model is the fundamental assumption that the potential profile across the molecule(s) between the electrodes is linear (section 2.1). We show, for three quite different experimental data sets, that this so-called “linear model” is not applicable to real systems (4.1). We suggest that this problem can be overcome by limiting the analysis to low bias voltages, because at those voltages such a linear potential profile is mostly a reasonable approximation (2.3.1). This approach has as further benefit that a much simpler polynomial I-V relations result (2.3.2) along the lines originally suggested by Simmons24 and used later by Brinkman et al.28 Both the scaled analysis and the polynomial approximation converge to the same characteristic magnitude and shape parameters (2.3.3 and 2.3.4). The proposed analysis is tested for analysis of I-V data measured across junctions of alkyl thiols at three very different contacting configurations (3). Alkyl chains are by far the cleanest and most reproducible of all molecular junctions and thus serve as a basic model system for molecular electronics studies. In addition, tunneling is expected to dominate charge transport across alkyl chains.14 Extracting characteristic parameters from actual data is tested (4.2) and compared to traditional numerical fits (4.3). Finally, it is shown that a model which allows for a slightly nonlinear potential profile provides an excellent fit over the full bias range (4.4). The discussion section reviews the possible merits in using the lumped parameters under typical experimental conditions where the junction area29-31 and width are not well-known. In cases of a roughly known contact area (within an order of magnitude), the product β0l can be extracted from G0, and the ratio (β0l/F) provides the barrier height. Alternatively, if lengthdependent data (I-V-l) are available, we show that the intercept of ln(G0·F) vs F plot can be used to extract experimentally the contact area, A, of a molecular junction. 2. Revisiting the Simmons Model: Power Expansion, Linear Presentation, and Scaling This section starts with the formal definitions of the Simmons model and its characteristic equilibrium parameters (2.1). Then the full model is reduced to its commonly used linear form (2.2), and, finally, a polynomial approximation to the Simmons model is presented (2.3). 2.1. Basic Underlying Assumptions. The Simmons model describes inelastic, nonresonant tunneling through an insulating barrier and was formulated by J.G. Simmons in three successive papers in 1963.22-24 This model has found wide acceptance and use in the analysis of oxide tunnel barriers as well as barriers

Analyzing Molecular Current-Voltage Characteristics

J. Phys. Chem. C, Vol. 111, No. 11, 2007 4433 expression of asymmetry in the barrier height. The price is a significantly more complicated mathematical formulation. 2.1.2. Linearity of the Simmons Relations. The widely used form of the Simmons relations (equation a3 in the Appendix) represents only a special case of the general Simmons model.22 This form is based on two very stringent assumptions: 1. The potential profile is spatially linear, meaning that it has either a rectangular or a trapezoidal shape, for symmetric (Figure 1a) and asymmetric (Figure 1b) electrodes respectively. 2. The applied bias adds linearly to the potential profile. The formal form of the linearity requirement is written as follows:

ψ h)

Figure 1. Schematic illustration for different barrier types and characteristic parameters. Barrier within the linear Simmons model include the symmetric rectangular barrier (a) and the asymmetric, trapezoidal barrier (b) where the variation of potential profile is linear following eq 1 (dark and light gray for equilibrium and biased potential profiles). Equation 1 is incorrect for nonlinear cases such as image force lowering (c) or traps within the barrier (d). The averaged potential is given by φ, and its equilibrium value (φ0) for the linear cases (a and b) equals to the average of the difference in work function at the two interfaces (eq 1). Panel e illustrates the characteristics parameters of a rectangular barrier: dimensionless thickness (βl) and the equilibrium barrier height (φ0). High and low shape factors (F ) β0l/φ0) are demonstrated on the two rows of panel e (constant φ0 at each column). The width ratio (R ≡ L/l) is shown on the top-right barrier: L is the width of the potential barrier above the Fermi level, which can be narrower than the molecular length (l) due to packing tilt or more fundamentally due to, e.g., image force lowering (c) or contact resistance.

formed by molecular films, since the very early days of this field32,33 and also recently.14,17,18 The underlying assumptions of the model are essential for defining the range of conditions over which it can be applied and its limitations, viz. • the potential is spatially averaged (2.1.1) and • the potential varies linearly with space and applied voltage (2.1.2). We discuss both these assumptions next. 2.1.1. AVeraged Potential. The Simmons model is based on averaging the potential-thickness profile of the tunnel barrier, ψ(V,x) resulting in a single characteristic parameter, which is the averaged barrier potential, ψ h (V). The “Simmons view” of a potential barrier for tunneling is shown schematically in Figure 1a,b where the dark and light gray depict the barrier at equilibrium (V ≡ 0) and under arbitrary applied bias, respectively. The averaged potentials for the two cases (φ0 and φ(V), respectively) are also shown. One consequence of the averaged potential assumption is that asymmetry between the height of the barrier at the two ends (Figure 1b) almost does not play a role (as long as the bias does not exceed the smaller of the two barriers). To account for asymmetry, Brinkman et al. revised the Simmons model by integrating over transmission probabilities rather than over the potential profile,28 allowing a direct

φleft + φright (Wleft + Wright)/e + V V ) ) φ0 + (1) 2 2 2

where ψ h is the averaged potential, φ is the potential difference (barrier height) at each electrode (left and right), W is the work function, e is the electron charge, and the subscript 0 emphasizes that this is an equilibrium (V ≡ 0) property. Thus, the so-called Simmons model really is the “linear Simmons model” as it is based on the assumption that the potential profile is linear in space and responds linearly to applied bias. A nonlinear potential profile was originally included by Simmons for describing the effect of the image force.22,23 Recently Hansen and Brandbyge used a parabolic potential profile,34 as shown schematically in Figure 1c, to improve modeling of experimental data. The first assumption of averaging could be wrong for cases where traps contribute dominantly to the transport, as considered by Probst35 and schematically shown in Figure 1d. Despite these limitations, the Simmons model has the great merit of simplicity, whereas the alternatives mentioned above have generally much more complicated mathematical forms. Keeping this in mind, we will show how we can use the Simmons model to analyze I-V characteristics to yield unique parameters, while at the same time maintaining simplicity. To achieve this, we will explore two independent approaches, that of scaling (2.2) and linearization (2.3). 2.2. Scaled Nature of the Simmons Relations. 2.2.1. Modified Simmons Relations. The relation describing the tunneling current (I) as a function of the applied bias (V) originally defined by Simmons22,23 (see equation a3, Appendix) includes five unknown barrier characteristics, namely: • equilibrium barrier height (φ0) • nominal barrier width (l) • contact area (A) • effective mass (m*) • ratio of the barrier width over the molecular length (R) However, mathematically there are only three distinguishable parameters. One of these is obviously the barrier height (φ0), which scales the voltage. Several alternatives can be conceived of to define the other two parameters. As will be justified below, we suggest rewriting the original Simmons relations (see equations a3 to a5 in the Appendix) in the following manner:

I)

{(

) [ ( x )] ) [ ( x )]}

2G0φ0 V exp Fφ0 1 ‚ 1Fφ0 - 2 2φ0

(

1+

V exp Fφ0 1 2φ0

1-

V 2φ0

1+

V 2φ0

-

(2)

for |V | e φ0, where the three independent parameters are • φ0, the averaged, equilibrium (V ≡ 0) barrier height (eq 1 above)

4434 J. Phys. Chem. C, Vol. 111, No. 11, 2007

Vilan

• F, the equilibrium shape factor, defined as follows:

F)

x

4π Rl h

2m*m0e ) 1.025Rl φ0

x

m* φ0

(3)

which (see equation a1 in the appendix) implies F ≡ β0l/φ0. • G0 the equilibrium conductance, defined as follows:

G0 )

e2 A 4π Rl x2m*m0eφ0 - 2 × 4πh R2l2 h 4π exp Rl 2m*m0eφ0 (4a) h x

(

)

(

G0 ) 310

)

A (Fφ0 - 2) exp(-Fφ0) R2l2

G0 ) 325Am*

(Fφ0 - 2) F2φ0

(4b)

exp(-Fφ0) ≈ 325

Am* exp(-Fφ0) (4c) F

where h is Planck’s constant and e, m0, and m* are the electron charge, mass, and effective mass, respectively. R is the ratio between the true junction width (L) and the experimentally known molecular length (l), as schematically shown in the top right of Figure 1e (see section 2.2.4). The numerical coefficients in the short forms have the following units: 310 in Ω-1 [Å/µm2], 325 in [Ω-1 V-1 µm-2], and 1.025 (eq 3) in [Å-1 V-0.5]. Note that both F and G0 (eqs 3 and 4) are by definition equilibrium parameters as no bias term is included. Because, by definition, Fφ0 ≡ β0l, eq 2 is equivalent to equation a4 in the Appendix, the two products will be occasionally exchanged, as (β0l) is a more familiar expression than (Fφ0). 2.2.2. Shape Factor, F. The importance of the shape factor, F,36 is demonstrated schematically in Figure 1e, as the ratio between the width and the height of the barrier.37 Figure 2a illustrates the effect of varying F on simulated I-V curves for junctions with constant β0l and G0. The dotted line (labeled G0) reflects the simple Ohmic relation (I ) G0V), which provides the current “magnitude”, and all curves converge to this curve at very low bias. At higher bias, the tunneling current rises faster than the Ohmic one, and the higher F is, the steeper the curve is. The problem of using β0l and φ0 as parameters is illustrated in the inset of Figure 2a, where we see a blown-up section of three I-V curves for junctions with constant F and G0. Hence, the inset shows that for a set of curves where β0l and φ0 differ but their ratio (F) is constant, the I-V curves cannot be distinguished visually, even at strong magnification (a range of 0.1 V is shown, at intermediate bias). Varying F has a minor effect on the total G0, but this is compensated for by having the contact area (A) vary between 250 and 700 µm2 (Figure 2a main panel) However, keeping F constant and varying β0l (φ0) affects G0 exponentially. Thus, the three curves in the inset to Figure 2a represent 13 orders of magnitude change in contact area, A (this is equivalent to keeping A constant and normalizing the curves by a 13 orders of magnitude change in G0). Figure 2a also allows us to illustrate visually the effect of the term in braces of eq 2. That effect is expressed as the ratio between any of curves I-IV and the dotted line, the asymptotic Ohmic relation (G0). This behavior differs from that of the common version of the Simmons relations (equations a3 and a4 in the appendix) where the term in braces includes also a

Figure 2. Illustration of the scaled nature of the Simmons I-V relation. (a) The effect of the shape factor, F on the shape of Simmons’ tunneling I-V curves is demonstrated on simulated I-V curves of identical dimensionless thickness (β0l ) 30) and varying shape factor, F and, accordingly, varying barrier height, φ0 (numbers in brackets), i.e., F (φ0): (I) 7.5 (4); (II) 10 (3); (III) 15 (2); and (IV) 20 (1.5). For all curves, a constant equilibrium conductance was used, G0 ) 10-9 Ω-1 (meaning that the nominal contact area varied between 250 and 700 µm2). The dotted line marked as G0 gives the Ohmic relation (I ) G0V). The inset shows the strong similarity between curves for systems with constant shape factor (F ) 15) and varying barrier height (β0l). Shown are a blow-up of curve III from the main panel and curves for two systems with identical F and varying φ0, i.e., φ0 (β0l): (III) 2 (30); (V) 1(15); and (VI) 3 (45). In this case, the contact area varied between 170 nm2 to 17 cm2 to be able to keep the same G0 ) 10-9 Ω-1 (b and c) Maps showing the sensitivity of numerical fitting process to the choice of free variables. The gray level of a pixel represents the MSE (mean square error, eq a9) between a “target” I-V curve and the I-V curve, based on the parameters given in the X-Y coordinates for that pixel. The third parameter is either the correct G0 (b) or the correct A/L2 (c). The varying fitting parameters were the barrier height (Y coordinate) and either the shape factor or dimensionless thickness (X coordinate) for panels b and c, respectively. The choice of specific coordinates is explained in the text. Maps of other combinations of coordinates appear in the Suppl. Mater. (S2).50 The “correct” minimal point is the central pixel. The gray levels are linear with log(MSE) and normalized to the span of MSE, with white and black equal to maximal and minimal MSE respectively. The fitting is derived by the strength of the gradient along the variables’ coordinates. Thus, a direction of low gradient (constant gray level) means large tolerance for variation in this parameter. “Target” I-V data were calculated from equation a4 (Simmons’ model) using φ0 ) 2 V, β0l ) 15, and A/L2 ) 106.

Analyzing Molecular Current-Voltage Characteristics strong contribution from the current magnitude, because the argument of the exponent in eq 2 will be between 0 and 1, in contrast to that of eq a4 which is close to β0l (and similarly for a3). Thus, the name “shape factor” reflects both the “aspect ratio” of the barrier (cf. Figure 1e) and the shape of the I-V curve (cf. Figure 2.a). Furthermore, considering that to a first approximation only the potential barrier is affected by the applied bias and the kinetic energy of the electrons is rather constant, it becomes clear that F, instead of the more common β0 parameter, is the natural parameter for describing the biasdependent current transport. In contrast, β0 is the characteristic length-dependent parameter, dictating primarily the magnitude of the current. 2.2.3. Equilibrium Conductance, G0. The equilibrium conductance (G0) is the first derivative of the current with respect to voltage near 0 V (G0 ) dI/dV|V f0, see section 2.3.2). As such it includes information on the current magnitude, as demonstrated by the dotted asymptotic line in Figure 2a (simulated I-V’s) and in Figure 3 (experimentally measured I-V’s, see section 4 for details). The use of G0 as a characteristic parameter is fairly common,12,18,21,26,38 though it was not related to analyses within the Simmons model. The uniqueness of the G0 concept within the Simmons model is that, in principle, it allows reducing the number of numerically fitted parameters, by extracting G0 (from a linear fit of the low-bias range) and then using G0 as a fixed parameter in a subsequent numerical fit over the full bias range. This idea was tested by simulations shown in Figure 2b,c, which are maps illustrating numerical curve fitting to the Simmons model (eqs 2 or a4). First, a fictitious I-V data set was constructed from eqs 2 or a4 with arbitrary junction parameters. We then calculated the current (for the same voltage array) for a large matrix of φ0-F values (with eq 2, assuming that G0 is known, Figure 2b) or of φ0-βl values (with eq a4, assuming that the contact area, A, is known Figure 2c), and the mean squared error (MSE; see appendix A4) between the new current array and the target current array was calculated. These MSE values are shown as gray levels, where the φ0-F (φ0-βl) matrices are used as X-Y coordinates. Figure 2b shows that for a fixed G0 the numerical fit is very sensitive to the shape factor parameter (the dark valley along F ) 7.5 in Figure 2b) and rather insensitive to the barrier height (the weak vertical contrast along the F ) 7.5 valley in Figure 2b; vertical contrast reflects variation in φ0. Figure 2c illustrates the common fitting approach,14,17,18,21 where the contact area (A) is used as a known quantity. In that case, the convergence of the numerical search (Figure 2c) is much sharper than for fixed G0 (Figure 2b), with a shallow valley along a constant β0l value.39 The reason for this difference between parts b and c of Figure 2 is that in Figure 2c we optimize two properties: the current magnitude (dictated by β0l) and the shape of the curve (dictated by φ0). In contrast, in Figure 2b, the magnitude of the current is almost completely accounted for by the known G0, and the numerical search looks mostly for the shape of the curve (up to the very fine differences as in the inset to Figure 2a). As can be seen, F is by far the dominant factor, and for each φ0, a complementary A value can be found to give the same (constant) G0 (see eqs 4b and 4c, A varies exponentially with φ0, if F is fixed). The choice of the independent parameter on the X axis for the two maps (Figure 2b, c) comes from this insight, and was also verified by constructing maps for other possible permutations (see S2 of the Supporting Information), producing a

J. Phys. Chem. C, Vol. 111, No. 11, 2007 4435 diagonal valley (cf. orthogonal one in Figure 2b,c), which is an indication for correlation between the two fitting parameters. To summarize, Figure 2b shows that, although F is a very robust characteristic of the Simmons I-V relations (for known G0), φ0 becomes well defined only if we know the contact area (A). Naturally, the question that remains is how well we know the true contact area (see section 5.3).29-31 2.2.4. Width Ratio, R, and EffectiVe Mass, m*. Historically, there was a direct translation from β0 (or F in our present terminology) to the barrier height, as β0 included only universal constants and the barrier height. However, attempts to fit molecular I-V data with the Simmons model always encountered the problem that the experimentally measured currents were larger than those expected theoretically.17,18,21 As can be seen from substituting eq 4 into eq 2, the current scales exponentially as well as inversely with the product Fφ0 (≡β0l), so that the high observed currents correspond to a reduced value for β0l. However, φ0 cannot vary much as it scales the bias (see eq 2) or dictates the shape of the curve. To compensate between these contradicting requirements to fit to the experimental magnitude and shape of the I-V curve, further parameters were introduced in β0 (F). Holmlin et al. used an asymmetry parameter (R) to yield a reduced decay parameter.17 However, the work of Simmons22,23 does not provide justification for such asymmetry effect (see eq 1 above). An effective mass, smaller than the free electron mass that appears in the original Simmons model, was introduced by Cui et al.18 They extracted a value of m* ∼ 0.16, in reasonable agreement with theoretical calculations that were based on the spectrum of the electronic states of the molecular wire.40 Another possible cause for a reduced F value is the uncertainty regarding the actual width of the potential profile (L, shown schematically in Figure 1c,e). Although L is the relevant tunneling property, it can differ from the experimentally known parameter, which is the molecular length (l), possibly modified by the tilt (Figure 1e), if that value is actually known. This is why we have interpreted R as a ratio of widths, R ) L/l. Such reduction in effective width can be due to image force narrowing of the potential22,23 (Figure 1c), to high contact resistance, which considerably alters the shape of the potential profile,7 or to significant molecule-electrode orbital hybridization at the interfaces.41 Through-space, rather than through-bond, transport can also reduce the actual width (L), compared to the nominal molecular length (l) for tilted molecules10 (Figure 1e). The presence of oxide layer on the electrode20 or of residual water (a distinct possibility if the top contact is evaporated on a cooled substrate, in the absence of a separate cold finger in the evaporation chamber that is at a temperature, well below that of the cooled sample) could lead to an actual thickness larger than the nominal one. Because of their different physical meanings, both R and m* are used here, but it is clear that they are coupled and lumped together in F. The two parameters contribute differently only to the pre-exponent in G0 (eq 4b vs eq 4c), which is expected to be negligible compared to the exponential term. Thus, whether to interpret the data in terms of R or m* is a matter of choice, as they cannot be distinguished based on analyses of experimental I-V data in the frame of the Simmons model. In summary, introducing the G0 and F as parameters allows to fit experimental I-V data in such a way that one can separate the magnitude of the current from the shape of the current-voltage curve. The importance of this parameter couple emerges also from the linearized analysis, as is discussed next.

4436 J. Phys. Chem. C, Vol. 111, No. 11, 2007

Vilan

2.3. Linearization. 2.3.1. Assuming a Linear Potential Profile near Equilibrium is Justified. Section 2.1.2 above emphasized that the common “Simmons model” refers only to linear potential profiles (eq 1, Figure 1a,b). Similarly, using β2-V plots (for β extracted from thickness dependent data, see, e.g., ref 18 and equation a2) relies on the same linearity assumption (eq 1). It is, however, very reasonable to assume that molecules (especially incompletely conjugated molecules, the ones of most interest for analysis with tunneling models) develop nonlinear potential profiles across them (such as those of, e.g., Figure 1c,d), if only because of their varying atomic composition or the spatial localization of the molecular orbitals. Nevertheless, even nonlinear profiles can be approximated with the linear model over small intervals, such as near equilibrium. Namely, we assume the potential varies linearly with the applied bias (eq 1) and with negligible bias effect on the net shape of the profile. Thus, linear (common) Simmons model should properly be used only near equilibrium where, for the small biases, quasiequilibrium conditions can be assumed, instead of forcing it over a large bias range, where its applicability will, to say the least, be questionable. Accepting this logic, we can, with the small signal assumption, replace the highly nonlinear eq 2 with its power expansion.24,28 2.3.2. Polynomial Approximation to the Linear Simmons Model. In a paper24 lesser known than his later ones, Simmons showed that eq a3 can be rewritten using a polynomial expansion (with V/2φ0 as the expansion variable), to give

I = G0(V + CV3 + PV5 + ‚‚‚)

(5)

Equation 5 can be further rearranged (ignoring the fifth order term) to get a linearized presentation of experimental I-V data

I = G0 (1 + CV2) V

(6)

Thus, plotting I/V (the integral conductance) vs the square of the bias will yield two functional parameters, G0 and C (see eq 8.b below). Equation 6 is very similar to the parabolic expression developed by Brinkman et al.28 but with two significant differences, namely that eq 6 gives the integral conductance (in contrast to the differential one in ref 28) and the absence of a linear term in eq 6.42 The I/V vs V2 relations (eq 6) are expected to be linear, as long as (1) the error in the power expansion is small; (2) the averaged potential decreases as V/2 (linearity in bias, eq 1); (3) the (low) bias has a negligible effect on the width of the barrier (linear profile, see section 4.4); and (4) the asymmetry between the barriers at the two sides of the junction is negligible.28 2.3.3. Relating Polynomial Coefficients to Barrier Parameters. The coefficients of eq 5 are directly related to those of eq 2. The equilibrium conductance, G0, is identical (eqs 4). The cubic coefficient (C) is closely related to the shape factor, F

C)

(

)

β0l + 3 F2 F2 ‚ 1≈ 96 96 β0l(β0l - 2)

(7a)

where, for convenience, the dimensionless thickness, β0l, is used, instead of Fφ0. The expansion was continued until the fifth order (P) coefficient, attempting to extract more independent information (see section 4.3)

P)

-1 -2 -3 F4 β0l/15 - (β0l) - 3(β0l) - 3(β0l) ‚ ≈ 0.3C2 2048 β 0l - 2 (7b)

However, practically, it was very difficult to extract a robust P coefficient even from simulations. First, it is almost a direct function of C (see approximated term in eq 7b, right-hand side), and second, at small bias, the fifth order contribution is very small, and at large bias, it is incorrect to approximate the original expression (eq 2) by a power series (eq 5). The approximated relations (right-hand side of eqs 7) follow Simmons24 and Brinkman et al.,28 where it was assumed that β0l . 2, which allows one to neglect all powers of β0l apart from the highest one. Similarly, with this approximation also the (Fφ0 - 2) term in the G0 expression (eq 4) reduces to Fφ0. The β0l . 2 approximation is questionable for molecules, where ∼0.5 < β0 < ∼1 Å-1 and ∼15 < l < ∼30 Å are typical values, which yield a β0l product that is not much larger than 2. Nevertheless, the β0l . 2 approximation has a major merit that it allows to extract F without knowing β0l (i.e., directly from C), something that requires knowing the contact area, A (see eqs a6 and a7). We will show later (section 4.2) that practically the β0l . 2 approximation does not introduce large errors. Equation 5 can be solved analytically, in contrast to eqs 2 or a3-a5 that can be solved only numerically. However, the accuracy of the approximation (eq 5) deteriorates as V f φ0, and the fifth order parameter (P) is rather useless. Therefore, from here onward, we shall consider only the graphical extraction of G0 and C from plots following eq 6. The polynomial coefficients can then be translated into junction parameters (e.g., barrier height and effective mass or width ratio) by combining eqs 4, 7a, and 7b (for a full glossary see S1 of the Supporting Information). 2.3.4. Linearized Approach Combined with Thickness Variation. Even though the main purpose of our analysis is to find ways to analyze I-V data from experiments without thicknessdependent information, we will now briefly consider how the linearized analyses discussed here can also be useful for data sets that include length information. The terms for G0 can be related to plots of G0 vs length (l).18,21,38 Substituting for Fφ0 ≡ β0l in eq 4, neglecting the “-2” term and rearranging we get

(

ln(G0l) ≈ ln 310

)

Aβ - βl R2

ln(G0F) ≈ ln(325Am*) - φ0F

(8a) (8b)

where, as noted above, the units of the numerical factor 310 is [Ω-1(Å/µm)2] and that of 325 is [Ω-1 V-1 µm-2]. Thus, for a system of varying barrier thickness, the decay parameter (β0) can be extracted from the slope of a plot of ln(G0l) vs the molecular length, l, as is well-known. Yet eq 8 provides also an interpretation for the pre-exponential term (or the intercept) of the semilogarithmic fit, which can be used to extract the effective electrical contact area for transport, A scaled by either β0/R2 or m*. In eq 8b, the independent parameter is now the extracted one, F, rather than the assumed length, l. A linear ln(G0F) vs F plot indicates a constant barrier height that can be extracted from the slope. The effective electrical contact area can be extracted from the intercept as above. If the plot is not linear, it could indicate that the barrier height varies between the different junctions. Clearly eq 8b could be useful for cases where there is uncertainty about the exact barrier width (Rl). It also allows

Analyzing Molecular Current-Voltage Characteristics

J. Phys. Chem. C, Vol. 111, No. 11, 2007 4437 test case since these molecular junctions are well established and the results of our test analyses can then be compared with those obtained in a more conventional manner, from the length dependence. The first type of molecular junction was formed by a single dithioloctane (HS-(CH2)8-SH), bound chemically on one side to a Au substrate and on the other side to a Au nanoparticle. The molecule is embedded in a monolayer of monothioloctane, i.e., all of the other molecules are bound only to the Au substrate. Transport was measured with a conductive probe atomic force microscope (CP-AFM; Figure 2a of ref 18). The second type of junction contained a monolayer of a several thousands of monothiol dodecane molecules, (HS(CH2)12H), chemically bound to a Au electrode, in a nanopore configuration, with no chemical bond to the other Au electrode, that was evaporated onto the nanopore structure (Figure 5 of ref 14). The last junction was formed by a bilayer of ∼108 monothiol docosane molecules (HS(CH2)22H), with the bottom layer bound chemically to lithographically patterned gold electrodes and the top layer bound chemically to a micrometer-sized floated gold flake.27 Further details on that system are given in the Supporting Information (S4) and will also be part of a separate publication. Graphical fitting, numerical searches, and simulations were done using MATLAB software. The specific functions employed are given in Appendix A4. The voltage array for simulated I-V data ranged from -φ0 to +φ0, with steps of 1 and 10 mV below and above φ0/10, respectively.

Figure 3. Curve fitting to typical molecular I-V data (a) single alkyldithiol;18 (b) alkylthiol monolayer;14 and alkylthiol bilayer.27 Dots represent measured data (diluted for clarity) and curves are different fits to the Simmons model, including (I) full-range numerical fitting with two free variables (gray line) (eq 2); (II) quasi-equilibrium approximation (dotted line, eq 6; see Figure 4); (III) bias dependent barrier width (green line, eq 9; see Figure 5); and (0) linear (Ohmic) fit over 10 mV (dashed line). Because the measured interval for (b) was 20 mV, the linear fit was done over the (50 mV range. The inset to each panel shows a blow-up of the (30 mV range, except for (b) where the inset covers the (200 mV range. The inset tables show the fitting parameters for each curve: equilibrium conductance (G0 in nS); barrier height (BH in V); shape factor (SF, in V-1); and decay parameter (beta in Å-1). The width ratio (R) was unity for curves I and II, and equals the following cubic polynomial for curve III: (a) R ) 1.00 0.03|V| + 0.23V2; (b) R ) 1.00 - 0.05|V| + 0.23V2; (c) RV>0 ) 1.00 - 0.03V + 0.31V2, and RV 1 can also be interpreted as R > 1, meaning a thicker barrier than what was assumed in the analysis (although a narrower barrier than what is assumed from physical measurements and modeling is the more common scenario). To overcome the need to know the contact area, we can try to extract more information from the data (i.e., extract three parameters), using a numerical fit to the original eq 2, as described next. 4.3. Equilibrium Approach, Compared to Numerical Fits over a Large Bias Range. The gray curves (I) of Figure 3 represent the common fitting approach,14,17 based on a (presumably) known contact area (A), with φ0 and β0l (or m*)39 as two free fitting parameters. This procedure also uses the full measurement range for fitting.44 It gives a β0 value that is rather similar to the one that is extracted at quasiequilibrium, reflecting mainly the magnitude of the current. However, the barrier height, extracted from fitting the data over the large bias range, is larger than that extracted from quasi equilibrium data, expressing the much slower increase in current at high bias (gray curves, I) than that expected, based on the quasiequilibrium fit (dotted curves, II). The small change in β0 between the low and full measured bias ranges suggests that β0 can be considered as a constant, dictated by the current magnitude (G0, provided that A is known). Thus, the fitting is mainly one of the barrier height, and, instead of letting β0 vary freely it can be fixed to its value extracted from G0 (extracted graphically from the data in a preceding step). We tested such fitting procedure (see Tables s2-s4 of the Supporting Information) and found only minor differences between the parameters, extracted by optimizing solely on φ0 or on both φ0 and β0l. We conclude that there is no need for a numerical fitting approach to find β0 as it is a direct function of G0 and A (cf. eqs a6 and a7).39 Alternatively, a known G0 can be introduced into the fitting search, combined with φ0 and A as free fitting parameters (where β0l ) f(G0, φ0, A)), to extract the contact area from the I-V data. However, as shown in Figure 2b, in such case, the fitting becomes very insensitive to the barrier height, which translates into changes of orders of magnitude in the (now) extracted contact area. For the single molecule junction such procedure yielded a contact area of only 4 Å2 (cf. nominal alkyl thiol footprint of 22 Å2, see Table s3 of the Supporting Information),

J. Phys. Chem. C, Vol. 111, No. 11, 2007 4439 which we consider as a reasonable accuracy. However, for the monolayer junction, the extracted contact area was ∼3 nm2 (cf. nominal value of 1600 nm2), which could be taken as an indication for defect-mediated transport. For the bilayer junction, the numerical search completely diverged. We take this result to reflect the inability of this fitting procedure, in this case, to match the shape of the curve, if the current magnitude is dictated by the value of G0, derived from the low bias fit. Interestingly, despite this huge divergence, it was still possible to extract a reasonable value for the shape factor, which again emerges as a robust parameter. However, we find that, by limiting the bias range over which we fit the data, we can extract a set of parameters with a contact area of the order of a few defects in the layer. Low A values are accompanied by low β0 and φ0 values, because G0 cannot change much since it reflects the current magnitude (see eq 4). Considering which of the various possible extracted parameters is the most characteristic one, it clearly emerges that G0 is the most robust one. The variation in F increases if high bias ranges are included, but considerably improves if only low bias ranges are taken. The reproducibility in β0 was similar to that in G0, but the conventional characteristic junction parameters (i.e., the nonlumped parameters φ0, m*, A) showed much poorer reproducibility. As noted before, extending the voltage range toward high bias could produce erroneous parameters simply because the linearity assumption does not hold there. The next section considers the possibility of recovering the bias dependence of the potential profile from experimental I-V characteristics. 4.4. Potential Profile across a Molecular Junction is Nonlinear. Information on the bias variation of the barrier characteristics is of major interest18 and commonly expressed as bias variation of the decay parameter, β (see eq a2, in contrast to the equilibrium β0 of eq a1). Engelkes et al. found a bias-independent β parameter but a strongly bias-dependent residual resistance (the reciprocal of the pre-exponential term in eq 4).26 As described in section 2.1, the commonly used version of the Simmons model assumes that the shape of the barrier is linear, instead of a more realistic parabolic shape22,23,34 (see Figure 1, panels a and b vs c and d). In such a case, both the averaged potential (ψ h ) and the barrier thickness (R(V)l) are expected to vary with the applied bias (see eq a5). This variation will be different from the linear dependence of ψ h predicted by eq 1. We will show now that accounting for the bias variation of the barrier width (R(V)l) can produce a good fit to the data over the full measured bias range, as shown by the green curves in Figure 3 (an equivalent approach is to use ψ h (V) as a biasdependent parameter). The potential width-bias (R-V) relations are extracted from the I-V data by assuming that all of the other parameters of eq 2 are independent of bias and equal their equilibrium values as extracted from Figure 4. This is done by defining a bias-dependent shape factor F(V) ) R(V)F(0), where,for simplicity, we assume R(0) ≡ 1. Substituting F(V) for F into eq 2 yields

I)

2G0φ0 R(V) (Fφ0 - 2) 2

[ (

{(

‚ 1-

)

V × 2φ0

x

exp Fφ0 1 - R(V)

1-

[ (

V 2φ0

)] (

)

- 1+

V × 2φ0

x

V 2φ0

exp Fφ0 1 - R(V)

1+

)]}

(9)

4440 J. Phys. Chem. C, Vol. 111, No. 11, 2007

Vilan Similar parabolic R-V curves were found for the other two data sets, confirming the general behavior of a bias-dependent potential profile. The green curves (III) in Figure 3 were calculated from eq 9 using G0 and F values of curve II. However, the low barrier height extracted for curves II could not be used for fitting the whole bias range (note that both eqs 2 and 9 are limited to V e φ0). In view of the insensitivity to the exact barrier height, shown in Figure 5a, we arbitrarily set φ0 ) 1V to allow fitting the whole bias range. Figure 5 is thus a phenomenological proof that fitting molecular I-V data over large bias ranges44 must consider nonlinear potential profiles. Naturally, a meaningful parameter extraction should rely on a physical model rather than use the empirical approach employed in eq 9. A physical model could be the detailed image force model by Simmons22,23 or the parabolic model by Hansen and Brandbyge.34 Adding G0 as a fitting constraint (see section 4.3, above) can be useful in decreasing the large number of fitting variables in these advanced models. We will now discuss how we can use the quasiequilibrium Simmons model, rather than the commonly used version, for analyzing molecular I-V characteristics. 5. Discussion

Figure 5. Bias variation of the width ratio (R) plotted for the monolayer junction14 for different equilibrium barrier height (φ0 ) 1, 1.6, and 2 eV, see curve labels) and constant shape factor (F ) 10.1) (a); or constant dimensionless thickness (β0l ) 16.16) (b). At each I-V point, the width ratio was extracted by numerically solving eq 9 using G0 ) 4.67 nS, and other parameters specified in the figure.

Thus, at each I-V data point, R(V) can be extracted by solving (numerically) eq 9 for R(V). The results for such calculation for the monolayer junction (Figure 3b) are shown in Figure 5. The extraction procedure was repeated for several sets of equilibrium parameters (reflecting a possible uncertainty in β0l due to the fact that the actual contact area is not known, see eq a7). Figure 5 compares the shape of the R-V relations obtained for three different barrier height values and either the same shape factor (Figure 5a) or the same dimensionless thickness (Figure 5b). The figure shows that the R parameter varies by ∼6% with the applied bias, which can be fit very well to a quadratic parabola (on each side of |V| ) 0). Interestingly, this 6% variation is equivalent to the error that Simmons predicted in his model (eq a5) if the integration correction factor (termed there β, yet equivalent to R in the symbols used here and nowadays) is set equal to 1.22 Varying φ0 with a constant F, changes only the scaling of the curve but not its shape (Figure 5a). In contrast, varying φ0 for fixed β0l drastically changes the shape of the R-V curve (Figure 5b). The strong increase of R away from equilibrium (for φ0 ) 1; Figure 5b), which reflects a significant widening of the barrier, appears to be an artifact due to an incorrect choice of β0l and φ0. This type of sensitivity to F (≡β0l/φ0) further demonstrates the importance of F as the characteristic parameter in bias-dependent analyses.

5.1. Why Use the Simmons Model? We conclude from the results presented in section 4 that formally the (linear) Simmons model cannot describe real molecular current-voltage characteristics. As reasons we note the lack of robustness of the extracted barrier parameters and the inability to find biasindependent parameters to cover all bias ranges, the low, nearequilibrium, and the high bias one, far from equilibrium44 (Figure 3). This bias dependence (Figure 5) explains the inadequacy of the linear Simmons model in describing charge transfer via molecular junctions. This conclusion is a fairly general one, as it was found for three very different molecular junctions. Because, at present, there is no convenient adequate alternative, we need to find ways to use the model, while being aware of its limitation, rather than forego using it altogether and we will do so in the following sections: • The deviations from linearity at high bias can be taken care of largely by limiting the range of fitting as will be discussed in the next section (5.2). • The critical issue of contact area will be addressed in the following section (5.3). • In the last two sections of this chapter, we suggest ways to analyze the experimental transport data within the linear Simmons model by using G0 and F as phenomenological parameters (5.4 and 5.5). 5.2. Simplifying Parameter Extraction by Linear Presentation of the Experimental Data. Presenting molecular I-V data in the form of Figure 4, rather than as numerical fits to eq 2 (cf. curves I in Figure 3), enables a simple visible evaluation of the relevance of the model to a given experimental data set. This is in clear contrast to the multiple possible solutions to eq 2 shown in Figure 3. Hence, presenting molecular I-V data as I/V vs V2 allows determining the bias range where the Simmons model is relevant or identifying a change in conduction behavior, as was evident for the bilayer junction (Figure 4.c). Recently, Beebe et al. used a linearized presentation for FowlerNordheim field emission transport to evaluate the onset of this transport mechanism.43 The linear I-V presentation, suggested here, completes this approach by covering the intermediate range, between the linear behavior at very low bias and the

Analyzing Molecular Current-Voltage Characteristics Fowler-Nordheim relations at very high bias, as used earlier by Aswal et al.20 A transition to Fowler-Nordheim tunneling will appear as curving down of the linear trend in a presentation such as that shown in Figure 4. Applying this presentation to the experimental data-sets tested here, reveals no transition to field emission for the monolayer14 and bilayer27 junctions. The single molecule junction18 does show a possible turning point at V ∼ 0.4 V, which is just outside the scale of Figure 4a (V2 ∼ 0.16). Indeed the data of Figure 4a start deviating below the linear trend line at ∼0.06 V2. The mathematical approximation (eq 5) is scientifically justified as it follows directly from the assumption of a linear potential profile (eq 1). If we limit the voltage range that we use to extract parameters to low bias quasiequilibrium conditions,44 we can compare consistently and in a standardized fashion between different experiments. However, such limitation leads to loss of important information regarding the bias variation of the potential barrier. At the same time, extending the linear assumption to bias ranges where the barrier depends nonlinearly on the bias merely increases the inaccuracy of extracted parameters (even if it appears to provide a visually good fit; cf. curves I in Figure 3, see main figures vs insets). Although the analysis of section 4.4 and Figure 5 is a qualitative one, it shows clearly that the potential profile over a molecule or over molecular layers is nonlinear, and not only the barrier height but also the barrier shape varies with applied bias. Therefore, analysis over a bias range beyond the nearequilibrium range44 must account for nonlinearity, by either using existing models or developing new ones. Such models are naturally more complicated than the one expressed by eq 2 and require extraction of the parameters numerically. The concept of separating the magnitude issue (using G0) from the shape of the curve (using numerical search), as suggested here can simplify the application of these models. 5.3. Do We Really Know the Contact Area? The ambiguity regarding whether the charge transport is attenuated by the length of the molecules or by the reduced contact area for transport (e.g., via defects in the monolayer) has been a worry for molecular (or monolayer) electronics since its early days29 and is relevant also for other fields that involve charge tunneling, such as spintronics.30,31 Remarkably, though, the issue is not considered in most molecular electronics publications, which use the nominal (geometrical) area of the electrical contact. Recently, Engelkes et al. tried to identify the actual contact area in conducting AFM measurements of Au/decanethiol/Au junctions. They found that the measured zero resistance (equivalent to reciprocal G0, see eq 4) was insensitive to the diameter of the AFM tip, suggesting that the current is confined to the narrowest path between a rough protrusion of the substrate electrode and the very end of the tip.46 The critical effect of the contact area on the extracted barrier parameters is clear from Figure 2c, which shows the sharp minimum if the data are fitted, using a known contact area, in contrast to the shallow valley that we get if we fit, using a known equilibrium conductance (Figure 2b). Allowing the contact area to be a free fitting parameter for the monolayer and bilayer (Tables s3 and s4 of the Supporting Information) results mostly in values that are much smaller than the nominal one. The view of a junction that is controlled by pinholes/defects is consistent with the low barrier height (approximately tens of mV) that we extract for fits that yield small contact areas, as such small barriers are reasonable for two metals at very near proximity (note, however, that the magnitude of φ0 and A decrease

J. Phys. Chem. C, Vol. 111, No. 11, 2007 4441 (increase) together if G0 is kept constant, as explained in section 2.2.2 and inset to Figure 2a). This suggests that the result obtained, using the contact area as a free parameter, has physical value. Nevertheless, our confidence in the “pinholes” results is not better than in the “nominal contact area” results, as there is not much difference in the accuracy of the two numerical fitting approaches. This is in contrast to the work of Dorneles et al., on Al/Al2O3/Al junctions where the accuracy of the fit improved drastically if the contact area was allowed to vary (i.e., three free fitting variables).31 The difference in achieving a good fit could be explained by that the behavior of the inorganic insulator fits the Simmons model closer than does that of molecular barriers. Possible deviations from the Simmons model are the major problem, if one wants to distinguish between transport via pinholes and tunneling across all of the area occupied by the molecules, solely based on a single I-V data set. Thus, it is critical to have an independent determination of the actual contact area, which can be done in cases were also lengthdependent data are available, as described in section 5.5 below. Even in cases where we have direct evidence for tunneling from temperature-independent currents or inelastic electron tunneling spectroscopy (IETS) that shows characteristic molecular vibronic modes involved in the charge transport,15,25,47 we cannot assume that the geometric (nominal) contact area is the effective one for electronic transport. This was demonstrated clearly for Al/Al2O3/Al junctions where the effective area for electronic transport was found to be 5 orders of magnitude smaller than the nominal contact area. This was explained by arguing that the current flux was limited to “hot-spots” in the insulator film. Transport across a “hot-spots” is still by tunneling but the tunneling length is now much shorter, resulting in exponentially higher conductance, and thus, most of the current will flow there.31 5.4. Use of G0 and F as Phenomenological Parameters. Figure 2, panels a and b, show that G0 and F, rather than β and φ0, are the fundamental, independent properties characterizing the Simmons tunneling model in terms of scaled relations (magnitude and shape). Moreover, these parameters have a natural presentation of data (Figure 4) allowing their linear (analytical) extraction, without knowing any other parameters. Extraction of the approximated F (i.e., F ) x96C, from eq a4, with C defined in eqs 5-7) bypasses the contact area ambiguity, discussed in the previous section (as it is independent of β0l, which necessarily depends on A). Wold and Frisbie12 already suggested G0 as a phenomenological characteristic parameter for comparing net current magnitude. We propose that the shape factor, F can be used similarly for comparative purposes, to indicate how fast the current increases away from equilibrium. The physical meaning of the equilibrium conductance (G0) is straightforward as it scales exponentially with the effective tunneling length or dimensionless thickness (β0l). The historical origin for allowing F36 as a free parameter (equivalent to R or m* as a free parameter) instead of being a fixed function of φ0 and l is that the measured currents were higher than expected or the barrier was narrower than expected.14,17,18 Clearly, transport across a barrier is affected by both the barrier’s height and width (see Figure 1e). The “width” of the barrier can be related to its geometrical width (l) possibly narrowed by image forces22,23 or wave function hybridization41 (R) and to coupling between the molecular orbitals and the electrode (m*). The barrier width scales the current exponentially, with only minor bias dependence (slightly stronger

4442 J. Phys. Chem. C, Vol. 111, No. 11, 2007 dependence if a nonlinear potential barrier is assumed, see Figure 5 and eq 9). This is in contrast to the effect of the barrier height which directly varies with bias (eq 1). F reflects the ratio of these two contributions, in contrast to the decay parameter, β0 which is related to their product. Thus, the “heavier” the charge carrier (m*) is or the wider the barrier (Rl) is, the larger F is and the steeper the I-V curve is away from equilibrium (Figure 2a). Increasing the height of the tunneling barrier will have the opposite effect (smaller F and shallower, more linear I-V). This can serve as a partial answer to the introductory example of how to interpret a set of systematic changes in a molecular junction. If a systematic structural change leads only to an increase in barrier height, it will be expressed as a decrease in both F and G0 (increase in β0). However, if either the net junction width or the effective mass increases, then F is expected to increase while G0 should decrease (β0l increases). Naturally, a combination of these effects will be more difficult to interpret. If, by comparing a series of junctions, we find a constant F (i.e., β0l is also constant) and a changing G0, this can indicate that pinholes dominate conduction and that their density changes. In addition, the linearized presentation of Figure 4 is a straight forward tool for evaluating the bias range where the (linear) Simmons model is applicable and for a possible change in mechanism as indicated by the kink in Figure 4c. The numerical values extracted here and presented in the insert Tables to Figure 3 (and in Tables s2-s4 in the Supporting Information) cannot provide an answer to the question what are the junction parameters. Rather, they stress the ambiguity in the extraction of these parameters. Nevertheless, we can characterize junctions in an acceptable manner as described next. 5.5. Suggested Routine for Parameter Extraction. In view of the discussion in the previous sections, I suggest the following guidelines for analyzing molecular I-V characteristics: (a) Predetermine the contact area by testing the specific contacting method on a known molecular system of varying length (e.g., alkyls), using eq 8a. (b) Measure as many data at low bias (eφ0/3) as possible. Logarithmic bias intervals are convenient for this.18 Naturally, this requires a low-level meter capable of measuring the expected low currents. (c) Use the presentation of Figure 4 and eq 6 to extract the characteristic G0 and C parameters and to validate if the Simmons model is at all relevant for analyzing the data. A complementary linearized presentation (e.g., Fowler-Nordheim, see ref 43) can be added to map the high bias range of the measured data. (d) In the case of deviation from the model (nonlinear I/V vs V2 presentation), only G0 can be regarded as a characteristic parameter and can be translated to β0 ) f(G0, l, A) (eq a7). Only if the I/V vs V2 plots are linear, can one use the following recommendations: (e) Extract approximate F value from C (eq a8). (f) In the case of high confidence in the effective electrical contact area (see point a above), translate G0 into β0l (using C and A, eq a6) and calculate F more accurately (F ) f(C, β0l), eq a8). (g) Combining F and β0l (eqs 3 and a1) yields φ0 ) β0l/F and lRxm* ) 1.025(Fβ0l)1/2, where the last product can be further interpreted depending on the experimental framework (e.g., systematic change in molecular length, conjugation or electrode work function). (h) Achieving a linear plot in any of the versions of eq 8 is required for extracting believable absolute numbers. Still, even

Vilan if linearity is not achieved, such a presentation can be very informative about qualitative differences between different junctions. 6. Summary The possibility for extracting information on transport characteristics from analyses of molecular current-voltage (cf. current-length) data was considered. For this purpose, the accepted Simmons model was reconsidered. Major failures are (1) the strong dependence of extracted parameters on the assumed contact area for transport, and (2) the strong deviation of experimental current-voltage data from the assumptions that underlie the Simmons model, viz., the barrier profile is linear with distance, and the barrier profile varies linearly with bias. The seemingly good graphical correlation between raw data and numerical fits to the Simmons model probably results from the scaled nature of the model, without proving its validity or accuracy. The applicability of the Simmons tunneling model to a specific I-V should be tested using the linearized version (I/V vs V2) rather than using the direct form (I vs V). The relevant range for the Simmons model is near equilibrium, and based on this, an alternative parametrization procedure was suggested. This procedure emphasizes the equilibrium conductance and shape factor as robust characteristics of the junctions, independent of any further information about the junction. These two parameters can be related to more common barrier characteristics, such as the decay parameter and barrier height but only if the contact area is known or from data where the molecular length (barrier width) is varied systematically. Thus, the proposed approach is complementary to the common thickness-dependent analysis, to which it can be added to allow separating the contribution of barrier height from that of an effective mass parameter. The quasiequilibrium I-V analysis will be especially useful for analyses of chemically unique molecular junctions, the thickness of which cannot be varied. Despite the gross oversimplification of the Simmons model, the justification for using it originates from the need for quantifying imperfect molecular junctions. Such quantification is crucial for achieving more controlled, reliable, and reproducible experiments, as well as acquiring a predictive power toward a more realistic description of transport across molecular junctions and more adequate interpretation for extracted averaged characteristic parameters. Acknowledgment. I thank the Kreitman foundation for a fellowship at Ben Gurion University of the Negev; M. Reed and S. Lindsay for kindly providing their I-V data; Y. Selzer and A. Solomon for discussions; and N. Kalisman for helpful remarks on numerical analyses. I am grateful to D. Cahen for numerous discussions, especially on contact area issues, and critical readings of, and constructive remarks on, the manuscript. Appendix A1. Glossary of Terms A ) contact area [µm2] C ) cubic polynomial coefficient [V-2] I ) current [A] G0 ) equilibrium conductance [Ω-1] L ) (true) junction width [Å] P ) fifth order polynomial coefficient [V-4] V ) bias [V] W ) work function [eV] e ) electron charge [C] h ) Plank’s constant [Js]

Analyzing Molecular Current-Voltage Characteristics

J. Phys. Chem. C, Vol. 111, No. 11, 2007 4443

l ) molecular length [Å]

A3. Conversion Equations. The two quasi-equilibrium coefficients can be combined to extract β0l, under a priori known effective mass, m*, and contact area, A, by combining eqs 4c and 7a

m0 ) free electron mass [kg] m* ) electron effective mass

(

R ) width ratio ()L/l)

0

φ0 ) barrier height [V]

0

ln(2πm0e3/h3) - ln

F ) shape factor [V-1] ψ ) potential [V] A2. Other Forms of the Simmons Model. a. Equilibrium (V ≡ 0) expression for the decay parameter, β0

β0 )

4π R 2m*m0eφ0 ) 1.025Rxm*φ0 h x

(a1)

The numerical constant is valid for β0 in [Å-1], φ0 in [V], and dimensionless R and m* values (1.025 is in units of [Å-1 V-0.5]). b. Nonequilibrium (V * 0) expression for the decay parameter, β0

β)

4π R 2m*m0e(φ0 - V/2) ) β0x1 - V/2φ0 (a2) h x

{(

[

)

x ( )] ( [ x ( )]}

4πRl V V 2m*m0e φ0 - φ0 + × h 2 2 4πRl V 2m*m0e φ0 + (|V| e φ0) (a3) exp h 2

)

exp -

The same relations, expressed in terms of the equilibrium decay parameter, β0, are

I)

{(

) ( x ) ( ) ( x )}

V e2 Aφ0 1exp - β0l 2πh (Rl)2 2φ0 1+

1-

V exp - β0l 2φ0

V 2φ0 1+

V 2φ0

(a4)

d. The general Simmons relation including a bias-dependent potential, ψ(V), and profile (width ratio), R(V), indicated by a subscript V

I)(

4π e2 A ‚ψ 2m*m0eψ h exp -lRV hV 2πh R 2l2 V h x V

[

(

(

(ψ h V + |V|) exp -lRV

)

4π 2m*m0e(ψ h V + |V|) h x

)

G0x6C (a6) Am*

β0l - ln(β0l - 2) ) ln(e2/4πh) - ln(G0R2l2/A)

(a7)

Although the logarithmic term on the left cannot be ignored (∼3), it is still possible, though, to calculate β0l without considering that term, to arrive at an initial guess. Rearranging eq 7a provides an explicit term for the shape factor, F, or the barrier height, φ0, as a function of the cubic coefficient, C, and the dimensionless thickness, β0l, extracted from G0 (eq a6)

F)

x x

1 - 2/β0l 6C ≈ 4x6C (a8) 1 - 3/β0l - 3/β02l2

β0l ) 4‚ φ0

φ0 ) V e2 A φ0 - × 2πh (Rl)2 2

(

Note that second term on the left ∼0 and can be neglected to get an explicit expression for β0l. In addition, any uncertainty in the effective mass, m*, has little effect on β0l, as these parameters are logarithmically related. In case of extracting G0 only, it can be translated to β0l, if the junction width (Rl) and real contact area, A, are known by rearranging eq 4b and substituting β0l ≡ Fφ

c. The widely accepted form of the linear Simmons relations is

I)

)

3 1 5 6 β0l - ln 1 + + ) 2 β0l β 2l2 β 3l3

β0 ) decay parameter [Å-1]

1 4

2 β 0l 1 (β0l) - 3/β0l - 3 ≈ 6C 1 - 2/β0l 4x6C

(a9)

where the approximated equality refers to the condition of β0l . 2. A4. Mathematical Tools. Mean square error (MSE) fitting was done by a MATLAB function named “fminsearch” which uses the Nelder-Mead simplex48 direct searching algorithm to find a parameter array, p, that minimizes the MSE defined as49

MSE )

1

N

∑ [f (Vi, p) - Ii]2

N i)1

(a10)

where p ) (G0, F, φ0), [Vi] and [Ii] are N long arrays of voltage and current, and f represents eq 2 (or any variation of it). Implicit equations (e.g., eq a6) were solved using the MATLAB function “fzero” or “fminsearch”. Supporting Information Available: Glossary of translation functions from polynomial coefficients into tunneling parameters (S1); additional maps of mean square error as function of fitting variables (S2); summary tables comparing the tunneling parameters extracted by various fitting procedures (S3); and a detailed experimental section for the bilayer junction (S4). This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes

)]

(a5)

(1) Collier, C. P.; Mattersteig, G.; Wong, E. W.; Luo, Y.; Beverly, K.; Sampaio, J.; Raymo, F. M.; Stoddart, J. F.; Heath, J. R. Science 2000, 289, 1172.

4444 J. Phys. Chem. C, Vol. 111, No. 11, 2007 (2) Chen, J.; Reed, M. A.; Rawlett, A. M.; Tour, J. M. Science 1999, 286, 1550. (3) For a comprehensive review of the field, see: Kagan and Ratner,4 and other articles in this special issue of MRS bulletin. (4) Kagan, C. R.; Ratner, M. A. MRS Bull. 2004, June, 376. (5) Tomfohr, J. K.; Sankey, O. F. Phys. Rev. B 2002, 65, 245105. (6) Magoga, M.; Joachim, C. Phys. ReV. B 1997, 56, 4722. (7) Galperin, M.; Nitzan, A.; Sek, S.; Majda, M. J. Electroanal. Chem. 2003, 550, 337. (8) Zahid, F.; Paulsson, M.; Polizzi, E.; Ghosh, A. W.; Siddiqui, L.; Datta, S. J. Chem. Phys. 2005, 123, 064707. (9) Berlin, Y. A.; Ratner, M. A. Radiat. Phys. Chem. 2005, 74, 124. (10) Slowinski, K.; Chamberlain, R. V.; Miller, C. J.; Majda, M. J. Am. Chem. Soc. 1997, 119, 11910. (11) York, R. L.; Nguyen, P. T.; Slowinski, K. J. Am. Chem. Soc. 2003, 125, 5948. (12) Wold, D. J.; Frisbie, C. D. J. Am. Chem. Soc. 2001, 123, 5549. (13) Definition of the decay length, β varies. Most formally β should scale the width of the potential, L; however as this property is hard to measure experimentally, β is defined clearer as scaling of the molecular length, l. Note that in latter case, β includes an additional contribution from the ratio R ) L/l. It is also common to express β in terms of methylene units or per carbon atom. Here we shall follow the 1/molecular length (Å-1) convention. (14) Wang, W.; Lee, T.; Reed, M. A. Phys. ReV. B 2003, 68, 035416. (15) Wang, W.; Lee, T.; Reed, M. A. Rep. Prog. Phys. 2005, 68, 523544. (16) Parikh, A. N.; Allara, D. L.; Azouz, I. B.; Rondelez, F. J. Phys. Chem. 1994, 98, 7577. (17) Holmlin, R. E.; Ismagilov, R. F.; Haag, R.; Mujica, V.; Ratner, M. A.; Rampi, M. A.; Whitesides, G. M. Angew. Chem.-Int. Ed. 2001, 40, 2316. (18) Cui, X. D.; Primak, A.; Zarate, X.; Tomfohr, J.; Sankey, O.; Moore, A. L.; Moore, T. A.; Gust, D.; Nagahara, L. A.; Lindsay, S. M. J. Phys. Chem. B 2002, 106, 8609. (19) Adams, D. M.; Brus, L.; Chidsey, C. E. D.; Creager, S.; Creutz, C.; Kagan, C. R.; Kamat, P. V.; Lieberman, M.; Lindsay, S.; Marcus, R. A.; Metzger, R. M.; Michel-Beyerle, M. E.; Miller, J. R.; Newton, M. D.; Rolison, D. R.; Sankey, O.; Schanze, K. S.; Yardley, J.; Zhu, X. Y. J. Phys. Chem. B 2003, 107, 6668. (20) Aswal, D. K.; Lenfant, S.; Guerin, D.; Yakhmi, J. V.; Vuillaume, D. Small 2005, 1, 725. (21) York, R. L.; Slowinski, K. J. Electroanal. Chem. 2003, 550, 327. (22) Simmons, J. G. J. Appl. Phys. 1963, 34, 1793. (23) Simmons, J. G. J. Appl. Phys. 1963, 34, 2581. (24) Simmons, J. G. J. Appl. Phys. 1963, 34, 238. (25) Selzer, Y.; Cai, L.; Cabassi, M. A.; Yao, Y.; Tour, J. M.; Mayer, T. S.; Allara, D. L. Nano Lett. 2005, 5, 61. (26) Engelkes, V. B.; Beebe, J. M.; Frisbie, C. D. J. Am. Chem. Soc. 2004, 126, 14287. (27) Vilan, A.; Hikmet, R. A. M. Manuscript in preparation 2007. On preparation of gold flakes, see: patent No. WO 2005085389; Hikmet, R. A. M., Koninklijke Philips Electronics N.V., The Netherlands, 2005

Vilan (28) Brinkman, W. F.; Dynes, R. C.; Rowell, J. M. J. Appl. Phys. 1970, 41, 1915. (29) Miles, J. L.; McMahon, H. O. J. Appl. Phys. 1961, 32, 1176. (30) Zhang, Z.-S.; Rabson, D. A. J. Appl. Phys. 2004, 95, 199. (31) Dorneles, L. S.; Schaefer, D. M.; Carara, M.; Schelp, L. F. Appl. Phys. Lett. 2003, 82, 2832. (32) Mann, B.; Khun, H. J. Appl. Phys. 1971, 42, 4398. (33) Polymeropoulos, E. E.; Sagiv, J. J. Chem. Phys. 1978, 69, 1836. (34) Hansen, K.; Brandbyge, M. J. Appl. Phys. 2004, 95, 3582. (35) Probst, O. M. Am. J. Phys. 2002, 70, 1110. (36) A very similar expression (up to a factor of 4/π) was used by Hansen et al.34 as the exponential coefficient in the WKB transmission function. (37) The shape factor can be viewed in analogy to the Reynolds number in fluid mechanics (the ratio of inertial forces to viscous forces, see: http:// en.wikipedia.org/wiki/Reynolds_number) which determines the flow pattern (i.e., laminar or turbulent) regardless of the absolute flow speed. Similarly, the shape factor dictates the shape of the tunneling I-V characteristics, regardless of the absolute current magnitude. (38) Wold, D. J.; Haag, R.; Rampi, M. A.; Frisbie, C. D. J. Phys. Chem. B 2002, 106, 2813. (39) Fitting with β0l as a free parameter is equivalent to use of m* or R, but β0l is the general lumped parameter, regardless of any assumptions. (40) Joachim, C.; Magoga, M. Chem. Phys. 2002, 281, 347-352. (41) Segev, L.; Salomon, A.; Natan, A.; Cahen, D.; Kronik, L.; Ami, F.; Chan, C. K.; Kahn, A. Phys. ReV. B 2006, 74, 165323. (42) The linear term in the analysis of Brinkman et al.28 expresses asymmetry between the barrier heights at the two electrodes. It emerges from integration over probabilities in contrast to the integration over potential used by Simmons22-24 and followed here. (43) Beebe, J. M.; Kim, B.; Gadzuk, J. W.; Frisbie, C. D.; Kushmerick, J. G. Phys. ReV. Lett. 2006, 97, 026801. (44) Full or high bias range refers to V < φ0, which is the high limit on eq 2. For V . φ0 the tunneling can be described by the Fowler-Nordheim relations.22,23,43 By the term near- or quasiequilibrium, we refer to the range where the I/V vs V2 plot is linear. (45) Nesher, G.; Vilan, A.; Cohen, H.; Cahen, D.; Amy, F.; Chan, C.; Hwang, J.; Kahn, A. J. Phys. Chem. B 2006, 110, 14363. (46) Engelkes, V. B.; Beebe, J. M.; Frisbie, C. D. J. Phys. Chem. B 2005, 109, 16801. (47) Kushmerick, J. G.; Lazorcik, J.; Patterson, C. H.; Shashidhar, R.; Seferos, D. S.; Bazan, G. C. Nano Lett. 2004, 4, 639. (48) Lagarias, J. C.; Reeds, J. A.; Wright, M. H.; Wright, P. E. SIAM J. Optim. 1998, 9, 112. (49) In comparison, the Microcal Origin nonlinear search algorithm relies on derivatives rather than on direct search. Yet the MSE function is similar to the χ2 function of Microcal Origin for a single independent variable (Vi). (50) Similar maps were produced also for the experimental data providing fairly irregular contours, thus not presented. A comparable example can be found in Figure 6 of ref 15.