Electrode-molecule energy level offsets in a gold-benzene diamine-gold single molecule tunnel junction.

One means for describing electron transport across single molecule tunnel junctions (MTJs) is to use density functional theory (DFT) in conjunction with a nonequilibrium Green's function formalism. This description relies on interpreting solutions to the Kohn-Sham (KS) equations used to solve the DFT problem as quasiparticle (QP) states. Many practical DFT implementations suffer from electron self-interaction errors and an inability to treat charge image potentials for molecules near metal surfaces. For MTJs, the overall effect of these errors is typically manifested as an overestimation of electronic currents. Correcting KS energies for self-interaction and image potential errors results in MTJ current-voltage characteristics in close agreement with measured currents. An alternative transport approach foregoes a QP picture and solves for a many-electron wavefunction on the MTJ subject to open system boundary conditions. It is demonstrated that this many-electron method provides similar results to the corrected QP picture for electronic current. The analysis of these two distinct approaches is related through corrections to a junction's electronic structure beyond the KS energies for the case of a benzene diamine molecule bonded between two gold electrodes.


I. INTRODUCTION
Determining the conductance for single molecule tunnel junctions (MTJs) has been a subject of experimental and theoretical investigation for over two decades [1,2]. To reconcile theoretical predictions with experimental observations for the conductance across a MTJ is often challenging due to the difficulties in predicting or measuring energy level spacings and offsets for occupied and unoccupied states in a MTJ, as well as the difficulties in determining the geometries for molecules bonded between two metallic electrodes. Combining theoretical and experimental efforts has nonetheless led to greater understanding for charge transport mechanisms in MTJs.
Theoretical calculations for electronic currents often require empirical corrections or highly accurate electronic structure computations, or in instances both approaches are required in tandem. From a simplified viewpoint, a molecule in a MTJ may be thought of as a tunnelling barrier with the magnitude of the potential barrier given by the energy difference between the quasi-Fermi levels in the electrodes and the highest or lowest occupied molecular orbital (HOMO/LUMO) energy levels. At low applied voltages and assuming non-resonant tunnelling, a potential barrier model provides a reasonable estimate of tunnelling currents.
In this picture, the barrier height and width are the key physical quantities for predicting charge transport. Of course, the barrier model as described is an oversimplification of the true electronic structure of a MTJ. More accurate electronic structure descriptions include effects due to charge transfer and geometry variations accompanying the bonding of a molecule between electrodes. To describe these effects, the treatment of hybridization between molecular and electrode electronic states is required.
To provide a description of the energy levels within a MTJ to predict transport properties, a common strategy is to combine the non-equilibrium Green's functions (NEGF) [3][4][5] formalism together with density functional theory (DFT) [6], an approach sometimes designated NEGF+DFT [7]. Significant understanding has been achieved in relating measured current-voltage (IV) characteristics to the electronic structure of MTJs within this framework. There remain however shortcomings in obtaining quantitative agreement to reconcile experiment within this approach that typically lead to a systematic over-estimation of electronic current. This deficiency is primarily attributed to two causes. The first of these is the well-documented band gap underestimation if using Kohn-Sham orbital energies as quasiparticle (QP) energies obtained from typical approximations to the DFT exchangecorrelation (XC) functional [8]. The 'band gap problem' results from spurious electron self-interactions errors [9] that are not explicitly cancelled or compensated for within many XC approximations. In terms of a simple barrier model, the band gap problem manifests as an underestimation of the potential barrier and hence an over-estimation of tunnelling transmission. A second difficulty facing many approximate DFT approaches in describing a MTJ is a poor description of image charge potentials for molecules absorbed at or bonded to metal surfaces [10]. An approach to overcome these shortcomings is to add self-energy corrections to the Kohn-Sham energy levels that compensate electron self-interactions and then to add corrections that also describe the image charge potential. This and related approaches have been shown to yield results in quantitative agreement with experiments [11][12][13][14][15][16]. Corrections for the electron self-interaction typically rely on computationally demanding GW calculations. Suggestions for alternative approaches using less demanding computations of the self-interaction correction and classical electrostatics have been demonstrated to improve the description of the image potential relative to uncorrected NEGF+DFT [10,11,15,16].
The methods described up to this point retain a QP description of the junction electronic structure. Alternative approaches focusing on calculation of the electronic current from a direct determination of a density matrix (DM) forego a QP picture [17][18][19][20][21][22][23]. These methods apply different methods but they share in common that non-equilibrium boundary conditions appropriate to electronic transport [24] are imposed directly to the calculation of DMs.
In contrast to NEGF+DFT methods based on the Kohn-Sham equations, direct calculation of a DM or determining a one electron reduced DM (1e-RDM) from a many-electron wavefunction does not need to invoke the use of QP energies. Hence relating electronic currents calculated from a determination of a DM or a RDM to the potential barrier model is not immediate.
Herein many-electron correlated scattering (MECS) [17][18][19] is applied which relies on a many-electron wavefunction method, configuration interaction (CI), to calculate a 1e-RDM subject to open system boundary conditions. Previously, MECS has been applied to single MTJs consisting of benzene dithiol [17,18], alkane dithiol [25], alkane diamine [26], silane dithiol/diamine [27], and point contacts [28]. These calculations demonstrate that MECS predicts electronic currents consistent with measured molecular conductance. The agreement with experiment suggests that deficiencies associated with the NEGF+DFT approach do not arise or are less severe when using a many-electron wavefunction treatment of a MTJ's electronic structure. To explore if improved molecular conductance predicted from MECS calculations can be be directly related to the improvements resulting from correcting the QP energy spectrum by self-interaction and image potential correction is the subject of this study. The experimentally well-characterized BDA MTJ has been chosen as the test case to compare between the two distinct approaches to electronic charge transport.

II. METHODS
Previously, Quek et al. showed with correction terms to QP energies that theoretical predictions can be made in agreement with low voltage conductance statistics measured from a series of traces for a gold break junction in a BDA molecular solution [11]. A gold electrode-BDA-gold electrode model based upon DFT calculations is also described by the authors. Their NEGF+DFT model is used as a reference point for calculation of IV characteristics which they subsequently used correction terms to improve. Based on various tip geometries built from semi-infinite electrodes with a bonding tip on an otherwise planar Au (111) surface or from pyramidal point contact models, it is found for molecules with amine linkers bonding to Au electrodes that IV characteristics are not sensitive to the exact bonding geometry of the linker to the electrode 'tip' [11,26]. To ensure that a common starting point between the previous calculations and this work, a key consideration in addition to the junction's explicit atomic structure is that similar electronic structure treatments for the uncorrected NEGF+DFT calculations are used in this study to enable a comparison to the previous calculations reported in ref. [11]. Hydrogen. The junction axis is taken to run between the two gold electrode tips in the diagram.
The DFT calculations undertaken by Quek et al. are obtained using localized atomic orbitals and periodic boundary conditions for treating the electrodes. Their basis makes use of numerical split valence polarized atomic orbitals (PAO), details of which are available in their supplementary material [11]. As the goal in this work is to apply a CI calculation to determine the 1e-RDM, it is advantageous to use a DFT program that does not apply periodic boundary conditions motivating the choice of the Turbomole quantum chemistry programs [29,30] for the DFT calculations. Turbomole applies Gaussian atomic basis sets and a double zeta basis including polarization is selected for the molecular atoms (C, N, H). However, care is taken in selecting basis sets to produce results comparable to those of Quek et al. This is accomplished by comparing the Kohn-Sham orbital energies from DFT calculations on isolated atoms and for the molecular orbitals of BDA. To describe Au atoms for the electrode model, a basis with 11 explicit electrons is applied by Quek et al. In contrast to their approach, CI calculations are required in the following as a first step to obtain corrections beyond the NEGF+DFT model. Explicit electrodes must also be included in the CI calculations and this is done using pyramidal electrodes with a tip geometry similar to that used in ref. [11]. Similar electrode structures have been studied in detail for the calculation of the conductance quantum in point contacts using the MECS approach [28]. For MECS calculations, a CI expansion is used to determine a many-electron wavefunction. The length of a CI expansion increases combinatorially with the number of electrons and single electron basis functions. Treating electrodes with a large number of electrons or basis functions becomes impractical due to the resulting rapid increase in the length of the CI expansion. To limit the number of explicit electrons, the electrode Au atoms are split into two groups. The four atoms closest to the molecule at each electrode (one Au apex atom and three Au atoms in the second layer of the electrode tip) are treated with an 11-electron basis set. For the remaining Au atoms, a basis with only one explicit electron and four basis functions, representing each Au atom's 6s orbital, with the addition of three diffuse p-type functions, are used. The limitations of the finite cluster approximation and extending the electrodes is considered in detail in ref. [28]. Using the same bond lengths and bond angles for the electrode tip and linker molecules as reported by Quek et al., the gold-BDA-gold MTJ is constructed and taken as a starting geometry. The structured is then allowed to relax (minimization of energy with respect to atomic positions) using the Perdew-Burke-Ernzerhof (PBE) XC functional which implements a generalize gradient approximation (GGA) for the charge density. The atomic positions are relaxed using the same force threshold as used in. ref. [11] The final junction geometry is shown in fig. 1.
For the NEGF+DFT calculations, electrode self-energies are computed from repeat units consisting of the two back planes of the explicit electrode. It is found that six layers of Au atoms (3 repeat units) form an electrode cell that is sufficiently large to make nextnearest-neighbour interactions between cells negligible. This model of the junction and semi-infinite electrodes are then used with the TiMeS quantum transport program [31] to calculate electron transmission. The TiMeS program uses semi-analytical expressions for the Green's functions in the NEGF formalism for which the limit for the artificial broadening parameter approaching zero has been taken. This is in contrast to most transport programs that rely on numerical determination of the Green's function resulting in a broadening and smoothing of resonance peaks. To determine electron transmission within the NEGF+DFT approach, the relevant Fock matrices, interaction matrices, and orbital overlaps are extracted from the Kohn-Sham calculations to build a device Hamiltonian and electrode self-energies in a localized atomic orbital basis to allow partitioning into a device region, left electrode, and right electrode.
After establishing a QP basis that can reasonably reproduce the junction transmission at the NEGF+DFT level, a CI calculation is then used to improve the electronic structure description of the junction. The Kohn-Sham orbitals are then used to build spin-projected Slater determinants or configuration state functions (CSFs) which are used to form the CI expansion [32][33][34]. To provide an accurate picture of the electronic current in the MECS formalism, it is necessary to include several excited states to allow coupling to the ground state as a voltage is applied. The voltage generates an electric field between the electrodes that couples excited states to the ground state resulting in polarization of the junction. The coupling can be explained by the fact that the energy change caused by an applied electric field is proportional to dipole interaction matrix elements of the applied electric field. To lowest order, the coupling to ground state is Here, E is the magnitude of the electric field, the summation runs over all excited manybody states, and the direction of the field is taken along the transport axis selected to be the x-direction. The electric field operator is anti-symmetric with respect to inversion symmetry in the junction. Excited states with the same symmetry with respect to inversion as the ground state will not couple in the summation for the polarizability energy. For the case of the junction studied in this work, which has C 2h symmetry, there are four irreducible representations (irreps) denoted a 1 , a 2 , b 1 , and b 2 . Application of the electric field breaks symmetry along the transport direction, a 1 states therefore couple to b 2 states, and likewise a 2 states will couple to b 1 states. The many-electron ground state has A 1 symmetry. As the junction polarizes, the current-carrying state to lowest order can be described as a superposition of non-polarized A 1 and B 2 states. Excited states of opposing parity under inversion to the ground state increasingly couple with increasing electric field strength, and therefore become increasingly important in the CI expansion as a voltage is applied across the MTJ.
The previous observations suggest a procedure for pre-selecting important CSFs to describe the junction as a voltage is applied. To accurately capture the polarization of the junction and likewise the current induced by the open system boundary conditions, a set of CSFs to be included in the MECS calculation can be found from the lowest excited CSFs that couple to the unperturbed (zero voltage) wavefunction through the dipole matrix elements. Following this reasoning, the many-electron transport basis is performed using a [35,36]. To generate the excited CSFs of opposing parity, an electron from the highest occupied a 1 molecular orbital (MO) is transferred to the lowest unoccupied b 2 MOs.

Monte Carlo Configuration Interaction (MCCI) calculation with a CI coefficient threshold
The three lowest energy B 2 CSFs provide approximations to the first three B 2 excited states.
Additionally, singly excited CSFs of A 1 symmetry are included in the CI expansion as Kohn-Sham orbitals are used so that Brillouin's theorem does not exclude the single excitations.
The role of the A 1 single excitations may be understood from Thouless' theorem [37] which relates two non-orthogonal Slater determinants |Ψ occupied with N electrons selected from M spin orbitals by which can for small coefficients c ij be approximated as where a † i , a j are creation and annihilation operators, respectively, |Ψ , |Ψ represents Slater determinants, and N is the normalization. Hence for weak electric fields across the MTJ, the effect of including the A 1 singles allows a 're-optimization' of the the reference state, as well the A 1 single excitations couple to the excited B 2 states as the electric field is introduced.
This procedure results in a CI vector including only 1890 CSFs. The resulting CI expansion provides a compact description of electron-electron correlations and polarizability of the MTJ.
The MECS method relies on the use of the Wigner phase-space distribution function [38,39] to constrain the 1e-RDM to open boundary conditions [17][18][19]. Constraints are applied to ensure the momentum flowing from the electrodes into the molecule remain unchanged as a voltage is applied -consistent with the condition that the electrodes are locally in equilibrium. To the contrary, electrons entering the left and right electrodes from the molecular region do not fulfil detailed balance as voltage difference is applied across the junction resulting in a net electronic current flow. It is assumed the electrons scattered into the electrodes rapidly re-equilibrate to ensure the electrodes remain locally in equilibrium.
To apply the boundary conditions for the electrodes, the Wigner function f (q, p) is calculated from the 1e-RDM obtained from the CI wavefunction allowing for a semiclassical picture of the momentum p at a given position q to be described by the quasi-probability distribution function f . Due the uncertainty principle, the Wigner phase-space distribution is not a true probability distribution function as it is not everywhere positive nor is it a unique phase space representation for quantum mechanics [24,38,39]. Nonetheless, in many cases the Wigner function is known to provide a semi-classical picture for electron distributions and this has been previously demonstrated for semiconductor devices and for the metal electrodes used in break junction experiments [28]. For MECS calculations, a plane normal to the current flow in each electrode is chosen and the Wigner function is averaged over this plane to give a one-dimensional representation of the momentum flowing across the plane. At zero voltage bias, the distributions in each electrode are symmetric indicating that the net momentum flow in and out of the MTJ balances to zero at each electrode, there is no current. As an electric field is applied across the molecular region, the electrode equilibrium condition requires the inward flow to remain the same as at zero voltage. However the potential difference across the MTJ gives rise to asymmetric scattering between left and right incoming electrons. MECS calculations apply a constrained energy minimization procedure to fix the momentum flow from the electrodes to their equilibrium (zero applied voltage) value as an electric field is applied. The momentum flow out of the molecular region into either electrode is unconstrained. The method of Lagrangian multipliers is applied to constrain the Wigner function, and additional Lagrangian multipliers are used to force the current magnitude to be constant between points at which current constraints are applied.
The latter constraints are introduced to mitigate the violation of current continuity due to the variational nature of the calculation and finite basis effects.
In order for the assumption that the incoming electron momentum distributions at the Wigner planes remain in equilibrium with the fictional electron reservoirs, the Wigner planes need to be positioned within the finite electrode model at a position that minimizes edge effects [28]. Consideration of the placement of the Wigner planes focuses on the region in the electrodes between 7.5 and 9. the most effect on electronic currents [28].
The CI variational problem is typically solved through the associated matrix eigenvalue problem. However in the MECS formulation additional constraints as those arising from imposition of the scattering boundary conditions [24] and current continuity as are required.
The associated eigenvalue becomes non-linear due to the constraint conditions and standard methods for solving the CI problem are not applicable. Treating the variational problem as a constrained minimization problem allows for the determination of a many-body wavefunction satisfying open system boundary conditions for the current-carrying state induced by the applied electric field E . The Lagrange multiplier penalty method function is applied to the minimization problem by defining the augmented Lagrangian L by where z is chosen as the transport axis. The first term is the unconstrained many-electron energy for the junction with application of an external electric field represented as a CI expansion, the second term is the sum of Wigner function and current continuity constraint violations with their respective Lagrangian multipliers λ i , and the third term is a penalty term introduced to aid in the numerical convergence and stability of the minimization procedure. Further details can be found in refs. [17][18][19].
The steps involved in the MECS approach may be summarized as: The choice of the Wigner function to impose open system boundary conditions is a convenient choice [24]. However, the MECS method is not restricted to use of the Wigner function, any distribution that defines the momentum flowing into and out of each electrode can be used to constrain the 1e-RDM. Also, the integration of the Wigner momentum across a plane in the electrodes is also not strictly required, constraints at several points across the plane in the electrodes can be used. However, for computational efficiency it is found that constraining the total momentum crossing the plane for the Wigner momenta p is sufficient to treat the electrodes as electron reservoirs in equilibrium.
The current-voltage measurements reported in ref. [11]   However in their work resonances from DFT are lower in energy relative to the Fermi level, possibly due to differences in basis sets, XC potential, and the details of the bonding geometry to the tip. The broader peaks found in the calculation of ref. [11] are due in part to their planar electrode structure relative to the pyramidal electrode model used in this work. The planar electrode model is a better approximation to the wide band limit expected to well approximate gold electrodes, whereas in this work the electrode is a nanowire structure leading to the gold tips as described in sect. II leading to a less bulk-like density of states. Also, as previously mentioned, the TiMeS program uses a semi-analytic version of the Green's functions and no artificial broadening is introduced. However, the overall hybridization of the molecular states to the tips of the gold electrodes is seen to be comparable between the two models. Additional differences likely arise from the specific bonding geometries and the atomic orbital basis set differences between the two calculations. However, the behavior and descriptions are in overall strong agreement. These NEGF+DFT results are taken as the zero th order reference between subsequent calculations that include corrections beyond the two DFT models used to generate fig. 2.
The MECS based transport calculations are next considered using the CI expansion as described in sect. II. The MECS transport calculations enables a determination of the molecular conductance which is then compared to the experimental and theoretical estimates found in ref. [11] to evaluate electronic structure corrections beyond the NEGF+DFT approximation. The electronic current density in the MECS approach is calculated directly from the 1e-RDM as where ρ( r, r ) is the 1e-RDM , e is the electron charge, is Planck's constant divided by 2π, and m 0 is the electron mass. The current density is integrated over planes normal to the transport direction. The MECS method relies on minimizing the many-electron energy over the explicit MTJ subject to the application of the open-system boundary conditions [18,24], as well as the inclusion of current constraints. The energy is an integrated quantity and the minimization procedure does not ensure that the RDM calculated from a finite CI expansions is uniformly accurate at each position. Hence the calculation of the current can vary from plane-to-plane. This motivates adding local current continuity conditions to the constrained minimization procedure. These additional constraints are imposed by requiring the current calculated on two subsequent planes normal to the transport direction to be within a specified tolerance. Fig. 3 shows The current density calculated from the 1e-RDM matrix can be shown in a finite basis set approximation to have a maximum least squares error, and that the maximum error decreases with increasing basis set size when using orbitals that diagonalize the 1e-RDM or natural orbitals. However, the error can be lower than this maximum and convergence need not be monotonic as the expansion is increased. As well, the error varies point-to-point [40].
A second source of error is that the molecular orbitals used to construct the CI wavefunction are calculated from a finite cluster model and are not open system (current-carrying) orbitals.
The ability of the CSFs constructed from these orbitals to carry a current relies on the use of complex CI coefficients within a MECS calculation. This leads to an additional error leading to point-point variations when constructing the 1e-RDM from these orbitals. This can be understood in analogy to approximating a momentum eigenstate from particle-in-abox (PiB) wavefuntions. For a rectangular potential well of length L, the (unnormalized) eigenfunctions defined on [−L/2, +L/2] are given by To approximate a momentum eigenstate, the approximation ψ(x) = ψ n (x) + iψ n+1 (x) is used. Assume n is odd, the current density in this approximation is proportional to which is at larger values of n is approximately equal to the current from a momentum eigenstate k n = nπ/L multiplied by a broad envelope function that is maximum at L = 0.
The second term represents a low amplitude, rapid oscillation in current with an envelope function that is a minimum at L = 0. Note the difference of a factor of 2 in the wavevector  fig. 4c with seventeen constraints, the larger number of constraints results in a slight reduction in violation of current continuity but also produces a significant reduction in the current through the junction.
The three examples in fig. 4 highlight the effects of the current constraints. In fig. 4a, the violation of current continuity is large but with the average current providing a reasonable approximation to measured currents. Fig. 4b provides a much improved current profile over the explicit junction region and the current magnitude remains in good agreement for a wide range of simulation parameter seems to indicate that the energy levels as corrected in ref. [11] are in some sense reproduced by the CI treatment of the molecular junction.
Indeed, if it is assumed the molecule gives rise to approximately a rectangular potential barrier, the width of which is determined by the geometry of the junction, then the fact that similar currents are obtained in both approaches implies the effective potential barriers must be of comparable energies.
To determine to what level these expectations are met, a comparison of the junction energy level alignments between the two methods is considered.
The level alignment between the HOMO of the BDA molecule with respect to the metal's level to the experimental ionization energy (IP) [41]. The difference in energy between the experimental and Kohn-Sham IPs is then used to compensate for the unphysical electron self-interaction. The energy difference is determined to be roughly -3.0 eV and this value is also taken to estimate the correction required when the molecule is bonded between two metallic electrodes. The second correction term applied is to include the effects of the image charge potential as described in for example refs. [10,11,16]. This correction term is  transfer processes, and to continue with the comparison between MECS and the corrected NEGF+DFT, it is necessary to consider how to define similar energies from the manyelectron wavefunction treatment of the MTJ. Pickup and Gościński show that ionization potentials and electron affinities can be calculated perturbatively from a many-body wave function using the one-body Green's function [42]. Their analysis relies on Hartree-Fock orbitals in which the first order corrections vanish but which can be generalized for the use of Kohn-Sham orbitals and allows for a direct comparison to the QP picture. Here an alternative approach relating the CI energies and wavefunctions directly to the QP energies is used to define the potential barrier height.
To identify the molecular states in the junction that are hybridized with the electrode states, a procedure analogous to a Mulliken population analysis is performed. The total charge for a normalized molecular orbital is unity and can be calculated for the i th MO φ i expanded in a basis of atomic orbitals ξ j as eV, corresponding to the onset of the transport resonance in the NEGF+DFT calculations.
Within this energy range, two MO levels display appreciable charge density in the molecular region indicating hybridization to the metal electrode states.
To determine the alignment of these molecular frontier states within the MTJ to the metal electrode work function after the introduction of electron correlations, the vertical IPs are calculated for the two molecular states with the highest charge densities in the energy range of interest. To determine the IP estimates requires two CI calculations to be performed, the first is for the neutral N electron model of the MTJ and the second is for a vertical ionization process from the molecular hybridized states for N − 1 electron models of the MTJ. The ionization from the molecular level is approximated from the equilibrium or N electron ground state CI expansion by operating on it by the destruction operator a i removing the i th spin orbital to create a hole state. The CI energies are calculated for the N and N − 1 many-electron states with the resulting IP estimate given as where the subscript denotes the single electron state for which a hole has been introduced.
The difference in energies between the ionization potential and the electrode work function then provides an estimate of the energy level alignment between the molecular hybridized states and the electrode work function. As before, the energy difference is then taken to approximate the potential barrier height in the junction. Thus by comparing the energy offsets predicted from the many-electron wavefunction calculations and the corrected QP treatments in relation to the the magnitude of the calculated MTJ conductances, the comparison between the two methods is achieved.    [13] in which an explicit GW calculation for the gold-BDA-gold junction has been performed. From fig. 3 of ref. [13], the energy difference between the molecular HOMO and the gold Fermi energy is on the order of 2 eV, and the barrier height from their DFT calculation is on the order of 1 eV. Hence the GW calculation does not predict as large of an increase in the barrier height relative to the DFT estimate compared to this work and ref. [11]. However the results do agree in that a more accurate electronic structure description of the tunnel junction results in an increase barrier height leading to a reduction in the conductance. It should also be noted that the barrier heights obtained from the three methods used to account for the effects of electron correlation beyond approximate DFT produce conductance values that are within the range of experimental error. Hence the ability of the QP corrected transport levels and the MECS calculation to produce similar current magnitudes can be related to improved estimates of the potential barrier height in the simple tunnel model of a MTJ.
Notably, it appears that a modest CI expansion can significantly improved the electronic structure description of the junction beyond a DFT/GGA calculation, and this description does not suffer from electron self-interaction errors due to its absence from the many-electron and providing self-interaction and image potential corrections can reproduce experimentally measured conductance values. In this study, an approach that explicitly includes correlations on the molecular junction through a many-electron wavefunction approach produces comparable results to the corrected K-S QP approach as well as the reported experimental measurements. The analysis reveals that a simple tunnel barrier model for the MTJ is treated comparably between two approaches for improving the electronic structure description beyond the reference NEGF+DFT calculation.
To perform the comparison, the MECS method was studied with respect to applying constraint conditions to the 1e-RDM obtained from the CI wavefunction for correlated electron transport. In particular, the sensitivity of the calculations to application of open system boundary conditions and limitations to current conservation were explored. Their effects based on reasonable physical considerations are found to be relatively moderate relative to experimental distributions for the BDA MTJ conductance. It should be noted that although the CI expansion used is defined with the use of Kohn-Sham orbitals, that the (non-relativistic) molecular Coulomb Hamiltonian defines the electronic structure for the MTJ for the MECS calculation. Effects from using approximations for the XC potential in DFT are absent, and hence improving a single determinant approximation by including a moderate level of electron correlation and polarizabilty can compare quite well to systematic corrections to the Kohn-Sham energy levels. Although conceptually and formally quite different, this work demonstrates that in the low voltage bias regime that a single particle picture of electron transport and a many-electron scattering approach both provide predictions comparable to well defined experiments for the conductance of a MTJ consisting of a benzene diamine molecule bonded between two gold electrodes. The calculations reinforce the importance of predicting electronegativity in a MTJ to accurately treat charge transport [43]. In this regard, the ability to optimize ionization potentials and electron affinities in correlated orbital theories [44] suggests an attractive path forward if combined with NEGF methods to treat charge transport in molecular junctions.