Molecular electronics in junctions with energy disorder

We investigate transport through molecular wires whose energy levels are affected by environmental fluctuations. We assume that the relevant fluctuations are so slow that they, within a tight-binding description, can be described by disordered, Gaussian distributed onsite energies. For long wires, we find that the corresponding current distribution can be rather broad even for a small energy variance. Moreover, we analyse with a Floquet master equation the interplay of laser excitations and static disorder. Then the disorder leads to spatial asymmetries such that the laser diving can induce a ratchet current.


Introduction
Chemical adsorption of sulfur atoms on gold surface allows a stable bond between gold tips and thiol groups of molecules. This has been exploited for measuring the conductance and the current-voltage characteristics of gold-molecule-gold junctions [1][2][3][4][5][6]. Repeated measurements even with the same sample, however, reveal small but noticeable differences which possibly stem from environmental fluctuations that impact upon the effective molecule parameters. Moreover, the particular form of the gold tip can have a significant influence on the transport properties [7].
A present line of experimental research is the measurement of molecular conductance when the electrons are excited by electromagnetic waves. There one expects various phenomena ranging from photo-assisted transport [8,9] to ratchet or non-adiabatic pump effects, i.e. the induction of dc currents by ac fields even in the absence of any voltage bias [10,11]. Moreover, it has been predicted that properly taylored laser pulses can give rise to short current pulses [12][13][14][15]. Since a dc current flows into one particular direction, a ratchet effect can occur only in "sufficiently asymmetric" systems [9]. In that respect, a static disorder is sufficient to break the reflection symmetry of an individual realization and, thus, may support a ratchet effect.
The quantitative prediction of the current through a molecule is still a great challenge despite the significant progress achieved in recent years [16][17][18][19]. For a more qualitative understanding of the mechanisms involved in molecular transport, it is thus advantageous to employ for the molecule a rather generic tight-binding model [8,9,[20][21][22][23][24][25]. Then a flexible method for the computation of transport properties is provided by master equations of the Bloch-Redfield type which allow one to include electron-electron and electron-phonon interactions, as well as time-dependent fields [9,26]. Similar methods have also been used for describing incoherent transport [27,28].
Here, we explore the role of slow fluctuations or static disorder for molecular conductance. Thereby we will assume that the relevant environmental fluctuations are so slow that they can be described as static disorder which defines an ensemble of wire Hamiltonians. Then a natural quantity of interest is the corresponding distribution of stationary currents. A setup for which this current distribution is also directly relevant is an array of molecular junctions that conduct in parallel. We employ a tight-binding model for the molecule and treat it with the Floquet master equation formalism derived in Ref. [26] which we review briefly in Section 2. In Section 3, we present results for a static model with a large voltage bias, while in Section 4, we investigate pumping effects caused by an interplay of ac driving fields and disorder. The analytical derivation of the current distribution for a wire with two sites is deferred to the Appendix. Figure 1. Bridged molecular wire model consisting of N = 5 sites with internal tunnelling matrix elements ∆ and effective wire-lead coupling strengths Γ L/R .

Wire-lead model and master equation
The system of the driven molecular wire, the leads, and the coupling between the molecule and the leads, as sketched in Fig. 1, is described by the Hamiltonian The wire is modelled by N tight-binding orbitals |n , n = 1, . . . , N, such that with the tunnel matrix element ∆ and the capacitive energy U . Each onsite energy E n (t) + ξ n contains a random contribution ξ n that subsumes the influence of environmental fluctuations. We assume that these fluctuations are Gaussian distributed and so slow that we can treat them as static disorder. Thus, the probability that the onsite energy of orbital n lies in an interval of size dξ around E n (t) + ξ n reads where the variance σ 2 is assumed to be position-independent. This implies that the energy fluctuations are spatially uncorrelated, such that ξ n ξ n ′ = σ 2 δ nn ′ . The onsite energies E n (t) = E n + Ax n cos(Ωt) are modulated by a harmonically time-dependent dipole force, where A denotes the electrical field amplitude multiplied by the electron charge and the distance between neighbouring sites, with x n = 1 2 (N + 1 − 2n) the scaled position of site |n . Our goal will be to compute for many realizations of the wire Hamiltonian the resulting dc current which provides the current distribution P(I). The last term in Eq. (2) captures the electron-electron interaction within a capacitor model and the operator N = ∑ n c † n c n describes the number of excess electrons residing on the molecule. Below we shall assume that U is so large that only states with zero or one excess electron play a role. The first and the last site of the molecule, |1 and |N , couple via the tunnelling Hamiltonian to the respective lead. The operator c † Lq (c † Rq ) creates an electron in the left (right) lead in state |Lq which is orthogonal to all wire states. The influence of the tunnelling Hamiltonian is fully characterised by the spectral density Γ ℓ (ε) = 2π ∑ q |V ℓq | 2 δ (ε − ε q ). If the lead states are dense and located at the centre of the conduction band, the spectral densities can be replaced by a constant, i.e. we assume Γ ℓ (ε) = Γ for both leads.
The leads are modelled as ideal Fermi gases which are initially at thermal equilibrium with the chemical potential µ L/R and, thus, are described by the density operator ρ leads,eq ∝ exp qℓ c qℓ is the electron number operator for lead ℓ = L, R. Since a typical metal screens all electric fields with a frequency below the plasma frequency, we assume that the bulk properties of the leads are not affected by the laser irradiation.

Perturbation theory
The derivation of a master equation starts from the Liouville-von Neumann equation ihρ(t) = [H (t), ρ(t)] for the total density operator ρ(t) for which one obtains by standard techniques the approximate equation of motion [9,[29][30][31][32] Here the first term corresponds to the coherent dynamics of both the wire electrons and the lead electrons, while the second term describes resonant electron tunnelling between the leads and the wire. The tilde denotes operators in the interaction picture with respect to the molecule and the lead Hamiltonian without the molecule-lead coupling, , where U 0 is the propagator without the coupling. The net (incoming minus outgoing) electrical current through the left contact is given by minus the time-derivative of the electron number in the left lead multiplied by the electron charge −e. From Eq. (6) follows for the current in the wide-band limit the expression where f ℓ is the Fermi function of the respective lead andf ℓ = 1 − f ℓ . Furthermore, · · · = tr wire ρ wire · · · denotes the expectation value with respect to the wire density operator. We emphasise that the expectation values in Eq. (7) depend directly on the dynamics of the isolated wire and are thus influenced by the driving.

Floquet theory
An important feature of the current formula (7) is its dependence on solely the wire operators while all lead operators have been eliminated. Therefore it is sufficient to consider the reduced density operator ρ wire = tr leads ρ of the wire electrons. The effort necessary to calculate ρ wire can be reduced significantly by exploiting the fact that the master equation (6) inherited from the total Hamiltonian H (t) a periodic time-dependence, which opens the way for a Floquet treatment.
One possibility for such a treatment is to use the Floquet states of the central system, i.e. the driven wire, as a basis. Thereby we also use the fact that in the wire Hamiltonian (2), the single-particle contribution commutes with the interaction term and, thus, these two Hamiltonians possess a complete set of common eigenstates. In analogy to the quasimomenta in Bloch theory for spatially periodic potentials, the quasienergies ε α come in classes ε α,k = ε α +khΩ, k ∈ Z, of which all members represent the same physical solution of the Schrödinger equation. Thus we can restrict ourselves to states within one Brillouin zone like for example 0 ≤ ε α <hΩ.
For the computation of the current it is convenient to have an explicit expression for the interaction picture representation of the wire operators. It can be obtained from the (fermionic) Floquet creation and annihilation operators defined via the transformation c α (t) = ∑ n ϕ α (t)|n c n , which reads in the interaction picturec with the important feature that the time difference t − t ′ enters only via the exponential prefactor [9].

Master equation and current formula
In the following, we assume the interaction U to be the dominant energy scale in the system, therefore we can restrict the wire Hilbert space to the N + 1 dimensional subspace of states with zero and one electron, such that a basis for the decomposition of the reduced operator is {|0 , c † α (t) |0 }, where |0 denotes the wire state in the absence of excess electrons. Moreover, it can be shown [26] that at large times, the density operator of the wire becomes diagonal in the electron number N . Therefore a proper ansatz reads Note that we keep terms with α = β , which means that we work beyond a rotating-wave approximation. Following our evaluation of the master equation [26], we arrive at a set of N 2 coupled equations of motion for ρ αβ (t) which in Fourier representation read + same terms with the replacement 1, L → N, R.
In an analogous manner we obtain for the dc current the expression The results of this section allow us to numerically compute the dc current through a driven conductor as well as studying the undriven limit. The current distributions discussed below are obtained by computing the dc current for typically 10 4 realizations of the wire Hamiltonian (2). Then these values are taken for a histogram with 150 bins which finally will be scaled such that we obtain a normalised probability density.

Electron transport with slowly fluctuating energies
We first address an undriven wire in the configuration sketched in Fig. 1 where the distribution of all wire levels is centred at energy E n = 0. The transport voltage is so large that all eigenenergies lie well within the voltage window and, consequently, the transport is unidirectional. Then in the absence of onsite energy fluctuations (σ = 0), the current can be evaluated analytically within a rotating-wave approximation and reads I max = eΓ/h(N + 1), i.e. it decays with increasing wire length [33]. The index "max" refers to the fact that any modification of the onsite energies can only reduce the current-which is confirmed by our simulations. The physical reason for this is that for equal onsite energies, solely the kinetic energy determines the eigenstates which, consequently, are delocalised. Different onsite energies, by contrast, tend to "localise" the eigenstates. Thus in the limit of small disorder, the current distribution P(I) is expect to possess a clear peak at I = I max and some minor contribution for lower values of I. Figure 2 shows the simulated current distributions for two different variances. For a small variance (panel a), the distributions for short wires show the expected behaviour. With an increasing wire length, the peak at I = I max disappears and is eventually replaced by an apparently parabolic distribution. This length dependence can be understood in the following way: For a short wire, the probability that a level is out of resonance is rather low and, thus, most realizations of the wire Hamiltonian will allow resonant inter-site tunnelling. With an increasing number of levels, however, the probability for having at least one misaligned level increases and a current significantly smaller than I max becomes more likely. The precise values will depend on the details and, consequently, we expect a broad distribution. This means that whenever a large number of levels plays a role, the transport through a molecule is extremely sensitive to even small disorder induced by environmental fluctuations.
With a larger variance, this scenario becomes even more pronounced as can be seen in Fig. 2b: Then the peak at I max is rather small and noticeable only for 2 and 3 sites. The most likely realization is a completely disordered wire with an accordingly low conductance. For N > 3, the distributions even possesse a significant peak at I = 0 which corresponds to isolating behaviour. A closer inspection of the numerical data reveals that the crossover between conducting and isolating behaviour occurs when the effective disorder √ Nσ exceeds the tunnel matrix element ∆.
Interestingly enough, for N = 2, 3 the distribution turns out to be even non-monotonic, which means that most one most likely finds either a current close to the theoretical maximum or a significantly smaller lower value. The non-monotonic behaviour for N = 2 can also be found analytically. The derivation of the corresponding current distribution (A.3) can be found in the Appendix. The excellent agreement of this analytical solution and the simulated distributions emphasis that the simulation with approximately 10 4 realisations ensures good convergence.

AC-driven disordered junctions
In order to investigate the influence of an AC driving, we employ the same model as above, but now with an additional dipole driving modelled by time-dependent onsite energies E n (t) = Ax n cos(Ωt) as discussed in Sect. 2. The driving frequency Ω = 2∆/h is chosen such that it matches the average splitting of the wire energies, while the amplitude A = ∆ corresponds to intermediately strong driving. The solid line in Fig. 3 shows the current distribution in the absence of a voltage bias, V = 0. The reflection symmetry of the ensemble relates to the symmetric shape of the distribution, which implies that the current vanishes in the ensemble average. An individual realization of the wire, however, generally does not possess reflection symmetry because the random energy shifts are spatially uncorrelated. This asymmetry in combination with the non-adiabatic driving induces a coherent ratchet current, i.e. a dc current even in the absence of any net voltage bias. In the present case, the ratchet current is of the order of 10-20 percent of the current observed above in the large bias limit. This order of magnitude is typical when the driving frequency or a multiple of the driving frequency lies close to an internal resonance, while the intensity is moderate [26]. In addition to the broad distribution of ratchet currents, P(I) exhibits a peak at I = 0. This stems from realizations for which the driving is well out of resonance. For a bias voltage V > 0, the ensemble no longer possesses reflection symmetry and the current distribution is shifted towards positive values (see Fig. 3). For sufficiently small voltages V ∆/e, non-adiabatic pumping against the voltage bias is still possible.
Rather surprisingly, the peak at zero current remains. It now corresponds to realizations for which on the one hand, the driving is off-resonant while on the other hand, both levels lie outside the voltage window. With an increasing bias voltage, the second condition is less frequently fulfilled, and eventually the distribution assumes a form similar to that found for a large voltage in the absence of driving. Already for V ≈ 4∆/e, the distribution is hardly distinguishable from the one shown in Fig. 2a for a wire with N = 2 sites in the absence of driving.

Conclusions
We have investigated the current through a molecular wire with disordered onsite energies. Such a disorder can stem from the interaction with slow fluctuations of background charges in the substrate. In particular, we computed the resulting current distribution for two typical cases, namely an "open transport channel" and a driven molecular wire for which random energy shifts break reflection symmetry and, thus, the driving can induce a ratchet current. The open transport channel is characterised by tight-binding levels with equal onsite energies, such that any misalignment stems from the disorder. Its main consequence is that as soon as the standard deviation of the onsite energies exceeds the tunnel matrix elements, the current distribution no longer peaks only at a finite value, but also at zero. For longer wires, only the peak at zero current remains. This isolating behaviour resembles Anderson localisation which is found for electrons in a one-dimensional disordered lattice [34]. Note however, that we here considered short wires far from the scaling limit in which Anderson localisation is usually studied.
Since the random energy shifts break reflection symmetry, driving the molecular wire with a laser field induces ratchet currents for which we found a relatively broad distribution. If the driving frequency is far from any molecular excitation energy, the ratchet current will be rather small, and we indeed found that this is the case for very many wire realizations. It has the consequence that the corresponding distribution possesses a spike at zero current. This means that non-adiabatic pumping of electrons against a voltage bias becomes generally impossible whenever the relevant wire energy levels lie well within the voltage window.
In conclusion, our results reveal that slow fluctuations or a static disorder can have a significant effect on molecular conduction. In various cases, the current distribution emerges to be rather flat, which means that even the magnitude of the current depends sensitively on environmental influences.

Appendix A. Analytical solution for two levels
A wire model with N = 2 sites represents an analytically solvable case for which one observes a non-monotonic current distribution and which can serve as test case for numerical implementations. Here we consider a two-level system with on-site energies ξ 1,2 , i.e. with a bias 2η = ξ 1 − ξ 2 . Since the random energy shifts ξ n are Gaussian distributed with variance σ 2 , the bias 2η is also Gaussian distributed but with variance 2σ 2 , i.e. its distribution function For the computation of the current, we restrict ourselves to the limit of a large transport voltage such that both eigenenergies of the two-level system lie within the voltage window. Then, the Fermi functions of the left and the right lead effectively become f L = 1 and f R = 0. In this case, transport can be described within rotating-wave approximation (RWA) which practically means that the reduced density operator of the wire is diagonal in energy representation [33]. Within RWA thus follows from the master equation (9) the occupation probability ρ αα = w 1 α /w 2 α and, thus, ρ 00 = 1 − ∑ α w 1 α /w 2 α . The coefficients w n α = | φ α |n | 2 denote the overlap between the eigenstate |φ α and the localised state |n (Note that in the undriven case, all non-vanishing contributions have sideband index k = 0, such that here the sideband index k can be omitted). Inserting this solution into the current formula (10), we obtain I = eΓ/h(1 + ∑ α w 1 α /w 2 α ). The remaining task is now to diagonalise the single-particle Hamiltonian which provides the coefficients w n α . For bias 2η and tunnelling matrix element ∆, the Hamiltonian in pseudospin notation reads H = ησ z + ∆σ x and possesses the eigenenergies ±δ = ±(η 2 + ∆ 2 ) 1/2 . The corresponding eigenvectors φ α are proportional to (δ +η, ∆) and (δ −η, ∆), respectively, such that w 1 α /w 2 α = (δ ± η) 2 /∆ 2 . Then we obtain for the current the expression I(η) = eΓ h 1 3 + 4η 2 /∆ 2 = I max 1 + 4η 2 /3∆ 2 , (A.1) which assumes its maximum I max = eΓ/3h in the unbiased limit η = 0. The probability distribution for the current relates to w(η) via where the summation considers all values of η that fulfil the condition I = I(η). After some straightforward algebra, we obtain by evaluating expression (A.2) the current distribution P(I) =