Tutorial
Theory of Electron, Phonon and Spin Transport in Nanoscale Quantum Devices
Hatef Sadeghi
At the level of fundamental science, it was recently demonstrated that molecular wires can mediate longrange phasecoherent 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 graphenelike twodimensional 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 (T_{e}) and phonon transmission coefficient (T_{ph}) through them. The aim of this tutorial article is to outline the basic theoretical concepts and review the stateoftheart 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
Download this tutorial: PDF format
A molecular junction consists of Bipyridine molecule connected to gold electrodes. Bottom right: Electron transmission coefficient (T_{e}). Bottom left: Phonon transmission coefficient (T_{ph}) . Electron (phonon) transmission coefficient describes the transmission probability of electrons (phonons) with energy E (ℏω) from left to right electrodes through the molecule.
I. Introduction: Molecular electronics
The siliconbased 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 selfassembled monolayers [4] and singlemolecule junctions [5] are of interest not only for their potential to deliver logic gates [6], sensors [7], and memories [8] with ultralow power requirements and sub10nm device footprints, but also for their ability to probe roomtemperature quantum properties at the molecular scale such as quantum interference [9] and thermoelectricity [10, 11]. There are five main area of research in molecularscale 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 singlemolecule 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 waveguides 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 semiclassical methods [12–14], 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.
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 timedependent Schrödinger equation and introduce density functional theory (DFT) and a tight binding description of quantum system. The scattering theory and nonequilibrium Green’s function method are discussed and different transport regimes (on and off resonances) are considered. One dimensional system and a more general multichannel 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 FranckCondon blockade regimes. We follow with a discussion of the physical interpretation of quantum systems including charge, spin and thermal currents, phonon thermal conductance, electronphonon 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 nonrelativistic 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 timeindependent 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 timedependent 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 timedependent 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 timeindependent. Of course equation 5 is a particular solution of the timedependent 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 timeindependent 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 timedependent 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 BornOppenheimer 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 tightbinding Hamiltonian using Hückel parameters or use DFT to construct the material specific meanfield 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 timeindependent Schrödinger equation could be written as a matrix equation:

(8) 
where

(9) 
and

(10) 
The evaluation of these integrals is the most timeconsuming step, but once [H] and [S] are obtained, the eigenvalues E_{n} and eigenvectors ϕ_{n} are easily calculated. If ⟨i and j⟩ are orthogonal then S_{ij} = δ_{ij} where δ_{ij} is the Kronecker delta (δ_{ij} = 1 if i = j and δ_{ij} = 0 if i≠j). Note that a system with the Hamiltonian H and overlap matrix S obtained using nonorthogonal basis could be transformed to a new Hamiltonian = S^{1∕2} × H × S^{1∕2} with orthogonal basis ( = 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 BornOppenheimer 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 groundstate 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 BornOppenheimer 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 groundstate density and by obtaining the groundstate density, one can in principle calculate the groundstate 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 noninteracting particles in an effective external potential, which has the same groundstate density as the original system as illustrated in figure 2 a. The difference between the energy of the noninteracting and interacting system is referred to the exchange correlation functional (figure 2 a).
Fig. 2. From manybody problem to density functional theory DFT. (a) BornOppenheimer approximation, HohenbergKohn theorem and KohnSham ansatz, (b) Schematic of the DFT selfconsistency 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 semilocal. 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 KohnSham 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 HartreeFock. One of the latest and most universal functionals, the Van der Waals density functional (vdWDF [24]), contains nonlocal 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 planewave basis set is natural since it is, by itself, periodic. However, to construct a tightbinding 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 cutoff radius, and are constructed from the orbitals of the atoms.
Meanfield Hamiltonian from DFT: To obtain the ground state meanfield 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 exchangecorrelation 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 kgrid 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 selfconsistent 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 multigrid or fast Fouriertransform method. Then the KohnSham equations are solved and a new density is obtained. This selfconsistent 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 selfconsistent 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. TightBinding Model
By expanding wave function over a finite set of atomic orbitals, Hamiltonian of a system can be written in a tightbinding 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 planewave DFT code is used, a LCAOlike 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 c′ are the neighboring identical cells containing states α and

(12) 
and

(13) 
Equation 11 could be written as

(14) 
where

(15) 
and

(16) 
More generally, the singleparticle tightbinding Hamiltonian in the Hilbert space formed by R_{α}⟩ could be written as:

(17) 
where ε_{α} is the corresponding onsite 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 onsite 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 meanfield DFT or a simple HMO Hamiltonians (described in this section) were obtained, electron and spin transport properties of a junction can be calculated.
As a simplest example, consider a close system of two singleorbital sites with onsite 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 antibonding 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 onsite energies ⟨jHj⟩ = ε_{0} and the hopping matrix elements ⟨jHj ± 1⟩ = ⟨j ± 1Hj⟩ = γ could be written as:

(24) 
Therefore the Schrödinger equation reads

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

(26) 
and

(27) 
where k = a_{0} is the dimensionless wave vector and π∕a_{0} < < π∕a_{0} in the first Brillouin zone. Equation 27 is called a dispersion relation (E  k) or electronic bandstructure 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 E_{min} (figure 3 c), by Taylor expansion of equation 27 , a parabolic band structure is obtained: E(k) ≈ ε_{0}  2γ + γk^{2}. Comparing this with a free electron parabolic band structure E(k) = ε_{0}  2γ + ℏ^{2}∕2m^{*}k^{2}, γ = ℏ^{2}∕2m^{*} 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_{x},k_{y},k_{z} where ∑ _{kp} →∫ _{∞}^{+∞}dk_{p}∕2π. From the dispersion relation (equation 27 ), DOS of a 1D chain is obtained as: D(E) = dk∕πdE = 1∕2πγ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 springconstant 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ödingerlike equation could be written as:

(29) 
Similar to what was discussed above, using x_{n}(t) = Ae^{i(knωt)}, equation 29 reads mω^{2} = K[2  e^{ik}  e^{ik}] 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 offdiagonal terms of 1D chain TB Hamiltonian which make sense to satisfy translational invariance. Therefore, Schrödinger equationlike 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} = e^{ikj} + 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 IIBmbox 2 , the bandstructure and density of states of a 1D chain were calculated. Now let’s consider two most used 2D lattices: a square lattice where the unitcell 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 bandstructure could be calculated [12–14] using equation 17 and the Bloch wave function of the form Ae^{ikxj+ikyl} (figure 5 ). The Schrödinger equation for two dimensional square lattice (figure 5 a) with onsite energies ε_{0} and hopping integrals γ could be written as:

(37) 
Using the Bloch function ϕ_{j,l} = Ae^{ikxj+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 timedependent Schrödinger equation 1 ,

(40) 
By expanding ψ_{t}⟩ over orthogonal basis j⟩ equation 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 I_{l} = dρ_{t}^{l}∕dt 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 I_{l+1} → l and left moving electrons I_{l1} → l. The corresponding current to a Bloch state ψ_{j}(t) = e^{ikjiE(k)t∕ℏ} are:

(48) 
and

(49) 
where v_{k} = ∂E(k)∕ℏ∂k = 2γsin(k)∕ℏ is the group velocity. Note that this defines the time that a wavepackage need to move from one site to another e.g. l + 1 → l, so the actual group velocity is v_{k} × a, where a is the spacing between the sites. Although the individual currents are nonzero and proportional to the group velocity, the total current I = I_{l+1→l} + I_{l1→l} 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 e^{ikj} 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 e^{ikj}∕ in later derivations. Furthermore, one important consequence of equations 46 and 47 is that if ψ_{j} = Ae^{ikj} + 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^{2}sink B^{2}sink due to the Ae^{ikl} and Be^{ikl}.
Initial states are usually assumed to be stationary. However, if a nonstationary initial state was prepared in a closed (isolated) system such as a finite 1D chain of N atom, the charge density would be timedependent (oscillatory) and therefore, current could be defined. As an example, for a system of two atoms coupled to each other by γ, Hamiltonian reads . If the initial states are and which are nonstationary states, the final state is obtained which is not stationary. For such a closed system, the current could be obtained from equations 44 and 45 .
D. ParrPariserPople (PPP) Hamiltonian
Equation 17 described a noninteracting 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, c_{i,s}^{†} is the electron creation operator that inserts an electron in state i and c_{j,s} is the electron annihilation operator which takes an electron out of state i [12–14]. n_{i,s} = c_{i,s}^{†}c_{j,s} is the electron number operator with n_{i} = ∑ _{s}n_{i,s}, ε_{i} is the energy of the orbital relative to the chiral symmetry point and γ_{i,j} are hopping integrals. Electronelectron interactions are missing from the noninteracting Hamiltonian (equation 50 ). In order to take the Coulomb electronelectron interaction into account, the ParrPariserPople (PPP) model can be used to write the interacting tightbinding Hamiltonian as

(51) 
where U_{ii} and U_{ij} are the onsite and long range Coulomb interactions, respectively given by the Ohno parametrization,

(52) 
and for i≠j

(53) 
where d_{ij} is the distance between sites i and j and U_{0} is the interaction amplitude e.g. U_{0} 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 selfconsistent field (SCF) regime in which the thermal broadening k_{B}T and coupling Γ to electrodes are comparable to Coulomb energy U_{0}. The SCF method (single electron picture) implemented within nonequilibrium Green’s function (NEGF) could be used to describe transport in this regime as discussed in sections mbox IIIAmbox to mbox IIIDmbox .
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 roomtemperature experiments suggest applicability of this method. A simplified BreitWigner formula ( mbox IIIAmbox 5 ) derived from this method also could be used to model onresonance 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 onresonance properties of system.
(2) The Coulomb blockade (CB) regime in which Coulomb energy U_{0} is much higher than both thermal broadening k_{B}T and coupling Γ where the SCF method is not adequate and multielectron master equation should be used to describe properties of the system (figure 6 ) as discussed in section mbox IIIEmbox . 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 U_{0} is comparable to the larger of the thermal broadening k_{B}T 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. Onsite 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 e^{ikLj}∕ traveling from the left to right (figure 7 ).
If the wave function in the left and right leads and scattering region are ψ_{j} = e^{ikLj}∕ + re^{ikLj}∕, ϕ_{j} = te^{ikRj}∕ and f_{j}, 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 (bandstructure) in the left and right leads are obtained as:

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

(60) 
where s_{a} = α_{L}ψ_{0} and s_{b} = α_{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) 

(64) 
Since s_{a} = α_{L}(1 + r)∕ and s_{b} = α_{R}t∕, 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 0^{L} and 0^{R} (figure 7 ) and

(67) 
where Σ_{L,R} = α_{L,R}^{2}g_{L,R} are called selfenergies due to the left and right contacts. The Green’s function in the surface of a semiinfinite 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 = 0^{L} due to a source at site l = 0^{L} (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 = 0^{L} due to a source at site l = 0^{L} is g_{L,R} = e^{ikL,R}∕γ_{L,R} (equation 66 ). Assuming two identical leads (k_{L} = k_{R} = 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 = (g_{aa}α_{L} + g_{bb}α_{R})∕γ and B = α_{L}^{2}α_{R}^{2}(g_{aa}g_{bb}  g_{ab}g_{ba})∕γ^{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 multichannel leads are obtained from

(72) 
and

(73) 
t_{i,j} (r_{i,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} = Sψ_{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 wavefunctions 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 g_{jl} satisfy the Green’s function equation (E  H)g_{jl} = δ_{jl}. We note that g_{jl} can be written as

(80) 
where

(81) 
and since any wavefunction 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)g_{jl} = δ_{jl},

(83) 
The first and second terms are zero from equation 82. For j = l  1, the third term γg_{l1,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 wavefunction

(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. G_{ij},

(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 wavefunction evolution from one point to the other in a system [12–14, 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 semiinfinite 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 (G_{10}) or 0 and 0 (G_{00}); 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 noninteracting parts g and Hamiltonian that connects them h. As shown in figure 8 , using the surface Green’s functions of the decoupled two semiinfinite leads 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 incident with a barrier, the wave is transmitted with the amplitude of t () and reflected with the amplitude of r (). Using the surface Green’s function of the leads (g_{00} and g_{11}), 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 semiinfinite 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 onsite 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 IIBmbox 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 g_{oooree}^{benzene} = 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 g_{oe}^{benzene}≠0 which is called constructive quantum interference. This implies nonzero transmission between any oe sites.
The Green’s function of the ring of 6 sites with onsite energies ε = 0 and hopping integrals γ can also be obtained by substituting its wavefunction (equation 35 ) into equation 91 . The eigenenergies of such system are: E_{n} = 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 HOMOLUMO gap E = E_{HL} where E_{HL} = (E_{H}  E_{L})∕2 = 0 and E_{H} (E_{L}) 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 g_{1}^{ring} = 0, g_{2}^{ring} = 0, g_{3}^{ring} = 0 and g_{oooree}^{ring} = 0 (the destructive quantum interference). This is a radical behaviour since contribution from all pares of HOMO – LUMO, HOMO1 – LUMO+1 and HOMO2 – 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 nonzero values g_{oe}^{ring}≠0 (the constructive quantum interference).
Consider a molecule which possesses only a HOMO ψ^{H}(l) of energy E_{H} and a LUMO ψ^{L}(l) of energy E_{L}, 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} = (g_{lm}(E))^{2} a destructive interference feature occurs at an energy E given by g_{lm}(E) = 0, or equivalently

(99) 
If the energy E at which the destructive interference feature occurs lies within the HOMOLUMO gap, then E  E_{H} > 0 and E_{L}  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 HOMOLUMO 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 HOMOLUMO gap. In the most symmetric case, where ψ^{H}(l)ψ^{*H}(m) = ψ^{L}(l)ψ^{*L}(m), this yields E = (E_{L}  E_{H})∕2 and therefore the interference dip occurs at the middle of the HOMOLUMO gap. On the other hand, if ψ^{H}(l)ψ^{*H}(m) << ψ^{L}(l)ψ^{*L}(m), then E ≃ E_{H} 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 HOMO1 should also be considered.
It is apparent from equation 98 that manipulating antiresonances (e.g. due to the destructive quantum interference) is easier than resonances. To manipulate a resonance, redox 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 antiresonances.
Figure 9 shows two examples of polycyclic aromatic hydrocarbons, benzene and pyridine. Six molecular orbitals and corresponding eigenenergies due to six p_{z} 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  x_{0}) can be defined as the limit of a function which exhibits a sharp peak about x_{0} 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 ψ_{n}^{i}ψ_{n}^{j}^{*} = δ_{ij}, we can find the expression for the trace of the Green’s function when η → 0 as trace(G_{i}) = ∑ _{n} = π ∑ _{n}δ(E  λ_{n}). Therefore DOS is obtained,

(104) 
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 g_{ba} 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 y_{a,b} = f_{a,b}^{m} and f^{m}⟩ is 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 n≠m are much smaller than in n = m. This yields to onresonance 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} = y_{a,b}^{2}cos(k_{L,R}) are the real part of the selfenergies. Γ_{L,R} = y_{a,b}^{2}sin(k_{L,R}) are the imaginary part of the selfenergies (Σ = σ + 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, y_{a}^{2} and y_{b}^{2} 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 onresonance (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 antiresonance 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 antiresonates 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 selfenergies 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 = te^{iθ} (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 N_{y} whereas the lattice is infinite in the x direction. A normalized wave function and the bandstructure of such structure are calculated as ψ_{kx}^{m} = sin(mπl∕(N_{y} + 1))e^{ikxj} and E(k_{x}) = ε_{0}  2γcos(mπ∕(N_{y} + 1)) + 2γcos(k_{x}), respectively. Similar to the onedimensional case (section mbox IICmbox ), current is associated to each ψ_{kx}^{m} since every miniband corresponds to a Bloch state. ψ_{kx}^{m} 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 E_{F} is given by G(E_{F}) = 2e^{2}M(E_{F})∕h where M(E_{F}) is the number of open conduction channels at E_{F}. In a onedimensional 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 N_{y} atomic sites, the Green’s function in the sites l and j due to the source in the sites l′ and j′ could be written as:

(109) 
where k_{x}^{m} is longitudinal momentum and v_{x}^{m} = ∂E(k_{x}^{m})∕ℏ∂k_{x}^{m} is the group velocity of channel m. Equation 109 could be rewritten as:

(110) 
where ϕ_{l}^{m} = sin(l). g_{lj,l′j′} consists of the sum of all allowed longitudinal modes e^{ikxmj } weighted by the corresponding transverse components ϕ_{l}^{m}. Note that for a given E, k_{x}^{m} could be both real and imaginary. If k_{x}^{m} value of an eigenstate has no imaginary part Im(k_{x}^{m}) = 0, this state defined to be open, or propagating since a complex k_{x}^{m} will only occur if the wave is tunnelling, or decaying. The sign of the imaginary part of k_{x}^{m} (the group velocity) can be used to define the direction of a decaying wave (a propagating wave) as summarized in the table I .
Four classes of possible scattering channels
left  right  
Decaying:  Im(k_{x}^{m}) > 0  Im(k_{x}^{m}) < 0 
Propagating:  Im(k_{x}^{m}) = 0,v_{x}^{m} < 0  Im(k_{x}^{m}) = 0,v_{x}^{m} > 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 N_{y} 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 T_{nn′} = s_{n,n′}(E,H)^{2} between two scattering channels n,n′ of 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 crosssection, labelled L = 1,2,…. Consider two vector spaces A and B, spanned by a countable set of basis functions. In what follows, the subspace B represents the structure of interest and subspace 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. X_{L} denotes the surface of A connected to B. Subspace B includes some open subspaces connected to reservoirs shown by dotes.
The Hamiltonian is H = H_{A} + H_{B} + H_{J}, where H_{J} allows transitions between the subspaces. Since H_{J} 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 H_{J} = 0, G reduces to the Green’s function g of the decoupled system, where

(113) 
and Dyson’s equation G = g(1  H_{J}g)^{1} yields

(114) 
Rewriting the above result for G_{AA} in the form

(115) 
yields

(116) 

(117) 
and

(118) 
These demonstrate that once G_{BB} is known, all other quantities are determined. To obtain an expression for transmission coefficients, it is convenient to introduce a set of states {⟩}, which span the subspace A and write g_{A} = ∑ _{nm}⟩g_{nm}⟨. Since part of A consists of a number of ideal, straight, normal leads of constant crosssection, described by a real Hamiltonian, it is convenient to associate a subset of the states {⟩} with open channels of these leads. For these states, the notation ⟩ = n,x⟩ is introduced, where n is a discrete label identifying the lead, quasiparticle 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 ⟩,⟩ orthogonal to the open channels (the close channels), and

(120) 
is the Green’s function of the semiinfinite lead between any position point x and x′ in the transport direction terminated at x = x_{L} and vanishes at x = x_{L} + a [34]. k_{x}^{n} is the longitudinal wave vector of channel n. If the lead belonging to channel n terminates at x = x_{L}, then on the surface of the lead, the Green’s function g_{n}(x,x′) takes the form g_{n}(x_{L},x_{L}) = g_{n}, where g_{n} = a_{n} + ib_{n} with a_{n} real and b_{n} equal to π times the density of states per unit length of channel n. Moreover, if v_{n} is the group velocity for a wave packet traveling along channel n, then ℏv_{n} = 2b_{n}∕g_{n}^{2}. It is interesting to note that if x and x′ are positions located between x_{L} and some point x_{n},

(121) 
If x_{n} 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 H_{J}, then the transmission amplitude t and the transmission coefficient T from channel n′ to channel n (n≠n′) are

(122) 
and

(123) 
and since

(124) 
one obtains

(125) 
Let’s introduce eigenstates of H_{B}, satisfying H_{B}f_{ν}⟩ = ϵ_{ν}f_{ν}⟩ and write

(126) 
From the expression for G_{BB} 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 g_{nm} is zero, Img_{nm} = 0 and g_{nm} = g_{nm}^{*}. 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 (G_{BB}^{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 selfenergy matrices and g_{B} is the retarded Green’s function of subspace B when H_{J} = 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, G_{BB} describes a quantum structure with welldefined energy levels, shifted by the self energy σ′ arising from contact with closed channels. Clearly no quasiparticle 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 lifetime matrix Γ.
For a system with nonorthogonal 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 selfenergy. Furthermore, figure 11 shows a slightly different approach to calculate the transmission (reflection) amplitude t (r) in a two terminal system with nonorthogonal basis set [28].
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 semiinfinite 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 H_{0} connected to the first nearest neighbour with matrix elements H_{1}. 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 semiinfinite 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 IIIAmbox , equation 68 and the section mbox IIIBmbox .
In the recursive methods, the following nonlinear equation 138 is solved iteratively using different algorithms such as a fixedpoint iterative or Newton’s scheme.

(138) 
where g_{s} is the surface Green’s function and η is a small number to avoid divergence of the Green’s function. In the fixedpoint iterative scheme, equation 138 takes the following form:

(139) 
where n indicates the iteration number and g_{s}^{1} = ((E iη)I H_{0})^{1}. The convergence of the fixedpoint 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 FermiDirac distribution function f(E μ) = 1∕(1 + e^{(Eμ)∕kBT}) associated with the electrochemical potential μ, k_{B} 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) 
G_{0} is called the ”Conductance Quantum”. In other words, current associated with a single Bloch state v_{k}∕L and generated by the electrochemical potential gradient is I = e(v_{k}∕L)DΔμ where D = ∂n∕∂E = L∕hv_{k}. 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) LandauerButtiker formula for multiterminal structures
Conductance measurements are often performed using a fourprobe structure to minimize the contact resistance effect. Also multiprobe structures are widely used to describe the Halleffect or in sensing applications. Based on the Landauer approach for two terminal system, Buttiker [37] suggested a formula to model multiprobe currents for structures with multiple terminals as:

(144) 
where I_{i} is current at ith terminal and T_{ij} is the transmission probability from terminal j to i. In a multiterminal 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 N_{i} is number of open conduction channels in lead i. In the four probe structure, if probe 3 and 4 are the outer voltage probes (I_{3} = I_{4} = 0) and probe 1 and 2 are the inner current probes, the four probe conductance is

(146) 
2) Equilibrium vs. nonequilibrium IV
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 nonlinear 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 ▽^{2}U = ρ∕ε where U is potential profile due to charge distribution ρ and ε is permittivity which could vary spatially [12–14]. In the nonequilibrium 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 G_{diff}(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. Nonequilibrium Green’s function formalism
If ESψ⟩ = Hψ⟩ describes the properties of the closed system H with nonorthogonal 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 nonequilibrium condition could be written [12–14]:

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

(149) 
where G^{R} = [ES  H  Σ]^{1} is retarded Green’s function (G^{A} = [G^{R}]^{†}), Σ = Σ_{1} + Σ_{2} + Σ_{0} is sum of selfenergies due to the electrodes Σ_{1}, Σ_{2}, and surroundings Σ_{0} such as dephasing contact or inelastic scattering e.g. electronphonon 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 [12–14, 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 selfenergies 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 H_{1,2} are the coupling matrices between electrode 1 (2) and scattering region and g_{1,2} are the surface Green’s function of the electrodes, Σ_{1,2} = H_{1,2}^{†}g_{1,2}H_{1,2}. Furthermore, current could be calculated as

(150) 
where Γ_{1} = i(Σ_{1}  Σ_{1}^{†}) is the imaginary part of the selfenergy, G^{n} is the density matrix, ∑ _{1}^{in} = s_{1}s_{1}^{†}∕2π is related to the source s⟩ and A = G^{R}ΓG^{A} is the spectral function as shown in figure 12 . Note that the definition of selfenergy 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,2}^{in} = Γ_{1,2}f_{1,2}(E)). However, in dephasing contact, Σ_{0}^{in} is not described by any Fermi function and since inflow and outflow should be equal Trace[∑ _{0}^{in}A] = Trace[Γ_{0}G^{n}]. Figure 12 summarize the basic nonequilibrium 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 rewritten as [12–14]:

(151) 
where T_{ij}(E) = Trace[Γ_{i}(E)G^{R}(E)Γ_{j}(E)G^{A}(E)] is the transmission coefficient for electrons with energy E passing from lead i to lead j.
It is worth mentioning that in the past decade, several numerical implementations of the scattering theory and nonequilibrium Green’s functions approach such as Gollum [29], SMEAGOL [28], TransSIESTA [39] and TURBOMOLE [40] have been developed.
Transport through onelevel system: Consider two identical 1D leads with onsite 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)G^{R}(E)Γ_{R}(E)G^{A}(E) where the retarded Green’s function is G^{R}(E) = (E  ε_{1}  Σ)^{1}, the selfenergies Σ = Σ_{L} + Σ_{R} is obtained from Σ_{L} = α^{2}e^{ik}∕γ and Σ_{R} = β^{2}e^{ik}∕γ and the broadening due to the left and right leads are Γ_{L} = i(Σ_{L}  Σ_{L}^{†}) = 2α^{2}sin(k)∕γ and Γ_{R} = i(Σ_{R}  Σ_{R}^{†}) = 2β^{2}sin(k)∕γ. Note that for a onelevel system, this equation can be exactly rewritten in the form of BreitWigner formula (equation 106 ). However, for the two or more level systems, BWF is a good approximation only to describe onresonances transport and if the level spacing is less than level broadening (see section mbox IIIAmbox 5 ).
Mapping quantum to semiclassical model: Current in the semiclassical 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 ingoing S_{1}D_{0} and outgoing ν_{1}N currents where ν_{1} is coupling strength to the contact 1, S_{1} = ν_{1}f_{1} is source and D_{0} = ∫ D(E)dE is total density of state in the scattering region. Therefore, current in the contact 1 is written as:

(152) 
where f_{1} and f_{2} are the Fermi distribution in the contact 1 and 2, respectively.
Fig. 14. Semiclassical method to calculate current [12–14].
In order to build a physical intuition of NEGF equations (figure 12 ), we compare the semiclassical picture (figure 14 ) with the quantum model (figure 12 ),

(153) 
A and G^{n} 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(G^{n}(E))∕2π. Clearly, the semiclassical 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 onresonance 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 multielectron master equation. In the multielectron picture, the overall Nelectron system has different probabilities P_{α} of being in one of the 2^{N} possible states α. Furthermore all of these probabilities P_{α} must add up to one. The individual probabilities could be calculated under steadystate 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 multielectron master equation [12–14]. 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 k_{B} is Boltzmann’s constant and T is temperature.
Oneelectron energy levels represent differences between energy levels corresponding to states that differ by one electron. If E(N) is the energy associated with the Nelectron 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 spindegenerate one electron system with energy ε where there are only two possibilities, either the state is full 1⟩ or 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 f_{1}(E) and f_{2}(E) Fermi functions.
Two electron system: When the number of electrons N increases, the number of possibilities increases by 2^{N}. In two electron system, there are four possibilities (2^{2}), both states empty 00⟩ or full 11⟩ and either one of them full and the other empty (01⟩ and 10⟩).
Figure 16 shows how to calculate current in two electron system [12–14]. 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 multielectron system was considered. In this case, the number of calculations rapidly increases, therefore efficient numerical algorithms are needed to solve the multielectron master equation.
Coulomb and FranckCondon 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 electronvibron coupling leads to suppression of tunnel current at low bias called FranckCondon (FC) blockade regimes.
To describe transport in this regime, a minimal model (the AndersonHolstein 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, spindegenerate electronic level, and (3) one vibron taken into account within the harmonic approximation. In this case, the AndersonHolstein Hamiltonian reads H = H_{mol} + H_{leads} + H_{T} with

(157) 
describing the electronic and vibrational degrees of freedom of the molecule,

(158) 
the noninteracting leads, and

(159) 
the tunnelling between the leads and molecule. Here, Coulomb blockade is taken into account via the charging energy U where eV,k_{B}T << U. The operator d_{σ} (d_{σ}^{†}) annihilates (creates) an electron with spin projection σ on the molecule, n_{d} = ∑ _{σ}d_{σ}d_{σ}^{†} denotes the corresponding occupationnumber operator. Similarly, c_{apσ} (c_{apσ}^{†}) 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 ~ n_{d}(b^{†} + b), which can be eliminated by a canonical transformation, leading to a renormalisation of the parameters ε and U, and of the leadmolecule coupling t_{a} → t_{a}e^{λ(b†+b) }. The master equations determining the molecular occupation probabilities P_{q}^{n} for charge state n and vibrons q is:

(160) 
P_{q}^{eq} denotes the equilibrium vibron distribution with a relaxation time τ and W_{q→q′}^{n→n′} denotes the total rate for a transition from n,q⟩ to n′,q′⟩.

(161) 
where f_{a} is the Fermi function and the transition rates Γ are calculated from Fermi’s golden rule.

(162) 
Here, ρ_{a} denotes DOS in lead a, M_{q→q′;a}^{n→n±1} denotes the FC matrix elements and s^{n→m} the spin factor [45] such that for sequential tunnelling and assuming twofold degeneracy they are s^{1→0} = s^{1→2} = 1,s^{0→1} = s^{2→1} = 2. The matrix elements M_{q→q′;a}^{n→n±1} defined for vibrations are

(163) 
where q_{1} = min{q,q′} and q_{2} = 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 phasecoherent 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 selfinteractions. The bias and gate voltage depended differential conductance is calculated by differentiating biasdependent current with respect to the bias voltage G_{diff} = 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 = (T_{1} + T_{2})∕2 where T_{1} (T_{2}) 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 FermiDirac distribution function f(E). Therefore, a simplified and approximated form of G_{diff} 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.
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 wavefunctions. 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 realtime 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 wellknown physical phenomena and the methods to model them theoretically.
A. Virtual leads versus physical leads
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 k_{i}, _{i}, 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 wavevectors k_{αi}, _{αi} as virtual leads and treating them on the same footing as physical leads [29].
This viewpoint is particularly useful when the Hamiltonians H_{0}^{i}, H_{1}^{i} describing the principle layers PLs (the identical periodic unit cells H_{0}^{i} connected to each other by H_{1}^{i}) of the physical lead i are block diagonal with respect to the quantum numbers associated with k_{αi}, _{α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, H_{0}^{i}, H_{1}^{i} are block diagonal and it is convenient to identify virtual leads α_{i} with each block, because we can compute the channels k_{αi}, _{α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 H_{0}^{i}, H_{1}^{i} 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 quasiparticles of type α_{i} is Q_{αi}(E), then the flux of Q from reservoir i is:

(166) 
T_{αi,βj}^{i,j} in this expression is transmission coefficient of quasiparticles 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 Q_{p} carried by particles is e, whereas the charge Q_{h} carried by holes is +e. In the presence of noncollinear 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 multiterminal junction having collinear magnetic moments, α_{i} = σ for all i and since there is no spinflip 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 _{σ}^{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 _{σ}^{j}(E) can be approximated by differentials, therefore by Taylor expanding f_{j} = f(E  μ_{j}), f_{j}  f = (μ_{j}  μ) ()(T_{j}  T) is obtained. Using this expression for f_{j}  f, equations 169 and 170 could be written as:

(171) 
Since n^{th} moment of probability distribution P(x) define as: ⟨x^{n}⟩ = ∫ dxP(x)x^{n}, the spindependent moment L_{ij,σ} in the presence of collinear magnetism is obtained as

(172) 
where f(E,T) = (1 + e^{(EEF)∕kBT})^{1} is the FermiDirac distribution function, T is temperature and k_{B} is Boltzmanns constant. Therefore, in the linearresponse 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 spindependent lowvoltage electrical conductance G(T,E_{F}), the Seebeck coefficient (somtimes called thermopower) S(T,E_{F}) = ΔV∕ΔT , the Peltier coefficient Π(T,E_{F}) and the thermal conductance due to the electrons κ_{el}(T,E_{F}) as a function of Fermi energy E_{F} 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 E_{F}, T_{σ,σ}^{ij}(E) ≈ T_{σ,σ}^{ij}(E_{F}) + _{E=EF}(E E_{F}), the low temperature approximation of L_{ij,σ}^{1} and L_{ij,σ}^{0} where T_{σ,σ}^{ij}(E) varies approximately linearly with E on the scale of k_{B}T could be written as:

(175) 
and

(176) 
and

(177) 
where α = k_{B}^{2}π^{2}∕3e^{2} WΩK^{2} is the Lorenz number. Therefore, Seebeck coefficient could be rewritten 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,E_{F}) ≈ G_{0}T(E_{F}) and κ_{el}(T,E_{F}) ≈ αTG(T,E_{F}). The latter is also called the WiedemannFranz 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 T_{h} and T_{c} are temperatures of hot and coldsides, respectively, ΔT = T_{h} T_{c} and T_{avg} = (T_{h} + T_{c})∕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 = S^{2}G∕κ, 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.T_{avg}) 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 twoterminal 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 T_{ph} for different vibrational modes as

(182) 
where f_{BE}(ω,T) = (e^{ℏω∕kBT}  1)^{1} is BoseEinstein distribution function and ℏ is reduced Plancks constant and k_{B} 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 δq′ and δq′ in x, y and z directions and the forces F_{i}^{q} = (F_{i}^{x},F_{i}^{y},F_{i}^{z}) in three directions q_{i} = (x_{i},y_{i},z_{i}) 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 K_{ij}^{qq′} for i≠j 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 K_{ii} = ∑ _{i≠j}K_{ij}. Once the dynamical matrix is constructed the Green’s function method as described in section mbox IIIBmbox could be used to calculate the phononic transmission coefficient T_{ph}.
It is worth to mention that in order to compute the electronphonon coupling matrices W^{λ} for a given phononic mode λ, one could use the Hamiltonian matrices {{⟨iHj⟩}} for the displaced configurations [47].

(185) 
where i⟩ are basis orbitals. The W^{λ} matrices are then used [47] to calculate the appropriate selfenergies (see section mbox IIIDmbox ) to account for inelastic scattering due to electronphonon 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 r_{ij} 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 IVCmbox ). = V ^{T}H 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].
E. Spectral adjustment
Although DFT has been successful at predicting the trends, it usually underestimates the position of the Fermi energy E_{F}, the exact energy levels (KohnSham eigenvalues [49]) and therefore the position of the HOMO and LUMO and the energy gap. To compare meanfield 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 HOMOLUMO gap using the values measured experimentally. A phenomenological scheme that improves the agreement between theoretical simulations and experiments in, for example, singlemolecule 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 are energy shifts. is the density matrix where (no, nu) denote the occupied and unoccupied states, respectively and S_{M} is overlap matrix. If experimental HOMO and LUMO energies are available, _{}can be chosen to correct HOMO and LUMO obtained from the mean field Hamiltonian. Alternatively, in the simplest case, the shifts _{}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
(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 energygap E_{g} of a molecule (sometimes called additional energy) could be calculated from IP and EA as: E_{g} = IP  EA [12–14]. 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 = e^{2} ln2∕(8πϵ_{0}a) where a is the distance between the image plane and the point image charge and ϵ_{0} = 8.85 × 10^{12}F∕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 G_{c} is a preexponential 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 G_{a} is a constant preexponential factor for each chemical reaction, E_{A} is the hopping activation energy, T is the temperature and k_{B} 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 spinup and spindown, respectively. In the paramagnetic case, the site energies for spinup and spindown are equal ε_{0}. In the presence of spinorbit coupling, there are coupling between spinup and spindown sites.
G. Spin Hamiltonian
To take electrons spin into account in a Hamiltonian H (Eψ = Hψ), in the absence of spinorbit coupling, the Schrödinger equation can be written,

(192) 
where ψ and are spinup and spindown 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 spinorbit correction need to be considered using the Dirac equation.

(193) 
Spinorbit Hamiltonian is often written as σ.M where σ is Pauli spin matrices.

(194) 
If M points in the direction (θ,ϕ), M_{x} = sinθcosϕ, M_{y} = sinθsinϕ and M_{z} = cosθ, therefore

(195) 
The spin wavefunction, spinors thus is obtained

. If there is no spinorbit 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 paramagnetic system, one could obtain the total transmission by multiplying one spin transmission by 2. Figure 18 shows examples of ferromagnetic and antiferromagnetic 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 nearestneighbor tightbinding Hamiltonian, the hoping matrix element H_{ij} between site i and site j is replaced with the modified element,
where

(197) 
and r_{i} and r_{j} 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. Twoterminal 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 p_{1} and h_{1} (p_{2} and h_{2}) carrying particle and hole channels, respectively.
I. Superconducting systems
Figure 19 a shows a twoprobe normalsuperconductornormal 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 H_{N}, then in the presence of superconductivity within the extended scattering region, the new system is described by the Bogoliubovde Gennes Hamiltonian

(198) 
where the elements of the matrix Δ are nonzero only in the region occupied by a superconductor, as indicated in figure 19 b. Physically, H_{N} describes particle degrees of freedom, H_{N}^{*} describes hole degrees of freedom and Δ is the superconducting order parameter.
The multichannel scattering theory for such a normalsuperconductingnormal structure could be written as [29]:

(199) 
where I_{L} (I_{R}) 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 N_{L} (N_{R}) is the number of open channels in the left (right) lead, R_{o},T_{o} (R_{a},T_{a}) are normal (Andreev) reflection and transmission coefficients for quasiparticles emitted from the right lead, R_{o}′,T_{o}′ (R_{a}′,T_{a}′) 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 R_{o} + T_{o} + R_{a} + T_{a} = N_{L} and R_{o}′ + T_{o}′ + R_{a}′ + T_{a}′ = N_{R}.
The currentvoltage relation of equation 199 is fundamentally different from that encountered for normal systems, because unitarity of the smatrix 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 quasiparticle distributions). In practice such a dependence arises only after the chemical potential of the superconductor adjusts itself selfconsistently to ensure that the current from the left reservoir is equal to the current entering the right reservoir. Insisting that I_{L} = I_{R} = I, the twoprobe 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∕2e^{2})G = 2∕R_{a} + 2∕R_{a}′. 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, counterions, 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 T_{new} = T_{old} × 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 E_{SB}, (2) The total energy of the surroundings in the presence of the backbone ghost states E_{Sb}, (3) The total energy of backbone in the presence of the surroundings ghost states E_{sB}. The binding energy is then obtained by E_{SB}^{binding} = E_{SB}  (E_{Sb} + E_{sB}). 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 meanfield 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 T_{e} and phonons T_{ph} 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 electronelectron 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, “Twodimensional 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, “Selfassembled 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, “Singlemolecule 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 Graphenelike 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. GarciaSuarez, 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. RubioBollinger, 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] “nanohubu: Fundamentals of nanoelectronics  part a.” https://nanohub.org/courses/fon1.
[14] “nanohubu: 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, “Densityfunctional theory for timedependent 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 collesalvetti correlationenergy 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ánchezPortal, A. Garcia, J. Junquera, P. Ordejon, and D. SanchezPortal, “The SIESTA method for ab initio orderN 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 porphyrinbased singlemolecule transistors with electroburnt 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’sfunction formalism for transport calculations with spd Hamiltonians and giant magnetoresistance in Co and Nibased magnetic multilayers,” Phys. Rev. B, vol. 59, no. 18, pp. 936–948, 1999.
[29] J. Ferrer, C. J. Lambert, V. M. GarcíaSuárez, D. Z. Manrique, D. Visontai, L. Oroszlany, R. RodríguezFerradás, I. Grace, S. Bailey, K. Gillemot, et al., “Gollum: a nextgeneration 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 statemediated, longrange 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 semiinfinite 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. 23, 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, “Densityfunctional 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 heteroatomsubstituted graphenelike molecules,” Nanoscale, vol. 8, no. 27, pp. 13199–13205, 2016.
[42] S. Sangtarash, H. Sadeghi, and C. J. Lambert, “Connectivitydriven bithermoelectricity in heteroatomsubstituted 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, “FranckCondon 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, “Redoxdependent franck–condon blockade and avalanche transport in a graphene–fullerene singlemolecule 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, “Franckcondon blockade in a singlemolecule 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. MoroLagares, 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 densityfunctional theories,” The Journal of Chemical Physics, vol. 98, no. 2, pp. 1372–1377, 1993.
[51] L. Hedin, “New method for calculating the oneparticle green’s function with application to the electrongas problem,” Phys. Rev., vol. 139, pp. A796–A823, Aug 1965.
[52] M. Strange, C. Rostgaard, H. Häkkinen, and K. S. Thygesen, “Selfconsistent gw calculations of electronic transport in thioland aminelinked 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,4bipyridine singlemolecule 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 singlemolecule 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 tunnelingtohopping transition of singlemolecule conductances,” Chemistry of materials, vol. 25, no. 21, pp. 4340–4347, 2013.
[57] G. Sedghi, V. M. GarcíaSuárez, L. J. Esdaile, H. L. Anderson, C. J. Lambert, S. Martín, D. Bethell, S. J. Higgins, M. Elliott, N. Bennett, et al., “Longrange electron tunnelling in oligoporphyrin molecular wires,” Nature nanotechnology, vol. 6, no. 8, pp. 517–523, 2011.
[58] M. D. Yates, J. P. Golden, J. Roy, S. M. StrycharzGlaven, S. Tsoi, J. S. Erickson, M. Y. ElNaggar, 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 tunnelingtohopping transition of singlemolecule conductances,” Chemistry of Materials, vol. 25, no. 21, pp. 4340–4347, 2013.