First-principles molecular transport calculation for the benzenedithiolate molecule

TL;DR

This paper employs first-principles calculations combining DFT and NEGF to analyze how different electrodes and many-body effects influence electron transport in benzenedithiolate molecular junctions.

Contribution

It introduces a comprehensive first-principles methodology to study molecular transport, including symmetry analysis and many-body effects with an extended Hubbard model.

Findings

Conductance mechanisms differ between gold and platinum electrodes due to symmetry.

Many-body effects significantly influence the charge stability and transport properties.

Electrode material and many-body interactions are crucial for molecular device performance.

Abstract

A first-principles approach based on Density Functional Theory and Non-Equilibrium Green's functions is used to study the molecular transport system consisting of benzenedithiolate connected with monoatomic gold and platinum electrodes. Using symmetry arguments we explain why the conductance mechanism is different for gold and platinum electrodes. We present the charge stability diagram for the benzenedithiolate connected with monoatomic platinum electrodes including many-body effects in terms of an extended Hubbard Hamiltonian and discuss how the electrodes and the many-body effects influence the transport properties of the system.

Tables

Table 1. Table 1: Many-body eigenenergies of the central region in case of Pt-BDT-Pt.

Level	Energy [eV]	Spin [$\hslash$]	Level	Energy [eV]	Spin [$\hslash$]
${10}_g$	-203.08	0
${10}^{\prime }$	-200.81	1	$9_g$	-202.48	1/2
${10}^{\prime \prime }$	-200.24	1	$9^{\prime }$	-200.95	1/2
${10}^{\prime \prime \prime }$	-200.11	0	$9^{\prime \prime }$	-199.46	1/2

Equations

$$\Gamma \gg \Delta E,U,k_{B}T,$$

$$\Gamma <\Delta E,U$$

$${\mathop{\hat{H}}}=\left(\begin{array}[]{ccccc}{\mathop{\hat{H}}}_{L_{L}}&{\mathop{\hat{V}}}_{L_{L}T_{L}}&&&{\mathop{\hat{V}}}_{L_{L}L_{R}}\\ {\mathop{\hat{V}}}_{L_{L}T_{L}}^{\dagger}&{\mathop{\hat{H}}}_{T_{L}}&{\mathop{\hat{V}}}_{T_{L}C}&\\ &{\mathop{\hat{V}}}_{T_{L}C}^{\dagger}&{\mathop{\hat{H}}}_{C}&{\mathop{\hat{V}}}_{CT_{R}}&\\ &&{\mathop{\hat{V}}}_{CT_{R}}^{\dagger}&{\mathop{\hat{H}}}_{T_{R}}&{\mathop{\hat{V}}}_{T_{R}L_{R}}\\ {\mathop{\hat{V}}}_{L_{L}L_{R}}^{\dagger}&&&{\mathop{\hat{V}}}_{T_{R}L_{R}}^{\dagger}&{\mathop{\hat{H}}}_{L_{R}}\\ \end{array}\right)\;.$$

$${\hat{H}}_C={\hat{H}}_0+{\hat{H}}_1$$

$${\hat{H}}_0=\sum _{ij}^N{\tilde{t}}_{ij}{\hat{a}}_i^{†}{\hat{a}}_j.$$

$${\tilde{t}}_{ij}=\{QE{Q}^{†}{\}}_{\bar{i}\bar{j}},$$

$${\hat{H}}_1=\frac{1}{2}\sum _{i\ne j}^N{U}_{ij}{\hat{n}}_i{\hat{n}}_j.$$

$$U_{ij}=\int d^{3}r\int d^3{r}^{\prime }{\left|{\Psi }_i(\mathbf{r})\right|}^2{\left|{\Psi }_j({\mathbf{r}}^{\prime })\right|}^{2}U(\mathbf{r},{\mathbf{r}}^{\prime })$$

$$U(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{1}{4\pi {ϵ}_0\eta |\mathbf{r}-{\mathbf{r}}^{\prime }|}$$

$$t_{\bar{i}\bar{j}}^C={\tilde{t}}_{\bar{i}\bar{j}}-\Delta t_{\bar{i}\bar{i}}{\delta }_{\bar{i}\bar{j}}.$$

$$\Delta t_{\bar{i}\bar{i}}$$

$${\hat{H}}_C=\sum _{ij}^N{t}_{ij}^C{\hat{a}}_i^{†}{\hat{a}}_j+\frac{1}{2}\sum _{i\ne j}^N{U}_{ij}{\hat{n}}_i{\hat{n}}_j.$$

$$I=\frac{e}{h}\int_{-\infty}^{\infty}d\omega\sum_{\begin{array}[]{c}i\in T_{L}\\ j\in C\end{array}}\Re\left(V_{ij}\mathcal{G}^{k}_{ij}(\omega)\right)\;.$$

$$I=\frac{e}{h}{\int }_{-\infty }^{\infty }d\omega \left(f_R(\omega )-f_L(\omega )\right)T(\omega )$$

$$T(\omega )=\text{Tr}\left({\mathbf{\Gamma }}_L{\mathcal{G}}_{CC}^r{\mathbf{\Gamma }}_R{\mathcal{G}}_{CC}^a\right),$$

$${\mathbf{\Delta }}_{\nu }^r={\mathbf{R}}_{\nu }-\frac{i}{2}{\mathbf{\Gamma }}_{\nu }.$$

$$\displaystyle\big{(}\bm{{\mathcal{G}}}_{CC}^{r,a}\big{)}^{-1}$$

$$\big{(}{\bf g}^{r,a}_{CC}\big{)}^{-1}=\big{(}{\bf g}^{r,a}_{CC,0}\big{)}^{-1}-{\bm{\Sigma}}^{r,a}$$

$${\left\{{\mathbf{g}}_{CC}^r\right\}}_{ij}=$$

$${\left({\mathbf{G}}_{T_{\nu }{T}_{\nu }}\right)}^{-1}={\left({\mathbf{g}}_{T_{\nu }{T}_{\nu }}\right)}^{-1}-{\mathbf{V}}_{T_{\nu }{L}_{\nu }}{\mathbf{g}}_{L_{\nu }{L}_{\nu }}{\mathbf{V}}_{L_{\nu }{T}_{\nu }}.$$

$$A_{\nu \nu }(\omega )$$

$$G$$

$${\omega }_1-\frac{V_b}{2}$$

$${\omega }_2-\frac{V_b}{2}$$

$${\omega }_2-\frac{V_b}{2}=+\frac{V_b}{2}.$$

$${V}_{T_{L}C}=\left(\begin{array}[]{cccccccc}-1.11&0.08&0.01&-0.34&0.14&-0.02&-0.02&0.00\\ -0.85&0.11&0.00&0.08&-0.05&-0.02&0.01&0.00\end{array}\right)$$

$$\tilde{t}=\left(\begin{array}[]{cccccccc}-1.46&-2.12&0.21&0.17&0.07&0&0.05&-0.03\\ -2.12&-0.29&-2.58&-2.54&0.01&0.27&-0.23&0.05\\ 0.21&-2.58&0.16&0.01&-0.08&-2.77&0.02&0.07\\ 0.17&-2.54&0.01&0.53&-2.77&-0.24&0.27&0\\ 0.07&0.01&-0.08&-2.77&0.16&0.01&-2.58&0.21\\ 0&0.27&-2.77&-0.24&0.01&0.52&-2.54&0.18\\ 0.05&-0.23&0.02&0.27&-2.58&-2.54&-0.29&-2.12\\ -0.03&0.05&0.07&0&0.21&0.18&-2.12&-1.46\end{array}\right)$$

$${V}_{CT_{R}}^{\dagger}=\left(\begin{array}[]{cccccccc}0.00&0.02&-0.14&0.03&-0.01&0.34&-0.08&1.11\\ 0.00&-0.01&0.05&0.02&0.00&-0.08&-0.11&0.85\end{array}\right)$$

$$U=\left(\begin{array}[]{cccccccc}7.72&4.86&3.24&3.25&2.27&2.21&1.96&1.49\\ 4.86&8.52&5.06&4.97&3.34&3.23&2.85&1.96\\ 3.24&5.06&8.67&3.32&3.00&5.04&3.34&2.27\\ 3.25&4.97&3.32&7.93&5.03&2.86&3.23&2.22\\ 2.27&3.34&3.00&5.03&8.67&3.32&5.06&3.24\\ 2.21&3.23&5.04&2.86&3.32&7.94&4.97&3.26\\ 1.96&2.85&3.34&3.23&5.06&4.97&8.52&4.86\\ 1.49&1.96&2.27&2.22&3.24&3.26&4.86&7.72\\ \end{array}\right)\;.$$

Full text

1 Introduction

First-principles calculations in the field of molecular electronics represent a major challenge in computational physics. In order to build electronic devices based on organic molecules it is crucial to understand the underlying charge transport phenomena in detail to be able to simulate the actual behaviour [1, 2]. However, it is a challenging task to predict quantitatively the experimental characteristics of molecular junctions. First of all, in many cases the charge transport properties of a system depend on the detailed contact geometry [3, 4, 5], which is experimentally difficult to identify. But even if all atomic positions were known, non-equilibrium calculations are always based on approximations which affect the results. Different theoretical methods exist for different transport regimes. In the coherent regime electrons proceed elastically, without exchanging energy. This is the case if there is a strong coupling between the leads and the molecule, i.e. if

[TABLE]

here, $\Gamma$ is the coupling between the leads and the molecule, $\Delta E$ the level spacing, $U$ the Hubbard interaction parameter and $k_{B}T$ the temperature. In the coherent regime, the energy levels of the electrodes and the molecule are hybridized and charge quantization is suppressed and the Landauer formalism becomes exact. The combination of Density Functional Theory and non-equilibrium Green’s functions (DFT+NEGF) has become a standard tool for coherent ab-initio transport calculations but there are also other approaches like the Lippmann-Schwinger method [6] or the wave-function matching technique [7]. In the incoherent regime the time that an electron spends on the system is sufficiently long for it to interact with other particles (electrons, phonons, etc.). In the case of electron-electron interaction we have

[TABLE]

for the incoherent regime. If the leads are weakly coupled to the system then there is charge quantization and the electrons propagate by sequential tunnelling. This incoherent regime can be modelled by a Hubbard-Hamiltonian for describing the molecule and non-interacting tight-binding leads. Non-equilibrium cluster perturbation theory (n-CPT) [8, 9] or the variational cluster approach (VCA) [10, 9] are techniques for solving this model accurately in the coherent limit. In the case of weak coupling between the molecule and the leads master equation (ME) [11] approaches became useful. The auxiliary master equation approach (AMEA) [12, 13, 14] is in principle exact, also for strong coupling between the molecule and the leads, but is limited to systems with only a view interacting sites in the central region. In order to calculate transport properties from first principles one can combine these techniques with DFT and carry out a charge self-consistency procedure in non-equilibrium, which is often neglected for simplicity. In our work we apply DFT+NEGF and an extended version of DFT+NEGF, which we denote by DFT+CPT, for calculations including electron-electron interactions. The approximation in DFT+CPT is to approximate the self-energy of the system by the self-energy of the isolated central region. The method allows treating larger molecules in the central region and is accurate in the coherent limit. A more sophisticated but numerically challenging technique for calculating transport properties is the combination of DFT with GW [15, 16, 17]. DFT can also be used in combination with the dynamical mean field theory (DMFT) [18], with ME [19] if the system is in the weak coupling regime or even with ME+CPT [20].

Although the 1,4-benzenedithiolate (BDT) molecule is usually considered in the field of molecular electronics there are still severe discrepancies between experimental [21, 22, 23, 24, 25, 26] and theoretical results [27, 19, 28, 29, 30, 6]. The calculated conductances tend to be much higher than the measured ones. Proposed reasons are the sensitivity of the conductance to the detailed geometry of the junction [29, 6], hydration of the BDT molecule to benzenedithiol [31, 16] and correlation effects [3, 32].

The purpose of this work is to to study the transport properties of the BDT molecule connected to monoatomic metal chains as leads. This system allows to study the impact of the junction geometry on the conductance and how transport is affected by low dimensional leads. Many publications [28, 27, 30] explain the transport through the gold-BDT-gold system with a single-level model where the conductance is carried by a single transport channel corresponding to the HOMO orbital of the BDT molecule. Ryndyk et al. [19] performed DFT+ME calculations including strong electron correlations and got a multi-scale Coulomb blockade. We will also address the impact of strong correlations. Our results obtained for the single level model for the BDT molecule will be compared with the results for the model studied by Ryndyk et al. [19].

The paper is organized as follows. Section 2 introduces briefly the methodology and is divided in two parts. The first part describes how the parameters in the extended Hubbard Hamiltonian are obtained from first-principles, while the second part presents a technique based on NEGF and CPT for solving this model. In section 3 the method is applied to a gold-BDT-gold and a platinum-BDT-platinum system. The conclusions are given in section 4.

2 Methods

2.1 First-principles calculation of the model parameters

The system we are considering consists of a central region (C) attached to left and right leads, see figure 1. To determine the tight-binding (TB) parameters of the system we employ the pseudopotential plane-wave code Quantum Espresso [33]. The unit cell, red dashed line in figure 1, consists of the central region, left and right transition layers ($T_{L}$, $T_{R}$), also called buffer regions, and the surface part of the remaining leads ($L_{L}$,$L_{R}$). The size of the transition layers is chosen such that the electronic properties of the outermost atoms in the unit cell do not change anymore when the unit cell is increased. This ensures that the central region has negligible impact on the periodicity of the remaining leads [34, 35].

For all DFT calculations in this paper we have used the Perdew Zunger (PZ) exchange correlation functional within the local density approximation (LDA) and non-relativistic ultrasoft pseudopotentials from the Standard Solid-State Pseudopotentials (SSSP) library (PSlibrary 0.3.1) [37].

We have used Maximally Localized Wannier Functions (MLWF) [38, 39] to get localized orbitals and to set up the TB Hamiltonian [34, 35]. All Wannier transformations are done with wannier90 [40] which produces orthogonal MLWF. Therefore the derivation of the NEGF formalism in section 2.2 is demonstrated only for orthogonal bases. A derivation for non-orthogonal bases in the coherent regime can be found in [41]. In the localized basis the Hamilton operator of the unit cell is then given by

[TABLE]

The outermost atoms of the unit cell are periodically continued and form the non-interacting leads ($L_{L}$, $L_{R}$) in the TB model. In this way we determine all single-particle parameters of the Hamiltonian.

We describe the electron-electron interaction by an extended Hubbard Hamiltonian in the central region, which has the form

[TABLE]

where ${\mathop{\hat{H}}}_{0}$ is the single-particle part and ${\mathop{\hat{H}}}_{1}$ the interaction part of the model Hamiltonian. We neglect correlation effects in the leads ($L_{L}$,$T_{L}$,$L_{R}$,$T_{R}$) and we expect correlation effects to be important only in the central region, because on the molecule the wave functions are more localized. The single-particle part of the model Hamiltonian in the central region in second quantization is given by

[TABLE]

The indices $i=\{\bar{i},\sigma\}$ and $j=\{\bar{j},\sigma^{\prime}\}$ include the wave function indices $\bar{i}$ and $\bar{j}$ and the spins $\sigma$ and $\sigma^{\prime}$. The operator ${{\mathop{\hat{a}_{i}^{\dagger}}}}$ creates an electron on orbital $\Psi_{i}({\bf r})$ which we obtain by the MLWF construction. The hopping parameters $\tilde{t}_{ij}$ are

[TABLE]

where $Q$ is the matrix which diagonalizes the TB Hamiltonian and $E$ are the corresponding eigenenergies. The interaction part of the model Hamiltonian is given by

[TABLE]

The density-density term form of the Hubbard model is only justified if the non-density-density term is small which also requires that the orbitals are maximally localized. The interaction parameters $U_{ij}$ are calculated via

[TABLE]

with the screened Coulomb potential

[TABLE]

where $\eta$ is the screening. According to [19] we set the screening to be constant. The Random Phase Approximation (RPA) [42] would be a possibility to determine the screening consistently.

The electron correlations already included in DFT, have to be subtracted in order to avoid double counting. Here we choose the around mean-field approximation (AMF) [43], because it generally gives good results for weakly correlated systems. That is, we have subtracted the Hartree term $\Delta t_{\bar{i}\bar{i}}$ of the interaction part from the bare hopping $\tilde{t}_{ij}$. This approach is justified as follows. If we apply the LDA approximation to the model Hamiltonian we end up with the same band structure as that is obtained from the original band structure calculation. The double counting correction leads to modified hopping parameters $\tilde{t}_{ij}\to t^{C}_{ij}$ with

[TABLE]

The Hartree term is given by

[TABLE]

where $n_{\bar{i}}^{0}$ is the equilibrium population of the $i$-th orbital obtained in the DFT calculation. The model Hamiltonian for the central region has finally the form

[TABLE]

The full self-consistent non-equilibrium transport calculation could be done by feeding the non-equilibrium potential, obtained from the non-equilibrium densities in the transport calculation, back into DFT and repeat this procedure until self-consistency is reached [44, 45, 46].

2.2 Non-equilibrium Green’s functions

We have used the Keldysh formalism [47] to compute the current through the central region. In steady-state the two-time Green’s function depends only on the time difference and can be treated efficiently in frequency space. The steady-state current from the left transition layer to the central region can be calculated by integrating over the real part of the Keldysh Green’s function

[TABLE]

Here, $i$ and $j$ denote the orbitals within the clusters $T_{L}$ and $C$, see figure 1, for calculating the current between the left transition layer and the central region. In this paper we employ non-equilibrium cluster perturbation theory (n-CPT) [48, 9] to compute the Keldysh Green’s function $\mathcal{G}$ of the interacting non-equilibrium system required in (13). Within n-CPT, (13) reduces to the well-known Landauer-Büttiker formula [49].

[TABLE]

$T(\omega)$ denotes the transmission function

[TABLE]

where ${\bm{\Gamma}}_{\nu}$ with $\nu\in\{L,R\}$ is the anti-hermitian part of the hybridization function

[TABLE]

We denote matrices with bold faced letters and for simplicity of the notation, we have omitted the $\omega$ arguments. To this end, the entire system is decomposed into the device ($C$), e.g. the BDT molecule, and the two leads, which are further decomposed into $L_{L}$ plus $T_{L}$ on the left side and $L_{R}$ plus $T_{R}$ on the right. We calculate the full Green’s function for indices of the central region ${\bm{\mathcal{G}}}_{CC}$ within the CPT approximation, i.e. we use the following Dyson equation

[TABLE]

In CPT the Green’s function ${\bf g}^{r,a}_{CC}$ of the isolated central cluster is required. From the Dyson equation

[TABLE]

it is clear that in the interacting case, the CPT approximation is equivalent to approximating the self-energy by that of the isolated cluster. Therefore it is important to cut the system in a way that the coupling between the transition layers and the central region is weak compared to the hopping parameters within the central region. The Green’s function of the isolated central cluster can be calculated using the Lehmann representation

[TABLE]

If the system is not too far away from the equilibrium situation, it will still be a good approximation to consider only excitations from the ground state to higher states and hence the corresponding self-energy.

The Green’s function of the isolated left lead with indices restricted to the transition region ${\bf G}_{T_{L}T_{L}}$ is needed in (17). We will also need ${\bf G}_{T_{R}T_{R}}$. To calculate ${\bf G}_{T_{\nu}T_{\nu}}$ for $\nu\in\{L,R\}$ we once more introduce a cluster decomposition of the leads into part $L_{\nu}$ and $T_{\nu}$, introduced before. The exact Dyson equation for this decomposition yields

[TABLE]

where ${\bf g}_{T_{L}T_{L}}$ stands for the Green’s function of the isolated transition region, which can be computed easily, as the matrix dimension is small. ${\bf g}_{L_{\nu}L_{\nu}}$ on the other hand represents the Green’s function of the isolated semi-infinite part $L_{\nu}$ of lead $\nu$ that is composed of identical unit cells. The hopping elements ${\bf V}_{T_{\nu},L_{\nu}}$ couple only to the unit cell near to the transition region in each lead and therefore only the surface part of ${\bf g}_{L_{\nu}L_{\nu}}$ is required, which can be determined iteratively by using the method of Sancho et al. [50]. An alternative technique for calculating the hybridization function is to use complex absorbing potentials [51].

A bias voltage $V_{b}$ drives the system out of equilibrium and enters the calculation as a shift of the on-site energies and the chemical potential by $-V_{b}/2$ for the left lead and by $+V_{b}/2$ for the right lead. In the Landauer-Büttiker formula $V_{b}$ enters also in the Fermi functions. A gate voltage $V_{g}$ can be applied by shifting the on-site energies of $\tilde{t}_{\bar{i}{\bar{j}}}$.

The spectral density and the conductivity are important quantities for our analyses and they are defined as

[TABLE]

3 Applications

3.1 The Au-BDT-Au system

In a first step we have applied the DFT+NEGF formalism to a system consisting of a BDT molecule sandwiched between two gold chains, schematically drawn in figure 1. As shown below in the case of monoatomic gold chains it is possible to include all lead atoms in the Wannier transformation in a controlled way. In mechanically controllable break junctions (MCBJ) or with an scanning tunnelling microscope (STM) the formation of monoatomic gold chains was demonstrated under elongation of the junction by several authors [52, 44, 53, 54, 55]. In a first step we have performed DFT calculations 111We used a cut-off energy for the wavefunction of $E_{\rm cowf}=64~{}$Ry and the charge density of $E_{con}=512~{}$Ry, the Methfessel-Paxton smearing technique with the smearing parameter $E_{\rm tsm}=0.01~{}$Ry, the vacuum distance $d_{\rm vac}=12~{}$Å  and $32$ k-points. for the periodic chain of monoatomic gold atoms. All the parameters are converged with respect to the total energy within $5~{}$meV. In order to find the equilibrium geometry, we have calculated the total energy of the chain for different distances $d_{Au}$ between the gold atoms. We find an optimal distance of $d_{Au}=2.515~{}$Å  which is close to the LDA result $2.51~{}$Å  in [56].

The band structure of the Au chain in an interval around $E_{F}$ is shown in the left panel of figure 2. The $\Gamma$-$X$ axis is always the axis along the gold chain in reciprocal space. The 6p bands are at higher energies. The 6s and 5d bands coincide exactly with the Wannier bands, which demonstrates the properness of the Wannier transformation. The dotted lines represent the inner and the outer energy window needed to disentangle the 6s and 5d bands from the 6p band in the Wannier transformation. In figure 2 tight binding (TB1) means that only hopping processes to nearest neighbour gold atoms are taken into account. The 6s band and the 5d bands separate into the two bands corresponding to the $d_{xz}$ and the $d_{yz}$ orbital (Channel 1) and four bands corresponding to orbitals with inversion symmetry with respect to the $xy$-plane (Channel 2). The definition of the coordinate system is according to figure 5. Each symmetry channel will only couple to orbitals of the same symmetry and therefore transport takes place in two separated channels. The TB1 approximation turns out to be a good approximation for the channel 1 bands, which are mainly interesting for transport. The right panel of figure 2 is the calculated equilibrium ($V_{b}=0$) transmission function for TB1 and is similar to the results obtained in [56]. In the coherent regime the transmission function is proportional to the number of bands at a certain energy.

After having calculated the Hamiltonian of the monoatomic Au chain, we considered the full system, including BDT 222We used a cut-off energy for the wavefunction of $E_{\rm cowf}=100~{}$Ry and the charge density of $E_{con}=400~{}$Ry, the Methfessel-Paxton smearing technique with the smearing parameter $E_{\rm sm}=0.001~{}$Ry, the vacuum distance $d_{\rm vac}=18~{}$Å  and $1$ k-point.. The BDT molecule can assume different geometry configurations between the two semi-infinite gold chains. In the line configuration the BDT molecule is forced to be on the gold chain axis. In the atop configuration the BDT molecule slightly hops out of the gold chain axis, while it twists out of the axis in the twisted configuration as shown in figure 1. There are also mixed configurations between atop and twisted, called atop-twisted and twisted-twisted. In order to determine the energetically favourable configuration, we compute the total energy as function of the distance between the gold atoms near to the BDT molecule $d_{Au-Au}$, namely the gap between the two (left and right) semi-infinite gold chains, for different configurations, see figure 1 and 3. The distance between the gold atoms was fixed $d_{Au}$ during the geometry optimizations.

At $d_{Au-Au}\approx 8.1$ the BDT molecule rearranges into thienothiophene. It occurs that the twisted configuration of the BDT atom is energetically most favoured and, consequently, we have chosen this configuration for the ensuing transport calculations. The optimized distance is $d_{Au-Au}=9.16~{}$Å. We have chosen 6 gold atoms per unit cell, which ensures that the eigenvalues of the atoms at the outer edges of the unit cell (see figure 1) are converged. We have chosen the size such that the changes in eigenenergies are less than $0.05~{}$eV.

The left panel in figure 4 depicts a comparison of the band structure obtained with the full plane wave basis and with the Wannier basis. The Wannier transformation is performed only at the $\Gamma$-point. In the Wannier transformation we picked out 47 bands, 35 corresponding to the 6 gold atoms and there are 12 BDT bands within this energy range. At the $\Gamma$-point the PWSCF bands coincide exactly with the Wannier levels, which demonstrates again the properness of the Wannier transformation. The obtained Wannier basis functions of the gold atoms coincide with the Wannier basis functions of the pure gold chain, which shows that the disentanglement procedure of the 6s and 5d bands from the 6p bands is done correct. Projections of the wavefunctions onto orthogonalized atomic wavefunctions show whether the Wannier levels belong to the gold leads and the transition layers or to the BDT molecule. The right panel presents the eigenlevels of the decoupled BDT molecule, the eigenenergies of ${\mathop{\hat{H}}}_{C}$. As mentioned above and outlined in [6, 29] transport takes place in separated channels. Channel 1 (blue dashed) consists of $p_{z}$-like BDT orbitals and channel 2 (cyan) of $p_{xy}$-like BDT orbitals. Each channel couples only to the corresponding channel in the gold leads.

The eigenorbitals of channel 1 are shown in figure 5. The shape of the orbitals can be explained by the symmetry of the benzene molecule, which has the point group $D_{6h}$. LUMO+2 is the $B_{2g}$-orbital, LUMO+1 and LUMO are the $E_{2u}$-orbitals, HOMO-4 is the $A_{2u}$-orbital and the others are a mixture of the benzene $E_{1g}$- and the sulfur $p_{z}$-orbitals.

Due to the fact that only the HOMO level is close to the Fermi energy it is the level mainly contributing to transport which is in general the case in systems with thiol-anchoring groups [58]. The HOMO orbital is an unsaturated sulfur $p$-orbital and can also be compared to [59, 31]. Passivating the unsaturated sulfur $p$-orbitals would shift the HOMO level away from the Fermi energy and produce lower conductances [16]. Charge transport through channel 1 is blocked due to the fact that there are no empty states with $p_{z}$-like symmetry in the monoatomic gold chain (only a $sd$-hybridized orbital belonging to channel 2 is above the Fermi level in figure 2). The levels corresponding to channel 2 are more than $2~{}$eV below the Fermi energy. Therefore, transport in this channel is due to the broadening of the leads and the current-voltage characteristic nearly linear. The conductivity for the BDT molecule connected with monoatomic Au chains is almost constant and in the order of $0.01~{}G_{0}$ for biasvoltages bellow $4~{}$V which is in the order of the experimental values of $0.01~{}G_{0}$ [23, 24, 25, 26]. But using gold tips or bulk-like gold as leads [29, 6] produces empty states with $p_{z}$-like symmetry and activates channel 1 which raises the conductance. Therefore the theoretical LDA result for the Au-BDT-Au system with bulk-like leads is $0.28~{}G_{0}$ [16].

Since channel 1 contains the HOMO-orbital and the typical benzene $p_{z}$ orbitals it would be the most important channel if supported by the leads. The coupling between leads and the channels depends on the geometry of the junction and therefore as mentioned above gold tips or bulk-like gold as leads [29, 6] activate channel 1. Also the use of platinum instead of gold opens channel 1 because platinum has one electron less than gold and therefore also orbitals with $p_{z}$-like symmetry contribute to transport.

3.2 The Pt-BDT-Pt system

In order to activate channel 1 we have replaced the gold atoms by platinum ones in our calculations 333We used a cut-off energy for the wavefunction of $E_{cowf}=60~{}$Ry and the charge density of $E_{con}=200~{}$Ry, the Methfessel-Paxton smearing technique with the smearing parameter $E_{sm}=0.005~{}$Ry, the vacuum distance $d_{vac}=12~{}$Å  and $60$ k-points.. The optimized distance between the platinum atoms $d_{Pt}=2.328~{}$Å  is comparable to the value of $2.34~{}$Å  in [60]. The plane-wave pseudopotential band structure, Wannier bands and TB1 bands along with the transmission function are shown in figure 6. In the Pt chain not just channel 2 but also channel 1 bands contribute to the equilibrium transmission at the Fermi energy.

We now consider the Pt-BDT-Pt system 444We used a cut-off energy for the wavefunction of $E_{cowf}=80~{}$Ry and the charge density of $E_{con}=320~{}$Ry, the Methfessel-Paxton smearing technique with the smearing parameter $E_{sm}=0.005~{}$Ry, the vacuum distance $d_{vac}=15~{}$Å  and $1$ k-point..

The optimized distance between the Pt atoms near to the BDT molecule is $d_{Pt-Pt}=8.58~{}$Å. It is necessary to take at least 8 platinum atoms per unit cell to ensure that the influence of the BDT molecule on the lead platinum atoms is negligible. The comparison of the PW band structure with the Wannier bands is given in the left panel of figure 7. In the Wannier transformation we retained 62 bands, 47 corresponding to the 8 platinum atoms and there are 15 BDT bands within this energy range. The right panel presents the decoupled levels of the BDT molecule and is comparable to figure 4 disregarding 3 more levels in channel 2 needed to get all the platinum levels in the calculation. In the following we will only discuss channel 1 because due to the level alignment the current in channel 2 is an order of magnitude smaller. The orbitals of channel 1 look similar to the ones in the Au-BDT-Au calculation in figure 5. The orbitals in figure 8 are the localized Wannier orbitals of the BDT molecule in the Pt-BDT-Pt calculation. These are the 8 orbitals entering the TB model for the central region. They all have $p_{z}$ symmetry.

The matrix form of the TB Hamiltonian of the central region, in the basis of the localized orbitals in figure 8, and the coupling matrices to the leads that have been used in the calculations are given in section A. The Fermi energy in the transport Hamiltonian is adjusted such that the local DOS of channel 1 of the atom further away from the BDT molecule has the same filling as in the pure atom chain.

The plot in the middle of figure 9 shows the spectral density $A_{CC}(\omega)$ of the isolated central cluster, namely the isolated BDT molecule. The curves labelled by DFT+CPT represent the result for the correlated system, which will be explained in section 3.4. The peaks in the DFT+NEGF result are at $-0.4~{}$eV, $-2.0~{}$eV and $-2.3~{}$eV. As discussed above the HOMO level closer to the Fermi level ($-0.4~{}$eV) is mainly responsible for the transport.

The spectral density at the surface atom of the left and the right lead is shown in the left and the right panel of figure 9. The spectral density of a homogeneous multiorbital chain is a superposition of semi-circles, see $A_{LL}(\omega)$ and $A_{RR}(\omega)$. Coupling the homogeneous chains to the the transition layers changes the semi-circular structure of the spectral density. The spectral densities $A_{TLi}(\omega)$ and $A_{TRi}(\omega)$ are located $i$ atoms away from the last point of the transition region, i.e. points with larger $i$ are closer to the central regions. We nicely see, how the spectral function gradually changes from the multi-circular structure of the homogeneous chain to the structure of the spectral densities $A_{TL3}(\omega)$ and $A_{TR3}(\omega)$, which consists of peaks at $-1.6$ and $-0.5$ and a spike at about $-1.3$. The grey shaded area indicates the filling. The remarkable point is that the spectral function of the transition region, that enters the transmission function and therefore the current, has a complex structure that will even change when a voltage is applied. As demonstrated by Cuniberti et al. [61, 4] low dimensional leads can affect the conductance due to the finite band width and the structure in the DOS. If the DOS is nearly constant one can use the wide-band limit (WBL) approximation which works well for systems with bulk metal electrodes [62]. However the DOS of a one-dimensional chain has Van Hove singularities near the band edge and the WBL approximation does not work. The key message therefore is that the details of the leads in such a molecular device are of crucial importance for the transport properties. This adds to the observation that the symmetry of the molecular orbitals and the leads can lead to selection rules as far as transport channels are concerned.

3.3 Stability diagrams of the Pt-BDT-Pt system

Next we present various aspects of the charge stability of the Pt-BDT-Pt system. The stability diagram computed with the DFT+NEGF formalism is presented in figure 10. We have restricted the bias voltage to $-3~{}$V$<V_{b}<3~{}$V, because otherwise we would have to include platinum $p$-orbitals, and likewise we have restricted the gate voltage to the interval $-2~{}$V$<V_{g}<2~{}$V, otherwise channel 2 becomes important. All diagrams are calculated at an inverse lead temperature of $\beta=300~{}1/$eV, which enters in the Fermi functions of the leads. The conductance at $V_{b}=V_{g}=0$ is $0.52~{}G_{0}$ and is considerable higher than the experimental ($0.01~{}G_{0}$ [23, 24, 25, 26]) and theoretical LDA ($0.28~{}G_{0}$ [16]) values of Au-BDT-Au systems. Beside the influence of the different contact material the formation of monoatomic chains results in large increases in the conductance [31, 52, 63]. The distance $d_{Pt-Pt}$ is also small compared to experimental values, where the monoatomic chains arise under applying stress to the leads. This results in a higher coupling strength between the leads and the central region.

Most strikingly, the stability diagram in figure 10 does not show the usual structure of crossing straight lines, resulting in rhombic patterns. The only such structure is the pair of straight lines starting at the left edge of the diagram with $V_{g}=0.8$ and $V_{g}=1.3$, respectively. They correspond to the level at $-2.3~{}$eV. At $V_{b}=3$ the Fermi energy of the right leads is shifted to $\mu=-1.5$. Above the Fermi level the remaining DOS of unoccupied states has a width of $\Delta E_{empty}=0.5$ eV, see figure 9. Hence, the transport window of the right lead is $(-1.5,-1.0)$ eV. With $V_{g}=0.8$ the level, which was originally at $-2.3$ eV is shifted to the lower edge of the transport window and with $V_{g}=1.3$ it is shifted to the upper edge. This constant distance between the lines is an effect of the band edge according to which a positive conductance is always followed by a negative one. In our case the distance between the lines along the $V_{b}$-axis is $\Delta V_{b}=2\Delta E_{empty}=1~{}$eV. The factor two results from the fact that half of the bias voltage is applied to each lead. The slope of the lines is constant and equal to $\pm 1/2$. If the leads are shifted with $\pm V_{b}/2$, one needs $2V_{b}$ to compensate a shift in $V_{g}$. As a consequence of weak coupling of this level to the leads the effects can be seen so clearly.

The levels at $-2.0~{}$eV and $-0.4~{}$eV produce a similar structure but at different energies. The structure can be explained by a ”band edge effect” and a ”supporting effect”. To explain the ”supporting effect” we have projected out the level at $-0.4~{}$eV. The resulting charge stability diagram is depicted in figure 11. Comparison to the full calculation illustrates that in a broad range of $V_{g}$ it is enough to use a single-level model for describing the transport through the BDT molecule and suggests approaches based on a single-level model [64].

The structure in the stability diagram is a consequence of the strong coupling to the leads. In the present case the leads within channel 1 consist again of two nearly decoupled channels, the ”supporting channel” and the ”conducting channel”. The real (imaginary) part $R^{L}$ ($\Gamma^{L}$) of the retarded left lead hybridization function $\Delta^{L,r}$ for these two channels are shown in figure 12.

Only the ”conducting channel” has states above the Fermi level. The real part $R^{L}$ shifts the levels in the central region. $R^{L}$ of the ”supporting channel” has anti-resonances at $\omega_{1}=-0.5~{}$eV and $\omega_{2}=-1.6~{}$eV. The spike at $-1.3~{}$eV is due to the small coupling between the channels in the leads and its influence is negligible. The imaginary part $\Gamma^{L}$ broadens all the levels and has maxima at the anti-resonances in $R^{L}$. By the ”supporting effect” we mean that at every $V_{b}$, where the positive branch of an antiresonance of $R_{L}$ hits the transport window of the opposite lead, the energy level in the central region will be shifted up to higher energies and therefore one needs smaller $V_{g}$ to compensate this and the conductivity structure in the stability diagram bends down, see figure 11. The blue dashed lines show the case where the positive branches of the anti-resonances of $R_{L}$ hit the upper edge of the opposite lead,

[TABLE]

and act as limit where the effect is maximal. The same effect with opposite sign happens if a negative branch of an antiresonance of $R_{L}$ hits the transport window of the opposite lead, there the conductivity structure in the stability diagram bends up. The red dashed line in figure 11 shows the case where the negative branch of the anti-resonances of $R_{L}$ hits the lower edge of the transport window, the Fermi energy, of the opposite lead,

[TABLE]

There the structure in the stability diagram is maximally bent up. This effect explains the appearance of more than one maximum in the current voltage characteristic even if there is only a single level in the central region.

Between $-\Delta E_{empty}<V_{b}<\Delta E_{empty}$ is the usual structure of a stability diagram, but smeared out due to the coupling of the leads. Disabling the supporting channel produces the stability diagram in figure 13 where the structure can be explained just by the ”band edge effect”.

3.4 Many-body effects in the Pt-BDT-Pt system

After having discussed the properties of the single-particle part of the Hamiltonian, we now include many-body effects using the DFT+CPT method. The interaction parameters $U_{ij}$ entering the model Hamiltonian in (12) are determined by numerical integration of (8) and listed in section A. We take a constant screening parameter $\eta=1.5$ as proposed by Ryndyk et al. [19]. The values obtained in our calculation for the nearest neighbour hopping between the carbon atoms within benzene in (27) and the on-site Hubbard interaction in (29) are comparable to those reported in [65].

In table 1 the lowest many-body eigenenergies of the isolated central region are listed. The ground state ($10_{g}$) is in the $10$-particle sector and is a singlet state. The other low-lying states in the $10$-particle sector have singlet or triplet symmetry. The lowest eigenenergies in the $9$-particle sector are doublet states.

Including interaction in principle shifts levels and spectral weight away from the Fermi energy and produces additional peaks and therefore signatures in the current-voltage characteristic at higher voltages. Therefore the benzene band gap increases compared to the LDA calculation [66] and can cause a reduction of the current. As suggested by [3, 32] the gap also widens in calculations based on the self-interaction correction (SIC) in DFT. In our transport system adding the full interaction the HOMO level is shifted down to $-0.6~{}$eV, see DFT+CPT result in figure 9. The peaks obtained by this calculation can be identified with excitations obtained by a true many-body calculation for the effective Hubbard model. There is a peak at $-0.6~{}$eV corresponding to the excitation from the $10_{g}$ to the $9_{g}$ state. The next peak is at $-2.1~{}$eV and corresponds to the excitation from the $10_{g}$ to the $9^{\prime}$ state. This excitation needs more energy than we have taken into account in our Wannier basis and are therefore negligible.

The stability diagram obtained by including electron-electron correlations within CPT does not change qualitatively for $V_{g}<0.6~{}$eV, see figure 14. The drastic change at $V_{g}=0.6~{}$eV is an artefact of CPT. Equation (19) is solved once for each $V_{g}$ independent of the applied $V_{b}$. At $V_{g}=0.6~{}$eV the HOMO level gets depleted and therefore the particle sector changes. We expect that using more sophisticated non-equilibrium approaches like the ME+CPT calculation [20] the drastic change at $V_{g}=0.6~{}$eV should disappear. The conductance at $V_{b}=V_{g}=0$ of $0.86~{}G_{0}$ is higher than the result obtained with LDA even though the HOMO is shifted down in energy and has less spectral weight. In the DFT+CPT calculation the spectral weight on the HOMO level is lower but distributed in a way that there is more weight in the orbitals near to the leads which causes the higher conductance. Thygesen et al. [16] have obtained the same trend at the BDT molecule connected with gold tips and studied with LDA and GW.

In contrast to Ryndyk et al. [19] the lead coupling effects overcome the influence of strong electron correlations in our transport system because of stronger coupling between leads and BDT due to the geometry of the transport system.

4 Conclusion

We have performed first-principle calculations based on DFT, MLWFs, NEGF and CPT to study the molecular system consisting of benzenedithiolate connected to semi-infinite monoatomic metal electrodes. DFT, within the plane-wave pseudopotential method, is used to calculate the electronic band structure of the transport system. Transforming the Kohn-Sham eigenvalues and eigenfunctions to a real-space basis of Maximally Localized Wannier Functions allows to extract a tight-binding Hamiltonian to model the transport system. Non-equilibrium Green’s functions are used in turn to calculate the charge transport through the BDT molecule. In the case of gold electrodes the HOMO level, which provides the dominant contribution to transport properties does not contribute to transport due to symmetry reasons and therefore the conductance is small. Platinum electrodes, on the other hand, enable transport via the HOMO level. Strong electron correlations are included on the BDT molecule using an extended Hubbard Hamiltonian. Since the system is in the intermediate coupling regime the spectral properties of the leads overcome the influence of many-body effects in this system. We find that ”band edge effects” and ”supporting effects” are more relevant for the structure of the stability diagram in this regime.

5 Acknowledgements

We would like to thank M. Aichhorn for fruitful discussions. This work was partially supported by the Austrian Science Fund (FWF) within the SFB-ViCoM (F4103 and F4115) and P26508, and NAWI Graz. The DFT calculations were partly done on high performance computing resources of the ZID at TU Graz.

Appendix A Parameters of the model Hamiltonian

The matrices in (26), (27) and (28) are the single-particle parameters of the Pt-BDT-Pt system. $\tilde{t}$ are the hopping parameters of the of the central region (C) and ${V}_{T_{L}C}$ and ${V}_{CT_{R}}$ the coupling matrices between the central region and the transition layers. The row and column indices of $\tilde{t}$ correspond to the basis functions in figure 8.

[TABLE]

[TABLE]

[TABLE]

The interaction parameters $U_{ij}$ entering the model Hamiltonian determined by numerical integration of (8) are

[TABLE]

The relative integration error is within $5~{}$%.