# Tutorial

## Theory of Electron, Phonon and Spin Transport in Nanoscale Quantum Devices

At the level of fundamental science, it was recently demonstrated that molecular wires can mediate long-range phase-coherent tunnelling with remarkably low attenuation over a few nanometre even at room temperature. Furthermore, a large mean free path has been observed in graphene and other graphene-like two-dimensional materials. These create the possibility of using quantum and phonon interference to engineer electron and phonon transport through nanoscale junctions for wide range of applications such as molecular switches, sensors, piezoelectricity, thermoelectricity and thermal management. To understand transport properties of such devices, it is crucial to calculate electron transmission coefficient (Te) and phonon transmission coefficient (Tph) through them. The aim of this tutorial article is to outline the basic theoretical concepts and review the state-of-the-art theoretical and mathematical techniques to treat electron, phonon and spin transport in nanoscale molecular junctions. This helps not only to explain new phenomenon observed experimentally but also provides a vital design tool to develop novel nanoscale quantum devices.

Index Terms: Molecular electronics, Nanoelectronics, Theory and modelling, Quantum interference, Phonon interference

#### I. Introduction: Molecular electronics

The silicon-based semiconductor industry is facing a grave problem because of the performance limits imposed on semiconductor devices when miniaturized to nanoscale [1]. To counter this problem, new nanoscale materials such as carbon nanotubes, graphene and transition metal dichalcogenide monolayers have been proposed [2]. Alternatively, the idea of using single molecules as building blocks to design and fabricate molecular electronic components has been around for more than 40 years [3], but only recently has it attracted huge scientific interest to explore their unique properties and opportunities. Molecular electronics including self-assembled monolayers [4] and single-molecule junctions [5] are of interest not only for their potential to deliver logic gates [6], sensors [7], and memories [8] with ultra-low power requirements and sub-10-nm device footprints, but also for their ability to probe room-temperature quantum properties at the molecular scale such as quantum interference [9] and thermoelectricity [10, 11]. There are five main area of research in molecular-scale electronics [5] namely: molecular mechanics, molecular optoelectronics, molecular electronics, molecular spintronics and molecular thermoelectrics.

By studying electron and phonon transport across a junction consisting of two or more electrodes connected to a single or few hundred molecules, one could study the electrical and mechanical properties of nanoscale junctions such as molecular electronic building blocks, sensors, molecular spintronic, thermoelectric, piezoelectric and optoelectronic devices. For example, when a single molecule is attached to metallic electrodes, de Broglie waves of electrons entering the molecule from one electrode and leaving from the other form a complex interference pattern inside the molecule. These patterns could be utilized to optimize the single-molecule device performance [6]. In addition, the potential of such junctions for removing heat from nanoelectronic devices (thermal management) and thermoelectrically converting waste heat into electricity has recently been recognized [10]. Indeed, electrons passing through single molecules have been demonstrated to remain phase coherent, even at room temperature. In this tutorial, the aim is to outline the basic theoretical concepts and review theoretical and mathematical techniques to model electron, phonon and spin transport in nanoscale molecular junctions. This helps not only to understand the experimental observations but also provides a vital design tool to develop strategies for molecular electronic building blocks, thermoelectric devices and sensors.

Transport on the molecular scale: Any nanoscale device consists of two or more electrodes (leads) connected to a scattering region (figure 1 ). Electrodes are perfect wave-guides where electrons and phonons can transmit without any scattering. The main scattering occurs either at the interface to leads or inside scattering region. The goal is to understand the electrical and vibrational properties of nano and molecular junctions where a nanoscale scatter or a molecule is bonded to electrodes with strong or weak coupling in the absence or presence of surroundings, such as an electric field (e.g. gate and bias voltages or local charge), a magnetic field, a laser beam or a molecular environment (e.g. water, gases, biological spices, donors and acceptors). There are different approaches to study the electronic and vibrational properties of molecular junctions such as semi-classical methods [1214], kinetic theory of quantum transport [15], scattering theory [16] and master equation approach [17]. In this paper, our focus is mostly on the scattering theory based on the Green’s function formalism and the master equation approach.

Fig. 1.  A scattering region is connected to reservoirs trough ballistic leads. Reservoirs have slightly different electrochemical potentials to drive electrons from the left to right leads. All inelastic relaxation process take place in the reservoirs and transport in the leads are ballistic.

Here, we begin with the Schrödinger equation and relate it to the physical description of a matter at the nano and molecular scale. Then we will discuss the definition of the current using the time-dependent Schrödinger equation and introduce density functional theory (DFT) and a tight binding description of quantum system. The scattering theory and non-equilibrium Green’s function method are discussed and different transport regimes (on and off resonances) are considered. One dimensional system and a more general multi-channel method are derived to calculate transmission coefficient T(E) for electrons (phonons) with energy E (ω) traversing from one electrode to the other through a scattering region. We then briefly discuss the master equation method to model transport in the Coulomb and Franck-Condon blockade regimes. We follow with a discussion of the physical interpretation of quantum systems including charge, spin and thermal currents, phonon thermal conductance, electron-phonon interaction, piezoelectric response, inclusion of a Gauge field and superconducting systems. Furthermore, environmental effects and different techniques used to model the experiment are discussed.

#### II. Schrödinger equation

The most general Schrödinger equation [18] describes the evolution of physical properties of a system in time and was proposed by the Austrian physicist Erwin Schrödinger in 1926 as:

 (1)

where i = , is the reduced Planck constant (h∕2π), Ψ is the wave function of quantum system, r and t are the position vector and time respectively, and Ĥ is the Hamiltonian operator which characterizes the total energy of any given wave function. For a single particle moving in an electric field, the non-relativistic Schrödinger equation reads as:

 (2)

where m is the particle’s reduced mass, V is its potential energy and 2 is the Laplacian. If we assume that the Hamiltonian is time-independent and write the wavefunction as a product of spatial and temporal terms: Ψ(r,t) = ψ(r)θ(t), the Schrödinger equation become two ordinary differential equations:

 (3)

and

 (4)

where Ĥ = 2 + V (r). Note that this is not a solution if the Hamiltonian is time-dependent e.g. when a laser is shined to a system or an AC high frequency voltage is applied in which V (r) is varied by time and a time-dependent DFT should be considered [19]. The solution of equation 3 could be written as: θ(t) = e-iEt∕. The amplitude of θ(t) does not change with time and therefore the solutions θ(t) are purely oscillatory. The total wave function

 (5)

differs from ψ(r) only by a phase factor of constant magnitude and the expectation value |Ψ(r,t)|2 is time-independent. Of course equation 5 is a particular solution of the time-dependent Schrödinger equation. The most general solution is a linear combination of these particular solutions as:

 (6)

In time independent problems, the spatial part needs to be solved only, since the time dependent phase factor in equation 5 is always the same. Equation 4 is called the time-independent Schrödinger equation and it is an eigenvalue problem where E’s are eigenvalues of the Hamiltonian Ĥ. Since the Hamiltonian is a Hermitian operator, the eigenvalues E are real. ψ(r) describes the standing wave solutions of the time-dependent equation, which are states with definite energies called ”stationary states” or ”energy eigenstates” in physics and ”atomic orbitals” or ”molecular orbitals” in chemistry.

The Schrödinger equation must be solved subject to appropriate boundary conditions. Since electrons are fermions, the solution must satisfy the Pauli exclusion principle and wavefunctions ψ must be well behaved everywhere. The Schrödinger equation can be solved analytically for a few small systems e.g the hydrogen atom. However, the Schrödinger equation is too complex in most cases to be solved even with the best supercomputers available today. Therefore, some approximations are needed [20] such as the Born-Oppenheimer approximation to decouple the movement of electrons and nuclei; density functional theory (DFT) to describe electron - electron interactions and pseudopotentials to treat nuclei and core electrons except those in the valence band. We will briefly discuss these approximations in the next section. A more detailed discussion can be found elsewhere [20]. To describe transport through molecules or nanoscale matters, one needs to build a simple tight-binding Hamiltonian using Hückel parameters or use DFT to construct the material specific mean-field Hamiltonian.

To reduce the size of Hamiltonian, it is appropriate to define basis functions where

 (7)

The wavefunction can then be represented by a column vector |ϕconsisting of the expansion coefficients ϕi. The time-independent Schrödinger equation could be written as a matrix equation:

 (8)

where

 (9)

and

 (10)

The evaluation of these integrals is the most time-consuming step, but once [H] and [S] are obtained, the eigenvalues En and eigenvectors ϕn are easily calculated. If i| and |jare orthogonal then Sij = δij where δij is the Kronecker delta (δij = 1 if i = j and δij = 0 if ij). Note that a system with the Hamiltonian H and overlap matrix S obtained using non-orthogonal basis could be transformed to a new Hamiltonian H = S-12 × H × S-12 with orthogonal basis (S = I where I is the Identity matrix.)

##### A. Density functional theory (DFT)

In order to understand the behavior of molecular electronic devices, it is necessary to possess a reliable source of structural and electronic information. A solution to the many body problem has been sought by many generations of physicists. The task is to find the eigenvalues and eigenstates of the full Hamiltonian operator of a system consisting of nuclei and electrons as shown in figure 2 . Since this is not practically possible for the systems bigger than a few particles, some approximations are needed. The atomic masses are roughly three orders of magnitudes bigger than the electron mass, hence the Born-Oppenheimer approximation [20] can be employed to decouple the electronic wave function and the motion of the nuclei. In other words, we solve the Schrödinger equation for the electronic degrees of freedom only. Once we know the electronic structure of a system, we can calculate classical forces on the nuclei and minimize these forces to find the ground-state geometry (figure 2 a).

Once the Schrödinger equation was solved, the wavefunction is known and all physical quantities of intereste could be calculated. Although the Born-Oppenheimer approximation decouple the electronic wave function and the motion of the nuclei, the electronic part of the problem has reduced to many interacting particles problem which even for modest system sizes i.e. a couple of atoms, its diagonalization is practically impossible even on a modern supercomputer. The virtue of density functional theory DFT [20, 21] is that it expresses the physical quantities in terms of the ground-state density and by obtaining the ground-state density, one can in principle calculate the ground-state energy. However, the exact form of the functional is not known. The kinetic term and internal energies of the interacting particles cannot generally be expressed as functionals of the density. The solution is introduced by Kohn and Sham in 1965. According to Kohn and Sham, the original Hamiltonian of the many body interacting system can be replaced by an effective Hamiltonian of non-interacting particles in an effective external potential, which has the same ground-state density as the original system as illustrated in figure 2 a. The difference between the energy of the non-interacting and interacting system is referred to the exchange correlation functional (figure 2 a).

Fig. 2. From many-body problem to density functional theory DFT. (a) Born-Oppenheimer approximation, Hohenberg-Kohn theorem and Kohn-Sham ansatz, (b) Schematic of the DFT self-consistency process.

Exchange and correlation energy: There are numerous proposed forms for the exchange and correlation energy V xc in the literature [20, 21]. The first successful - and yet simple - form was the Local Density Approximation (LDA) [21], which depends only on the density and is therefore a local functional. Then the next step was the Generalized Gradient Approximation (GGA) [21], including the derivative of the density. It also contains information about the neighborhood and therefore is semi-local. LDA and GGA are the two most commonly used approximations to the exchange and correlation energies in density functional theory. There are also several other functionals, which go beyond LDA and GGA. Some of these functionals are tailored to fit specific needs of basis sets used in solving the Kohn-Sham equations and a large category are the so called hybrid functionals (eg. B3LYP [22], HSE [23] and Meta hybrid GGA [22]), which include exact exchange terms from Hartree-Fock. One of the latest and most universal functionals, the Van der Waals density functional (vdW-DF [24]), contains non-local terms and has proven to be very accurate in systems where dispersion forces are important.

Pseudopotentials: Despite all simplifications shown in figure 2 , in typical systems of molecules which contain many atoms, the calculation is still very large and has the potential to be computationally expensive. In order to reduce the number of electrons, one can introduce pseudopotentials which effectively remove the core electrons from an atom. The electrons in an atom can be split into two types: core and valence, where core electrons lie within filled atomic shells and the valence electrons lie in partially filled shells. Together with the fact that core electrons are spatially localized about the nucleus, only valence electron states overlap when atoms are brought together so that in most systems only valence electrons contribute to the formation of molecular orbitals. This allows the core electrons to be removed and replaced by a pseudopotential such that the valence electrons still feel the same screened nucleon charge as if the core electrons were still present. This reduces the number of electrons in a system dramatically and in turn reduces the time and memory required to calculate properties of molecules that contain a large number of electrons. Another benefit of pseudopotentials is that they are smooth, leading to greater numerical stability.

Basis Sets: In order to turn the partial differential equations (e.g. the Schrödinger equation 1 ) into algebraic equations suitable for efficient implementation on a computer, a set of functions (called basis functions) is used to represent the electronic wave function. For a periodic system, the plane-wave basis set is natural since it is, by itself, periodic. However, to construct a tight-binding Hamiltonian, we need to use localised basis sets discussed in the next section, which are not implicitly periodic. An example is a Linear Combination of Atomic Orbital (LCAO) basis set which are constrained to be zero after some defined cut-off radius, and are constructed from the orbitals of the atoms.

Mean-field Hamiltonian from DFT: To obtain the ground state mean-field Hamiltonian of a system from DFT, the calculation is started by constructing the initial atomic configuration of the system. Depending on the applied DFT implementation, the appropriate pseudopotentials for each element which can be different for every exchange-correlation functional might be needed. Furthermore, a suitable choice of the basis set has to be made for each element present in the calculation. The larger the basis set, the more accurate our calculation - and, of course, the longer it will take. With a couple of test calculations we can optimize the accuracy and computational cost. Other input parameters are also needed that set the accuracy of the calculation such as the fineness and density of the k-grid points used to evaluate the integral([21, 25]). Then an initial charge density assuming no interaction between atoms is calculated. Since the pseudopotentials are known, this step is simple and the total charge density will be the sum of the atomic densities.

The self-consistent calculation [21] (figure 2 b) starts by calculating the Hartree and exchange correlation potentials. Since the density is represented in real space, the Hartree potential is obtained by solving the Poisson equation with the multi-grid or fast Fourier-transform method. Then the Kohn-Sham equations are solved and a new density is obtained. This self-consistent iterations end when the necessary convergence criteria are reached such as density matrix tolerance. Once the initial electronic structure of a system was obtained, the forces on the nucleis is calculated and a new atomic configuration to minimize these forces obtained. New atomic configuration is the new initial coordinate for the next self-consistent calculation. This structural optimization is controlled by the conjugate gradient method for finding the minimal ground state energy and the corresponding atomic configuration [21]. From the obtained ground state geometry of the system, the ground state electronic properties of the system such as the total energy, the binding energies between different part of the system, the density of states, the local density of states, the forces could be calculated. It is apparent that the DFT could potentially provide an accurate description of the ground state properties of a system such as the total energy, the binding energy and the geometrical structures. However, all electronic properties related to excited states are less accurate within DFT.

##### B. Tight-Binding Model

By expanding wave function over a finite set of atomic orbitals, Hamiltonian of a system can be written in a tight-binding model. The main idea is to represent the wave function of a particle as a linear combination of some known localized states. A typical choice is to consider a linear combination of atomic orbitals (LCAO). If the LCAO basis is used within DFT, the Hamiltonian H and overlap S matrices used within the scattering calculation (section III ) could be directly extracted. However, if a plane-wave DFT code is used, a LCAO-like based Hamiltonian could be constructed using Wannier functions. For a periodic system where the wave function is described by a Bloch function, equation 8 could be written as

 (11)

where c and care the neighboring identical cells containing states α and

 (12)

and

 (13)

Equation 11 could be written as

 (14)

where

 (15)

and

 (16)

More generally, the single-particle tight-binding Hamiltonian in the Hilbert space formed by |Rαcould be written as:

 (17)

where εα is the corresponding on-site energy of the state |α, V α is the electrical potential and γαβ are the hopping matrix elements between states |αand |β.

Simple TB Hamiltonian: For conjugated hydrocarbons, the energy of molecular orbitals associated with the π electrons could be determined by a very simple LCAO molecular orbital method called Hückel molecular orbital method (HMO). Therefore, a simple TB description of system could be constructed by assigning a Hückel parameter to on-site energy εα of each atom in the molecule connected to the nearest neighbours with a single Hückel parameter (hopping matrix element) γαβ. Obviously, more complex TB models could be made using HMO by taking the second, third, forth or more nearest neighbours hopping matrix elements into account.

It is worth mentioning that once a material specific LCAO mean-field DFT or a simple HMO Hamiltonians (described in this section) were obtained, electron and spin transport properties of a junction can be calculated.

1) Two level system

As a simplest example, consider a close system of two single-orbital sites with on-site energies ε and -ε coupled to each other by the hoping integral γ. The Hamiltonian of such system is written as: H = , so the Schrödinger equation reads:

 (18)

The eigenvalues E are calculated by solving det(H - EI) = 0 where I = is the identity matrix.

 (19)

E- and E+ are called the bonding and anti-bonding states. There must be two orthogonal eigenvectors and corresponding to each eigenvalue. By substituting equation 19 into equation 18 ,

 (20)

If ε = 0 and E = ±γ, simplest normalised eigenstates could be written as:

 (21)

If γ = 0 and E = ±ε, the wave functions are fully localised on each site:

 (22)

Effective Hamiltonian: The solution of the two level system can be used to obtain solutions of a larger system. For a given N ×N Hamiltonian H, a new effective energy dependent (N - 1) × (N - 1) Hamiltonian Ĥ(E) could be obtained by decimating a given site p,

 (23)

Therefore, any N × N Hamiltonian H of an arbitrary system with N site could be reduced to 2 × 2 energy dependent effective Hamiltonian Ĥ(E) by decimating all sites but two sites of interest using equation 23 .

2) One dimensional (1D) infinite chain

Consider an infinite linear chain of hydrogen atoms shown in figure 3 a. A single orthogonal orbital nearest neighbor tight binding Hamiltonian of such system with on-site energies j|H|j= ε0 and the hopping matrix elements j|H|j ± 1= j ± 1|H|j= -γ could be written as:

 (24)

 (25)

where -∞ < j < +. The solution of this equation could be obtained using the Bloch function as

 (26)

and

 (27)

where k = ka0 is the dimensionless wave vector and -π∕a0 < k < π∕a0 in the first Brillouin zone. Equation 27 is called a dispersion relation (E - k) or electronic band-structure of a 1D chain. Since -1 < cos(k) < 1, hence ε0 - 2γ < E < ε0 + 2γ; therefore the bandwidth is 4γ. In the bottom of the band around Emin (figure 3 c), by Taylor expansion of equation 27 , a parabolic band structure is obtained: E(k) ε0 - 2γ + γk2. Comparing this with a free electron parabolic band structure E(k) = ε0 - 2γ + 22m*k2, γ = 22m* is inversely proportional to the effective mass m*. This implies that electrons in a system with smaller (larger) bandwidth are heavier (lighter).

Fig. 3. One dimensional (1D) infinite chain. (a) hydrogen atoms in an infinite chain with one orbital per atom, (b) 1D balls and springs, (c,d) electronic and phononic band structures and (e,f) density of states (DOS) of a and b.

The band structure of a perfect 1D chain (equation 27 ) is a continuous function and does not have any energy band gap. Therefore, electrons injected to a 1D chain transmit perfectly within the bandwidth (metal). However, 1D metallic chain does not exist in reality. The reason is, if a band gap was opened (e.g. by mechanical force), all energy levels below valance band move down and therefore, total energy of the system decreases leading to a more stable structure (Peierls distortion). Unless gained electric energy is higher than lost mechanical energy needed to distort the system, this will continue. Therefore, atoms in a 1D crystal oscillate, so that the perfect order is broken.

Density of States: From the dispersion relation (e.g. equation 27 ), one can obtain the density of states (DOS) using:

 (28)

where λk are the eigenvalues of the system, δ is Kronecker delta and k kx,ky,kz where kp -∞+dkp2π. From the dispersion relation (equation 27 ), DOS of a 1D chain is obtained as: D(E) = dk∕πdE = 12πγsin(k). At k = ±π where E = ε± 2γ (close to the band edges), D →∞. This singularity in the DOS called Van Hove singularity. Figure 3 a shows the band structure and the DOS of a infinite 1D chain.

1D balls and springs: We have yet discussed the electronic properties of a quantum system e.g. 1D chain. Now consider a chain of atoms with mass m connected to the next nearest neighbours with the springs with spring-constant K = -γ as shown in figure 3 . On the one hand, the derivative of the energy U with respect to the position x describes forces on a particular atom (F = -U). Provided a local minimum is reached, one can expand the potential energy about this minimum in terms of the powers of atomic displacements. Since all forces on all atoms are zero, the Taylor expansion does not have linear terms. By neglecting higher orders (harmonic approximation), F = -x = -Kx. Note that being at a local minimum, the matrix of second derivatives must be positive definite, and thus will have only positive eigenvalues. On the other hand, from Newton’s second law F = -m. Therefore, the Schrödinger-like equation could be written as:

 (29)

Similar to what was discussed above, using xn(t) = Aei(kn-ωt), equation 29 reads -2 = -K[2 - e-ik - eik] and therefore the phononic dispersion relation of a 1D chain of ball and springs (figure 3 b) is obtained as

 (30)

Comparing equations 27 and 30 , it is apparent that equation 30 is obtained by changing E ω2 and ε0 2γ∕m in equation 27 . ε0 = 2γ∕m is the negative of sum of all off-diagonal terms of 1D chain TB Hamiltonian which make sense to satisfy translational invariance. Therefore, Schrödinger equation-like relation for phonons could be written as

 (31)

where D = -K∕M is the dynamical matrix, M is the mass matrix and K is Hessian matrix calculated from the force matrix.

3) One dimensional (1D) finite chain and ring

To analyse the effect of boundary conditions in the solution of the Schrödinger equation, consider three examples shown in figure 4 . As the first example, consider a 1D finite chain of N atoms. As a consequence of introducing boundary conditions at the two ends of the chain, the energy levels and states are no longer continuous in the range of ε0 - 2γ < E < ε0 + 2γ; instead there are discrete energy levels with corresponding states in this range. This is obtained by writing the Schrödinger equation in 1 < j < N (equation 25 ) and at the boundaries j = 1 and j = N. At j = 1 the Schrödinger equation reads

 (32)

and at j = N

 (33)

using the Bloch function ϕj = eikj + ce-ikj, the solution for the 1D finite chain problem is obtained as:

 (34)

where n [1,...,N]. Similarly, solutions of 1D finite ring of N atom could be obtained (figure 4 ) as:

 (35)

where n [0,...,N - 1]. Clearly, the allowed energy levels of 1D finite chain is different than 1D ring. This demonstrates that a small change in a molecular system may significantly affect the energy levels and corresponding orbitals. This becomes more important when few number of atoms was investigated e.g. the molecules, so two very similar molecules may show different electronic properties.

Figure 4 also shows a solution for a phononic toy model consist of N ball connected to each other by springs with spring constant -γ.

 (36)

where n [0,...,N - 1], A = 1 for n = 0 and A = 1, otherwise. Note that to satisfy translational invariance condition, the diagonal terms in dynamical matrix are +2γ except in the boundaries (j = 1 and j = N) where they are +γ.

Fig. 4. 1D finite chain and ring. The energy levels and corresponding wave functions or orbitals for 1D finite chain and ring. The phononic mode for a finite chain of balls and springs with mass m.

4) Two dimensional (2D) square and hexagonal lattices

In section mbox II-Bmbox 2 , the band-structure and density of states of a 1D chain were calculated. Now let’s consider two most used 2D lattices: a square lattice where the unit-cell consist of one atom is connected to the first nearest neighbour in two dimensions (figure 5 a) and a hexagonal lattice where a unit cell consist of two atoms is connected to the neighbouring cells in which first (second) atom in a cell is only connected to the second (first) atom in any first nearest neighbour cell (figure 5 b). The TB Hamiltonian and corresponding band-structure could be calculated [1214] using equation 17 and the Bloch wave function of the form Aeikxj+ikyl (figure 5 ). The Schrödinger equation for two dimensional square lattice (figure 5 a) with on-site energies ε0 and hopping integrals γ could be written as:

 (37)

Using the Bloch function ϕj,l = Aeikxj+ikyl, the band structure of the 2D square lattice is obtained:

 (38)

Using similar approach, the band structure of a hexagonal lattice (e.g. graphene) could be obtained as:

 (39)

Figures 5 b,c,f,g show the bandstructure of the square and hexagonal lattices. The number of conduction channels (figures 5 d,h) could be calculated from the band structure of a crystalline system as we will discuss later in section III . Graphene is a semimetal and can be used as electrodes to probe molecules [7, 9, 26, 27].

Fig. 5. Two dimensional square and hexagonal lattices. Lattice geometry of (a) square and (e) hexagonal lattices, the bandstructure of (b,c) square and (f,g) hexagonal lattices and the number of conduction channels in (d) square and (h) hexagonal lattices.

##### C. Current carried by a Bloch function

The time evolution of density matrix ρt = |ψt⟩⟨ψt| allows us to obtain current associated with a particular quantum state |ψt. Using the time-dependent Schrödinger equation 1 ,

 (40)

By expanding |ψtover orthogonal basis |jequation 40 could be written as:

 (41)

Current carried by a Bloch function in a 1D chain: For a 1D infinite chain with the Hamiltonian of the form of equation 24 , the rate of change of charge Il = tldt at site l could be obtained by calculating the expectation value of both side of equation 41 over the state |l

 (42)

which could be simplified as

 (43)

where

 (44)

and

 (45)

These equations could be rewritten as:

 (46)

and

 (47)

The charge density is changing at atom site l as a result of two currents: right moving electrons Il+1 l and left moving electrons Il-1 l. The corresponding current to a Bloch state ψj(t) = eikj-iE(k)t∕ are:

 (48)

and

 (49)

where vk = ∂E(k)∂k = 2γsin(k)is the group velocity. Note that this defines the time that a wave-package need to move from one site to another e.g. l + 1 l, so the actual group velocity is vk × a, where a is the spacing between the sites. Although the individual currents are non-zero and proportional to the group velocity, the total current I = Il+1l + Il-1l for a pure Bloch state is zero due to an exact balance between left and right going currents. It is worth to mention that to simplify the notation, a Bloch state eikj is often normalized with its current flux 1 calculated from equations 48 and 49 to obtain a unitary current. Hence we will mostly use a normalized Bloch state eikj in later derivations. Furthermore, one important consequence of equations 46 and 47 is that if ψj = Aeikj + Be-ikj, although the charge density ρj = |ψj|2 is oscillating by j, the current is not oscillating by j and it is equal to the sum of the individual currents at site l only, Im(ψl*ψl+1) = |A|2sink -|B|2sink due to the Aeikl and Be-ikl.

Initial states are usually assumed to be stationary. However, if a non-stationary initial state was prepared in a closed (isolated) system such as a finite 1D chain of N atom, the charge density would be time-dependent (oscillatory) and therefore, current could be defined. As an example, for a system of two atoms coupled to each other by -γ, Hamiltonian reads $H=\bigl(\begin{smallmatrix} 0 &-\gamma\\ -\gamma& 0 \end{smallmatrix}\bigr)$. If the initial states are $\psi_1(t)=\bigl(\begin{smallmatrix} 1\\0 \end{smallmatrix}\bigr)$ and $\psi_2(t)=\bigl(\begin{smallmatrix} 0\\1 \end{smallmatrix}\bigr)$ which are non-stationary states, the final state is obtained $\psi(t)=\psi_1(t)+\psi_2(t)=\bigl(\begin{smallmatrix} cos(\gamma t/\hbar)\\i sin(\gamma t/\hbar) \end{smallmatrix}\bigr)$ which is not stationary. For such a closed system, the current could be obtained from equations 44 and 45 .

##### D. Parr-Pariser-Pople (PPP) Hamiltonian

Equation 17 described a non-interacting Hamiltonian which could also be written in the form of

 (50)

where i and j run over the orbitals centred on each site. For the orbital centred on site i and with spin s, ci,s is the electron creation operator that inserts an electron in state i and cj,s is the electron annihilation operator which takes an electron out of state i [1214]. ni,s = ci,scj,s is the electron number operator with ni = sni,s, εi is the energy of the orbital relative to the chiral symmetry point and γi,j are hopping integrals. Electron-electron interactions are missing from the non-interacting Hamiltonian (equation 50 ). In order to take the Coulomb electron-electron interaction into account, the Parr-Pariser-Pople (PPP) model can be used to write the interacting tight-binding Hamiltonian as

 (51)

where Uii and Uij are the on-site and long range Coulomb interactions, respectively given by the Ohno parametrization,

 (52)

and for ij

 (53)

where dij is the distance between sites i and j and U0 is the interaction amplitude e.g. U0 is equal to 11.26eV for gold and carbon and 9.95eV for sulfur.

### III. Nanoscale transport

Nanoscale transport can be described by three regimes:

(1) The self-consistent field (SCF) regime in which the thermal broadening kBT and coupling Γ to electrodes are comparable to Coulomb energy U0. The SCF method (single electron picture) implemented within non-equilibrium Green’s function (NEGF) could be used to describe transport in this regime as discussed in sections mbox III-Ambox to mbox III-Dmbox .

Fig. 6. Transport regimes. Transmission coefficient T of electrons with energy E traversing from one electrode to the other through a scattering region. Transport mechanism in a molecular junction could be either in tunnelling (off-resonance) regime where electrons tunnel through the molecule that is usually modelled with NEGF, or on-resonance where electrons transmit with high rate through an energy level that is modelled using master equation. The intermediate state (cross-over) between on and off resonance regimes are difficult to interpret either with NEGF or master equation.

In molecular junctions smaller than ~ 3 - 4nm, it is shown that transport remain elastic and phase coherent at room temperature. Therefore, using SCF models to describe the properties of these molecular junctions is well accepted in the mesoscopic community. Based on a single electron picture, the NEGF method coupled to the SCF Hamiltonian describes properties of the system on and off resonances (figure 6 ). Good agreement between these models and many room-temperature experiments suggest applicability of this method. A simplified Breit-Wigner formula ( mbox III-Ambox 5 ) derived from this method also could be used to model on-resonance transport through a device provided the level spacing is larger than resonances width. However, in those cases where the Coulomb energy has larger contribution, this method cannot describe the on-resonance properties of system.

(2) The Coulomb blockade (CB) regime in which Coulomb energy U0 is much higher than both thermal broadening kBT and coupling Γ where the SCF method is not adequate and multi-electron master equation should be used to describe properties of the system (figure 6 ) as discussed in section mbox III-Embox . This is usually needed to model the properties of molecular junctions at low temperature where an electrostatic gate voltage is applied through back gate.

(3) The intermediate regime (figure 6 ) in which the Coulomb energy U0 is comparable to the larger of the thermal broadening kBT and coupling Γ. There is no simple approach to model this regime. Neither the SCF method nor master equation could be used to well describe the transport in this regime because SCF method does not do justice to the charging, while the master equation does not do justice to the broadening.

##### A. Transport through an arbitrary scattering region

Consider the nanoscale junction of figure 7 where an arbitrary scattering region with Hamiltonian H connected to two single channel electrodes. On-site energies and coupling in the left (right) lead L(R) are εL (εR) and -γL (-γR), respectively. The leads are connected to the site a and b of the scattering region with the couplings -αL and -αR. The aim is to find the transmission t and reflection r amplitudes for a Bloch wave normalized with its current flux eikLj traveling from the left to right (figure 7 ).

If the wave function in the left and right leads and scattering region are ψj = eikLj + re-ikLj, ϕj = teikRj and fj, respectively; the Schrödinger equation in the left and right leads, the scattering region and connection points could be written as:

 (54)

 (55)

 (56)

 (57)

 (58)

From equations 54 and 58 , the E - k relations (band-structure) in the left and right leads are obtained as:

 (59)

Equation 56 could be re-written as = g where g = (E -H)-1 is Green’s function and called source which is a zero vector with non-zero elements in the connection points only (at site j = a and j = b). For the junction in figure 7 , has only two non-zero elements due to the source,

 (60)

where sa = -αLψ0 and sb = -αRϕ0. Furthermore, from equations 55 and 57 , the recurrence relation implies that:

 (61)

and

 (62)

Hence, by substituting ψ1 and ϕ-1 in equation 61

 (63)

From equation 60 and 63 ,

 (64)

Since sa = -αL(1 + r) and sb = -αRt∕, transmission t and reflection r amplitudes could be obtained.

 (65)

where

 (66)

is the surface Green’s function in the left and right leads at sites 0L and 0R (figure 7 ) and

 (67)

where ΣL,R = αL,R2gL,R are called self-energies due to the left and right contacts. The Green’s function in the surface of a semi-infinite lead (equation 66 ) can be obtained from the Green’s function of a doubly infinite crystalline lead. For example for the left electrode in figure 7 , the Green’s function of a doubly infinite crystalline chain is:

 (68)

To calculate the Green’s function at site j = 0L due to a source at site l = 0L (the surface Green’s function), equation 68 should vanish at site a (figure 7 ). This can be achieved by adding an appropriate wave function to equation 68 ,

 (69)

Hence, the Green’s function at site j = 0L due to a source at site l = 0L is gL,R = -eikL,RγL,R (equation 66 ). Assuming two identical leads (kL = kR = k and γL = γR = γ), equation 65 could be written as:

 (70)

where d = 1 + Δ1 + iΔ2 and Δ1 = Acos(k) + Bcos(2k), Δ2 = Asin(k) + Bsin(2k), A = (gaaαL + gbbαR)∕γ and B = αL2αR2(gaagbb - gabgba)∕γ2. From equation 70 , the transmission amplitude at E = 0 (e.g. k = π∕2) is

 (71)

Finally, transmission probability T = tt. More generally, the total transmission T and reflection R probabilities for multi-channel leads are obtained from

 (72)

and

 (73)

ti,j (ri,j) is the transmission (reflection) amplitude describing scattering from the jth channel of the left lead to the ith channel of the right (same) lead. Scattering matrix S is defined from ψOUT = IN and could be written by combining reflection and transmission amplitudes as:

 (74)

The S matrix is a central object of scattering theory and charge conservation implies that the S matrix is unitary: SS = I.

1) Transmission and reflection amplitudes from total Green’s function

As demonstrated in equation 65 , if the total Green’s function of a junction consisting two or few electrodes connected to an arbitrary scattering region is known, the transmission amplitude t (and transmission probability T ) for electrons traversing from one lead to the other could be calculated. The main task now is to find a method to calculate the Green’s function of whole system including crystalline leads (equation 68 ) connected to an arbitrary scattering region. Consider the nanoscale junction shown in figure 7 . The wave functions including and can be multiplied by any arbitrary amplitude

 (75)

without affecting the transport. Note that A does not depend on j. Using this amplitude, the wave-functions reads

 (76)

This equation looks like the Green’s function (equation 68 ) for j l. If we show that for j l,

 (77)

then ψj is the Green’s function of whole system at site j due to a source at site l and therefore transmission coefficient from any point to any other point can be obtained (equation 65 ). To demonstrate that equation 77 is valid for j l, consider

 (78)

where

 (79)

We shall show that gjl satisfy the Green’s function equation (E - H)gjl = δjl. We note that gjl can be written as

 (80)

where

 (81)

and since any wave-function can be added or subtracted from the Green’s function and the result is still a Green’s function, by subtracting from g, the new Green’s function ĝ is obtained:

 (82)

Substituting this into the Green’s function equation (E - H)gjl = δjl,

 (83)

The first and second terms are zero from equation 82. For j = l - 1, the third term -γgl-1,l = 1. Therefore

 (84)

is the Green’s function of whole system and describes the wave function at any site j due to a source at site l. Similarly, the wave-function

 (85)

is the Green’s function for a source in the left lead where j l. Therefore, the Green’s function of whole system e.g. Gij,

 (86)

the transmission amplitude t at j = 1 due to a source at l = 0 and the reflection amplitude r at j = 0 due to a source at l = 0 could be calculated:

 (87)

The transmission T and reflection R coefficients can then be obtained from equations 72 and 73 .

2) Scattering theory and Green’s function

Green’s function method has been widely used in the literature to model electron and phonon transport in nano and molecular scale devices and has been successful to predict and explain different physical properties. Green’s function is a wave function in a specific point of the system due to an impulse source in another point. In other words, Green’s function is the impulse response of the Schrödinger equation. Therefore, as shown in previous section, Green’s function naturally carries all information about the wave-function evolution from one point to the other in a system [1214, 28, 29].

The Green’s function G of a system with N site described by Hamiltonian H is defined as:

 (88)

where I is the identity matrix. Using the completeness condition,

 (89)

The Green’s function could be written in terms of eigenstates ψn and eigenenergies λn of H,

 (90)

and therefore the Green’s function element between point a and b is,

 (91)

Figure 8 shows how Green’s function could be used to calculate the transmission and reflection amplitudes in a simplest one dimensional system where two semi-infinite crystalline 1D leads are connected to each other through coupling β (representing the scattering region). The main question is what the amplitudes of the transmitted and reflected waves are? There are two main steps, first to calculate the total Green’s function matrix elements between sites 0 and 1 (G10) or 0 and 0 (G00); and secondly project these to the wavefunction to calculate transmission t and reflection r amplitudes (equation 87 ). For this example, the transmission and reflection probabilities are obtained from T = tt and R = rr.

Dyson’s equation describes the exact Green’s function of a system G = (g-1 -h)-1 in terms of the Green’s function of non-interacting parts g and Hamiltonian that connects them h. As shown in figure 8 , using the surface Green’s functions of the decoupled two semi-infinite leads $g=\bigl(\begin{smallmatrix} g_{00} & 0 \\ 0 & g_{11} \end{smallmatrix}\bigr)$ and the Hamiltonian that couples them h, the total Green’s function could be obtained from Dyson’s equation (first step). The second step is to use equation 87 to calculate t and r from the total green’s function G. It is worth to mention that Dyson’s equation could take different equivalent forms such as G = (g-1 - h)-1, G = g + ghG or G = g + gV g where V = (g-1 - h)-1.

Fig. 8.  Transport through a scatter connected to two 1D leads. For a Bloch wave $e^{ikj}/\sqrt{v_k}$ incident with a barrier, the wave is transmitted with the amplitude of t ($te^{ikj}/\sqrt{v_k}$) and reflected with the amplitude of r ($re^{-ikj}/\sqrt{v_k}$). Using the surface Green’s function of the leads (g00 and g11), the Hamiltonian of the scattering region in which bridge two leads h and Dyson’s equation, the total Green’s function G could be calculated. The Green’s function could then be used to calculate the transmission t and reflection r amplitudes.

3) Green’s function of N site finite chain and ring

Similar to the Green’s function of a semi-infinite chain (equation 69 ), from the Green’s function of the doubly infinite crystalline chain (equation 68 ), the Green’s function of N site finite chain and ring shown in figure 4 could be obtained [9] using appropriate boundary conditions as

 (92)

and

 (93)

These are useful equations to remember because they can help to understand quantum interference effects in simple molecules. As an example, consider a ring of 6 sites (N = 6) e.g. benzene with on-site energies ε = 0 and hopping integrals γ = -1. In the middle of the energy band e.g. E = 0, from the dispersion relation of a 1D chain (equation mbox II-Bmbox 2 ), the wave vector is k = π∕2. The Green’s function of benzene ring between any site i and j at the middle of the band is obtained from equation 93 as:

 (94)

For any odd to odd (oo) or even to even (ee) connectivities (fig. 4 ), (3 -|j - l|) is an odd number and therefore goooreebenzene = 0. This is called destructive quantum interference because the transmission between oo or ee sites are zero. In contrast, for any odd to even (oe) connectivity, (3 -|j -l|) is an even number and therefore goebenzene0 which is called constructive quantum interference. This implies non-zero transmission between any oe sites.

The Green’s function of the ring of 6 sites with on-site energies ε = 0 and hopping integrals γ can also be obtained by substituting its wavefunction (equation 35 ) into equation 91 . The eigenenergies of such system are: En = 2γcos(2nπ∕N) = [-2γ,-γ,-γ,+γ,+γ,+2γ]. The highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) levels are degenerate. At the middle of HOMO-LUMO gap E = EHL where EHL = (EH - EL)2 = 0 and EH (EL) is the energy of HOMO (LUMO) level, the Green’s function is obtained from equation 91 ,

 (95)

It is convenient to introduce the notation,

 (96)

where

 (97)

and A = [-1,1,1] for x = [1,2,3], respectively. For any oo or ee connectivities, j - l is an even number leading to a phase shift of 2π in the second term of equation 97 , that does not change the sign of second term. Therefore the magnitude of the first and second terms in equation 97 are equal with opposite sign and therefore g1ring = 0, g2ring = 0, g3ring = 0 and goooreering = 0 (the destructive quantum interference). This is a radical behaviour since contribution from all pares of HOMO – LUMO, HOMO-1 – LUMO+1 and HOMO-2 – LUMO+2 states to the Green’s function are zero.

In contrast, for any oe connectivity, j - l is an odd number leading to a phase shift of π in the second term of equation 97 , that changes the sign of second term. Therefore, the magnitude and sign of the first and second terms in equation 97 are equal, leading to non-zero values goering0 (the constructive quantum interference).

Consider a molecule which possesses only a HOMO ψH(l) of energy EH and a LUMO ψL(l) of energy EL, whose Green’s function from equation 91 is given by

 (98)

In this equation, ψH(l) and ψL(l) are the amplitudes of the HOMO and LUMO orbitals on connection site l, while ψH(m) and ψL(m) are the amplitudes of the HOMO and LUMO orbitals on connection site m. Since the core transmission coefficient for connectivity lm is given by τlm = (glm(E))2 a destructive interference feature occurs at an energy E given by glm(E) = 0, or equivalently

 (99)

If the energy E at which the destructive interference feature occurs lies within the HOMO-LUMO gap, then E - EH > 0 and EL - E > 0. This can only occur if the left hand side of equation 99 is positive and therefore the condition that a destructive interference feature occurs within the HOMO-LUMO gap is that the orbital products must have the same sign. Conversely, if they have opposite signs, there will be no destructive interference dip within the HOMO-LUMO gap. In the most symmetric case, where ψH(l)ψ*H(m) = ψL(l)ψ*L(m), this yields E = (EL - EH)2 and therefore the interference dip occurs at the middle of the HOMO-LUMO gap. On the other hand, if |ψH(l)ψ*H(m)| << |ψL(l)ψ*L(m)|, then E EH and the dip is close to the HOMO. In this case, for a real molecule with many orbitals, the approximation of retaining only the LUMO and HOMO breaks down, and the effect of the HOMO-1 should also be considered.

It is apparent from equation 98 that manipulating anti-resonances (e.g. due to the destructive quantum interference) is easier than resonances. To manipulate a resonance, red-ox state of a molecule should change whereas small environmental effects such as an inhomogeneous charge distribution or nearby ions could lead to a significant change in the position of anti-resonances.

Figure 9 shows two examples of polycyclic aromatic hydrocarbons, benzene and pyridine. Six molecular orbitals and corresponding eigenenergies due to six pz orbitals are shown in figure 9 a,b. HOMO and LUMO states are degenerate in benzene. These degeneracies are lifted in pyridine due to the presence of heteroatom (nitrogen). Consequently, the molecular orbitals are also affected.

Since Green’s function is the wave function due to a given source, we can visualise Green’s function just like molecular orbitals for given electron injection point and energy. Examples of Green’s function visualisation are illustrated in figure 9 c,d. Similar to wave functions visualisations, the radius and color of circles represent the amplitude and sign of Green’s function matrix elements, respectively due to a source in site shown by the green arrows. Figures 9 e,f show transmission coefficients of para and meta connectivities of benzene and a para connectivity of pyridine connected to two 1D leads through a weak coupling. Corresponding Green’s function illustrations are provided in the inset of figure 9 e,f at two different energies. Clearly, from the size of circles the main features of transmission is predictable. Therefore, one could use the Green’s function illustration of a molecule to predict transport intuitively.

4) Density of states from Green’s Function

The density of states for a system with eigenvalues λ is obtained from equation 28 . However, having calculated the Green’s function G from equation 90 , the density of states can be calculated. We note that a delta function δ(x - x0) can be defined as the limit of a function which exhibits a sharp peak about x0 and whose integral over space is 1. For Instance

 (100)

To prevent Green’s function to diverge at E = λ, equation 90 can be written as

 (101)

where η is a small number and

 (102)

and

 (103)

Since the eigenstates are orthonormal ψniψnj* = δij, we can find the expression for the trace of the Green’s function when η 0 as trace(Gi) = - n = -π nδ(E - λn). Therefore DOS is obtained,

 (104)

5) Breit-Wigner formula (BWF)

In the SCF regime, provided the coupling to electrodes was weak enough, level broadening on resonances due to electrodes was small enough and level spacing (differences between the eigenenergies of a quantum system) was large enough, the Green’s function gba in equation 65 for a system described by Hamiltonian H (figure 7 ) and energies close to an eigenvalue λm of H, is approximately

 (105)

where ya,b = fa,bm and |fmis the eigenvectors of H. This is a good approximation for E’s close to an eigenvalue λm if above mentioned conditions satisfied because the Green’s function g = n terms in nm are much smaller than in n = m. This yields to on-resonance transmission T for electrons with energy E passing through a molecule described by a Lorentzian like transmission function called BWF:

 (106)

where εn = λ - σL - σR, σL,R = ya,b2cos(kL,R) are the real part of the self-energies. ΓL,R = ya,b2sin(kL,R) are the imaginary part of the self-energies (Σ = σ + iΓ) which describe broadening due to the coupling of a molecular orbital to the electrodes. λ is the eigenenergy of the molecular orbital shifted slightly by an amount σ = σL + σR due to the coupling of the orbital to the electrodes. In this expression, ya2 and yb2 are the local DOS on the scattering region at the contact point. This formula shows that when the electron resonates with the molecular orbital (e.g. when E = εn), electron transmission is a maximum. The formula is valid when the energy E of electron is close to an eigenenergy λ of the isolated molecule, and if the level spacing of the isolated molecule is larger than (ΓL + ΓR). If ΓL = ΓR (a symmetric molecule attached symmetrically to the identical leads), T(E) = 1 on-resonance (E = εn).

If a bound state (e.g. a pendant group εp) is coupled (by coupling integral α) to a continuum of states, Fano resonances could occur [30, 31]. Fano resonance contains an anti-resonance followed by a resonance with an asymmetric line profile in between. Fano resonance originates from a close coexistence of resonant transmission and resonant reflection. This could be modeled by considering εn = λ-σL -σR + α2(E -εp) in BWF. At E = εp, electron transmission is destroyed (the electron anti-resonates with the pendant orbital) and at E = εn, the electron transmission is resonated by εn. The level spacing between this resonance and antiresonance is proportional to α.

Two levels BWF: As we discussed above, if a level broadening is smaller than level spacing between resonances, BWF can be used in the weak coupling regime. In case of two degenerate states, since these resonances can be close such that their level spacing is smaller than broadening, it is useful to drive a new form of BWF for two levels system. For a two levels system, the Hamiltonian H in figure 7 is given by

 (107)

where ε1 and ε2 are the energy levels coupled to each other by V d. If ε1 (ε2) is weakly bonded to the left (right) lead, the transmission coefficient T(E) could be obtained form equation 65 as [32, 33]:

 (108)

where σL,R (ΓL,R) are the real (imaginary) part of the self-energies due to the left L and right R leads.

Wigner delay time: Wigner delay time is the measure of time spent by an electron to pass from a scattering region of an open system. If the transmission amplitude of a given system t = |t|e (106 ) is defined by its magnitude |t| and phase θ, the Wigner delay time describes the phase difference between a scattered wave and a freely propagating one. Therefore, the Wigner delay time τw = dθ∕dE.

6) Open and close channels in leads

To calculate the number of open conduction channels in a 3D arbitrary crystalline lead, it is useful to consider a simple 2D cubic lattice with one orbital per site where each site is connected to its first nearest neighbouring sites as shown in the figure 5 . For simplicity, consider a finite system in the y direction Ny whereas the lattice is infinite in the x direction. A normalized wave function and the band-structure of such structure are calculated as ψkxm = sin(mπl∕(Ny + 1))eikxj and E(kx) = ε0 - 2γcos(mπ∕(Ny + 1)) + 2γcos(kx), respectively. Similar to the one-dimensional case (section mbox II-Cmbox ), current is associated to each ψkxm since every mini-band corresponds to a Bloch state. ψkxm are called channels. If we assume that the injected electrons from each lead to any individual channel are uncorrelated, the conductance at a given Fermi energy EF is given by G(EF) = 2e2M(EF)∕h where M(EF) is the number of open conduction channels at EF. In a one-dimensional lead with one orbital per site, there are either one open conduction channel or it is closed.

For the above quasi one dimensional system where x is the transport direction and y is the transverse direction with Ny atomic sites, the Green’s function in the sites l and j due to the source in the sites land jcould be written as:

 (109)

where kxm is longitudinal momentum and vxm = ∂E(kxm)∂kxm is the group velocity of channel m. Equation 109 could be re-written as:

 (110)

where ϕlm = sin(l). glj,lj consists of the sum of all allowed longitudinal modes eikxmj weighted by the corresponding transverse components ϕlm. Note that for a given E, kxm could be both real and imaginary. If kxm value of an eigenstate has no imaginary part Im(kxm) = 0, this state defined to be open, or propagating since a complex kxm will only occur if the wave is tunnelling, or decaying. The sign of the imaginary part of kxm (the group velocity) can be used to define the direction of a decaying wave (a propagating wave) as summarized in the table I .

TABLE I
Four classes of possible scattering channels
 left right Decaying: Im(kxm) > 0 Im(kxm) < 0 Propagating: Im(kxm) = 0,vxm < 0 Im(kxm) = 0,vxm > 0

Propagating (decaying) channels are conventionally called open (close) channels. It is worth to mention that retarded Green’s function of equation 110 is obtained by summing up all Ny scattering channels some of which are open channels and some others are closed channels. Figures 5 d,h show two examples of the number of conduction channels for the square and hexagonal lattices. The number of channels has a maximum in the middle of the band for the square lattice, whereas for the hexagonal lattice, there are fewer open channels (e.g. only two for graphene) in the middle of the band.

##### B. Generalized model to calculate transmission coefficient

In this section, an expression for the transmission coefficient Tnn = |sn,n(E,H)|2 between two scattering channels n,nof an open vector space A, in contact with a closed subspace B is obtained. The result is very general and makes no assumptions about the presence or otherwise of resonances. More precisely, we describe a quantum structure connected to ideal, normal leads of constant cross-section, labelled L = 1,2,. Consider two vector spaces A and B, spanned by a countable set of basis functions. In what follows, the sub-space B represents the structure of interest and sub-space A the normal leads, as shown in figure 10 .

Fig. 10. A sketch of a closed subspace B, in contact with subspace A through couplings W. XL denotes the surface of A connected to B. Subspace B includes some open subspaces connected to reservoirs shown by dotes.

The Hamiltonian is H = HA + HB + HJ, where HJ allows transitions between the subspaces. Since HJ can be written

 (111)

the Green’s function G for the combined space A B has the form

 (112)

To derive the more general formula where degenerate states can simultaneously resonate, note that when HJ = 0, G reduces to the Green’s function g of the decoupled system, where

 (113)

and Dyson’s equation G = g(1 - HJg)-1 yields

 (114)

Rewriting the above result for GAA in the form

 (115)

yields

 (116)

 (117)

and

 (118)

These demonstrate that once GBB is known, all other quantities are determined. To obtain an expression for transmission coefficients, it is convenient to introduce a set of states {|n⟩}, which span the subspace A and write gA = nm|ngnmm|. Since part of A consists of a number of ideal, straight, normal leads of constant cross-section, described by a real Hamiltonian, it is convenient to associate a sub-set of the states {|n⟩} with open channels of these leads. For these states, the notation |n= |n,xis introduced, where n is a discrete label identifying the lead, quasi-particle type, transverse kinetic energy and any other quantum numbers of an open channel and x is a position coordinate parallel to the lead. With this notation,

 (119)

where the prime indicates a sum over states |n,|morthogonal to the open channels (the close channels), and

 (120)

is the Green’s function of the semi-infinite lead between any position point x and xin the transport direction terminated at x = xL and vanishes at x = xL + a [34]. kxn is the longitudinal wave vector of channel n. If the lead belonging to channel n terminates at x = xL, then on the surface of the lead, the Green’s function gn(x,x) takes the form gn(xL,xL) = gn, where gn = an + ibn with an real and bn equal to π times the density of states per unit length of channel n. Moreover, if vn is the group velocity for a wave packet traveling along channel n, then vn = 2bn|gn|2. It is interesting to note that if x and xare positions located between xL and some point xn,

 (121)

If xn is some asymptotic position far from the end of the lead belonging to channel n and far from the scattering region (e.g. contact) defined by HJ, then the transmission amplitude t and the transmission coefficient T from channel nto channel n (nn) are

 (122)

and

 (123)

and since

 (124)

one obtains

 (125)

Let’s introduce eigenstates of HB, satisfying HB|fν= ϵν|fνand write

 (126)

From the expression for GBB given in equation 114 , this yields

 (127)

where

 (128)

Combining this with equation 119 yields

 (129)

In general, since the energy E lies in a region where the contribution to the density of states from gnm is zero, Imgnm = 0 and gnm = gnm*. For this reason it is convenient to introduce the notation

 (130)

 (131)

 (132)

 (133)

 (134)

and

 (135)

Clearly the matrices σ, σ(n) and Γ(n) are Hermitian. With this notation (GBB-1)μν = (E - ϵν)δμν - Σμν + iΓμν. Furthermore equation 125 becomes

 (136)

where the trace is over all internal levels of B and

 (137)

In these expressions, Γ(n) is a Hermitian matrix of inverse lifetimes, Γ = nΓ(n), σ and σ are Hermitian self-energy matrices and gB is the retarded Green’s function of subspace B when HJ = 0.

The form of equations 136 and 137 highlights the essential difference between open and closed channels. In the absence of open channels, σ and Γ are identically zero and if the subspace B is closed, GBB describes a quantum structure with well-defined energy levels, shifted by the self energy σ arising from contact with closed channels. Clearly no quasi-particle transport is possible through such a structure. When contact is made with open channels, the levels are further shifted by the self energy matrix σ and more crucially are broadened by the life-time matrix Γ.

For a system with non-orthogonal basis sets, in equation 128 , δμν should be replaced with the overlap matrix Sμν = fμ|fν. It is interesting to note that the vector spaces A representing the normal leads include both crystalline structures connected to the outside world and any close system coupled to the vector spaces B representing the structure of interest. In the latter case, the only effect of the closed part of the vector spaces A is to contribute in the scattering by its self-energy. Furthermore, figure 11 shows a slightly different approach to calculate the transmission (reflection) amplitude t (r) in a two terminal system with non-orthogonal basis set [28].

Fig. 11. Generalized transport model using Green’s function method. Generalized transport model using equilibrium Green’s function method [28] and its equivalent model for a simple 1D problem.

Surface Green’s function: An important step to calculate the total Green’s function of a system is to calculate the Green’s function at the surface of a semi-infinite lead. A lead is a perfect wave guide connected to a reservoir in the infinity from one side and to the scattering region from another side. It consists of identical slices described by Hamiltonian H0 connected to the first nearest neighbour with matrix elements H1. Note that for any lead the periodic slices could be chosen large enough to avoid the second and higher nearest neighbours interactions. The Green’s function of the semi-infinite lead at the point of contact to the scattering region is called the surface Green’s function. There are two main methods to calculate the surface Green’s function: analytic and recursive methods. In the analytic methods, first Green’s function of a doubly infinite lead is calculated. Then, a proper wave function is added to the Green’s function of a doubly infinite lead such that the Green’s function is vanished at the site next to the surface of the lead. This was discussed in the section mbox III-Ambox , equation 68 and the section mbox III-Bmbox .

In the recursive methods, the following non-linear equation 138 is solved iteratively using different algorithms such as a fixed-point iterative or Newton’s scheme.

 (138)

where gs is the surface Green’s function and η is a small number to avoid divergence of the Green’s function. In the fixed-point iterative scheme, equation 138 takes the following form:

 (139)

where n indicates the iteration number and gs1 = ((E -)I -H0)-1. The convergence of the fixed-point method can be quite poor, so more sophisticated methods such as the Newton’s scheme may be used [35].

##### C. The Landauer formula

Landauer used the scattering theory of transport as a conceptual framework to describe the electrical conductance and wrote ”Conductance is transmission” [36]. In the Landauer approach a mesoscopic scatterer is connected to two ballistic leads (see figure 1 ). The leads are connected to the reservoirs where all inelastic relaxation processes take place. The reservoirs have slightly different electrochemical potentials μL - μR 0 to drive electrons from the left to the right lead. Current therefore could be written as:

 (140)

where e = |e| is the electronic charge, T(E) is the transmission coefficient and f is Fermi-Dirac distribution function f(E -μ) = 1(1 + e(E-μ)∕kBT) associated with the electrochemical potential μ, kB is Boltzmann’s constant and T is temperature. The Fermi functions can be Taylor expanded over the range eV ,

 (141)

where μL - μR = eV . By including the spin, the electrical conductance G = I∕V reads:

 (142)

At T = 0K, - = δ(μ) where δ(μ) is the Kronecker delta. For an ideal periodic chain where T(E) = 1 at T = 0K, the Landauer formula becomes:

 (143)

G0 is called the ”Conductance Quantum”. In other words, current associated with a single Bloch state vk∕L and generated by the electrochemical potential gradient is I = e(vk∕L)DΔμ where D = ∂n∕∂E = L∕hvk. It is worth mentioning that the Landauer formula 141 describes the linear response conductance, hence it only holds for small bias voltages, δV 0.

1) Landauer-Buttiker formula for multi-terminal structures

Conductance measurements are often performed using a four-probe structure to minimize the contact resistance effect. Also multi-probe structures are widely used to describe the Hall-effect or in sensing applications. Based on the Landauer approach for two terminal system, Buttiker [37] suggested a formula to model multi-probe currents for structures with multiple terminals as:

 (144)

where Ii is current at ith terminal and Tij is the transmission probability from terminal j to i. In a multi-terminal system, it is convenient to assume one of the probes as reference voltage V ref = 0 and write currents based on that. As an example, for a four probe structure, current in each probe can be written as:

 (145)

where Ni is number of open conduction channels in lead i. In the four probe structure, if probe 3 and 4 are the outer voltage probes (I3 = I4 = 0) and probe 1 and 2 are the inner current probes, the four probe conductance is

 (146)

2) Equilibrium vs. non-equilibrium I-V

The Landauer formula only holds in the linear response regime. A transmission coefficient T(E) of particle with energy E from one electrode to the other is obtained in the steady state condition where the junction is assumed to be close to equilibrium (δV 0). In this regime, transmission coefficient is assumed to be not voltage dependent and current is calculated using equation 140 .

However, if transmission coefficient changes by applied bias voltage (the non-linear regime), bias voltage dependent transmission coefficient T(E,V b) should be calculated. In order to take the effect of electric field on T(E) into account, new Hamiltonian for a given field is calculated. This is obtained by calculating the potential profile applied to the junction due to the given electric field (e.g. bias voltage) using Poisson’s equation 2U = -ρ∕ε where U is potential profile due to charge distribution ρ and ε is permittivity which could vary spatially [1214]. In the non-equilibrium condition, the Landauer formula (equation 140 ) takes the form,

 (147)

where V b and V g are bias and gate voltages, respectively and T(E,V b,V g) is transmission coefficient calculated at each bias and gate voltages. It is worth mentioning that in some experiments, the measured conductance G = I∕V b is noisy. Therefore, the differential conductance map Gdiff(V b,V g) = dI(V b,V g)∕dV b is plotted. This could be calculated by differentiation of equation 147 with respect to the bias voltage V b.

##### D. Non-equilibrium Green’s function formalism

If ES|ψ= H|ψdescribes the properties of the closed system H with non-orthogonal basis set S, then once it connects to the outside world and becomes an open system (see figure 12 ), a modified Schrödinger equation in the non-equilibrium condition could be written [1214]:

 (148)

where the terms Σ|ψand |sdescribe the outflow and inflow (e.g. see equation 60 ), respectively arises from the boundary conditions. Equation 148 could be rewritten as

 (149)

where GR = [ES - H - Σ]-1 is retarded Green’s function (GA = [GR]), Σ = Σ1 + Σ2 + Σ0 is sum of self-energies due to the electrodes Σ1, Σ2, and surroundings Σ0 such as dephasing contact or inelastic scattering e.g. electron-phonon coupling, emission and absorption. Dephasing contact terms could be described by the SCF method whereas for inelastic processes one needs to use for instance Fermi’s golden rule to describe these self energies.

There are a few proposals in the literature to treat incoherent and inelastic processes [1214, 38]. Buttiker [38] suggests to treat the inelastic and incoherent scattering by introducing a new probe to the original coherent system. This could be seen as assigning new self-energies associated with any inelastic or incoherent processes. However, if the incoherent and inelastic effects are treated by introducing an extra electrode, corresponding distribution function e.g. Fermi function for electrons is assigned which in general may not be the case. More generally, any incoherence and/or inelastic processes are introduced by appropriate self energies which not necessarily described by equivalent Fermi functions in the contact.

For a normal, coherent elastic junction if H1,2 are the coupling matrices between electrode 1 (2) and scattering region and g1,2 are the surface Green’s function of the electrodes, Σ1,2 = H1,2g1,2H1,2. Furthermore, current could be calculated as

 (150)

where Γ1 = i1 - Σ1) is the imaginary part of the self-energy, Gn is the density matrix, 1in = s1s12π is related to the source |sand A = GRΓGA is the spectral function as shown in figure 12 . Note that the definition of self-energy here is a bit different than what was discussed in section 3.2. From the basic law of equilibrium, in a special situation where we have only one contact connected; the ratio of the number of electrons to the number of states must be equal to the Fermi function in the contact ( 1,2in = Γ1,2f1,2(E)). However, in dephasing contact, Σ0in is not described by any Fermi function and since inflow and outflow should be equal Trace[ 0inA] = Trace0Gn]. Figure 12 summarize the basic non-equilibrium Green’s function (NEGF) equations to calculate current in a most general junction where surroundings are present. In the absence of surroundings, current (equation 150 ) in lead i could be re-written as [1214]:

 (151)

where Tij(E) = Tracei(E)GR(Ej(E)GA(E)] is the transmission coefficient for electrons with energy E passing from lead i to lead j.

Fig. 13. Two terminal system with two 1D leads connected to a scattering region ε1.

It is worth mentioning that in the past decade, several numerical implementations of the scattering theory and non-equilibrium Green’s functions approach such as Gollum [29], SMEAGOL [28], TransSIESTA [39] and TURBOMOLE [40] have been developed.

Transport through one-level system: Consider two identical 1D leads with on-site energies ε0 and hoping integrals γ connected to a scattering region ε1 with coupling integrals α and β as shown in figure 13 . The transmission coefficient T for electrons with energy E traversing from the left to right leads can be calculated as T(E) = ΓL(E)GR(ER(E)GA(E) where the retarded Green’s function is GR(E) = (E - ε1 - Σ)-1, the self-energies Σ = ΣL + ΣR is obtained from ΣL = α2eik∕γ and ΣR = β2eik∕γ and the broadening due to the left and right leads are ΓL = iL - ΣL) = -2α2sin(k)∕γ and ΓR = iR - ΣR) = -2β2sin(k)∕γ. Note that for a one-level system, this equation can be exactly re-written in the form of Breit-Wigner formula (equation 106 ). However, for the two or more level systems, BWF is a good approximation only to describe on-resonances transport and if the level spacing is less than level broadening (see section mbox III-Ambox 5 ).

Mapping quantum to semi-classical model: Current in the semi-classical method is defined as the flow of charge N in the unit time t as I = . As shown in figure 14 , in the contact one, there are two in-going S1D0 and outgoing ν1N currents where ν1 is coupling strength to the contact 1, S1 = ν1f1 is source and D0 = D(E)dE is total density of state in the scattering region. Therefore, current in the contact 1 is written as:

 (152)

where f1 and f2 are the Fermi distribution in the contact 1 and 2, respectively.

Fig. 14. Semi-classical method to calculate current [1214].

In order to build a physical intuition of NEGF equations (figure 12 ), we compare the semi-classical picture (figure 14 ) with the quantum model (figure 12 ),

 (153)

A and Gn look like matrix version of the density of states and the total number of electrons, respectively. However, the exact relations are: D = -∞+dEtrace(A(E))2π and N = -∞+dEtrace(Gn(E))2π. Clearly, the semi-classical picture is missing most of the interesting effects such as quantum interference and connectivity dependence of conductance [6, 41, 42].

##### E. Master equation

As discussed in section III , if the Coulomb energy is much higher than thermal broadening and coupling to electrodes (the on-resonance transport regime) in a molecular junction, the SCF method is not adequate to describe properties of the junction. For example, when Fermi energy is moved from one side of a resonance to the other using bias or gate voltages, the redox state of molecule may change by gaining or losing an electron. In order to account for this effect, the occupation probability of each state should be calculated. This is described by the multi-electron master equation. In the multi-electron picture, the overall N-electron system has different probabilities Pα of being in one of the 2N possible states α. Furthermore all of these probabilities Pα must add up to one. The individual probabilities could be calculated under steady-state conditions where there is no net flow into or out of any state (see figures 15 and 16 )

 (154)

where R(α β) is the rate constants obtained by assuming a specific model for the interaction with the surroundings. In a system that electrons can only enter or exit from the source and drain contacts, these rates are given in figures 15 and 16 for one and two electron systems. Equation 154 is called a multi-electron master equation [1214]. The general principle of equilibrium statistical mechanics could be used to calculate the probability Pα that the system is in state α with energy Eα and Nα electrons as

 (155)

where kB is Boltzmann’s constant and T is temperature.

Fig. 15. Spin-degenerate one electron system with energy ε. There are two possibilities either the state is full |1or empty |0. The rate to move an electron into the state (|0to |1) is a product of the rates that electron can go in and the Fermi function in the leads (R(|0⟩→|1) = γ1f1 + γ2f2) [1214]

One-electron energy levels represent differences between energy levels corresponding to states that differ by one electron. If E(N) is the energy associated with the N-electron state, the energy associated with the addition of one electron is called affinity energy EA = E(N) - E(N + 1). Similarly the energy needed to remove one electron is called ionization energy IP = E(N - 1) - E(N).

One electron system: Figure 15 shows the master equation for spin-degenerate one electron system with energy ε where there are only two possibilities, either the state is full |1or empty |0. Current is then calculated as:

 (156)

where γ1 and γ2 are the rates that electron can go in and out from the left and right electrodes with f1(E) and f2(E) Fermi functions.

Two electron system: When the number of electrons N increases, the number of possibilities increases by 2N. In two electron system, there are four possibilities (22), both states empty |00or full |11and either one of them full and the other empty (|01and |10).

Fig. 16. Two electron system. There are four possibilities, both states empty |00or full |11and either one of them full and the other empty (|01and |10) [1214].

Figure 16 shows how to calculate current in two electron system [1214]. Note that as soon as a state become full, there need an additional energy (Coulomb repulsion energy) to have second electron in the other state in addition to the level spacing between the two energy levels. Furthermore, it is not correct to assume one Fermi function for all transitions. Due to the Coulomb blockade energy, each level needs certain electrochemical potential to overcome the barrier and current flow. Clearly, the expression for current in two electron system is more complicated than one electron system. This becomes even more complicated when multi-electron system was considered. In this case, the number of calculations rapidly increases, therefore efficient numerical algorithms are needed to solve the multi-electron master equation.

Coulomb and Franck-Condon blockade regimes: The electronic properties of weakly coupled molecules are dominated by Coulomb interactions and spatial confinement at low temperatures [43]. This could lead to Coulomb blockade (CB) regimes in which the transport is blocked due to the presence of an electron trapped in the junction. In addition, charge transfer can excite vibrational modes or vibrons, and strong electron-vibron coupling leads to suppression of tunnel current at low bias called Franck-Condon (FC) blockade regimes.

To describe transport in this regime, a minimal model (the Anderson-Holstein Hamiltonian) can be used [44] that captures the CB, FC and the Kondo effect if three assumptions are made: (1) the relaxation in the leads assumed to be sufficiently fast leading to Fermi functions for the distribution of the electrons in thermal equilibrium at all times; (2) the transport through the molecule is dominated by tunnelling through a single, spin-degenerate electronic level, and (3) one vibron taken into account within the harmonic approximation. In this case, the Anderson-Holstein Hamiltonian reads H = Hmol + Hleads + HT with

 (157)

describing the electronic and vibrational degrees of freedom of the molecule,

 (158)

 (159)

the tunnelling between the leads and molecule. Here, Coulomb blockade is taken into account via the charging energy U where eV,kBT << U. The operator dσ (dσ) annihilates (creates) an electron with spin projection σ on the molecule, nd = σdσdσ denotes the corresponding occupation-number operator. Similarly, capσ (capσ) annihilates (creates) an electron in lead a (a = L,R) with momentum p and spin projection σ. Vibrational excitations are annihilated (created) by b (b). They couple to the electric charge on the molecule by the term ~ nd(b + b), which can be eliminated by a canonical transformation, leading to a renormalisation of the parameters ε and U, and of the lead-molecule coupling ta tae-λ(b+b) . The master equations determining the molecular occupation probabilities Pqn for charge state n and vibrons q is:

 (160)

Pqeq denotes the equilibrium vibron distribution with a relaxation time τ and Wqqnn denotes the total rate for a transition from |n,qto |n,q′⟩.

 (161)

where fa is the Fermi function and the transition rates Γ are calculated from Fermi’s golden rule.

 (162)

Here, ρa denotes DOS in lead a, Mqq;ann±1 denotes the FC matrix elements and snm the spin factor [45] such that for sequential tunnelling and assuming twofold degeneracy they are s10 = s12 = 1,s01 = s21 = 2. The matrix elements Mqq;ann±1 defined for vibrations are

 (163)

where q1 = min{q,q′} and q2 = max{q,q′}. This minimal model captures the main features of resonant tunnelling due to Coulomb energy and vibrons effect in low temperature [45, 46].

##### F. Minimal phase-coherent CB description

Since transport through molecular junctions is usually phase coherent even at room temperature and at low temperatures exhibits Coulomb blockade, we can employ a minimal model by including the effect of the Coulomb energy and computing the differential conductance as a function of the bias and gate voltages. This model preserves phase coherence, implements a derivative discontinuity to describe Coulomb blockade and avoids self-interactions. The bias and gate voltage depended differential conductance is calculated by differentiating bias-dependent current with respect to the bias voltage Gdiff = dI(V b,V g)∕dV b using the Landauer formula:

 (164)

Note that a factor of 2 in this equation accounts for spin so that T = (T1 + T2)2 where T1 (T2) is transmission function for majority (minority) spins. This method can lead to long calculation times, which can be reduced by considering only the first term of the Taylor expansion of the Fermi-Dirac distribution function f(E). Therefore, a simplified and approximated form of Gdiff is obtained in the limit of low temperatures

 (165)

To include the effect of additional or Coulomb energy, for each new bias, we can start from the ’bare’ tight binding Hamiltonian and then we (1) add the gate voltage to each site energy, (2) diagonalise the Hamiltonian (3) if an eigenvalue drops below the drain voltage and the level becomes occupied we add the Coulomb energy to all other eigenvalues and then transform back to the original basis. (4) T(E,V b,V g) and the corresponding differential conductance is computed using this new Hamiltonian and steps 1 to 4 are repeated. Figure 17 shows the stability diagram obtained by applying this minimal model calculation with and without additional energy.

Fig. 17.  Stability diagram.(a) a tight-binding model consists of two one-dimensional leads with hopping elements of γ = -2 (which sets the energy scale) connected to a scattering region containing 10 sites with hopping elements γi = -[0.25, 0.125, 0.25, 0.25, 0.075, 0.2, 0.2, 0.25, 0.25, 0.4]. All on site energies (εi and ε0) are set to zero. The scattering region is coupled to the left and right electrodes by hopping elements α = 0.1γ. (b,d) The differential conductance dI(V b = 0,V g)∕dV b. (c,e) The stability diagram dI(V b,V g)∕dV b is obtained by applying minimal model calculation. In b,c the the Coulomb energy U is zero whereas in d,e U = 0.125γ. Chessboard pattern due to quantum interference is obtained in c. By including additional energy U, the size of Coulomb diamonds (blockade region) is increased accordingly and excited states do not penetrate to the blockade region.

#### IV. Modelling the experiment

So far we have discussed, different transport regimes and the methods to model electron and phonon through nanoscale junctions. However, all of these tools are useful if they can explain new observed physical phenomenon or predict new properties of a future physical system and suggest new design principles. Experiments in the field of molecular electronics either study new junction physical properties such as conductance and current or focus on using well characterized junctions for future applications. It is important to understand the limits in theory and experiment. For example, there are certain quantities that are not directly measurable but can be calculated such as wave-functions. In contrast, there are quantities that are not predictable but accessible in experiment such as the position of the Fermi energy, the overall effect of the inhomogeneous broadening on the transport, or screening effects which is related to the exact junction configuration in the real-time experiment.

Furthermore, more reliable predictions can be made if a series of junctions were studied and the overall trends compared between theory and experiment. The bottom line is theory and experiment are not two isolated endeavours. They need to communicate with each other in order to discover new phenomena. Those quantities that cannot be computed reliably, but for which experimental data is available, can be used to correct and refine theoretical models. Usually to explain new phenomena, one needs to build a model based on a working hypothesis. To make an initial hypothesis, a theorist needs to know how to take into account different physical phenomena such as the effect of environment, presence of an electric or magnetic field. In the following, the aim is to make few bridges between the well-known physical phenomena and the methods to model them theoretically.

Let’s start by considering the differences between a lead and a channel. From a mathematical viewpoint, channels connect an extended scattering region to a reservoir and the role of lead i is simply to label those channels ki,qi, which connect to a particular reservoir i. Conceptually, this means that from the point of view of solving a scattering problem at energy E, a single lead with N(E) incoming channels can be regarded as N(E) virtual leads, each with a single channel. We could take advantage of this equivalence by regarding the above groups of channels with wave-vectors kαi,qαi as virtual leads and treating them on the same footing as physical leads [29].

This viewpoint is particularly useful when the Hamiltonians H0i, H1i describing the principle layers PLs (the identical periodic unit cells H0i connected to each other by H1i) of the physical lead i are block diagonal with respect to the quantum numbers associated with kαi,qαi. For example, this occurs when the leads possess a uniform magnetization, in which case the lead Hamiltonian is block diagonal with respect to the local magnetization axis of the lead and α represents the spin degree of freedom σ. This occurs also when the leads are normal metals, but the scattering region contains one or more superconductors, in which case the lead Hamiltonian is block diagonal with respect to particle and hole degrees of freedom and α represents either particles p or holes h. More generally, in the presence of both magnetism and superconductivity, α would represent combinations of spin and particles and holes degrees of freedom.

In all of these cases, H0i, H1i are block diagonal and it is convenient to identify virtual leads αi with each block, because we can compute the channels kαi,qαi belonging to each block in separate calculations and therefore guarantees that all such channels can be separately identified. This is advantageous, because if all channels of H0i, H1i were calculated simultaneously, then in the case of degeneracies, arbitrary superposition of channels with different quantum numbers could result and therefore it would be necessary to implement a separate unitary transformation to sort channels into the chosen quantum numbers. By treating each block as a virtual lead, this problem is avoided.

##### B. Charge, spin and thermal currents

When comparing theory with experiment, we are usually interested in computing the flux of some quantity Q from a particular reservoir. If the amount of Q carried by quasi-particles of type αi is Qαi(E), then the flux of Q from reservoir i is:

 (166)

Tαiji,j in this expression is transmission coefficient of quasi-particles of type αi. In the simplest case of a normal conductor, choosing Qαi = -e , independent of αi, this equation yields the electrical current from lead i. αi may represent spin, and in the presence of superconductivity it may represent hole (αi = h) or particle (αi = p) degrees of freedom. In the latter case, the charge Qp carried by particles is -e, whereas the charge Qh carried by holes is +e. In the presence of non-collinear magnetic moments, provided the lead Hamiltonians are block diagonal in spin indices, choosing αi = σi and Qαi = -e in Eq. (166 ) yields for the total electrical current

 (167)

Note that in general it is necessary to retain the subscripts i,j associated with σi or σj, because the leads may possess different magnetic axes.

Similarly the thermal energy carried by the electrons from reservoir i per unit time is

 (168)

For the special case of a normal multi-terminal junction having collinear magnetic moments, αi = σ for all i and since there is no spin-flip scattering, Tσ,σi,j = Tσ,σi,jδσ,σ [29]. In this case, the total Hamiltonian of the whole system is block diagonal in spin indices and the scattering matrix can be obtained from separate calculations for each spin. we assume that initially the junction is in thermodynamic equilibrium, so that all reservoirs possess the same chemical potential μ0. Subsequently, we apply to each reservoir i a different voltage V i, so that its chemical potential is μi = μ0 - eV i. Then from equation (166 ), the charge per unit time per spin entering the scatterer from each lead can be written as

 (169)

and the thermal energy per spin per unit time is

 (170)

where e = |e| and fσi(E) = f(E - μi) - f(E - μ) is the deviation in Fermi distribution of lead i from the reference distribution f(E - μ). In the limit of small potential differences or small differences in reservoir temperatures, the deviations in the distributions from the reference distribution fσj(E) can be approximated by differentials, therefore by Taylor expanding fj = f(E - μj), fj - f = -(μj - μ) -()(Tj - T) is obtained. Using this expression for fj - f, equations 169 and 170 could be written as:

 (171)

Since nth moment of probability distribution P(x) define as: xn= dxP(x)xn, the spin-dependent moment Lij,σ in the presence of collinear magnetism is obtained as

 (172)

where f(E,T) = (1 + e(E-EF)∕kBT)-1 is the Fermi-Dirac distribution function, T is temperature and kB is Boltzmanns constant. Therefore, in the linear-response regime, the electric current I and heat current passing through a device is related to the voltage difference ΔV and temperature difference ΔT by

 (173)

where electrical conductance G (thermal conductance κel) is the ability of the device to conduct electricity (heat) and the Seebeck coefficient S (Peltier Π) is a measure of generated voltage (temperature) due to a temperature (voltage) differences between two sides of the device. In the presence of two leads labeled i = 1,2, the spin-dependent low-voltage electrical conductance G(T,EF), the Seebeck coefficient (somtimes called thermopower) S(T,EF) = -ΔV∕ΔT , the Peltier coefficient Π(T,EF) and the thermal conductance due to the electrons κel(T,EF) as a function of Fermi energy EF and temperature T can be obtained as

Note that the thermal conductance is guaranteed to be positive, because the expectation value of the square of a variable is greater than or equal to the square of the expectation value. Furthermore, by Taylor expanding Tσ,σij(E) around EF, Tσ,σij(E) Tσ,σij(EF) + |E=EF(E -EF), the low temperature approximation of Lij,σ1 and Lij,σ0 where Tσ,σij(E) varies approximately linearly with E on the scale of kBT could be written as:

 (175)

and

 (176)

and

 (177)

where α = kB2π23e2 WΩK-2 is the Lorenz number. Therefore, Seebeck coefficient could be re-written as:

 (178)

Equation 178 is also called the Mott Formula. From equations 175 -177 , the electrical conductance and thermal conductance due to electrons are obtained as G(T,EF) G0T(EF) and κel(T,EF) αTG(T,EF). The latter is also called the Wiedemann-Franz law and shows that thermal conductance due to electrons is proportional to the electrical conductance.

Efficiency of a thermoelectric martial η is defined as the ratio between the work done per unit time against the chemical potential difference (between two hot and cold reservoir) and the heat extracted from the hot reservoir per unit time. The maximum efficiency ηmax could be written as:

 (179)

where Th and Tc are temperatures of hot- and cold-sides, respectively, ΔT = Th -Tc and Tavg = (Th + Tc)2. The thermoelectric conversion efficiency (equation 179 ) is the product of the Carnot efficiency () and a reduction factor as a function of the material’s figure of merit Z = S2G∕κ, where S, G, and κ = κel + κph are the Seebeck coefficient, electrical conductance, and thermal conductance due to both electrons and phonons, respectively. More commonly a dimensionless figure of merit (ZT = Z.Tavg) is used to account for the efficiency of the thermoelectric materials. The thermoelectric figure of merit could be written as

 (180)

where the electronic thermoelectric figure of merit for a two-terminal system is

 (181)

To calculate the total ZT , not only the thermal conductance due to the electrons are needed but also it is absolutely crucial to take the phonons contribution to the thermal conductance (κph) into account as described in the next section.

##### C. Phonon thermal conductance

To calculate the heat flux through a molecular junction carried by the phonons, equation 166 could be used where the thermal conductance due to the phonons κph could be obtained [10] by calculating the phononic transmission Tph for different vibrational modes as

 (182)

where fBE(ω,T) = (eω∕kBT - 1)-1 is Bose-Einstein distribution function and is reduced Plancks constant and kB is Boltzmanns constant. To calculate the vibrational modes of a system, we use the harmonic approximation method to construct the dynamical matrix D. From the ground state relaxed xyz coordinate of the system, each atom is displaced from its equilibrium position by δqand -δq in x, y and z directions and the forces Fiq = (Fix,Fiy,Fiz) in three directions qi = (xi,yi,zi) on each atom are calculated. For 3n degrees of freedom (n = number of atoms), the 3n × 3n dynamical matrix D is constructed from Hessian matrix K

 (183)

where Kijqq for ij are obtained from finite differences

 (184)

and the mass matrix M = . To satisfy momentum conservation, the Ks for i = j (diagonal terms) are calculated from Kii = - ijKij. Once the dynamical matrix is constructed the Green’s function method as described in section mbox III-Bmbox could be used to calculate the phononic transmission coefficient Tph.

It is worth to mention that in order to compute the electron-phonon coupling matrices Wλ for a given phononic mode λ, one could use the Hamiltonian matrices {{⟨i|H|j⟩}} for the displaced configurations [47].

 (185)

where |iare basis orbitals. The Wλ matrices are then used [47] to calculate the appropriate self-energies (see section mbox III-Dmbox ) to account for inelastic scattering due to electron-phonon interactions.

##### D. Piezoelectric response

In piezoelectric materials, electric field is generated due to conformational changes. The piezoelectric response between atoms i and j of a molecule is given by:

 (186)

where rij is the distance between atom i and j in the equilibrium geometry. Displacement derivatives is given by

 (187)

where W and V are normal modes and vectors of Hessian matrix K = , respectively (section mbox IV-Cmbox ). H = V TH where H = is dipole derivatives matrix. Note that just like equation 184 , these derivatives can be calculated using finite differences. Clearly, larger piezoelectric response is expected for molecules with larger dipole moments. Recently, a large converse piezoelectric effect was measured in helicene single molecules [48].

Although DFT has been successful at predicting the trends, it usually underestimates the position of the Fermi energy EF, the exact energy levels (Kohn-Sham eigenvalues [49]) and therefore the position of the HOMO and LUMO and the energy gap. To compare mean-field theory with experiment, some corrections are needed. One way is to use hybrid functionals e.g. B3LYP [50] or many body calculations e.g. GW approximation [51]. The latter is computationally expensive [52] and can only be used for molecules with few atoms. An alternative way is to correct the HOMO-LUMO gap using the values measured experimentally. A phenomenological scheme that improves the agreement between theoretical simulations and experiments in, for example, single-molecule electronics consists of shifting the occupied and unoccupied levels of the molecule (M) downwards and upwards, respectively to increase the energy gap. The procedure is conveniently called spectral adjustment in nanoscale transport [29, 53] and has been widely used in the literature to correct theoretical transport gap [54, 55]. The Hamiltonian K = H - ES of a given M region could be modified as:

 (188)

where $\Delta_\mathrm{o,u}$ are energy shifts. $\rho_\mathrm{M} = \sum_{no}\,|\Psi_{no}\rangle\langle\Psi_{no}|\nonumber$ is the density matrix where (no, nu) denote the occupied and unoccupied states, respectively and SM is overlap matrix. If experimental HOMO and LUMO energies are available, $\Delta_\mathrm{o,u}$ can be chosen to correct HOMO and LUMO obtained from the mean field Hamiltonian. Alternatively, in the simplest case, the shifts $\Delta_\mathrm{o,u}$ are chosen to align the highest occupied and lowest unoccupied molecular orbitals (ie the HOMO and LUMO) with (minus) the ionization potential (IP) and electron affinity (EA) of the gas phase molecule

 $\Delta_\mathrm{o}^0=-\varepsilon_\mathrm{HOMO}-IP+Z$ $\Delta_\mathrm{u}^0=-\varepsilon_\mathrm{LUMO}-EA-Z$ (189)

The ionization potential (IP = E(+e) -E(0)) and electron affinity (EA = E(0) -E(-e)) are calculated from total energy calculation of the gas phase molecule. E(0) is the total energy of the neutral molecule, E(+e) is the energy of the molecule with one electron removed (i.e. positively charged), and E(-e) is the total energy of the molecule with one added electron. The energy-gap Eg of a molecule (sometimes called additional energy) could be calculated from IP and EA as: Eg = IP - EA [1214]. The important conceptual point is that the electrochemical potential μ should lie between the affinity levels (above μ) and ionization levels (below μ) in the ground state. Note that IP and EA are traditionally defined positive energies below the vacuum level, whereas HOMO and LUMO level are negative, if they are below the vacuum level. The Coulomb interactions in the isolated molecule are screened if it is placed in close proximity to a metallic surface. Due to this image charge interactions, the occupied energy levels shift up whereas the unoccupied states move down resulting in shrinking of energy gap. This could be taken into account by using a simple image charge model, where the molecule is replaced by a point charge located at the middle point of the molecule and where the image planes are placed ~ 1 Å above the electrodes’ surfaces. Then the shifts are corrected by screening effects Z = e2 ln2(8πϵ0a) where a is the distance between the image plane and the point image charge and ϵ0 = 8.85 × 10-12F∕m is the vacuum permittivity. F. Hopping versus tunnelling Charge transport mechanism through the molecular scale junctions are either coherent transport via tunnelling (called also superexchange) or incoherent thermally activated hopping. Transport through the short molecules with the length of < 3 - 4nm has been demonstrated [56, 57] to be coherent tunnelling. This is characterized by an exponential dependence of conductance G with length given by

 (190)

where Gc is a pre-exponential factor dependent on junction contacts and nature of metallic leads, β is the tunnelling decay constant (called also attenuation or β factor), and L is the length of the molecule. The coherent process is also characterized by temperature independence. In contrast, incoherent hopping is believed to be the charge transport mechanism along longer molecules. In incoherent hopping transport regime, the conductance follows an Arrhenius relation

 (191)

where Ga is a constant pre-exponential factor for each chemical reaction, EA is the hopping activation energy, T is the temperature and kB is the Boltzmann’s constant. In this regime, the conductance varies as the inverse of molecular length and also decays exponentially with temperature. Therefore, the length and temperature dependent conductance measurements are widely used in the literature to distinguish between different transport regimes [58, 59].

Fig. 18.  Magnetisim. Schematic represents a scattering region with two 1D chain for each spin. ε1 and ε2 are site energies for spin-up and spin-down, respectively. In the para-magnetic case, the site energies for spin-up and spin-down are equal ε0. In the presence of spin-orbit coupling, there are coupling between spin-up and spin-down sites.

##### G. Spin Hamiltonian

To take electrons spin into account in a Hamiltonian H (= ), in the absence of spin-orbit coupling, the Schrödinger equation can be written,

 (192)

where ψ and ψ are spin-up and spin-down component of wavefunction. If electrons travel with high velocities, relativistic effects can become significant. This is not usually a case in solids but near the nuclei of atoms where the electric fields are high, weak relativistic effects might be expected. Therefore, a spin-orbit correction need to be considered using the Dirac equation.

 (193)

Spin-orbit Hamiltonian is often written as σ.M where σ is Pauli spin matrices.

 (194)

If M points in the direction (θ,ϕ), Mx = sinθcosϕ, My = sinθsinϕ and Mz = cosθ, therefore

 (195)

The spin wavefunction, spinors thus is obtained

 ψ = , ψ =

. If there is no spin-orbit coupling, the Hamiltonian is block diagonal for spin up and spin down. Therefore, transmission coefficient for spin up and down could be treated independently. If an scattering region cause spin flip, the new Hamiltonian is twice bigger than Hamiltonian without spin and to calculate transmission coefficient, one need to consider whole Hamiltonian. For a para-magnetic system, one could obtain the total transmission by multiplying one spin transmission by 2. Figure 18 shows examples of ferromagnetic and anti-ferromagnetic systems.

##### H. Inclusion of a Gauge field

For a scattering region of area a, if a magnetic field B is applied the magnetic flux ϕ = B ×a. To compute transport properties in the presence of a magnetic field, a Peierls substitution is considered by changing the phase factors of the coupling elements between atomic orbitals. For example in the case of a nearest-neighbor tight-binding Hamiltonian, the hoping matrix element Hij between site i and site j is replaced with the modified element,

where

 (197)

and ri and rj are the positions of site i and j and A is the vector potential. The gauge should be chosen such that the principal layers of the leads remain translationally invariant after the substitution.

Fig. 19. Two-terminal device consists of two physical leads connected to a scattering region containing two superconductors with order parameters 1 and 2. The left (right) physical lead consists of two virtual leads p1 and h1 (p2 and h2) carrying particle and hole channels, respectively.

##### I. Superconducting systems

Figure 19 a shows a two-probe normal-superconductor-normal device with left and right normal reservoirs connected to a scattering region containing one or more superconductors. If the complete Hamiltonian describing a normal system is HN, then in the presence of superconductivity within the extended scattering region, the new system is described by the Bogoliubov-de Gennes Hamiltonian

 (198)

where the elements of the matrix Δ are non-zero only in the region occupied by a superconductor, as indicated in figure 19 b. Physically, HN describes particle degrees of freedom, -HN* describes hole degrees of freedom and Δ is the superconducting order parameter.

The multi-channel scattering theory for such a normal-superconducting-normal structure could be written as [29]:

 (199)

where IL (IR) is the current from the left (right) reservoir, μL -μ (μR -μ) is the difference between the chemical potential of the left (right) reservoir and the chemical potential μ of the superconducting condensate and the voltage difference between the left and right reservoirs is (μL - μR)∕e. In this equation,

 (200)

where NL (NR) is the number of open channels in the left (right) lead, Ro,To (Ra,Ta) are normal (Andreev) reflection and transmission coefficients for quasi-particles emitted from the right lead, Ro,To(Ra,Ta) are normal (Andreev) reflection and transmission coefficients from the left lead and all quantities are evaluated at the Fermi energy E = μ. As a consequence of unitarity of the scattering matrix, these satisfy Ro + To + Ra + Ta = NL and Ro+ To+ Ra+ Ta= NR.

The current-voltage relation of equation 199 is fundamentally different from that encountered for normal systems, because unitarity of the s-matrix does not imply that the sum of each row or column of the matrix a is zero. Consequently, the currents do not automatically depend solely on the applied voltage difference (μL - μR)∕e (or more generally on the differences between incoming quasi-particle distributions). In practice such a dependence arises only after the chemical potential of the superconductor adjusts itself self-consistently to ensure that the current from the left reservoir is equal to the current entering the right reservoir. Insisting that IL = -IR = I, the two-probe conductance G = I∕((μL - μR)∕e) takes the form of

 (201)

If a superconductor is disordered, then as its length L increases, all transmission coefficients vanish. Consequently, the superconductor possesses zero resistivity (equation 201 ). In this limit, the above equation reduces to (h∕2e2)G = 2∕Ra + 2∕Ra. In contrast with a normal scatterer, this shows that in the presence of Andreev scattering, as L tends to infinity, the resistance ( = 1/conductance) remains finite and therefore the resistivity (i.e. resistance per unit length) vanishes.

##### J. Environmental effects

Environmental effects such as presence of water, counter-ions, solvents, pH, molecules conformational change, dopants, nearby molecules, charge trap in a substrate can affect transport properties through a molecular junction. In order to model these effects, a statistical analysis needs to be carried out. Since a molecular junction in the presence of surrounding molecules is a dynamic object at room temperature, a molecular dynamics simulation is usually needed to understand the range of possible configurations of the system. A few configuration then should be extracted and full DFT calculations carried out to obtain the mean field Hamiltonian of the system in the presence of surrounding molecules. An alternative way to study the environmental effect is to create a series of more ideal configurations systematically in the presence of surrounding molecules (e.g. by moving surroundings in different directions). Then without geometry optimization, one could calculate the binding energy of surroundings to the molecule for each configuration and study transport properties of the confirmations with higher binding energies. Alternatively, Boltzmann distribution Tnew = Told × e-Ebinding∕k BT can be used to account for the weight of each confirmation on total ensemble transmission coefficient. Both of these methods are widely used in the literature to model environmental effects.

To calculate the binding energy between surroundings S and the backbone B, three total energy calculations need to be carried out. (1) The total energy of whole structure including surroundings and the backbone ESB, (2) The total energy of the surroundings in the presence of the backbone ghost states ESb, (3) The total energy of backbone in the presence of the surroundings ghost states EsB. The binding energy is then obtained by ESBbinding = ESB - (ESb + EsB). It is worth mentioning that since different effects such as physobrtion and charge transfer can play important role in these simulations, SCF methods need to be used to calculate transport from the mean-field DFT Hamiltonian. Recently, machine learning algorithms have been proposed to predict behavior of a junction in the absence/presence of surroundings. In these methods, machine is trained using a series of molecular junctions which are fully characterised using SCF methods. Machine then predicts the behavior of new junctions using its data base and pattern recognition algorithms.

#### V. Conclusion

To understand, predict and explain new phenomenon in the nanoscale junctions, one needs to study electron, phonon and spin transport through these junctions. The focus of this paper was on exploring theoretical methods to study electron, phonon and spin transports through molecules and low dimensional materials. Systematic study of such systems begins with understanding their vibrational and electronic structure. Then, one need to calculate the core quantities, transmission function of electrons Te and phonons Tph traversing through such systems from one electrode to the other. This is not only vital to understand their properties and related experiments but also to develop new concepts and design strategies for future applications. We mainly considered junctions where the transport is assumed to be elastic and coherent. In addition, we showed that how these methods could be extended to the incoherent and inelastic regimes. A good agreement between experiments and theories devolved using these methods [6, 9, 11] demonstrate that they are accurate enough to predict trends and develop new design strategies for future applications. However, in terms of predicting actual numbers, there are rooms for improvements. Apart from those quantities that depends on the actual shape of the contact (e.g. screening effects) and cannot be easily calculated; new computationally cheap approaches are needed to take into account electron-electron interactions more accurately. This would help to predict more accurate energy gap and level spacing leading to more realistic modells.

#### References

[1] “Itrs roadmap 2015,” ITRS, 2015.

[2] M. Chhowalla, D. Jena, and H. Zhang, “Two-dimensional semiconductors for transistors,” Nat. Rev. Mater., vol. 1, no. 11, p. 16052, 2016.

[3] “Visions for a molecular future,” Nature Nanotechnology, vol. 8, no. 6, pp. 385–389, 2013.

[4] J. L. Christopher, L. A. Estroff, J. K. Kriebel, R. G. Nuzzo, and G. M. Whitesides, “Self-assembled monolayers of thiolates on metals as a form of nanotechnology,” Chemical Reviews, vol. 105, no. 4, pp. 1103–1170, 2005. PMID: 15826011.

[5] S. V. Aradhya and L. Venkataraman, “Single-molecule junctions beyond electronic transport,” Nature Nanotechnology, vol. 8, no. 6, pp. 399–410, 2013.

[6] S. Sangtarash, C. Huang, H. Sadeghi, G. Sorohhov, J. Hauser, T. Wandlowski, W. Hong, S. Decurtins, S.-X. Liu, and C. J. Lambert, “Searching the Hearts of Graphene-like Molecules for Simplicity, Sensitivity, and Logic,” Journal of the American Chemical Society, vol. 137, no. 35, pp. 11425–11431, 2015.

[7] H. Sadeghi, L. Algaragholy, T. Pope, S. Bailey, D. Visontai, D. Manrique, J. Ferrer, V. Garcia-Suarez, S. Sangtarash, and C. J. Lambert, “Graphene sculpturene nanopores for DNA nucleobase sensing,” Journal of Physical Chemistry B, vol. 118, no. 24, pp. 6908–6914, 2014.

[8] T. Prodromakis, C. Toumazou, and L. Chua, “Two centuries of memristors,” Nature Materials, vol. 11, no. 6, pp. 478–481, 2012.

[9] H. Sadeghi, J. A. Mol, C. S. Lau, G. A. D. Briggs, J. Warner, and C. J. Lambert, “Conductance enlargement in picoscale electroburnt graphene nanojunctions,” Proceedings of the National Academy of Sciences, vol. 112, no. 9, pp. 2658–2663, 2015.

[10] H. Sadeghi, S. Sangtarash, and C. J. Lambert, “Oligoyne molecular junctions for efficient room temperature thermoelectric power generation,” Nano letters, vol. 15, no. 11, pp. 7467–7472, 2015.

[11] C. Evangeli, K. Gillemot, E. Leary, M. T. González, G. Rubio-Bollinger, C. J. Lambert, and N. Agraït, “Engineering the thermopower of C60 molecular junctions,” Nano Lett., vol. 13, no. 5, pp. 2141–2145, 2013.

[12] S. Datta, Quantum Transport : Atom to Transistor. Cambridge Univ Pr, 2005.

[13] “nanohub-u: Fundamentals of nanoelectronics - part a.” https://nanohub.org/courses/fon1.

[14] “nanohub-u: Fundamentals of nanoelectronics - part b.” https://nanohub.org/courses/fon2.

[15] R. Gebauer and R. Car, “Kinetic theory of quantum transport at the nanoscale,” Phys. Rev. B, vol. 70, p. 125324, Sep 2004.

[16] M. Büttiker, “Scattering theory of current and intensity noise correlations in conductors and wave guides,” Phys. Rev. B, vol. 46, pp. 12485–12507, Nov 1992.

[17] E. Bonet, M. M. Deshmukh, and D. C. Ralph, “Solving rate equations for electron tunneling via discrete quantum states,” Phys. Rev. B, vol. 65, p. 045317, Jan 2002.

[18] E. Schrödinger, “An Undulatory Theory of the Mechanics of Atoms and Molecules,” Phys. Rev., vol. 28, pp. 1049–1070, dec 1926.

[19] E. Runge and E. K. U. Gross, “Density-functional theory for time-dependent systems,” Phys. Rev. Lett., vol. 52, pp. 997–1000, Mar 1984.

[20] E. Engel and R. M. Dreizler, Density Functional Theory, vol. 2011. Springer Verlag, 2011.

[21] N. Harrison, “An introduction to density functional theory,” NATO Science Series Sub Series III: Computer and systems sciences, vol. 187, pp. 45–70, 2003.

[22] C. Lee, W. Yang, and R. G. Parr, “Development of the colle-salvetti correlation-energy formula into a functional of the electron density,” Phys. Rev. B, vol. 37, pp. 785–789, Jan 1988.

[23] J. Heyd, G. E. Scuseria, and M. Ernzerhof, “Hybrid functionals based on a screened coulomb potential,” The Journal of Chemical Physics, vol. 118, no. 18, pp. 8207–8215, 2003.

[24] M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, “Van der waals density functional for general geometries,” Phys. Rev. Lett., vol. 92, p. 246401, Jun 2004.

[25] J. M. J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, D. Sánchez-Portal, A. Garcia, J. Junquera, P. Ordejon, and D. Sanchez-Portal, “The SIESTA method for ab initio order-N materials simulation,” Journal of physics. Condensed matter :, vol. 2745, p. 22, mar 2002.

[26] H. Sadeghi, S. Sangtarash, and C. J. Lambert, “Electron and heat transport in porphyrin-based singlemolecule transistors with electro-burnt graphene electrodes,” Beilstein J. Nanotechnol., vol. 6, no. 1, pp. 1413–1420, 2015.

[27] H. Sadeghi, S. Sangtarash, and C. Lambert, “Robust Molecular Anchoring to Graphene Electrodes,” Nano Lett., vol. 17, no. 8, pp. 4611–4618, 2017.

[28] S. Sanvito, C. J. Lambert, J. H. Jefferson, and a. M. Bratkovsky, “General Green’s-function formalism for transport calculations with spd Hamiltonians and giant magnetoresistance in Co- and Ni-based magnetic multilayers,” Phys. Rev. B, vol. 59, no. 18, pp. 936–948, 1999.

[29] J. Ferrer, C. J. Lambert, V. M. García-Suárez, D. Z. Manrique, D. Visontai, L. Oroszlany, R. Rodríguez-Ferradás, I. Grace, S. Bailey, K. Gillemot, et al., “Gollum: a next-generation simulation tool for electron, thermal and spin transport,” New Journal of Physics, vol. 16, no. 9, p. 093029, 2014.

[30] A. E. Miroshnichenko, S. Flach, and Y. S. Kivshar, “Fano resonances in nanoscale structures,” Rev. Mod. Phys., vol. 82, pp. 2257–2298, Aug 2010.

[31] U. Fano, “Effects of configuration interaction on intensities and phase shifts,” Phys. Rev., vol. 124, pp. 1866–1878, Dec 1961.

[32] L. Oroszlny, A. Kormnyos, J. Koltai, J. Cserti, and C. J. Lambert, “Nonthermal broadening in the conductance of double quantum dot structures,” Physical Review B, vol. 76, p. 045318, 2007.

[33] S. Sangtarash, A. Vezzoli, H. Sadeghi, N. Ferri, H. M. O’Brien, I. Grace, L. Bouffier, S. J. Higgins, R. J. Nichols, and C. J. Lambert, “Gateway state-mediated, long-range tunnelling in molecular wires,” Nanoscale, vol. 10, pp. 3060–3067, 2018.

[34] H. Sadeghi, Theory of electron and phonon transport in nano and molecular quantum devices: design strategies for molecular electronics and thermoelectricity. PhD thesis, 2016.

[35] D. John and D. Pulfrey, “Green’s function calculations for semi-infinite carbon nanotubes,” physica status solidi (b), vol. 243, no. 2, pp. 442–448, 2006.

[36] R. Landauer, “Electrical transport in open and closed systems,” Z. Phys. B: Condens. Matter, vol. 68, no. 2-3, pp. 217–228, 1987.

[37] M. Buttiker, “Symmetry of electrical conduction,” IBM Journal of Research and Development, vol. 32, no. 3, pp. 317–334, 1988.

[38] M. Buttiker, “Coherent and sequential tunneling in series barriers,” IBM Journal of Research and Development, vol. 32, no. 1, pp. 63–75, 1988.

[39] M. Brandbyge, J.-L. Mozos, P. Ordejón, J. Taylor, and K. Stokbro, “Density-functional method for nonequilibrium electron transport,” Phys. Rev. B, vol. 65, p. 165401, Mar 2002.

[40] F. Furche, R. Ahlrichs, C. Httig, W. Klopper, M. Sierka, and F. Weigend, “Turbomole,” Wiley Interdisciplinary Reviews: Computational Molecular Science, vol. 4, no. 2, pp. 91–100.

[41] S. Sangtarash, H. Sadeghi, and C. J. Lambert, “Exploring quantum interference in heteroatom-substituted graphene-like molecules,” Nanoscale, vol. 8, no. 27, pp. 13199–13205, 2016.

[42] S. Sangtarash, H. Sadeghi, and C. J. Lambert, “Connectivity-driven bi-thermoelectricity in heteroatom-substituted molecular junctions,” Phys. Chem. Chem. Phys., vol. 20, no. 14, pp. 9630–9637, 2018.

[43] C. C. Escott, F. A. Zwanenburg, and A. Morello, “Resonant tunnelling features in quantum dots,” Nanotechnology, vol. 21, no. 27, p. 274018, 2010.

[44] J. Koch, F. von Oppen, F. V. Oppen, and F. von Oppen, “Franck-Condon Blockade and Giant Fano Factors in Transport through Single Molecules,” Physical Review Letters, vol. 94, no. 20, p. 206804, 2005.

[45] C. S. Lau, H. Sadeghi, G. Rogers, S. Sangtarash, P. Dallas, K. Porfyrakis, J. Warner, C. J. Lambert, G. A. D. Briggs, and J. A. Mol, “Redox-dependent franck–condon blockade and avalanche transport in a graphene–fullerene single-molecule transistor,” Nano letters, vol. 16, no. 1, pp. 170–176, 2015.

[46] E. Burzuri, Y. Yamamoto, M. Warnock, X. Zhong, K. Park, A. Cornia, and H. S. J. Van Der Zant, “Franck-condon blockade in a single-molecule transistor,” Nano Lett., vol. 14, no. 6, pp. 3191–3196, 2014.

[47] T. Frederiksen, M. Paulsson, M. Brandbyge, and A.-P. Jauho, “Inelastic transport theory from first principles: Methodology and application to nanoscale devices,” Physical Review B, vol. 75, no. 20, p. 205413, 2007.

[48] O. Stetsovych, P. Mutombo, M. vec, M. mal, J. Nejedl, I. Csaov, H. Vzquez, M. Moro-Lagares, J. Berger, J. Vacek, I. G. Star, I. Star, and P. Jelnek, “Large converse piezoelectric effect measured on a single molecule on a metallic surface,” Journal of the American Chemical Society, vol. 140, no. 3, pp. 940–946, 2018. PMID: 29275621.

[49] J. M. Seminario, “An introduction to density functional theory in chemistry,” Theoretical and Computational Chemistry, vol. 2, no. C, pp. 1–27, 1995.

[50] A. D. Becke, “A new mixing of hartree–fock and local density-functional theories,” The Journal of Chemical Physics, vol. 98, no. 2, pp. 1372–1377, 1993.

[51] L. Hedin, “New method for calculating the one-particle green’s function with application to the electron-gas problem,” Phys. Rev., vol. 139, pp. A796–A823, Aug 1965.

[52] M. Strange, C. Rostgaard, H. Häkkinen, and K. S. Thygesen, “Self-consistent gw calculations of electronic transport in thiol-and amine-linked molecular junctions,” Physical Review B, vol. 83, no. 11, p. 115108, 2011.

[53] T. Markussen, C. Jin, and K. S. Thygesen, “Quantitatively accurate calculations of conductance and thermopower of molecular junctions,” physica status solidi (b), vol. 250, no. 11, pp. 2394–2402, 2013.

[54] T. Kim, P. Darancet, J. R. Widawsky, M. Kotiuga, S. Y. Quek, J. B. Neaton, and L. Venkataraman, “Determination of energy level alignment and coupling strength in 4,4-bipyridine single-molecule junctions,” Nano Letters, vol. 14, no. 2, pp. 794–798, 2014. PMID: 24446585.

[55] R. Frisenda, S. Tarkuç, E. Galán, M. L. Perrin, R. Eelkema, F. C. Grozema, and H. S. van der Zant, “Electrical properties and mechanical stability of anchoring groups for single-molecule electronics,” Beilstein J. Nanotechnol., vol. 6, no. 1, pp. 1558–1567, 2015.

[56] X. Zhao, C. Huang, M. Gulcur, A. S. Batsanov, M. Baghernejad, W. Hong, M. R. Bryce, and T. Wandlowski, “Oligo (aryleneethynylene) s with terminal pyridyl groups: synthesis and length dependence of the tunneling-to-hopping transition of single-molecule conductances,” Chemistry of materials, vol. 25, no. 21, pp. 4340–4347, 2013.

[57] G. Sedghi, V. M. García-Suárez, L. J. Esdaile, H. L. Anderson, C. J. Lambert, S. Martín, D. Bethell, S. J. Higgins, M. Elliott, N. Bennett, et al., “Long-range electron tunnelling in oligo-porphyrin molecular wires,” Nature nanotechnology, vol. 6, no. 8, pp. 517–523, 2011.

[58] M. D. Yates, J. P. Golden, J. Roy, S. M. Strycharz-Glaven, S. Tsoi, J. S. Erickson, M. Y. El-Naggar, S. Calabrese Barton, and L. M. Tender, “Thermally activated long range electron transport in living biofilms,” Phys. Chem. Chem. Phys., vol. 17, pp. 32564–32570, 2015.

[59] X. Zhao, C. Huang, M. Gulcur, A. S. Batsanov, M. Baghernejad, W. Hong, M. R. Bryce, and T. Wandlowski, “Oligo(aryleneethynylene)s with terminal pyridyl groups: Synthesis and length dependence of the tunneling-to-hopping transition of single-molecule conductances,” Chemistry of Materials, vol. 25, no. 21, pp. 4340–4347, 2013.