Many-body localized molecular orbital approach to molecular transport

Many-body localized molecular orbital approach to molecular transport

Dmitry A. Ryndyk111New affiliation: Institute for Materials Science, Dresden University of Technology, D-01069 Dresden, Germany, Andrea Donarini, Milena Grifoni, and Klaus Richter Institute for Theoretical Physics, University of Regensburg, D-93040 Regensburg, Germany

An ab initio based theoretical approach to describe nonequilibrium many-body effects in molecular transport is developed. We introduce a basis of localized molecular orbitals and formulate the many-body model in this basis. In particular, the Hubbard-Anderson Hamiltonian is derived for single-molecule junctions with intermediate coupling to the leads. As an example we consider a benzenedithiol junction with gold electrodes. An effective few-level model is obtained, from which spectral and transport properties are computed and analyzed. Electron-electron interaction crucially affects transport and induces multiscale Coulomb blockade at low biases. At large bias, transport through asymmetrically coupled molecular edge states results in the occurrence of “anomalous” conductance features, i.e., of peaks with unexpectedly large/small height or even not located at the expected resonance energies.

73.63.-b, 85.65.+h

I Introduction

Experimental and theoretical investigations of electron transport through single molecules are in the focus of the rapidly growing field of molecular electronics Cuniberti et al. (2005); Cuevas and Scheer (2010); Song et al. (2011). Electron-electron interaction plays an important role, controlling the position of resonant levels and leading to phenomena such as Coulomb blockade and Kondo effects. Depending on the ratio between the energy scales associated with an effective charging energy and coupling to the leads, molecular junctions can be classified in three groups. In the case of very small coupling, the molecular orbitals are only weakly disturbed, strong charge quantization and Coulomb blockade take place and transport is mainly determined by sequential tunneling Park et al. (2002); Kubatkin et al. (2003); Osorio et al. (2007). In the opposite case of large coupling, the electronic molecular orbitals are hybridized with states in the leads, charge quantization is suppressed, transport is mainly coherent and the conductance is of the order of the conductance quantum  Smit et al. (2002); Tal et al. (2008, 2009). Finally, in the intermediate situation, the effective electronic spectrum of a molecule is determined essentially by the hybridization, and the interplay between charge quantization and coherent transport may be important Reed et al. (1997); Lörtscher et al. (2007); Kim et al. (2011). Let us consider as a commonly used benchmark example a gold-benzenedithiol-gold (Au-BDT-Au) molecular junction with the central molecule --. The experiments Reed et al. (1997); Lörtscher et al. (2007) show that, while it is difficult or impossible to observe the typical Coulomb blockade features in this type of molecular junctions, there are signatures of correlated transport, namely a conductance gap at small voltages and a complex shape of the conductance peaks. One can conclude that transport in the case of intermediate coupling can be correlated and partially incoherent.

In parallel with the experimental investigations, a number of theoretical methods were suggested to describe the structure and electronic properties of molecular junctions. In the case of coherent transport ab initio DFT+NGF methods Taylor et al. (2001); Damle et al. (2001); Xue and Ratner (2003); Th. Frauenheim et al. (2002); Brandbyge et al. (2002); Pecchia and Di Carlo (2004); Ke et al. (2004), combining density functional theory (DFT) and nonequilibrium Green function (NGF) techniques Kadanoff and Baym (1962); Keldysh (1964); Langreth (1976); Rammer and Smith (1986); Haug and Jauho (1996), have become standard Cuniberti et al. (2005); Cuevas and Scheer (2010). However, the use of DFT, which is a powerful tool to deal with ground state properties, may become questionable when applied to transport and nonequilibrium situations, especially when inelastic and interaction effects are essential. Indeed, the use of the solutions of the Kohn-Sham equation as physical quasiparticle states can not be rigorously established. Besides, the DFT+NGF is a mean-field theory and lacks to describe many-body effects in systems with strong electron interactions. Thus, being based on the basis of an atomistic modeling of molecular junctions, the rigorous extension of the DFT approach to describe transport phenomena remains a challenge.

Figure 1: (Color online) Schematic picture of a Au-BDT-Au molecular junction. The dashed box defines the extended molecule. It comprises the central molecular region and parts of the leads, also marked with boxes.

On the other hand, model-based approaches are particularly important to understand the physics of correlated transport. The models are solved usually within two main approaches: the NGF technique in the case of strong coupling to the leads, and the quantum master equation (QME) Blum (1996); Breuer and Petruccione (2002) in the weak-coupling limit. The quantum master equation is usually formulated in the basis of many-body eigenstates of the molecule. It gives a fairly complete description of sequential tunneling, the main features of Coulomb blockade and even can capture Kondo physics for temperatures of the order of or larger than the Kondo temperature König et al. (1996). Finally, for some nonperturbative effects covering the whole range of temperatures, e.g. the Kondo effect in the crossover regime, more sophisticated treatments are necessary, e.g. numerical renormalization group approaches Costi et al. (1994); Andergassen et al. (2010), other numerical methods Segal et al. (2010); Pletyukhov and Schoeller (2012) or Keldysh field theories Smirnov and Grifoni (2012).

The challenge for molecular transport theory is to combine an ab initio approach, required to take into account a realistic geometry and capable to provide the electronic structure of molecular junctions, with a many-body quantum transport technique, which is necessary to incorporate the correlation effects. In this direction a number of different ab initio based approaches were suggested Hettler et al. (2003); Delaney and Greer (2004); Ferretti et al. (2005); Mirjani and Thijssen (2011); Bergfield et al. (2011); Frederiksen et al. (2007); Thygesen and Rubio (2007, 2008); Korytár et al. (2010); Korytár and Lorente (2011); Greuling et al. (2011); Karolak et al. (2011); Tröster et al. (2012); Strange et al. (2011); Rangel et al. (2011). It should be noted that many-body calculations usually require sophisticated analytical and numerical methods and can be very time consuming. Hence the methods of transport theory can not usually be applied directly to large realistic systems; instead a combined approach is preferable, where an effective model takes into account only the states predominantly participating in transport. Progress in this direction was achieved recently in applications to Coulomb blockade phenomena Delaney and Greer (2004), correlated transport through atomic wires Ferretti et al. (2005); Mirjani and Thijssen (2011), many-body interference Bergfield et al. (2011) and Kondo effect Korytár et al. (2010); Korytár and Lorente (2011); Greuling et al. (2011); Karolak et al. (2011); Tröster et al. (2012). Besides, Coulomb blockade phenomena in benzene and benzenedithiol junctions were also considered on the basis of semiempirical atomistic models Begemann et al. (2008); Bergfield and Stafford (2009); Sobczyk et al. (2012). However, a systematic ab initio based many-body theory of molecular transport is still lacking.

Our aim is to further develop such a theory. The main problem to be solved on the way from an atomistic model to transport calculations is that a huge number of microscopic (single-atom) degrees of freedom must be reduced to an effective few-electron (or few many-body state) interacting model, as a prerequisite of successive many-body transport methods. The main building blocks of our approach are an effective electron model (an orthonormal basis of localized molecular orbitals and a many-body Hamiltonian in this basis) extracted from atomistic calculations, and a nonequilibrium quantum transport method, based on nonequilibrium Green functions or on the quantum master equation.

In the case of strong or intermediate coupling to the leads the electronic molecular levels should be considered together with the levels in the leads. A corresponding generic atomistic model is shown in Fig. 1. A so-called extended molecule (inside the dashed box) is placed between equilibrium electrodes. The extended molecule consists of the central region (roughly the molecule) and two leads (the regions of electrodes near the molecule). The size and boundaries of the central region are actually not fixed and should be determined in a way to include all relevant electronic states, as we will see below.

After having defined the appropriate size of the extended molecule, we proceed as follows:

i) We perform Hartree-Fock or DFT ab initio calculations within the extended molecule and obtain the molecular orbitals, which are orthogonal and can serve as a basis for a many-body Hamiltonian.

ii) The central region (e.g. an organic molecule) and the metallic leads have quite different physical properties, and the approximations required to describe interactions are also different. Thus, it is convenient to transform the basis of molecular orbitals obtained in (i) into a basis of localized molecular orbitals (LMOs), which can be spatially separated into a basis for the central region and a basis in the lead ends. Besides, the advantage of LMOs is that the Coulomb interaction in this basis can be simplified to the Hubbard form. This procedure enables us to obtain a many-body Hubbard-Anderson Hamiltonian in the central region with ab initio calculated parameters, including the coupling to the leads. The rest of the leads can be treated within some mean-field approximation.

iii) Using the Hubbard-Anderson Hamiltonian, many-body methods (nonequilibrium Green functions and the master equation in this paper) can be effectively applied to analyze spectral and transport properties of molecular junctions.

The paper is organized as follows. The ab initio basis of localized molecular orbitals is introduced in Sec. II. The electron-electron interaction and Hubbard approximation are discussed in Sec. III. The parameters of the effective Hubbard-Anderson Hamiltonian for the subspace of electronic states relevant for transport through the central region are derived in Sec. IV. In Sec. V we use nonequilibrium Green functions to analyze the coupling to the electrodes and spectral properties. The many-body spectrum of the central region is discussed in Sec. VI and the transport results are shown in Sec. VII. We give a short conclusion in Sec. VIII.

Ii Localized molecular orbitals

The first stage of our approach includes ab initio calculations of the optimized geometry and the relevant basis of electronic states. For the preliminary geometry optimization and molecular dynamics of whole structures we apply the DFT code Siesta Sie (). For the extended molecule only (without the full electrodes) the full-electron quantum chemistry code Firefly (former PC GAMESS) Granovsky () is used. The final geometry optimization is performed using a hybrid DFT method, usually B3LYP. Then, the calculation of the molecular orbitals (MOs) is performed inside the extended molecule (including both the central region and leads). Our test calculations show that a simple local density approximation (LDA) gives reasonable results for the considered systems. At this stage we obtain the MOs with associated energies . MOs have the advantage of being normalized and orthogonal. However, the canonical MOs of the extended molecule can not serve as a good basis for systems with interactions, because both electron-electron and electron-vibron interactions require physically different approximations in different, metal or molecular, parts of the junction. It is hence more convenient to use localized molecular orbitals (LMOs). The interactions between localized orbitals have a simple and transparent form and the appropriate approximations can be found. Besides, the number of required, relevant basis states is smaller and better controlled for LMOs than for MOs.

Figure 2: (Color online) The energy spectrum of the molecular orbitals of the extended molecule shown in Fig. 1. The occupied levels are marked red, the empty levels yellow.

To proceed, we divide all MOs into four groups (Fig. 2). Most relevant for transport are the transport levels near the Fermi energy of the electrodes. These levels are selected for the localization procedure and include both occupied and valence MOs in an appropriate energy range around the Fermi energy. This energy range should be larger than the energy scales of the external bias voltage and temperature. The other criterion is that the obtained LMOs should be localized strongly enough inside the central region and in the leads: only in this case it is possible to separate the system into parts and use different approximations for the central region and leads separately. Alternatively, the partial localization of only relevant MOs (for example only -type orbitals) is possible. In any case, these transport electron states play the main role. All other polarization states, which are further away from the Fermi energy or do not participate in transport for some other reasons, can still affect the interaction between transport electrons and, in particular, screen the Coulomb interaction. The polarization MOs can be localized in the same way as the transport orbitals. Finally, the core orbitals at low energies can be included into the effective (pseudo)potential and empty orbitals at high energies are neglected.

Figure 3: (Color online) Localized molecular orbitals in the central region of an Au-BDT-Au molecular junction.
Figure 4: (Color online) Localized molecular orbitals in the leads (several out of many are shown by different colors).

For benzenedithiol, considered in this paper, the transport window was chosen to be about 4 eV, and contains about 60 MOs, mainly due to the large density of states in the metal leads. The parameters of the obtained effective model are rather robust against the exact choice of this number.

The LMOs are obtained from the MOs by the unitary transformation


The indices with bars denote the states without the spin degree of freedom. We apply the Foster-Boys localization method Foster and Boys (1960), which minimizes the spatial extent of the orbitals and maximizes the distance between orbital centers. Thus, we obtain maximally localized orbitals. Out of the about 60 orbitals only 5 are localized in the central region. While the initial MOs spread across the extended molecule, the LMOs are spatially localized in the central region (Fig. 3) or in the leads (Fig. 4).

Due to the unitary transformation (1) the LMOs are still orthogonal and normalized, but the expression


is no longer diagonal. Moreover, the related Hamiltonian takes the form


where , , and are the Hamiltonians of the left lead, the central region, and the right lead separately. The direct coupling between left and right leads can be neglected in most cases, because the LMOs of different leads are only very weakly overlapping.

Iii Coulomb interaction and Hubbard approximation

Having the LMOs at hand we can calculate the Coulomb matrix elements. The electron-electron interactions are described by the Hamiltonian


where and denotes the spin. The matrix elements for the scalar Coulomb interaction are defined as


where is the Kronecker symbol. For the systems with localized wave functions , where the overlap between two different states is weak, the main matrix elements are those with and . We checked that, indeed, the overlap of 3 or 4 different orbitals can be neglected in comparison with the overlap of only 2 orbitals. In these cases it suffices to replace by the Hubbard Hamiltonian


describing only density-density interactions with and the Hubbard parameters defined as


As a further advantage of LMOs, the local nature of electron correlation is better described in the localized basis.

The interaction is the bare Coulomb interaction , if we include all electrons into the localization procedure. Actually we restrict the effective Hamiltonian to transport electronic states, on which we performed the Hubbard approximation. The remaining polarization electrons are included only through a screened Coulomb interaction. The screening can be described by the effective interaction potential in RPA (or GW) approximation, which is, however, energy dependent. To keep the simplicity of the Hubbard approximation, we use a screened Coulomb interaction with . This approximation gives reasonable values of for the benchmark -orbital model of benzene CH which we compared with the optimized semi-empirical PPP model Pariser and Parr (1953). In this way all coefficients are derived, but further semi-empirical corrections could be included for better agreement with experiments.

Iv Many-body model

The next important step is the derivation of the entire effective Hamiltonian in the basis of the LMOs for the central region. As we already explained, we divide the extended molecule into the central part (actually the molecule in our particular case) and the leads (Fig. 5). The full Hamiltonian is the sum of the noninteracting molecular Hamiltonian , the electron-electron interaction Hamiltonian , the Hamiltonians of the ends of the leads , and the tunneling Hamiltonian describing the molecule-to-lead coupling:


In our case the central Hamiltonian is replaced by the Hubbard cluster Hamiltonian:


where are the bare energies of the LMOs, including the shifts due to an external voltage.

Figure 5: The model of the extended molecule.

The noninteracting molecular Hamiltonian is obtained from the LMO Hamiltonian , Eq. (3). The zero-voltage energies are calculated from the HF or DFT MO energies , from which the contribution of the Hartree terms due to Hubbard interactions should be substracted:




In this expression denote the populations of the corresponding LMOs in the ab initio calculation. The sum is taken over all occupied molecular orbitals.

The coupling to the leads is described by the tunneling Hamiltonian


and the Hamiltonians of the left and right leads are


where is the index of a state, and is the spin. Note that in our case the states in the leads are not plane waves, but are represented by LMOs, calculated simultaneously with the LMOs in the active region. The leads are considered at the mean-field (DFT) level. The equilibrium electrodes, which can have different electrochemical potentials, determine the boundary conditions for the leads.

For the 5-level model (Fig. 3), which represents actually 10 levels with spin, we obtain the following parameters (for spin up or down, all numbers are in ):


The matrix of the Hubbard parameters calculated from expression (8) reads


The diagonal entries in this matrix correspond to the interaction of two electrons in the same LMO state but with different spin. The off-diagonal terms denote coupling beween two different LMOs spin of the electronic state. Of course, these expressions should be transformed into the 10-level basis before performing further calculations.

At finite bias voltage (defined by the left and right electrical potentials, ) the energies are shifted. In linear approximation these shifts are described by : , , where the parameters characterize the symmetry of the voltage drop across the junction, stands for the symmetric case. Note that the energies and other parameters can also be ab initio calculated at finite voltage, but that is very time-consuming.

The coupling matrix elements in Eq. (13) are obtained directly from the localization procedure. Indeed, the Hamiltonian of an extended molecule takes the form Eq. (3), and the off-diagonal terms describe the coupling to the leads. The number of states in the leads is many times larger than in the central region. Thus, to leading approximation, we can average over the lead level distributions and couplings (so-called wide-band limit). In this approximation the level-width function


is energy independent.

We now return to our 5-level model. The couplings of the first two states (localized at the left sulfur atom) and the last two states (localized at the right sulfur atom) are characterized by the level-width functions , , , . All other couplings are small and do not play an essential role. Thus, all parameters of the model Hamiltonian Eqs. (10)-(14) are well defined and we can proceed with transport calculations. In view of the experiments [Reed et al., 1997; Lörtscher et al., 2007], we will perform calculations below room temperature, , implying .

As we discussed in the introduction, transport at finite voltage can be described in the framework of nonequilibrium Green function or quantum master equation approaches implying numerical methods. For benzene-based junctions several methods were used, including coherent DFT based, the master equation approach in the sequential tunneling limit Hettler et al. (2003); Begemann et al. (2008), sophisticated approximations in the framework of the nonequilibrium Green function method Thygesen and Rubio (2008); Bergfield and Stafford (2009); Strange et al. (2011); Rangel et al. (2011), and other methods Delaney and Greer (2004). In this paper we use both NGF and QME methods, trying to attack the problem from both sides. It should be noted, however, that for our case, , both NGF and a QME with second order rates can only give a qualitative description of the transport problem. Very recently, a QME approach for a single-level junction able to describe the regime has been proposed Kern and Grifoni (2012). An extension of this method to a multilevel system will be subject of future investigations.

V Nonequilibrium Green function approach to spectral properties

We follow the formulation pioneered by Meir, Wingreen and Jauho  Meir and Wingreen (1992); Jauho et al. (1994); Jauho (2006). The lesser (retarded, advanced) Green function matrix of a nonequilibrium molecule can be found from the Dyson-Keldysh equations in the integral form


or from the corresponding equations in the differential form




is the total self-energy of the molecule composed of the interaction self-energy and the tunneling (coupling to the left (L) and right (R) lead) self-energies


where is the Green function of the leads.

The retarded tunneling self-energy can be represented as


where is the real part of the self-energy, which usually can be included in the single-particle Hamiltonian , and describes level broadening due to coupling to the leads. In the case of noninteracting leads with continuous energy spectra, the level-width function is determined by the expression (17).

For the corresponding lesser function of the noninteracting leads one finds




is the equilibrium Fermi distribution function with chemical potential ().

The expression for the interaction self-energy can not be obtained exactly. In the nonequilibrium Hartree-Fock approximation one has Ryndyk et al. (2009)


We do not consider more sophisticated cases here.

The current from the left () or right () lead into the central system is described by the expression (for spin-unpolarized leads)


Applying the NGF technique to our case, we should take into account that our system initially consists not of three, but of five parts: the large electrodes, quantum leads and the central region (Fig. 1). Accordingly, the full Hamiltonian has the form


where the central part is the same as before, Eq. (3), and the additional terms describe the electrodes and the coupling between the electrodes and the leads.

Figure 6: (Color online) The spectral function within the different approximations: initial DFT (black dashed), restricted HF (red), unrestricted HF (green) and equation-of-motion (dashed blue).

The solution of Eq. (18) for the central part is in this case


with the lead self-energies (see Eq. (23)) and level-width functions defined as


The lead Green functions for noninteracting leads (or for leads described by the effective mean-field Hamiltonians ) are defined correspondingly by the electrode self-energies


The calculation of the electrode self-energies is done by standard methods Cuniberti et al. (2005); Cuevas and Scheer (2010). The lesser Green functions are calculated in the same way. We assume additionally that the distribution function in the leads is the same as in the corresponding electrodes.

The equations are solved self-consistently within four approximations: initial DFT with the mean-field energies Eq. (2), the restricted HF approximation (RHF) with , the unrestricted HF approximation (23) and, finally, the equation of motion method (EOM) Song et al. (2007); Ryndyk et al. (2009).

First we analyze the equilibrium (zero bias) spectral function of the central region,


The first thing one can see (Fig. 6) is that the RHF approximation gives a spectral function similar to the ab initio (DFT) one. This is not surprising as the ab initio calculation is also RHF and we simply extracted the Hartree contribution when calculating the energies in Eq. (12). The other important point is that the results obtained in HF and EOM approximations are distinctly different and a gap is opened at the Fermi surface. The analysis of the populations of the single-particle states shows that two empty states are located at the left and one at the right side of the molecule (second and fifth states in Fig. 3). These two empty states have the same spin (in the HF approximation the ground state is spin polarized and degenerate, but quantum fluctuations can switch between different spin orientations), indicating that the true ground state, with quantum fluctuations taken into account, can be spin singlet or triplet. We discuss this point in the next section.

Vi The many-body spectrum:
exact diagonalization and the ground state properties

To gain insight into the many-body energy spectrum of the central system we perform an exact diagonalization of Eq. (10) obtaining the set of many-body eigenstates . Calculating the tunneling matrix elements we obtain from Eqs. (10) and (13) the Hamiltonian


First, we analyze the many-body spectrum. With 5 LMOs we get 1024 many-body states in the Fock space. The lowest 8-particle states consist of a series of alternating singlets and triplets (see Table 1). In particular, the ground state is a singlet, practically degenerate with a triplet (). It follows, at a distance of roughly a second pair of singlet-triplet quasi degenerate states. The 9-particle states are all doublets with a relatively regular distance of the order of . Finally, it is important to keep in mind that the energies of the lowest four 8-particle levels lie all below the one of the 9-particle ground state. The 7-particle states have much higher energies.

Level Energy [eV] Spin []
-91.1849 0
-91.1848 1
-90.8653 0
-90.8648 1
-90.7866 1/2
-90.3693 1/2
-90.0891 1/2
Table 1: The eigenenergies and the associated spins of the lowest 8 and 9 particle levels of an Au-BDT-Au molecular junction.

Note that in the sequential-tunneling master equation method the exact many-body states can be partially occupied at finite temperatures, but not at zero temperature, and the level broadening is not taken into account. This can give some noticeable difference in the position of the transport resonances compared to the NGF calculations, where the levels can be partially occupied even at low temperatures, and where the real part of the HF self-energy, Eq. (23), describes the energy shift of the single-particle levels.

The population probabilities are found from the master equation


where the tunneling rates are, in second order in the tunneling Hamiltonian,




This expression connects the tunneling rates to the level-width function; thus the calculated by the NGF method can be used, see Eqs. (33,34). In the wide-band limit one has , where is the density of states.

To check that the simple (diagonal) form of the master equation can be used, we have analyzed the many-body spectrum of the considered system and came to the conclusion that no coherences are needed for the description of the transport since the degeneracies are not of orbital but of spin nature (e.g. triplets for the 8-particle and doublets for the 9-particle states). However, there cannot be mixing of states with different total spin since otherwise the mixing will depend on the choice of the direction of the quantization axis. The solution of the Eqs. (40-42) is straightforward and can be obtained by direct numerical integrations in stationary and time-dependent cases.

As we discussed before, in the equilibrium state at zero voltage there are 8 electrons distributed due to thermal smearing between the states and , see Tab. 1. An equilibrium occupation with 8 electrons is in agreement with the HF calculations of Fig. 6. In Table 2 the composition of the many-body states in terms of the five localized molecular orbitals of Fig. 3 is quantified in terms of the average populations of the single particle states obtained from exact diagonalization of (see Eq. (10)). The composition of the states and is similar to the HF average populations. As discussed in Sec. VII and shown in Fig. 8, the LMO occupations change at finite bias.

NGF 1.9725 1.0934 1.9989 1.9294 1.0220
1.9801 1.1337 1.9973 1.8685 1.0204
1.9798 1.1339 1.9972 1.8685 1.0207
1.0390 1.8047 1.9997 1.1753 1.9815
1.0269 1.8159 1.9995 1.1764 1.9812
1.9990 1.4657 1.9992 1.5365 1.9996
1.9510 2.0000 1.9989 1.9999 1.0502
1.0772 1.9917 1.9998 1.9780 1.9532
Table 2: The average populations of the single-particle states at zero bias voltage. Calculations within the NGF method are shown in the second line. They agree with the composition of the states and obtained from exact diagonalization of (see Eq. (10)).

Vii Transport characteristics

We are now able to calculate and interpret the current-voltage characteristics of the benzenedithiol junction. The current at finite voltage, which is given by


is presented together with the differential conductance in Fig. 7. The curves are asymmetric with respect to a bias inversion because the junction geometry was chosen to be slightly asymmetric. As the main result we find a multi-scale Coulomb blockade. The large region of suppressed current is about 2 Volt wide. However, the current is completely blocked only in a much smaller region of bias voltage, as small steps in the current (peaks in the conductance) are present at lower biases.

Figure 7: (Color online) Current-voltage curve (black solid line) and differential conductance (red dashed line).

As a first step in the analysis of the current voltage characteristics we consider the average particle number in the central system presented in the left panel of Fig. 8. At low biases the average particle number is 8 corresponding to the neutral configuration of benzenedithiol. The many-body state with the minimal grand-canonical energy () is in fact the 8-particles ground state (see Fig. 10). When the bias drop is raised in the junction the average particle number takes values between 8 and 9 ensuring that the dominant transitions are negative ion resonances.

A further insight into the dynamics is obtained by monitoring the average occupation of the different localized molecular orbitals shown in the right panel of Fig. 8. At low biases the symmetric central orbital (the third orbital from the top in Fig. 3) is completely occupied, . Its occupation undergoes a sensible variation only at the voltages of the large current steps . Large variations in the population of the asymmetric LMOs centered around the molecule-lead interface (orbitals 1, 2, 4 and 5 in Fig. 3) are instead associated to the small current steps present at lower voltages. Interestingly, at a bias of the effective spatial symmetry of the system is recovered with the populations of the asymmetric states being all equal.

A deeper understanding of the dynamics of the system is obtained by the analysis of the occupation of the many-body states (Fig. 9), their energies (Tab. 1), and the transition rates among them (Tab. 3, schematically represented in Fig. 10). If the calculation of the current is performed taking into account hundreds of many-body states, the essential physics at the biases presented in Fig. 7 is captured by considering the lowest four 8-particle levels (for a total of 8 states) and the lowest three 9-particle ones (6 states).

Figure 8: Bias dependence of the average electron number (left) and average individual populations (right) on the molecular junction.

The tunnelling events from (to) the source or the drain connect these many-body states. The tunnelling rate is the product of a geometrical part ( of Eq. (42)) and an energetic contribution encoding the energy conservation in the tunnelling event and the Pauli exclusion principle (see Eq. (41)). The energetic contribution ensures that the rate changes (and correspondingly the current through the system) every time the resonant condition is fulfilled. With this argument it is already possible to assign a specific transition to most of the peaks in the conductance of Fig. 7. In particular the transitions are associated with the peak at the most negative bias and to the second peak from the left. The first small peak at positive bias is anomalous and we will return to it later. We only note that its position depends on the temperature and that it moves to the resonance at low temperatures. The rightmost conductance peaks are instead associated to the transitions and , respectively.

Figure 9: The occupation of the manybody states
Figure 10: The energies, tunneling rates and the associated populations of the most relevant states (see discussion in the text).

The approximate symmetries of molecular geometry are very important since they introduce selection rules which distinguish between transitions which are energetically equally allowed. In Tab. 3 we report the transition rates between the different many-body states. Here the values are given in and the spin is chosen to fulfill spin conservation in the tunnelling event. In the case of a doublet to triplet transition the value of the rate reported is the one involving the triplet state with maximum projection along the quantization axis. Except for the transition all transitions show a very pronounced left-right asymmetry. It is much easier for example for an electron to tunnel in (or out) of the molecule from (to) the left instead of the right lead when this tunnelling event involves the many-body eigenstates and . This asymmetry is essential to explain the dynamics of the system at low biases and can be understood in terms of the spatial distribution of the many-body eigenstates.

0.0017 0.1025 0.0024
0.0056 0.2039 0.004
0.1033 0.0003 0.0074
0.1676 0.0013 0.0488

0.0442 0.0127 0
0.0866 0.0278 0.0022
0.0044 0.0056 0.1157
0.0013 0.011 0.2185
Table 3: The transition rates between the different manybody states.

The left transition rate is larger than the right one when the transition from an 8 to a 9-particle state is associated to a larger variation of the density in the orbitals 1 or 2 than in the orbitals 4 or 5. Analogous arguments hold for the reverse situation.

Let us now return to the interpretation of the current voltage characteristics with the help of Fig. 10. By convention, to a positive bias voltage corresponds a stationary particle current flowing from right to left while the electrical current flows from left to right. We concentrate first on the negative bias. From an accurate analysis of the definition of the tunnelling rates (Eq. (41)) it is not difficult to prove that the first step in the current is due to the resonant condition between the and states at the left lead. Current flows since the system oscillates between the and states by receiving an electron from the left lead and by releasing it to the right one. The asymmetry between the transition rates, , ensures than even after the opening of the current channel the occupation of the (together with the almost degenerate ) is still the largest one. In the right panel of Fig. 10 we schematically represent the tunnelling rates and the associated populations of the most relevant levels for a bias just above (in absolute value) the first negative bias conductance peak. Starting from this population distribution it is then natural to observe the next visible current step related to the transition . Since this time the left tunnelling rate dominates, the population of the 8-particle states decreases substantially in favor of the 9-particle ones. Generally, a more uniform mixing of states with different particle number is associated with a larger fluctuation of the number of electrons in the central system and thus with a larger current.

The dynamics at positive bias is more complex. In particular the first conductance peak occurs at a bias at which even the ground to ground state transition is not yet open. This anomalous behaviour is understandable when taking into account the large left-right asymmetry of the rates. As schematically represented in the left panel of Fig. 10, even before the (right lead) resonance between the and the state opens a conventional current channel, the states and get strongly populated. The fundamental reason is the large probability to tunnel out of the system at the left lead through the transition which is also energetically favorable. Very soon the states and become the new effective ground states for the system (see Fig. 9). In this scheme it is thus not surprising that i) the first conductance peak is located at an ”average” between the resonance and the one; ii) the next two conductance peaks at positive bias occur at the and resonant conditions.

Viii Conclusions

In conclusion, we developed a many-body localized molecular orbital approach to transport through molecular junctions with the following protocol:

  1. Geometry optimization using DFT and hybrid DFT (usually B3LYP based) methods.

  2. Molecular vibrons can be calculated after the geometry optimization (not considered in this paper).

  3. Molecular orbitals of the extended molecule are obtained. Localized molecular orbitals (LMOs) are constructed and form the basis for all subsequent calculations.

  4. A Hubbard interaction is introduced for the LMOs in the central region: only density-density Coulomb integrals are taken into account.

  5. Electron-vibron interaction can be included in the central region (not considered in this paper).

  6. The leads are kept as effectively noninteracting (mean-field approximation). The interaction Hamiltonian between leads and central region yields the relevant tunneling couplings.

  7. A spectral analysis and transport calculations are performed on the basis of the ab initio based Hubbard-Anderson model.

Using the benchmark example of a benzene-dithiol molecular junction, we performed the full line of calculations in the framework of this approach. We determined the geometry of the junction, calculated molecular orbitals and transformed them into localized molecular orbitals. Upon using an energy range of about 4 eV around the Fermi energy of gold, we obtained a basis of 5 LMOs with energies . Then we calculated the Coulomb matrix elements for these orbitals and coupling matrix elements between the central region and the leads. Using the parameters , and , obtained from ab initio calculations, we calculated the spectral function in the framework of the nonequilibrium Green function approach (in the RHF, HF and NEOM approximations). Besides, the model was transformed into the many-body eigenstate basis, and the quantum master equation (applied in the sequential tunneling limit) was used to calculate the current. It is shown that transport through asymmetrically-coupled molecular edge states results in suppressed peaks of the differential conductance at small voltage and unexpectedly large peaks at higher voltages. The origin of these anomalies could be explained upon analyzing the occupation probabilities of the many-body states as well as their compositions in terms of LMOs. In general, we could qualitatively understand the equilibrium state and main transport properties of the considered molecular junction with strong electron-electron interaction and intermediate coupling to the leads.

Nevertheless, the further development of the theory is necessary with respect to both ab initio and quantum transport aspects. The results presented in this paper are only partially self-consistent because the parameters , and are calculated at zero voltage, but used at all voltages. Actually, it is possible to extend the theory to include the recalculation of the parameters at finite voltage and the influence of the nonequilibrium charge in the central region on the leads. A related issue is the effect of the external field on the LMOs energies, which we treat using a simplified linear approximation. The Hubbard interaction plays the main role, but the corrections due to non density-density interactions and polarization of the molecule can be important as well. Finally, we expect that the method proposed in Ref. Kern and Grifoni, 2012 could be of importance to treat the parameter regime typical for molecular junctions with intermediate coupling to the leads.


The ab initio calculations were done by the quantum chemistry code Firefly Granovsky () and partially by the DFT code Siesta Sie (). The results were analyzed and the LMOs were visualized with the help of MacMolPlt Bode and Gordon (1998) and Chemcraft. Many-body modeling and transport calculations were performed by our own codes. We thank Michael Hartung for his help with the local Linux cluster used for the numerical calculations.

This work was funded by Deutsche Forschungsgemeinschaft within the Priority Program SPP 1243 and Collaborative Research Center SFB 689.


  • Cuniberti et al. (2005) G. Cuniberti, G. Fagas, and K. Richter (Eds.), Introducing Molecular Electronics, vol. 680 of Lecture Notes in Physics (Springer, Berlin, 2005).
  • Cuevas and Scheer (2010) J. C. Cuevas and E. Scheer, Molecular electronics: An Introduction to Theory and Experiment (World Scientific, 2010).
  • Song et al. (2011) H. Song, M. A. Reed, and T. Lee, Advanced Materials 23, 1583 (2011).
  • Park et al. (2002) J. Park, A. N. Pasupathy, J. I. Goldsmit, C. Chang, Y. Yaish, J. R. Petta, M. Rinkoski, J. P. Sethna, H. D. Abruna, P. L. McEuen, D. C. Ralph, Nature 417, 722 (2002).
  • Kubatkin et al. (2003) S. Kubatkin, A. Danilov, M. Hjort, J. Cornil, J.-L. Brédas, N. Stuhr-Hansen, P. Hedegård, and T. Bjørnholm, Nature 425, 698 (2003).
  • Osorio et al. (2007) E. A. Osorio, K. O’Neill, N. Stuhr-Hansen, O. F. Nielsen, T. Bjørnholm, and H. S. J. van der Zant, Adv. Mater. 19, 281 (2007).
  • Smit et al. (2002) R. H. M. Smit, Y. Noat, C. Untiedt, N. D. Lang, M. C. van Hemert, and J. M. van Ruitenbeek, Nature 419, 906 (2002).
  • Tal et al. (2008) O. Tal, M. Krieger, B. Leerink, and J. M. van Ruitenbeek, Phys. Rev. Lett. 100, 196804 (2008).
  • Tal et al. (2009) O. Tal, M. Kiguchi, W. H. A. Thijssen, D. Djukic, C. Untiedt, R. H. M. Smit, and J. M. van Ruitenbeek, Phys. Rev. B 80, 085427 (2009).
  • Reed et al. (1997) M. A. Reed, C. Zhou, C. J. Muller, T. P. Burgin, and J. M. Tour, Science 278, 252 (1997).
  • Lörtscher et al. (2007) E. Lörtscher, H. B. Weber, and H. Riel, Phys. Rev. Lett. 98, 176807 (2007).
  • Kim et al. (2011) Y. Kim, T. Pietsch, A. Erbe, W. Belzig, and E. Scheer, Nano Letters 11, 3734 (2011).
  • Taylor et al. (2001) J. Taylor, H. Guo, and J. Wang, Phys. Rev. B 63, 121104 (2001).
  • Damle et al. (2001) P. S. Damle, A. W. Ghosh, and S. Datta, Phys. Rev. B 64, 201403 (2001).
  • Xue and Ratner (2003) Y. Xue and M. A. Ratner, Phys. Rev. B 68, 115406 (2003).
  • Th. Frauenheim et al. (2002) Th. Frauenheim, G. Seifert, M. Elstner, T. Niehaus, C. Köhler, M. Amkreutz, M. Sternberg, Z. Hajnal, A. Di Carlo, and S. Suhai, J. Phys.: Condens. Matter 14, 3015 (2002).
  • Brandbyge et al. (2002) M. Brandbyge, J.-L. Mozos, P. Ordejón, J. Taylor, and K. Stokbro, Phys. Rev. B 65, 165401 (2002).
  • Pecchia and Di Carlo (2004) A. Pecchia and A. Di Carlo, Rep. Prog. Phys. 67, 1497 (2004).
  • Ke et al. (2004) S.-H. Ke, H. U. Baranger, and W. Yang, Phys. Rev. B 70, 085410 (2004).
  • Kadanoff and Baym (1962) L. Kadanoff and G. Baym, Quantum Statistical Mechanics (Benjamin, New York, 1962).
  • Keldysh (1964) L. V. Keldysh, Zh. Eksp. Teor. Fiz. 47, 1515 (1964), [Sov. Phys. JETP 20, 1018 (1965)].
  • Langreth (1976) D. Langreth, in Linear and Nonlinear Electron Transport in Solids, edited by J. Devreese and E. van Doren (Plenum, New York, 1976).
  • Rammer and Smith (1986) J. Rammer and H. Smith, Rev. Mod. Phys. 58, 323 (1986).
  • Haug and Jauho (1996) H. Haug and A.-P. Jauho, Quantum Kinetics in Transport and Optics of Semiconductors, vol. 123 of Springer Series in Solid-State Physics (Springer, Berlin, 1996).
  • Blum (1996) K. Blum, Density Matrix Theory and Applications (Plenum Press, New York, 1996).
  • Breuer and Petruccione (2002) H. P. Breuer and F. Petruccione, The theory of open quantum systems (Oxford University Press, Oxford, 2002).
  • König et al. (1996) J. König, J. Schmid, H. Schoeller, and G. Schön, Phys. Rev. B 54, 16820 (1996).
  • Costi et al. (1994) T. A. Costi, A. C. Hewson, and V. Zlatic, J. Phys.: Cond. Matter 6, 2519 (1994).
  • Andergassen et al. (2010) S. Andergassen, V. Meden, H. Schoeller, J. Splettstoesser, and M. R. Wegewijs, Nanotechnology 21, 272001 (2010).
  • Segal et al. (2010) D. Segal, A. J. Millis, and D. R. Reichman, Phys. Rev. B 82, 205323 (2010).
  • Pletyukhov and Schoeller (2012) M. Pletyukhov and H. Schoeller, Phys. Rev. Lett. 108, 260601 (2012).
  • Smirnov and Grifoni (2012) S. Smirnov and M. Grifoni, arXiv:1203.4360 (2012).
  • Hettler et al. (2003) M. H. Hettler, W. Wenzel, M. R. Wegewijs, and H. Schoeller, Phys. Rev. Lett. 90, 076805 (2003).
  • Delaney and Greer (2004) P. Delaney and J. C. Greer, Phys. Rev. Lett. 93, 036805 (2004).
  • Ferretti et al. (2005) A. Ferretti, A. Calzolari, R. Di Felice, F. Manghi, M. J. Caldas, M. B. Nardelli, and E. Molinari, Phys. Rev. Lett. 94, 116802 (2005).
  • Mirjani and Thijssen (2011) F. Mirjani and J. M. Thijssen, Phys. Rev. B 83, 035415 (2011).
  • Bergfield et al. (2011) J. P. Bergfield, G. C. Solomon, C. A. Stafford, and M. A. Ratner, Nano Letters 11, 2759 (2011).
  • Frederiksen et al. (2007) T. Frederiksen, M. Paulsson, M. Brandbyge, and A.-P. Jauho, Phys. Rev. B 75, 205413 (2007).
  • Thygesen and Rubio (2007) K. S. Thygesen and A. Rubio, J. Chem. Phys. 126, 091101 (2007).
  • Thygesen and Rubio (2008) K. S. Thygesen and A. Rubio, Phys. Rev. B 77, 115333 (2008).
  • Korytár et al. (2010) R. Korytár, M. Pruneda, J. Junquera, P. Ordejón, and N. Lorente, J. Phys.: Condens. Matter 22, 385601 (2010).
  • Korytár and Lorente (2011) R. Korytár and N. Lorente, J. Phys.: Condens. Matter 23, 355009 (2011).
  • Greuling et al. (2011) A. Greuling, M. Rohlfing, R. Temirov, F. S. Tautz, and F. B. Anders, Phys. Rev. B 84, 125413 (2011).
  • Karolak et al. (2011) M. Karolak, D. Jacob, and A. I. Lichtenstein, Phys. Rev. Lett. 107, 146604 (2011).
  • Tröster et al. (2012) P. Tröster, P. Schmitteckert, and F. Evers, Phys. Rev. B 85, 115409 (2012).
  • Strange et al. (2011) M. Strange, C. Rostgaard, H. Häkkinen, and K. S. Thygesen, Phys. Rev. B 83, 115108 (2011).
  • Rangel et al. (2011) T. Rangel, A. Ferretti, P. E. Trevisanutto, V. Olevano, and G.-M. Rignanese, Phys. Rev. B 84, 045426 (2011).
  • Begemann et al. (2008) G. Begemann, D. Darau, A. Donarini, and M. Grifoni, Phys. Rev. B 77, 201406 (2008).
  • Bergfield and Stafford (2009) J. P. Bergfield and C. A. Stafford, Phys. Rev. B 79, 245125 (2009).
  • Sobczyk et al. (2012) S. Sobczyk, A. Donarini, and M. Grifoni, Phys. Rev. B 85, 205408 (2012).
  • (51) Siesta version 3.0,
  • (52) A. A. Granovsky, Firefly version 7.1.G,
  • Foster and Boys (1960) J. M. Foster and S. F. Boys, Rev. Mod. Phys. 32, 300 (1960).
  • Pariser and Parr (1953) R. Pariser and R. Parr, J. Chem. Phys. 21, 466 (1953), 21, 767 (1953).
  • Kern and Grifoni (2012) J. Kern and M. Grifoni, arXiv:1209.4995 (2012).
  • Meir and Wingreen (1992) Y. Meir and N. S. Wingreen, Phys. Rev. Lett. 68, 2512 (1992).
  • Jauho et al. (1994) A.-P. Jauho, N. S. Wingreen, and Y. Meir, Phys. Rev. B 50, 5528 (1994).
  • Jauho (2006) A.-P. Jauho, Journal of Physics: Conference Series 35, 313 (2006).
  • Ryndyk et al. (2009) D. A. Ryndyk, R. Gutiérrez, B. Song, and G. Cuniberti, Energy Flow Dynamics in Biomaterial Systems (Springer, Berlin, 2009), chap. Green function techniques in the treatment of quantum transport at the molecular scale, p. 213.
  • Song et al. (2007) B. Song, D. A. Ryndyk, and G. Cuniberti, Phys. Rev. B 76, 045408 (2007).
  • Bode and Gordon (1998) B. M. Bode and M. S. Gordon, J. Mol. Graphics Mod. 16, 133 (1998).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description