A Single-Level Tunnel Model to Account for Electrical Transport through Single Molecule- and Self-Assembled Monolayer-based Junctions

We present a theoretical analysis aimed at understanding electrical conduction in molecular tunnel junctions. We focus on discussing the validity of coherent versus incoherent theoretical formulations for single-level tunneling to explain experimental results obtained under a wide range of experimental conditions, including measurements in individual molecules connecting the leads of electromigrated single-electron transistors and junctions of self-assembled monolayers (SAM) of molecules sandwiched between two macroscopic contacts. We show that the restriction of transport through a single level in solid state junctions (no solvent) makes coherent and incoherent tunneling formalisms indistinguishable when only one level participates in transport. Similar to Marcus relaxation processes in wet electrochemistry, the thermal broadening of the Fermi distribution describing the electronic occupation energies in the electrodes accounts for the exponential dependence of the tunneling current on temperature. We demonstrate that a single-level tunnel model satisfactorily explains experimental results obtained in three different molecular junctions (both single-molecule and SAM-based) formed by ferrocene-based molecules. Among other things, we use the model to map the electrostatic potential profile in EGaIn-based SAM junctions in which the ferrocene unit is placed at different positions within the molecule, and we find that electrical screening gives rise to a strongly non-linear profile across the junction.


The rate equation formulation
Let M be the number of spin-split orbitals in the molecule. The number of possible electronic configurations is equal to n = 2 M (i.e., each orbital can be empty or occupied). Each configuration can be described by a set of occupation numbers c α (i), with α = 1, , n and i = 1, , M , where c α (i) = 0 or 1, depending on whether the i-th orbital is occupied or not. Call P α the probability of the α configuration, with P α ≥ 0 and α=1 n P α = 1.
Under certain approximations, which neglect memory effects and discard correlations between the electrodes and the molecule, the rate equation governing the change in the configuration probabilities over time can be expressed as [1,2,3] dP α dt = −P α β α where Γ α→ β is the rate of the α → β transition. Equation (2) is the well-known Pauli master equation [4]. We note here that there are several ways to derive this equation from fundamental theories, such as nonequilibrium quantum statistical mechanics [5,6,7], as well as from a semiclassical Boltzmann kinetic equation [1]. Here, rather than repeating those standard derivations, we focus on the use of the rate equation to compute transport proeperties in the context of the experiments reported in the main text. The main goal is to obtain expressions for the charge current across the molecule in different asymptotic regimes.
We can express Eq. (2) in matrix form, where Notice that, consistent with the normalization condition, we have α=1 implying that as it can be easily verified. On the other hand, Thus, if Γ β→α = Γ α→ β for all transitions, then the r.h.s. of this equation is equal to zero. This means that a trivial stationary solution exits where P α = 1/n for all α = 1, , n. However, when Γ β→α Γ α→ β , other nontrivial stationary solution exist as well.

Stationary solutions
In order to obtain a stationary solution to the rate equations set for all α = 1, , n. This implies Thus, to find the set of stationary probabilities {P α }, one needs to find the right eigenvector corresponding to the zero eigenvalue of the matrix Λ.

Transition rates
Let E α be the total energy of the molecule and N α be the total number of electrons in the α configuration.
When an electron hops from one of the leads and into the molecule, energy conservation requires where α(β) is the molecule's configuration before(after) the hopping and ε is the energy of the electronic state in the lead. For this transition to take place, the state with energy ε in the lead must have a finite occupation number, namely f ε − µl k B T > 0, where µ l is the lead's chemical potential (l = R, L) and T is the temperature. Here, f (x) denotes the Fermi-Dirac distribution, In the opposite case, when an electron hops from the molecule and into the lead, we have Now, the occupation number of the state with energy ε must be such that f Call γ R and γ L the level widths due to the coupling to the right and left leads, respectively. We can split the transition rate into two contributions, where [1,2,3] where d H (c β , c α ) is the Hamming distance between the binary sets c β and c α . Namely, only transitions where the number of electrons in the molecule changes by one are allowed. Notice that the bias voltage is equal to V = (µ L − µ R )/e, where e denotes the electron charge. Equation (17) can be derived in a number of ways, with the most standard being Fermi's Golden Rule or time-dependent perturbation theory on the level widths γ R and γ L .
The rate equation approach is valid when k B T ≫ γ R,L (since it assumes a perfectly sharp energy lever in the molecule), and when the tunneling through the molecule is sequentially incoherent. Below, we discuss the validity of the rate equation in more detail. Implicit in Eq. (17) is the assumption that the molecule's energy levels are sharp, such that γ R , γ L are much smaller than other energy scales of the problem, such as k B T , eV , and the separation of energy levels in the molecule.

Current
The current coming from the left lead is equal to where ∆N α→ β = N β − N α .

Total energy
The total energy in the molecule can be broken down as follows (constant charging energy model): where E c is the charging energy, V g is the gate voltage, and {ε i } i=1, ,M are the energies of the orbitals.

Single-level case
Let us apply this formulation compute the stationary current of a molecule with a single orbital (spinless case), in which case M = 1 and n = 2. P 0(1) corresponds to the probability of the empty(filled) state. The stationary problem is defined by the matrix where and with and The eigenvector of the Λ matrix with zero eigenvalue corresponds to and The current coming from the left lead is equal to where

Finite level width
When the broadening of the energy level is not negligible, we have to modify the calculations to account for the uncertainty in ǫ 1 . Let γ be the total level width and D 1 (ε) the density of state profile associated to the single-level configuration; for instance, consider the Lorentzian profile with dεD 1 (ε) = 1. The modified expressions for the transition rates are and where Notice that and Therefore, and Single-level case Going back to the expression defining the current through the left lead, we find This expression generalizes the previous one, Eq. (31), to include a finite level width γ. It is straightforward to check that Eq. (51) recoves Eq. (31) when γ → 0.
It is natural to assume that the total level width can be broken into three components, where γ 0 represents the broadening caused by effects other than the leakage of charge through the leads.

Adding spin
To add spin, we split the configuration where the molecule level is occupied into two (↑, ↓), resulting in a total of three configurations: i = 0, ↑, ↓ (we forbid double occupancy by assuming that the charging energy E c is a very large energy scale, namely, E c ≫ k B T , eV , |ε 1 |). Let us assume that the molecular level is spin degenerate. Then, the total current through the left lead is given by the expression The rate equations are Solving for the steady state yields and Assuming spin degeneracy in the leads, we find and Therefore, and where we introduced Notice that rates and probabilities do not depend on spin. Thus, we can recast the problem in terms of P 0 and P 1 = P ↑ + P ↓ . Then, we find and Plugging them into the expression for the current, we get where Notice that because of the term Γ 0→1 in the denominator of the prefactor in Eq. (72), the expression for the current in the presence of spin is not exactly equal to twice that for the spinless case. However, if we are only interested in linear response, we can set µ L = µ R = µ in Eq. (73), in which case we obtain Γ 0→1 ≈ (γ R + γ L )/ , provided that ε 1 − eV g < µ (namely, when the energy level is brought below the Fermi energy in the leads). Then, the factor of 2 is approximately cancelled and we recover the expression for the spinless current. The current for the spinfull case is only exactly equal to twice that for the spinless case when the charging energy in the molecule is zero (non-interacting limit), in which case conductance through the molecule is spin degenerate.

Exact solution of the single-level case (spinless)
It is possible to solve exactly the fully coherent single-level case by using the Keldysh non-equilibrium technique [7], or even scattering theory, since no many-body interactions are present [8]. The result is the following: the probability of the level to be occupied is equal to where γ = γ R + γ L (absence of any level broadening other than leakage of charge through the leads). The probability of the empty level configuration is P 0 = 1 − P 1 . The expressions for the probabilities are identical to those obtained with the rate equations after the broadening of the energy level is incorporated.
The fact that the coherent and incoherent formulations yield the same results for the probabilities is not surprising. For single channel leads and a single level in the molecule, interference plays no role since there is only one conduction path. When the molecule has multiple independent paths for electrons to hop in and out, then the coherent and incoherent predictions depart, since interference between paths can result in enhancement or depletion of certain configuration occupations.
An expression for the current was derived by Jauho, Wingreen, and Meir [9] using the Keldysh Green's function technique. Their result is Contrary to our previous derivation using rate equations, this expression fully takes into account coherence. Yet, Eq. (75) and Eq. (51) are identical, provided that we set γ = γ R + γ L (namely, no level broadening other than that due leakage through the leads). To some extend this should come as a surprise, as the coherent transport formulation contains incoherent, sequential regime as a limit.
Notice that for the spinfull case, one simply need to insert a factor of 2 on the right-hand-side of Eq. (75).

Asymptotic limits for the current
Notice that where eV b = µ L − µ R and E F = (µ L + µ R )/2. Defining ε ′ = ε − ε 1 + eV g , we can then rewrite Eq. (75) as where It is easy to show that Without loss of generality, we can set eV g = E F . Thus, ε 1 becomes the position of the energy level with respect the Fermi energy in the leads at zero bias. Then, Let us look at some asymptotic limits.
The current in this case shows an activation behavior, with the activation energy being the offset between the energy level in the molecule and the Fermi energy in the leads. Linear bias regime. On resonance, we find • γ ≪ |ǫ 1 | ≪ k B T ≪e |V b |: Weak broadening, intermediate temperature, large bias, nearly on resonance.
The current is approximatelly temperature and bias independent (non-linear bias regime).
The current is again approximatelly temperature and bias independent.
• γ ≪ k B T ≪ e|V b | < |ǫ 1 |: Similar to the previous case, but even more off-resonance.
The current shows activation behavior and is highly non-linear.
Exact solution of the single-level case (spinless) The current decreases with the inverse of the temperature (linear bias regime).
Notice that Eq. (93) is the starting point of a well-known theoretical description of electronic transport in "soft" molecular electronics [10,11]. The current is temperature independent and becomes linear with the bias voltage when e|V b | ≪ γ:

Conclusions
Given that for single-channel, single-level conductance both fully coherent and sequentially incoheren approaches lead to the same expression for the current, we can conclude that the most general expression (at low bias) is given by where we allow the total level width to include some broadening due to energy relaxation mechanisms other than leakage through the leads, namely, γ = γ R + γ L + γ 0 .