NEGF formalism

Device configuration

This section describes the most important aspects of the formalism in QuantumATK for simulating device systems using the non-equilibrium Green’s function (NEGF) method. For a more detailed description of the theory you may consult some of the many research papers on the subject [BMO+02], [STBG05], [HJ08], [Dat97].

A typical device configuration is illustrated in Fig. 107.

../_images/twoprobe_algorithm_1.png

Fig. 107 Geometry of a device configuration with two electrodes (two-probe configuration). In the left and right electrode regions the system is periodic in the transport direction (the C direction). The central region seamlessly connects the left and right electrode regions, and the part of the central region where the atom positions follow the periodic arrangement of the electrodes is called the left and right electrode extension. The background color symbolizes the perturbation of the electrode effective potential relative to the electrode self-consistent value. A sufficient fraction of the electrodes must be included in the central region to screen out the perturbation of the scattering region, such that the outermost part of the central region attains an effective potential similar to the bulk potential of the electrodes. The matching between these potentials forms the boundary condition for the electrostatic problem.

The system can be divided into three regions: left, central, right. The implementation relies on the so-called screening approximation, which assumes that the properties of the left and right regions, the electrodes, can be described by solving a bulk problem for the fully periodic electrode cell. The screening approximation will formally be fulfilled when the current through the system is sufficiently small that the electrode regions can be described by an equilibrium electron distribution.

Fig. 107 illustrates how the electrode regions must be extended into the central region, in order to screen out the perturbations from the scatterer (in this case the benzene molecule) inside the device. For metallic electrodes it is usually sufficient to extend the electrodes 5-10 Å into the central region.

Non-equilibrium electron distribution

The left and right regions are equilibrium systems with periodic boundary conditions, and the properties of these systems are obtained using a conventional electronic structure calculation. The challenge in calculating the properties of a device system lies in the calculation of the non-equilibrium electron distribution in the central region.

The assumption is that the system is in a steady state such that the electron density of the central region is constant in time. The electron density is given by the occupied eigenstates of the system. Since the chemical potential is different in the two electrodes, the contribution from each electrode to the total electron density in the central region must be calculated independently:

\[n({\bf r}) = n^L({\bf r}) + n^R({\bf r}).\]

The contribution from the left (\(n^L\)) and right (\(n^R\)) electrodes can be obtained by calculating the scattering states, which are the eigenstates of the system when scattering boundary conditions are used. Fig. 108 illustrates a left moving scattering state with origin in the right electrode.

../_images/electdist.png

Fig. 108 The electron distribution in a device configuration. The left and right regions have an equilibrium electron distribution with chemical potentials \(\mu_L\) and \(\mu_R\), related through the applied sample bias, \(\mu_R - \mu_L = e V_b\) . The electrons with energies in the bias window \(\mu_L \leq \varepsilon \leq \mu_R\) , give rise to a steady state electrical current. The figure illustrates a left moving scattering state with origin in the right electrode.

The left and right density are now calculated by summing up the occupied scattering states,

\[n^L({\bf r}) = \sum_{\alpha}|\psi_\alpha({\bf r})|^2 f\left(\frac{\varepsilon_\alpha-\mu_L}{k T}\right),\]
\[n^R({\bf r}) = \sum_{\alpha}|\psi_\alpha({\bf r})|^2 f\left(\frac{\varepsilon_\alpha-\mu_R}{k T}\right).\]

The scattering states of the system are calculated by first calculating the Bloch states in the electrodes and subsequently solving the Schrödinger equation of the central region using the Bloch states as matching boundary conditions.

The non-equilibrium Green’s function (NEGF) method

Instead of using the scattering states to calculate the non-equilibrium electron density, QuantumATK uses the non-equilibrium Green’s function (NEGF) method; the two approaches are formally equivalent and give identical results [BMO+02].

The electron density is given in terms of the electron density matrix, as described in section Electron density. We divide the density matrix into left and right contributions,

\[D = D^L + D^R.\]

The left density matrix contribution is calculated using the NEGF method as [BMO+02]

\[D^L = \int \rho^L(\varepsilon) f\left(\frac{\varepsilon-\mu_L}{k_B T_L}\right) d\varepsilon ,\]

where

\[\rho^L(\varepsilon) \equiv \frac{1}{2\pi} G(\varepsilon) \Gamma^L(\varepsilon) G^\dagger(\varepsilon)\]

is the spectral density matrix. Note that while there is a non-equilibrium electron distribution in the central region, the electron distribution in the electrode is described by a Fermi function \(f\) with an electron temperature \(T_L\) .

In this equation, \(G\) is the retarded Green’s function, and

\[\Gamma^L = \frac{1}{i}(\Sigma^L - (\Sigma^L)^\dagger)\]

is the broadening function of the left electrode, given in terms of the left electrode self energy, \(\Sigma^L\) . A similar equation exists for the right density matrix contribution.

The following section describes the calculation of \(G\) and \(\Sigma\) in more detail.

Retarded Green’s function

The key quantity to calculate is the retarded Green’s function matrix,

\[G(\varepsilon ) = \frac{1}{(\varepsilon +i\,\delta_+)S -H} ,\]

where \(\delta_+\) is an infinitesimal positive number. \(S\) and \(H\) are the overlap and Hamiltonian matrices, respectively, of the entire system.

The Green’s function is only required for the central region and can be calculated from the Hamiltonian of the central region by adding the electrode self energies:

\[G(\varepsilon) =[ (\varepsilon+i \delta_+) S -H - \Sigma^L(\varepsilon) -\Sigma^R(\varepsilon) ]^{-1}.\]

Therefore, the calculation of the Green’s function of the central region at a specific energy requires the inversion of the Hamiltonian matrix of the central region. In QuantumATK this matrix is stored in a sparse format, and the inversion is done using an \(\mathcal{O}(N)\) algorithm [PSrensenH+08].

Self energy

The self energy describe the effect of the electrode states on the electronic structure of the central region. The self energy can be calculated from the electrode Hamiltonian. QuantumATK provides 4 different methods for calculating the self energy:

Complex contour integration

The energy integral needed to obtain the density matrix within the NEGF framework is evaluated through a complex contour integration. It is divided into two parts; an integral over equilibrium states and an integral over non-equilibrium states.

The integral over the equilibrium states can be done using two different methods:

  • SemiCircleContour: Defines a semi-circular complex contour that gives high computational efficiency [BMO+02].
  • OzakiContour: Evaluation of the complex contour integral is based on the residue theorem and a continued-fraction representation of the Fermi–Dirac districution [Oza07][ONK10]. This is a highly stable method, but is also less efficient than SemiCircleContour.

The integral over the non-equilibrium states must be performed along the real axis, and uses the RealAxisContour method. Note that for large biases, this is often the most demanding part of the calculation.

Lastly, one of two different contour methods must be chosen (usually just left at the default value):

  • SingleContour uses a single contour for the calculation. This is appropriate for small biases, and is the fastest method.
  • DoubleContour uses a double contour for the calculation. This gives best stability and must be used for high biases. The method can also handle if there are bound states inside the bias window.

See also ContourParameters.

Spill-in terms

In terms of the density matrix, \(D\) , the electron density of the central region is given by

\[n({\bf r}) = \sum_{ij} D_{ij} \phi_i({\bf r}) \phi_j({\bf r}).\]

The Green’s function of the central region gives the density matrix of the central region, \(D_{CC}\) , however, to calculate the density correctly close to the central cell boundaries the terms involving, \(D_{LL},D_{LC},D_{CR},D_{RR}\) are also needed. These terms are denoted spill-in terms [SMB+16].

ATK includes all the spill-in terms, both for calculating the electron density \(n ({\bf r})\) and the Hamiltonian integral. This gives additional stability and well-behaved convergence in the device algorithm [SMB+16].

Effective potential

Once the non-equilibrium density is obtained the next step in the self-consistent calculation is the calculation of the effective potential as the sum of the exchange-correlation and electrostatic Hartree potentials. The calculation of the exchange-correlation potential is straightforward since it is a local or semi-local function of the density. However, the calculation of the electrostatic Hartree potential requires some additional consideration for a device system. The following describes the calculation of the Hartree potential in more detail.

The starting point is the calculation of the self-consistent Hartree potential in the left and right electrodes. The Hartree potential of a bulk system is defined up to an arbitrary constant. However, in a device setup the Hartree potentials of the two electrodes are aligned through their chemical potentials (i.e. their Fermi levels), since these are related by the applied bias:

\[\mu_R - \mu_L = e V_b.\]

The Hartree potential of the central region is obtained by solving the Poisson equation, using the bulk-like Hartree potentials of the electrodes as boundary conditions at the interfaces between the electrodes and the central region. QuantumATK offers 5 different methods for solving the Poisson equation of a device system:

  • FastFourier2DSolver uses a Fourier transform in the directions perpendicular to the transport direction, and a real space multigrid method in the transport direction. This is a fast and accurate method and it is the default in QuantumATK.
  • MultigridSolver uses a real space multigrid method in all directions. This method is slower and slightly less accurate than the FastFourier2DSolver. However, it is very flexible and allows for different types of boundary conditions, implicit solvents, and the use of metallic and dielectric spatial regions to simulate gates.
  • DirectSolver uses a direct real space solver method in all directions. This method is both accurate and fast, but requires significantly more memory than the other solver methods. Similar to the MultigridSolver, it is very flexible and allows for different types of boundary conditions, implicit solvents, and the use of metallic and dielectric spatial regions to simulate gates.
  • ParallelConjugateGradientSolver uses an iterative solver based on the conjugate gradient method. Similar to the MultigridSolver, ìt allows for different types of boundary conditions, implicit solvents, and the use of metallic and dielectric spatial regions to simulate gates. This method is similar to MultigridSolver in terms of accuracy and outperforms it when running calculations in parallel, while requiring much less memory than DirectSolver.
  • FastFourierSolver uses a Fourier transform in all directions. This is a very fast method and is still used in some other implementations of the NEGF method. It requires that the electrodes are identical, but even for such so-called homogeneous systems the method can in some cases be problematic, especially at finite bias, and it is recommended not to use this method for device configurations.

Total energy and forces

A device system is an open system where charge can flow in and out of the central region from the left and right reservoirs. Since the particle number is not conserved it is necessary to use a grand canonical potential to describe the energetics of the system [Tod98]:

\[\Omega[n] = \varepsilon[n]- e \, N_L \, \mu_L - e \, N_R \, \mu_R,\]

where \(N_{L/R}\) is the number of electrons contributed to the central region from the left/right electrode.

Due to the screening approximation, the central region will be charge neutral, and therefore \(N_L + N_R = N\) , where \(N\) is the ionic charge in the central region. Thus, at no applied bias, we ahve \(\mu_L = \mu_R\) and the particle terms in \(\Omega\) will be constant when atoms are moved in the central region. However, at finite bias, \(\mu_L \neq \mu_R\), and the particle terms in \(\Omega\) will be important.

The forces are given by

\[{\bf F}_i = - \frac{\partial \Omega[n]}{\partial \bf R_i}.\]

It can be shown that the calculation of this force is identical to the calculation of the equilibrium force, as described in Total energy and forces. In the non-equilibrium case it is however required that the density and energy density matrix is calculated with the NEGF framework [BMO+02], [Tod98], [ZRSH11].

See also TotalEnergy.

Transmission coefficient

When the self-consistent non-equilibrium density matrix has been obtained, it is possible to calculate various transport properties of the system. One of the most notable is the TransmissionSpectrum from which you can obtain the current and differential conductance.

By the transmission amplitude \(t_k\) we define the fraction of a scattering state \(k\) which propagates through a device. The transmission coefficient at energy \(\varepsilon\) is obtained by summing up the transmission amplitude from all the states at this energy,

\[T(\varepsilon) = \sum_k t_k^\dagger t_k \delta(\varepsilon-\varepsilon_k).\]

The transmission coefficient may also be obtained from the retarded Green’s function using

\[T(\varepsilon) = G(\varepsilon) \Gamma^L(\varepsilon) G^\dagger (\varepsilon) \Gamma^R(\varepsilon),\]

and this is how it is calculated in QuantumATK. The transmission amplitude of individual scattering states may be obtained through the TransmissionEigenvalues.

Electrical current

To calculate the current you must first calculate a TransmissionSpectrum, and then from this extract the electrical current using the current() method on the object. This approach has the advantage that once the TransmissionSpectrum is calculated, it is fast to calculate the current for different electrode temperatures [SMB+16].

See TransmissionSpectrum for more details.

References

[BMO+02](1, 2, 3, 4, 5) M. Brandbyge, J.-L. Mozos, P. Ordejón, J. Taylor, and K. Stokbro. Density-functional method for nonequilibrium electron transport. Phys. Rev. B, 65:165401, Mar 2002. doi:10.1103/PhysRevB.65.165401.
[Dat97]S. Datta. Electronic Transport in Mesoscopic Systems. Cambridge University Press, 1997. Cambridge Studies in Semiconductor Physics and Microelectronic Engineering.
[HJ08]H. Haug and A.-P. Jauho. Quantum Kinetics in Transport and Optics of Semiconductors. volume 123. Springer-Verlag Berlin Heidelberg, 2 edition, 2008. Springer Series in Solid-State Sciences. doi:10.1007/978-3-540-73564-9.
[LSLSR85]M. P. Lopez Sancho, J. M. Lopez Sancho, and J. Rubio. Highly convergent schemes for the calculation of bulk and surface Green functions. J. Phys. F: Metal Physics, 15(4):851, 1985. URL: http://stacks.iop.org/0305-4608/15/i=4/a=009.
[ONK10]T. Ozaki, K. Nishio, and H. Kino. Efficient implementation of the nonequilibrium green function method for electronic transport calculations. Phys. Rev. B, 81:035116, 2010. doi:10.1103/PhysRevB.81.035116.
[Oza07]Taisuke Ozaki. Continued fraction representation of the fermi-dirac function for large-scale electronic structure calculations. Phys. Rev. B, 75:035123, 2007. doi:10.1103/PhysRevB.75.035123.
[PSrensenH+08]D. E. Petersen, H. H. B. Sørensen, P. C. Hansen, S. Skelboe, and K. Stokbro. Block tridiagonal matrix inversion and fast transmission calculations. J. Comput. Phys., 227(6):3174–3190, 2008. doi:10.1016/j.jcp.2007.11.035.
[SLJB99]S. Sanvito, C. J. Lambert, J. H. Jefferson, and A. M. Bratkovsky. General green’s-function formalism for transport calculations with spd hamiltonians and giant magnetoresistance in co- and ni-based magnetic multilayers. Phys. Rev. B, 59:11936–11948, May 1999. doi:10.1103/PhysRevB.59.11936.
[STBG05]K. Stokbro, J. Taylor, M. Brandbyge, and H. Guo. Ab-initio based Non-equilibrium Green’s Function Formalism for Calculating Electron Transport in Molecular Devices., pages 117–151. Springer, 2005.
[SMB+16](1, 2, 3) D. Stradi, U. Martinez, A. Blom, M. Brandbyge, and K. Stokbro. General atomistic approach for modeling metal-semiconductor interfaces using density functional theory and nonequilibrium green’s function. Phys. Rev. B, 93:155302, Apr 2016. doi:10.1103/PhysRevB.93.155302.
[SrensenHP+08]H. H. B. Sørensen, P. C. Hansen, D. E. Petersen, S. Skelboe, and K. Stokbro. Krylov subspace method for evaluating the self-energy matrices in electron transport calculations. Phys. Rev. B, 77:155301, Apr 2008. doi:10.1103/PhysRevB.77.155301.
[SrensenHP+09]H. H. B. Sørensen, P. C. Hansen, D. E. Petersen, S. Skelboe, and K. Stokbro. Efficient wave-function matching approach for quantum transport calculations. Phys. Rev. B, 79:205322, May 2009. doi:10.1103/PhysRevB.79.205322.
[Tod98](1, 2) T. N. Todorov. Local heating in ballistic atomic-scale contacts. Philosophical Magazine Part B, 77(4):965–973, 1998. doi:10.1080/13642819808206398.
[ZRSH11]R. Zhang, I. Rungger, S. Sanvito, and S. Hou. Current-induced energy barrier suppression for electromigration from first principles. Phys. Rev. B, 84:085445, Aug 2011. doi:10.1103/PhysRevB.84.085445.