Comment on “A single level tunneling model for molecular junctions: evaluating the simulation methods” by E. M. Opodi, X. Song, X. Yu and W. Hu, Phys. Chem. Chem. Phys., 2022, 24, 11958”

Abstract

The present Comment demonstrates important flaws of the paper Opodi et al. Phys. Chem. Chem. Phys., 2022, 24, 11958 Their crown result (“applicability map”) aims at indicating parameter ranges wherein two approximate methods (called method 2 and 3) apply. My calculations reveal that the applicability map is a factual error. Deviations of I₂ from the exact current I₁ do not exceed 3% for model parameters where Opodi et al. claimed that method 2 is inapplicable. As for method 3, the parameter range of the applicability map is beyond its scope, as stated in papers cited by Opodi et al. themselves.

---

Comparing currents I₁, I₂, and I₃ through tunneling molecular junctions computed via three single level models (see below), Opodi et al.¹ claimed, e.g., that:

(i) The applicability of the method based on I₃,² which was previously validated against experiments on benchmark molecular junctions (e.g., ref. 3–5) is “quite limited”.

(ii) The “applicability map” (Fig. 5 of ref. 1) should be used in practice as guidance for the applicability of methods 2 and 3 (i.e., based on I₂ and I₃) because (ii₁) not only method 3 (ii₂) but also method 2 is drastically limited.

(iii) Model parameters for molecular junctions previously extracted from experimental I–V-data need revision.

Before demonstrating that these claims are incorrect, let me briefly summarize the relevant information available prior to ref. 1. Unless otherwise noted (e.g., the difference between Γ and expressed by eqn (5)), I use the same notations as ref. 1.

(a) As a particular case of a formula deduced earlier,⁶eqn (1) expresses the exact current in the coherent tunneling regime through a junction consisting of N molecules (set to N = 1 unless otherwise specified) mediated by a single level whose energy offset relative to the electrodes' unbiased Fermi energy is ε₀, coupled to wide, flat band electrodes (hence Lorentzian transmission). The effective level coupling to electrodes is expressed in terms of energy independent quantities Γ_s,t, representing the level couplings to the two electrodes—substrate (label s) and tip (label t)—, which also contribute to the finite level width Γ_a = (Γ_s + Γ_t)/2.

In the symmetric case assumed following Opodi et al.

and ε₀ is independent of bias.

In this Comment, I use the symbol —a quantity denoted by Γ in ref. 2 and in all studies on junctions fabricated with the conducting probe atomic force microscopy (CP-AFM) platform cited in ref. 1—in order to distinguish it from the quantity denoted by Γ by Opodi et al.¹ Comparison of the present eqn (2) and (3)—in ref. 2 these are eqn (3) and (4), respectively— with eqn (2) and (3) of ref. 1 makes it clear why: is one half of the quantity denoted by Γ = Γ_L + Γ_R = 2Γ_L = 2Γ_R by Opodi et al.¹

(b) Eqn (2) follows as an exact result from eqn (1) in the zero temperature limit (T → 0), when the Fermi distribution f(ε) ≡ 1/[1 + exp(ε/k_BT)] reduces to the step function.

This low temperature limit (expressed by eqn (6) and (7) below) assumes a negligible variation of the transmission function (which is controlled by ε₀ and ) within energy ranges of widths ∼k_BT around the electrodes' Fermi level wherein electron states switch between full (f ≈ 1) and empty (f ≈ 0) occupancies. Away from resonance, (usually a small value, , see eqn (11)) plays a negligible role and

|ε₀ ± eV/2| ≫ k_BT (6a)

is sufficient for the low temperature limit to apply.⁷ However, closer to resonance the aforementioned weakly energy-dependent transmission also implies a sufficiently large . This implies a relationship between the transmission width and the width (∼k_BT) of the range (∼(ε ± eV/2 − k_BT, ε ± eV/2 + k_BT), cf. 1) wherein the electrode Fermi functions rapidly vary. In practice, very close to resonance, loosely speaking, this means^8–10(c) Eqn (3) was analytically derived² from eqn (2) for sufficiently large arguments of the inverse trigonometric functionsThe bias range wherein eqn (3) holds within the above accuracy can be expressed as followswhereis the zero-temperature low bias conductance per molecule (G/N) in units of the universal conductance quantum G₀ = 2e²/h = 77.48 μS. Strict on-resonance (ε₀ ≡ 0) single-channel transport is characterized by g ≡ 1. In the vast majority of molecular junctions fabricated so far, tunneling transport is off-resonant (g ≪ 1), and g < g_max = 0.01 safely holds in all experimental situations of which I am aware, including all CP-AFM junctions considered in ref. 1. Imposing

g < 0.01 (10)

in eqn (9) yieldsand viaeqn (8)^3,11

|eV| < 1.4|ε₀|  or, equivalently (12a)

|eV| < 1.25eV_t  (cf. eqn (13)) (12b)

Along with the low-T limit assumed by eqn (2), (11) and (12a) are necessary conditions for eqn (3) to apply.

Aiming at aiding experimentalists interested in I–V data processing, who do not know ε₀a priori, in ref. 2 I rephrased eqn (7) by saying that eqn (3) holds for biases not much larger than the transition voltage V_t (eqn (12b)). V_t is a quantity that can be directly extracted from experiment without any theoretical assumption from the maximum of V²/|I| plotted vs. V.^3,12 The fact that eqn (3) should be applied only for biases compatible with eqn (13) has been steadily emphasized (e.g., ref. 13 and discussion related to Fig. 2, 3 and eqn (4) of ref. 11).

(d) Eqn (3) is particularly useful because it allows expression of the transition voltage V_t in terms of the level offset²

which can thus be easily estimated. ε₀ is a key quantity in discussing the structure–function relationship in molecular electronics. Thermal corrections to eqn (13), which are significant for small offsets (|ε₀| ≲ 0.4 eV) even at the room temperature (k_BT = 25 meV) assumed by Opodi et al. were also quantitatively analyzed (e.g., Fig. 4 in ref. 11).

To sum up, method 2 applies to situations compatible with eqn (6), and method 3 applies in situations compatible with eqn (6) and (12). This is a conclusion of a general theoretical analysis that needs no additional confirmation from numerical calculations like those of ref. 1.

Switching to the above claims, the following should be said:

To claim (i)

Fig. 2 of ref. 1 shows (along with |I_1,2| also) currents |I₃| computed for biases −1.5 V < V < +1.5 V at couplings and offsets ε₀ = {0.1;0.5;1} eV. The lower cusps visible there depict currents vanishing (I → 0, log|I| → −∞) at V = 0. The upper symmetric cusps (log|I| → +∞ at V → ±2ε₀) depicted by the blue lines were obtained by mathematical application of eqn (3) beyond the physically meaningful bias range of eqn (12a), for which eqn (3) was theoretically deduced.

To make clear this point, I corrected in the present Fig. 1 and 2 of ref. 1. The blue lines in the present Fig. 1 depict I₃ in the bias range e|V| < 1.4|ε₀| for which eqn (3) was theoretically deduced and for which it makes physical sense. The dashed orange lines are merely mathematical curves for I₃ computed using eqn (3) at e|V| > 1.4|ε₀|, where they have no physical meaning. “By definition”, these orange curves are beyond the scope of method 3.

Fig. 1 Currents I₁, I₂, and I₃ computed using the parameters of Fig. 2 of Opodi et al. Redrawing their figure emphasizes the difference between the current I₃ in the bias range for which eqn (3) was theoretically deduced² (blue curves) and I₃ computed outside of bias range (orange dashed curves), wherein eqn (3) is merely a mathematical formula without any physical sense. Notice that, throughout, the green (I₂) and blue (I₃) curves excellently agree with the exact red curves (I₁) precisely in the parameter ranges predicted by theory, i.e.eqn (6) and (12a). The tick symbols at V = 1.5 V depicted in all panels emphasize that method 2 is very accurate, invalidating thereby the “applicability map” shown by Opodi et al. in their Fig. 5 (also reproduced in the present Fig. 2a). Notice that the model parameter values in the panels (a) to (l) depicted here are exactly the model parameter values used by Opodi et al. in their Fig. 2A–L, respectively.

The presentation adopted in Fig. 2 by Opodi et al. also masks the fallacy of applying eqn (3) at biases |eV| > 2|ε₀|. There, the denominator in the RHS becomes negative and bias and current have opposite directions (i.e., I > 0 for V < 0 and I < 0 for V > 0). Visible at |eV| = 2|ε₀| are the nonphysical cusps of the blue (I₃) curves in Fig. 2 by ref. 1, same as the cusps of the orange curves of the present Fig. 1.

With this correction, the present Fig. 1 reveals what it should. Namely that, as long as the off-resonance condition of eqn (11) is satisfied—i.e., excepting for Fig. 1A—, in all other panels the blue and green curves (I₃ and I₂, respectively) practically coincide. Significant differences between the exact red curve (I₁) and the approximate green and blue (I₂ and I₃, respectively) are only visible in situations violating the low temperature condition (eqn (6)): in panels D, E, G, and H violating eqn (6b), and at biases incompatible with eqn (6a).

To claim (ii)

Refuting claim (ii₁) is straightforward. Based on their Fig. 5, Opodi et al. cannot make a statement on method 3: they consider parameters ε₀ < 1 eV at the bias V = 1.5 V(>1.4|ε₀|), which is incompatible with eqn (12a). Noteworthily, the condition expressed by eqn (12a) defied by ref. 1 was clearly stated in references that Opodi et al. have cited.

To reject claim (ii₂), I show in Fig. 2b–d deviations of the current I₂ from the exact value I₁ in snapshots taken horizontally (i.e., constant Γ) and vertically (i.e., constant ε₀) across the “applicability map” (cf.Fig. 2a). As visible in Fig. 2b and c, in all regions where Opodi et al. claimed that method 2 is invalid the contrary is true; the largest relative deviation I₂/I₁ − 1 does not exceed 3%.

Fig. 2 (a) Tick symbols depicting excellent agreement between I₂ and I₁ in all panels of Fig. 1 overimposed on the “applicability map” adapted after Fig. 5 (courtesy Xi Yu) of ref. 1 contradict the claim of Opodi et al. on the inaplicability of method 2. (b) and (c) Deviations of I₂ from the exact current I₁ reveal that method 2 excellently works in situations where Opodi et al. claimed the contrary. Importantly, showing parameter values |ε₀| < 1 eV at bias V = 1.5 V, panel a is beyond the scope of method 3 (cf.Eqn (12a)).

What the physical quantity is underlying the color code depicted in their Fig. 5 (reproduced here in Fig. 2a) is not explained by Opodi et al. Anyway, the conclusion of Opodi et al. summarized in their Fig. 5 contradicts their results shown in their Fig. 2; all panels of that figure reveal excellent agreement between the green (I₂) and red (exact I₁) curves at V = 1.5 V. For the reader's convenience, the thick symbols at V = 1.5 V in the present Fig. 1 overlapped on the “applicability map”of ref. 1 emphasize this aspect. Inspection of these symbols (indicating that method 2 is excellent) overimposed on Fig. 2a reveals that they (also) lie in regions where Opodi et al. claimed that method 2 fails. Once more, their “applicability map” is factually incorrect.

To claim (iii)

In their Fig. 3, 4A, B and Fig. S1–S4, S8 as well as in Table 2 Opodi et al. made unsuitable comparisons: the values for the CP-AFM junctions taken from their ref. 38, 39, 44 and 57 (present ref. 4, 5, 14 and 15) are values of , while those estimated by themselves are values of . Confusing and Γ, no wonder that they needed re-fitting of the original I–V data. If they had correctly re-fitted the CP-AFM data using I₃, with all the values of N given in the original works (namely, their ref. 38, 39, 44 and 57, the present ref. 4, 5, 14 and 15), up to minor inaccuracies inherently arising from digitizing the experimental I–V-curves, they would have reconfirmed the values of ε₀ reported in the original publications, and would have obtained values of two times larger than those originally reported for (cf.eqn (5)).

I still have to emphasize a difference of paramount importance between I–V data fitting based on eqn (1) and (2) on one side, and eqn (3) on the other side. Eqn (1) and (2) have three independent fitting parameters (ε₀,NΓ_g²,Γ_a) → (|ε₀|,NΓ²,Γ) while eqn (3) has only two independent fitting parameters (|ε₀|,NΓ_g²) → (|ε₀|,NΓ²).

All the narrative on the N–Γ-entanglement and wording on “twin sisters” used in Section 3.4 of the original article clearly reveal that ref. 1 overlooked that, when using I₃, N and Γ are two parameters whose values are impossible to separate; they build a unique fitting parameter . Data fitting using I₃ and three fitting parameters (ε₀,Γ,N) has an infinity number of solutions for Γ and N but they all yield a unique value of NΓ².

Were method 3 “quite limited” and the deviations of I₃ from I₁ or I₂ significant, Opodi et al. would have been able to determine three best fit parameters (|ε₀|,Γ,N); at least for the “most problematic” junctions where they claimed important departures of I₃-based estimates from those based on I₁ and I₂. If this is indeed the case, the value of N can be determined from data fitting.¹⁰ Their MATLAB code (additionally relevant details in the ESI†) clearly reveals how they arrived at showing such differences for real junctions considered. In that code, they keep N fixed and adjust ε₀ and Γ. As long it is reasonably realistic, an arbitrarily chosen value of N has no impact on directly measurable properties. It changes the value of Γ but neither NΓ² ∝ G₃ ≈ G nor the level offset ε₀ changes, because method 3 performs well in almost all real cases.

However, defying available values of N for the CP-AFM junctions to which they referred, Opodi et al. spoke of values up to N ∼ 10⁵. Employing such artificially large N's nonphysically reduces down to values incompatible with eqn (6b), arriving thereby at the idea that eqn (3) no longer applies.

In spot checks, I also interrogated curves shown by Opodi et al. for single-molecule mechanically controllable break junctions. I arrived so at the junction of 4,4′-bisnitrotolane (BNT),¹⁶ the real junction for which they claimed the most severe failure of method 3. If Fig. S7D–F and Fig. 4D of ref. 1 were correct, both ε₀ and Γ based on I₃ would be in error by a factor of two.

To reject this claim, in Fig. 3 I show curves for I₁, I₂, and I₃ computed with the values of ε₀ and Γ indicated by Opodi et al. in their Fig. S7D, E and F, respectively. They should coincide with the black curves of Fig. S7D, E and F, respectively if the latter were correct. According to Opodi et al., all these would represent fitting curves of the same experimental curve (red points in Fig. S7D–F, ESI†).

Fig. 3 (a) I–V curves for single-molecule junctions¹⁶ obtained using N = 1 and parameters taken from the figures of ref. 1 indicated in the legends. To convince himself or herself that the present curves for I₂ and I₃ are correct and different from those of Fig. S7E of ref. 1, the reader can easily generate the present green and blue curves by using the GNUPLOT script of the ESI.† (b) The curves for I₁, I₂, and I₃ computed with the parameter values (indicated in the inset) taken from Fig. S7D of ref. 1 do not support the failure of method 3 claimed by Opodi et al.

Provided that MATLAB is available, the reader can run the code “generateIVfitIV.m” included in the ESI,† to convince himself or herself that the red curve of Fig. 3a and not the black curve of Fig. S7D (ESI†) represents the exact current I₁ computed using eqn (1). Otherwise, running the GNUPLOT script also put in the ESI,† will at least convince the reader that the green and the blue lines of this figure and not the black curves in Fig. S7E and F, respectively do represent the currents I₂ and I₃ computed using eqn (2) and (3) for the parameters and the bias range indicated. The reader will realize that the three curves shown in Fig. 3 cannot represent best fits of the same experimental I–V-curve (red points in Fig. S7D–F of ref. 1).

If Opodi et al. had calculated I₂ and I₃ using the parameters indicated in their Fig. S7D–F (same as in the present Fig. 3a), they would not have obtained the black curves of their Fig. S7D–F but the red, green, and blue curves of Fig. 3a. The values ε₀ = 0.57 eV and Γ = 147 meV of Fig. S7F (ESI†) (so much different from ε₀ = 0.27 eV and Γ = 71 meV of Fig. S7D and ε₀ = 0.28 eV and Γ = 78.6 meV of Fig. S7E, ESI†) can by no means be substantiated from these calculations. All aforementioned values of ε₀ and Γ of Fig. S7D–F are exactly the same as the values shown in Table 2 and Fig. 4D of ref. 1, and used as argument against method 3.

As additional support, I also show (Fig. 3b) the curves for I₁, I₂, and I₃, all computed with the same parameters, namely those of Fig. S7D of ref. 1 (ε₀ = 0.27 eV and Γ = 71 meV). In accord with the general theoretical considerations presented under (b) and (c), the differences between I₂ and I₃ (blue and green curves) are negligible, while deviations of I₂ from I₁ are notable only for biases close ±2ε₀ that invalidate eqn (6a).

To sum up, the claim of Opodi et al. on the failure of method 3 for the specific case considered above is incorrect because is based on values incompatible with calculations.

As made clear under (c) above, eqn (12a) is a condition deduced analytically. If it holds, eqn (3) is within ∼1% as good as eqn (2). It makes little sense to check numerically a general condition deduced analytically, or even worse (as done in ref. 1) to claim that it does not apply for biases incompatible with eqn (12a).

The interested scholar needs not the “applicability map” (Fig. 5 of ref. 1, to be corrected elsewhere (I. Bâldea, to be submitted)). In the whole parameter ranges where Opodi et al. claimed the opposite, method 2 turned out to be extremely accurate (cf.Fig. 2b and c). Likewise, “by definition” (cf.Eqn (12a)), method 3 should not be applied at biases above eV > 2|ε₀| shown there, which makes the “applicability map” irrelevant for method 3.

Theory should clearly indicate the parameter ranges where an analytic formula is valid. This is a task accomplished in case of eqn (3). In publications also cited by Opodi et al.,^3,11 particular attention has been drawn on not to apply eqn (3) at biases violating eqn (12a)³ and/or for energy offsets (|ε₀| ≲ 0.5 eV) where thermal effects (k_BT ≃ 25 meV ≠ 0) matter.¹¹ It is experimentalists' responsibility not to apply it under conditions that defy the boundaries under which it was theoretically deduced.

Note added after submission:

(i) The “applicability map” of Opodi et al. has been corrected; see Fig. S15 in I. Bâldea, Phys. Chem. Chem. Phys., 2023, 25, 19750–19763, https://dx.doi.org/10.1039/D3CP00740E.¹⁷

(ii) A formula expressing the exact current I₁ in closed analytic form has been reported, which obviates the numerical integration in eqn (1); see eqn (7) in I. Bâldea, Phys. Chem. Chem. Phys., 2024, https://dx.doi.org/10.1039/D3CP05046G.¹⁸

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

Financial support from the German Research Foundation (DFG Grant No. BA 1799/3-2) in the initial stage of this work and computational support by the state of Baden-Württemberg through bwHPC and the German Research Foundation through Grant No. INST 40/575-1 FUGG (bwUniCluster 2.0, bwForCluster/MLS&WISO 2.0/HELIX, and JUSTUS 2.0 cluster) are gratefully acknowledged.

Notes and references

1. E. M. Opodi, X. Song, X. Yu and W. Hu, Phys. Chem. Chem. Phys., 2022, 24, 11958–11966.
2. I. Bâldea, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 035442.
3. I. Bâldea, Z. Xie and C. D. Frisbie, Nanoscale, 2015, 7, 10465–10471.
4. Z. Xie, I. Bâldea and C. D. Frisbie, J. Am. Chem. Soc., 2019, 141, 3670–3681.
5. Z. Xie, I. Bâldea and C. D. Frisbie, J. Am. Chem. Soc., 2019, 141, 18182–18192.
6. C. Caroli, R. Combescot, P. Nozieres and D. Saint-James, J. Phys. C: Solid State Phys., 1971, 4, 916.
7. G. Sedghi, V. M. Garcia-Suarez, L. J. Esdaile, H. L. Anderson, C. J. Lambert, S. Martin, D. Bethell, S. J. Higgins, M. Elliott, N. Bennett, J. E. Macdonald and R. J. Nichols, Nat. Nanotechnol., 2011, 6, 517–523.
8. I. Bâldea, Phys. Chem. Chem. Phys., 2017, 19, 11759–11770.
9. I. Bâldea, Adv. Theory Simul., 2022, 5, 202200158.
10. I. Bâldea, Int. J. Mol. Sci., 2022, 23, 14985.
11. I. Bâldea, Org. Electron., 2017, 49, 19–23.
12. J. M. Beebe, B. Kim, J. W. Gadzuk, C. D. Frisbie and J. G. Kushmerick, Phys. Rev. Lett., 2006, 97, 026801.
13. I. Bâldea, Chem. Phys., 2012, 400, 65–71.
14. C. E. Smith, Z. Xie, I. Bâldea and C. D. Frisbie, Nanoscale, 2018, 10, 964–975.
15. Q. V. Nguyen, Z. Xie and C. D. Frisbie, J. Phys. Chem. C, 2021, 125, 4292–4298.
16. L. A. Zotti, T. Kirchner, J.-C. Cuevas, F. Pauly, T. Huhn, E. Scheer and A. Erbe, Small, 2010, 6, 1529–1535.
17. I. Bâldea, Phys. Chem. Chem. Phys., 2023, 25, 19750–19763.
18. I. Bâldea, Phys. Chem. Chem. Phys., 2024 10.1039/D3CP05046G.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2cp05110a

---