Comment on “A single level tunneling model for molecular junctions: evaluating the simulation methods”

: The present Comment demonstrates important ﬂaws of the paper Phys. Chem. Chem. Phys. 2022, 24, 11958 by Opodi et al. 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 2 from the exact current I 1 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.

Opodi et al. 1 claimed, e.g., that: (i) The applicability of the method based on I 3 , 2 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 2 and I 3 ) because (ii 1 ) not only method 3 (ii 2 ) 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 6 , 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 ε 0 , coupled to wide, flat band electrodes (hence Lorentzian transmission).The effective level coupling to electrodes Γ g ≡ √ Γ s Γ t 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 ε 0 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. 1 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 (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 B T )] reduces to the step function.
This low temperature limit (expressed by eqn (6a) and (6b) below) assumes a negligible variation of the transmission function (which is controlled by ε 0 and Γ) within energy ranges of widths ∼ k B T 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, Γ ≪ |ε 0 |, see eqn (11)) plays a negligible role and is sufficient for the low temperature limit to apply. 7However, 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 B T ) of the range (∼ (ε ) wherein the electrode Fermi functions rapidly vary.In practice, very close to resonance, loosely speaking, this means [8][9][10] Γ ∼ k B T (6b) (c) Eqn (3) was analytically derived 2 from eqn (2) for sufficiently large arguments of the inverse trigonometric functions The bias range wherein eqn (3) holds within the above accuracy can be expressed as follows where is the zero-temperature low bias conductance per molecule (G/N) in units of the universal conductance quantum G 0 = 2e 2 /h = 77.48µS.Strict on-resonance (ε 0 ≡ 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) yields Γ < 0.1005 and via eqn (8) 3,11 |eV | < 1.4 |ε 0 | or, equivalently (12a) |eV | < 1.25 eV t (cf.eqn ( 13)) (12b) Along with the low-T limit assumed by eqn (2), eqn (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 ε 0 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 2 /|I| plotted vs V . 3,12The fact that eqn (3) should be applied only for biases compatible with eqn (12a) has been steadily emphasized (e.g., ref 13 and discussion related to Fig. 2 3) is particularly useful because it allows expression of the transition voltage V t in terms of the level offset 2 which can thus be easily estimated.ε 0 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 | < ∼ 0.4 eV) even at the room temperature (k B T = 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): The presentation adopted in Fig. 2  With this correction, the present Fig. 1 reveals what it should.Namely that, as long as the off-resonance condition of eqn ( 11) is satisfiedi.e., excepting for Fig. 1A-, in all other panels the blue and green curves (I 3 and I 2 , respectively) practically coincide.Significant differences between the exact red curve (I 1 ) and the approximate green and blue (I 2 and I 3 , 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 1 ) is straightforward.Based on their Fig. 5, Opodi et al. cannot make a statement on method 3: they consider parameters ε 0 < 1 eV at the bias V = 1.5 V(> 1.4 |ε 0 |), 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 2 ), I show in Fig. 2b, c, and d deviations of the current I 2 from the exact value I 1 in snapshots taken horizontally (i.e., constant Γ) and vertically (i.e., constant ε 0 ) 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 2 /I 1 − 1 does not exceed 3%.
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 2 ) and red (exact I 1 ) 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, 4B, S1 to S4, and 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,15) are values of Γ = Γ/2, while those estimated by themselves are values of Γ = 2 Γ.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 3 , 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,15), up to minor inaccuracies inherently arising from digitizing the experimental I-V -curves, they would have reconfirmed the values of ε 0 reported in the original publications, and would have obtained values of Γ = 2 Γ 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 ε 0 , NΓ 2 g , Γ a → |ε 0 | , NΓ 2 , Γ while eqn (3) has only two independent fitting parameters All the narrative on the N-Γ-entanglement and wording on "twin sisters" used in Sec.3.4 of the original article clearly reveal that ref. 1 overlooked that, when using I 3 , N and Γ are two parameters whose values are impossible to separate; they build a unique fitting parameter NΓ 2 ≡ 4N Γ2 .Data fitting using I 3 and three fitting parameters (ε 0 , Γ, N) has an infinity number of solutions for Γ and N but they all yield a unique value of NΓ 2 .
Were method 3 "quite limited" and the deviations of I 3 from I 1 or I 2 significant, Opodi et al. would have been able to determine three best fit parameters (|ε 0 | , Γ, N); at least for the "most problematic" junctions where they claimed important departures of I 3 -based estimates from those based on I 1 and I 2 .If this is indeed the case, the value of N can be determined from data fitting 10 .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 ε 0 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Γ 2 ∝ G 3 ≈ G nor the level offset ε 0 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 5 .Employing such artificially large N's nonphysically reduces Γ (NΓ 2 ≈ constant, Γ ∝ 1/ √ N) 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), 16 the real junction for which they claimed the most severe failure of method 3.If Fig. S7D to S7F and 4D of ref. 1 were correct, both ε 0 and Γ based on I 3 would be in error by a factor of two.
To reject this claim, in Fig. 3 I show curves for I 1 , I 2 , and I 3 computed with the values of ε 0 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, E, and F).
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 represents the exact current I 1 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 2 and I 3 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, E, and F of ref 1).
If Opodi et al. had calculated I 2 and I 3 using the parameters indicated in their Fig.S7D, E, and F (same as in the present Fig. 3a), they would not have obtained the black curves of their Fig.S7D, E, and F but the red, green, and blue curves of Fig. 3a.The values ε 0 = 0.57 eV and Γ = 147 meV of Fig. S7F (so much different from ε 0 = 0.27 eV and Γ = 71 meV of Fig. S7D and ε 0 = 0.28 eV and Γ = 78.6 meV of Fig. S7E) can by no means be substantiated from these calculations.All aforementioned values of ε 0 and Γ of Fig. S7D, E, and 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 1 , I 2 , and I 3 , all computed with the same parameters, namely those of Fig. S7D of ref 1 (ε 0 = 0.27 eV and Γ = 71 meV).In accord with the general theoretical considerations presented under (b) and (c), the differences between I 2 and I 3 (blue and green curves) are negligible, while deviations of I 2 from I 1 are notable only for biases close ±2ε 0 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 |ε 0 | 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) 3 and/or for energy offsets (|ε 0 | < ∼ 0.5 eV) where thermal effects (k B T ≃ 25 meV = 0) matter 11 .It is experimentalists' responsibility not to apply it under conditions that defy the boundaries under which it was theoretically deduced.
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    3) is merely a mathematical formula without any physical sense.Notice that, throughout, the green (I 2 ) and blue (I 3 ) curves excellently agree with the exact red curves (I 1 ) precisely in the parameter ranges predicted by theory, i.e. eqn (6) and eqn (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).
by Opodi et al. also masks the fallacy of applying eqn (3) at biases |eV | > 2 |ε 0 |.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 |ε 0 | are the nonphysical cusps of the blue (I 3 ) curves in Fig. 2 by ref. 1, same as the cusps of the orange curves of the present Fig. 1.

Fig. 1
Fig. 1 Currents I 1 , I 2 , and I 3 computed using the parameters of Fig. 2 of Opodi et al.. Redrawing their figure emphasizes the difference between the current I 3 in the bias range for which eqn (3) was theoretically deduced 2 (blue curves) and I 3 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 2 ) and blue (I 3 ) curves excellently agree with the exact red curves (I 1 ) precisely in the parameter ranges predicted by theory, i.e. eqn (6) and eqn (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).

1 Fig. 2 7 - 8 IFig. 3
Fig. 2 (a) Tick symbols depicting excellent agreement between I 2 and I 1 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,c) Deviations of I 2 from the exact current I 1 reveal that method 2 excellently works in situations where Opodi et al claimed the contrary.Importantly, showing parameter values |ε 0 | < 1 eV at bias V = 1.5 V, panel a is beyond the scope of method 3 (cf.eqn (12a)).