Manybody localized molecular orbital approach to molecular transport
Abstract
An ab initio based theoretical approach to describe nonequilibrium manybody effects in molecular transport is developed. We introduce a basis of localized molecular orbitals and formulate the manybody model in this basis. In particular, the HubbardAnderson Hamiltonian is derived for singlemolecule junctions with intermediate coupling to the leads. As an example we consider a benzenedithiol junction with gold electrodes. An effective fewlevel model is obtained, from which spectral and transport properties are computed and analyzed. Electronelectron 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.
pacs:
73.63.b, 85.65.+hI 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). Electronelectron 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 goldbenzenedithiolgold (AuBDTAu) 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 KohnSham equation as physical quasiparticle states can not be rigorously established. Besides, the DFT+NGF is a meanfield theory and lacks to describe manybody 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.
On the other hand, modelbased 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 weakcoupling limit. The quantum master equation is usually formulated in the basis of manybody 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 manybody 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 manybody 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), manybody 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 manybody 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 (singleatom) degrees of freedom must be reduced to an effective fewelectron (or few manybody state) interacting model, as a prerequisite of successive manybody transport methods. The main building blocks of our approach are an effective electron model (an orthonormal basis of localized molecular orbitals and a manybody 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 socalled 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 HartreeFock 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 manybody 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 manybody HubbardAnderson 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 meanfield approximation.
iii) Using the HubbardAnderson Hamiltonian, manybody 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 electronelectron interaction and Hubbard approximation are discussed in Sec. III. The parameters of the effective HubbardAnderson 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 manybody 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 fullelectron 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 electronelectron and electronvibron 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.
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.
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
(1) 
The indices with bars denote the states without the spin degree of freedom. We apply the FosterBoys 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
(2) 
is no longer diagonal. Moreover, the related Hamiltonian takes the form
(3) 
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 electronelectron interactions are described by the Hamiltonian
(4) 
where and denotes the spin. The matrix elements for the scalar Coulomb interaction are defined as
(5)  
(6) 
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
(7) 
describing only densitydensity interactions with and the Hubbard parameters defined as
(8) 
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 semiempirical PPP model Pariser and Parr (1953). In this way all coefficients are derived, but further semiempirical corrections could be included for better agreement with experiments.
Iv Manybody 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 electronelectron interaction Hamiltonian , the Hamiltonians of the ends of the leads , and the tunneling Hamiltonian describing the moleculetolead coupling:
(9) 
In our case the central Hamiltonian is replaced by the Hubbard cluster Hamiltonian:
(10) 
where are the bare energies of the LMOs, including the shifts due to an external voltage.
The noninteracting molecular Hamiltonian is obtained from the LMO Hamiltonian , Eq. (3). The zerovoltage 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:
(11) 
where
(12) 
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
(13) 
and the Hamiltonians of the left and right leads are
(14) 
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 meanfield (DFT) level. The equilibrium electrodes, which can have different electrochemical potentials, determine the boundary conditions for the leads.
For the 5level model (Fig. 3), which represents actually 10 levels with spin, we obtain the following parameters (for spin up or down, all numbers are in ):
(15) 
The matrix of the Hubbard parameters calculated from expression (8) reads
(16) 
The diagonal entries in this matrix correspond to the interaction of two electrons in the same LMO state but with different spin. The offdiagonal terms denote coupling beween two different LMOs spin of the electronic state. Of course, these expressions should be transformed into the 10level 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 timeconsuming.
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 offdiagonal 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 (socalled wideband limit). In this approximation the levelwidth function
(17) 
is energy independent.
We now return to our 5level 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 levelwidth 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 benzenebased 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 singlelevel 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 DysonKeldysh equations in the integral form
(18)  
(19) 
or from the corresponding equations in the differential form
(20) 
(21) 
Here
(22) 
is the total selfenergy of the molecule composed of the interaction selfenergy and the tunneling (coupling to the left (L) and right (R) lead) selfenergies
(23)  
(24) 
where is the Green function of the leads.
The retarded tunneling selfenergy can be represented as
(25) 
where is the real part of the selfenergy, which usually can be included in the singleparticle Hamiltonian , and describes level broadening due to coupling to the leads. In the case of noninteracting leads with continuous energy spectra, the levelwidth function is determined by the expression (17).
For the corresponding lesser function of the noninteracting leads one finds
(26) 
where
(27) 
is the equilibrium Fermi distribution function with chemical potential ().
The expression for the interaction selfenergy can not be obtained exactly. In the nonequilibrium HartreeFock approximation one has Ryndyk et al. (2009)
(28)  
(29) 
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 spinunpolarized leads)
(30) 
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
(31) 
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.
The solution of Eq. (18) for the central part is in this case
(32) 
with the lead selfenergies (see Eq. (23)) and levelwidth functions defined as
(33)  
(34) 
The lead Green functions for noninteracting leads (or for leads described by the effective meanfield Hamiltonians ) are defined correspondingly by the electrode selfenergies
(35)  
(36) 
The calculation of the electrode selfenergies 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 selfconsistently within four approximations: initial DFT with the meanfield 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,
(37) 
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 singleparticle 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 manybody spectrum:
exact diagonalization and the ground state properties
To gain insight into the manybody energy spectrum of the central system we perform an exact diagonalization of Eq. (10) obtaining the set of manybody eigenstates . Calculating the tunneling matrix elements we obtain from Eqs. (10) and (13) the Hamiltonian
(38)  
(39) 
First, we analyze the manybody spectrum. With 5 LMOs we get 1024 manybody states in the Fock space. The lowest 8particle 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 singlettriplet quasi degenerate states. The 9particle 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 8particle levels lie all below the one of the 9particle ground state. The 7particle 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 
Note that in the sequentialtunneling master equation method the exact manybody 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 selfenergy, Eq. (23), describes the energy shift of the singleparticle levels.
The population probabilities are found from the master equation
(40) 
where the tunneling rates are, in second order in the tunneling Hamiltonian,
(41) 
with
(42) 
This expression connects the tunneling rates to the levelwidth function; thus the calculated by the NGF method can be used, see Eqs. (33,34). In the wideband 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 manybody 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 8particle and doublets for the 9particle 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. (4042) is straightforward and can be obtained by direct numerical integrations in stationary and timedependent 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 manybody 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 
Vii Transport characteristics
We are now able to calculate and interpret the currentvoltage characteristics of the benzenedithiol junction. The current at finite voltage, which is given by
(43) 
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 multiscale 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.
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 manybody state with the minimal grandcanonical energy () is in fact the 8particles 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 moleculelead 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 manybody 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 manybody states, the essential physics at the biases presented in Fig. 7 is captured by considering the lowest four 8particle levels (for a total of 8 states) and the lowest three 9particle ones (6 states).
The tunnelling events from (to) the source or the drain connect these manybody 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.
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 manybody 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 leftright 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 manybody 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 manybody 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 
The left transition rate is larger than the right one when the transition from an 8 to a 9particle 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 8particle states decreases substantially in favor of the 9particle 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 leftright 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 manybody localized molecular orbital approach to transport through molecular junctions with the following protocol:

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

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

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

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

Electronvibron interaction can be included in the central region (not considered in this paper).

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

A spectral analysis and transport calculations are performed on the basis of the ab initio based HubbardAnderson model.
Using the benchmark example of a benzenedithiol 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 manybody 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 asymmetricallycoupled 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 manybody 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 electronelectron 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 selfconsistent 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 densitydensity 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.
Acknowledgments
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. Manybody 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.
References
 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. StuhrHansen, P. Hedegård, and T. Bjørnholm, Nature 425, 698 (2003).
 Osorio et al. (2007) E. A. Osorio, K. O’Neill, N. StuhrHansen, 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 SolidState 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, http://www.uam.es/siesta.
 (52) A. A. Granovsky, Firefly version 7.1.G, http://classic.chem.msu.su/gran/firefly/index.html.
 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).