Quantum Transport in a Silicon Nanowire FET Transistor: Hot Electrons and Local Power Dissipation

Abstract

A review and perspective is presented of the classical, semi-classical and fully quantum routes to the simulation of electro-thermal phenomena in ultra-scaled silicon nanowire field-effect transistors. It is shown that the physics of ultra-scaled devices requires at least a coupled electron quantum transport semi-classical heat equation model outlined here. The importance of the local density of states (LDOS) is discussed from classical to fully quantum versions. It is shown that the minimal quantum approach requires self-consistency with the Poisson equation and that the electronic LDOS must be determined within at least the self-consistent Born approximation. To bring in this description and to provide the energy resolved local carrier distributions it is necessary to adopt the non-equilibrium Green function (NEGF) formalism, briefly surveyed here. The NEGF approach describes quantum coherent and dissipative transport, Pauli exclusion and non-equilibrium conditions inside the device. There are two extremes of NEGF used in the community. The most fundamental is based on coupled equations for the Green functions electrons and phonons that are computed at the atomically resolved level within the nanowire channel and into the surrounding device structure using a tight binding Hamiltonian. It has the advantage of treating both the non-equilibrium heat flow within the electron and phonon systems even when the phonon energy distributions are not described by a temperature model. The disadvantage is the grand challenge level of computational complexity. The second approach, that we focus on here, is more useful for fast multiple simulations of devices important for TCAD (Technology Computer Aided Design). It retains the fundamental quantum transport model for the electrons but subsumes the description of the energy distribution of the local phonon sub-system statistics into a semi-classical Fourier heat equation that is sourced by the local heat dissipation from the electron system. It is shown that this self-consistent approach retains the salient features of the full-scale approach. For focus, we outline our electro-thermal simulations for a typical narrow Si nanowire gate all-around field-effect transistor. The self-consistent Born approximation is used to describe electron-phonon scattering as the source of heat dissipation to the lattice. We calculated the effect of the device self-heating on the current voltage characteristics. Our fast and simpler methodology closely reproduces the results of a more fundamental compute-intensive calculations in which the phonon system is treated on the same footing as the electron system. We computed the local power dissipation and “local lattice temperature” profiles. We compared the self-heating using hot electron heating and the Joule heating, i.e., assuming the electron system was in local equilibrium with the potential. Our simulations show that at low bias the source region of the device has a tendency to cool down for the case of the hot electron heating but not for the case of Joule heating. Our methodology opens the possibility of studying thermoelectricity at nano-scales in an accurate and computationally efficient way. At nano-scales, coherence and hot electrons play a major role. It was found that the overall behaviour of the electron system is dominated by the local density of states and the scattering rate. Electrons leaving the simulated drain region were found to be far from equilibrium.

1. Introduction

1.1. Context

Silicon metal-oxide semiconductor-field effect transistors (MOSFETs) have underpinned the last fifty years of the microelectronics revolution. The scalability of the basic planar single gate MOSFET has already approached the long-anticipated physical limits to device performance that require considerable architectural ingenuity to overcome. These limits include short-channel effects (SCE) such as: quantum mechanical tunnelling between source and drain and between drain and substrate; drain-induced barrier lowering (DIBL); threshold voltage roll-off. However, since the supply voltage has not been downscaled as smaller feature sizes have been achieved, the power densities have increased close to the limit of 150 W·cm⁻² beyond which air cannot cool the device temperature [1]. The move to a 3D FinFET architecture has given some respite [2,3]. Electronic devices at nanometre scales have already been fabricated and tested [4,5,6,7]: effective gate lengths of some of the recently produced transistors are well under 20 nm [8]. There is a consensus emerging that the most promising successor to the planar MOSFET is the gate-all-around nanowire-FET (NWFET) (a cross-sectional model is displayed in Figure 1 of Section 3); it has far superior electrostatic control to the FinFET [9]. Silicon NWFETs have been fabricated with circular cross sections with diameters in a range of 9 to 1 nm [7]. These devices show very good performances and ideal sub-threshold slopes (60 mv/dec). They also showed a substantial decrease in the on-current as the cross section decreases that cannot be explained by the reduction in area. Instead, it was suggested that increases in scattering rates or the enhancement of other scattering mechanism such as surface scattering or stresses in the fabrication might be responsible. In a recent study [10], we have shown that NWFETs with ultra-scaled square cross section and circular cross-section have very similar properties.

1.2. Electro Thermal Properties of NWFET

The gate-all-around NWFET has excellent SCE reduction at very small gate lengths (~10 nm) but requires small diameters (~5 nm) for best performance. At these scales, however, there are issues with atomistic effects (the effect of the random discrete nature of dopants and defects [11]) but the major future issue has been long-known: any conventional electronic logic is ultimately limited by heat dissipation [12]. Comparatively few studies have been made of the electro-thermal properties of NWFETs. It is known that silicon nanowire structures have a reduced thermal conductivity compared to bulk [13,14] due to boundary scattering of phonons at the oxide sheath-silicon channel interface, which further exacerbates self-heating and hence device performance. Such devices are expected to be prone to the appearance of self-heating effects and localised hot spots. Unfortunately, the small dimensions of ultra-scaled nanowires provides a serious challenge to measurements of self-heating effects and local hot spots (although recent progress in scanning thermometry is promising [15,16]).

1.3. Need for Efficient Modelling

Fabrication of nano-devices is of high lithographic precision and it is expensive. These factors make a comprehensive experimental investigation of a wide range of prototype devices prohibitive. From a design point of view, there is thus scope for good fast physics-based device modelling tools for NWFETs that assists the control of adverse electro-thermal effects. Good theoretical analysis of the electrical and carrier transport properties of proposed nanowire devices is required prior to optimised design and fabrication; in particular, the understanding of the interaction between the carrier flow and the electrostatic control.

The simulation of the thermal performance, power dissipation and carrier lattice interaction are also important in the design of self-cooling and heat management of nanostructures. High channel temperatures lead to raised temperatures in the interconnects at the silicon source extension-metal contact [17]. The scope for conduction cooling via the source/drain/interconnects [18,19,20] has been appreciated for some time [21]. The number of research programmes on nano-cooling is rising [22]. In addition, research into the impact of quantum effects on the thermodynamic efficiency of nanostructures [22] has been a long-standing concern.

1.4. Review of Modelling Approaches

Electronic transport and heat transport at nano scales have been studied by a variety of methodologies. Semi-classical techniques complemented with quantum corrections such as the Drift-Diffusion (DD) approach [23,24,25,26,27] and the Monte Carlo method [28,29,30,31,32,33,34,35,36,37] have been widely used.

The Drift-Diffusion approach is particularly computationally efficient but assumes that the electron charge fluid is both in local equilibrium [38] and moves due to electric fields and gradients of concentration. The equation to be solved is the charge density conservation law usually called the continuity equation. Generalizations of the DD approach and energy balance in concomitance with the heat equation have been carried out, a thorough investigation is presented in [39]. An extra equation can be added to the Drift-Diffusion model to treat the local changes in momentum or velocity of the charge fluid; the model that results is termed the hydrodynamic model. However, quantum effects such as confinement and tunneling need to be artificially added.

The Monte Carlo method [28] which is founded on a particle representation of the semi-classical Boltzmann equation is ideal to treat non-equilibrium processes semi-classically; the flows are energy-resolved but the carrier or electron transport is treated classically. It is assumed that the electron de Broglie wavelength is much smaller than the scale of variation of the electrical potential. This last assumption is violated in small nanostructures. Any quantum coherent transport is neglected in the semi-classical formalism. In addition, the scale of the thermalisation (as will be shown in this paper) is much larger that the dimension of the nanostructure.

A powerful quantum transport methodology, the Keldysh non-equilibrium Green function (NEGF) formalism, was applied from the 70’s in order to describe non-equilibrium quantum transport process beyond the linear response in nanostructures [40]. The formalism was used to study high electric field effects [41]. The original formalism was developed in the 60’s by Schwinger [42], Keldysh [43], Kadanov and Baym [44] to address highly non-equilibrium processes such as those in Masers. The formalism uses matrix equations containing a mix of dynamical and kinetic Green functions [43].

The transport methodology satisfies the second law of thermodynamic as the entropy generated by the transfer of electrons from one reservoir to another reservoir is always positive [22]. The non-equilibrium scattering rates are self-consistently dependent on the local electric field (or equivalently the local electrostatic potential), which is not the case for Boltzmann transport. This feature was pointed out using a super-resolvent non-linear time-dependent quantum kinetic equation based on the density matrix formalism [45,46,47], and confirmed in later studies by the NEGF formalism [41]. Electro-thermal simulations have been carried out using a variety of methodologies ranging from semiclassical [48], a blend of quantum transport with classical heat equation to a pure quantum mechanical [49] and atomistic methodologies [50]. All these methodologies have their weaknesses and strengths. The semi-classical methodologies are very efficient and handle micrometer size devices however they ignore coherence, tunneling and renormalization of the density of states and therefore lose accuracy for devices of few tens of nanometers. The pure atomistic methodologies (mostly based in the NEGF formalism) consider the wave nature of electrons, coherence and quantum effects that require substantial computational power, and it is computationally feasible only for very small dimensions. However, effects and boundary conditions that are far from the active device region could affect the electrical and heat performance of the device, for example, remote optical phonon scattering and remote coulomb interactions [51,52,53,54], phonon reflections or transmission through interfaces [55,56].

1.5. Content of This Paper

In the present paper, Section 2 provides a brief review and perspective on the key concepts and the transition from semi-classical to fully quantum mechanical modelling of electro-thermal effects applicable to quasi-one dimensional NWFETs. Section 2.1 discusses the basic concepts of electro-thermal modelling. Section 2.2 describes the crucial local density of states in energy and the need for a transition from semi-classical to full quantum transport for electro-thermal device modelling. Section 2.3 gives a brief overview of the non-equilibrium Green function (NEGF) methodology augmented by explicit technical details in an Appendix (Appendix A discusses the local density of states for nanowires and the multi-valley picture; Appendix B describes how contact self-energies are used to project the open many-body problem into a finite device simulation domain; Appendix C discusses details of the Keldysh NEGF formalism; Appendix D gives explicit details of the phonon Green functions and the generic electron-phonon self-energies; Appendix E demonstrates explicitly that the NEGF formalism shows no heat-dissipation for elastic acoustic phonon scattering; Appendix F describes detailed balance within the NEGF formalism; Appendix G describes the importance of the Kramers-Kronig relations for self-energies).

In Section 3, we present our simulation methodology for the silicon gate-all-around nanowire FET that uses the NEGF picture for the electrons coupled to a semi-classical Fourier heat equation. The device structure and the quasi-1D representation is introduced in Section 3.1. The NEGF methodology for the density and current density is outlined in Section 3.2. Section 3.3 discusses energy balance and power. Section 3.4 presents the results and discussion including comparison with other simulated cases. Section 3.5 discusses perspectives and challenges. Finally, the conclusions are presented in Section 4. In the following we exclude radiative emission and absorption as this is not relevant for conventional integrated circuits, although this is relevant to photo-excited hot electron phenomena.

2. Methodology and Models

2.1. Basic Concepts of Electro-Thermal Modelling

The key quantity in electro-thermal phenomena is the local energy current density J_Energy(r,t) synonymous with the heat current density. If we consider a finite volume Ω bounded by an enclosing surface S then the total power generated within the volume Ω is balanced by the net integrated energy current density flowing out of the surface S:

$${{\displaystyle \iiint }}_{\mathsf{\Omega }} P\left(\mathrm{r},t\right){\mathrm{d}}^3\mathrm{r}=-{{\displaystyle \iint }}_S{J}_{Energy}\left(\mathrm{r},t\right)·\mathrm{d}S$$

(1)

where $P\left(\mathrm{r},t\right)$ is the local power density, synonymous with the local rate of heat generation. Applying the divergence theorem we obtain:

$$P\left(\mathrm{r},t\right)=-\nabla ·J_{Energy}\left(\mathrm{r},t\right)$$

(2)

If Ω represents the active volume of an FET device there is no net heat generation in the device when it is switched off as it is in thermal equilibrium with the lattice and its environment. If the device is switched on (gate voltage applied) and a (chemical potential difference) voltage difference established between the source and drain contacts then an electrical current flows characterised by the current density J(r,t) = ej(r,t) where j(r,t) is the particle current density.

In a purely ballistic channel [57,58] for which the electron flow is quantum coherent (no dissipative scattering), the electron transport involves a net injection of electrons at high energy and momentum from the source followed by net extraction at the drain where the excess energy and momentum are dissipated: the flow is quantum transmission and the conductance G = 1/R (R:The electrical resistance) of the channel is quantised in units of e²/πħ. This extreme case (best studied for low temperatures) illustrates that quantum transport can lead to hot spots, in this case strong heating occurs in the drain contact. If a defect or impurity occurs in a ballistic channel such that its interaction with an electron is representable as a simple scalar potential perturbation, the flow then involves reflection and tunnelling transmission that may be described by T_c the transmission coefficient for the channel. Then G_min = T_ce²/πħ which is Landauer’s formula [59]. If T_c < 1 there is elastic (potential) scattering in the channel. Again, the dissipation of energy is in the contact region; the dissipation of momentum is partly within the channel via the defect/impurity. This simple picture ignores the electrical field profile E(r,t) = −▽ϕ(r,t) in the channel. In practice, it is neither abrupt nor uniform; instead it may be determined self-consistently from Poisson’s Equation:

$$\nabla ·{\epsilon }_d\nabla \varphi \left(\mathrm{r},t\right)=-\rho (\mathrm{r},t)$$

(3)

where ρ is the charge density in the channel and ε_d is the dielectric permittivity.

Silicon nanowire transistors do not operate in the purely ballistic regime although elements of the simple transmission/reflection flows may occur as we discuss later. In general, the applied electric field supplies energy and momentum to the carriers throughout the device and the electrostatic potential varies smoothly between source and drain. For dissipative transport the energy and momentum gained locally from the field will be transferred to the phonon system by scattering within the device involving a loss-gain process in which phonons may be emitted or absorbed by the electron assembly. Inelastic scattering gives rise to energy dissipation whereas the momentum dissipation occurs for elastic and inelastic processes.

More detailed non-phenomenological analysis requires an electron transport formalism and ultimately a coupled phonon transport picture. The deepest semi-classical picture is provided by the well-known Boltzmann transport equation for the electron probability distribution f[r,k,t]:

$$\frac{\partial f\left[\mathrm{r},\mathrm{k},t\right]}{\partial t}+\mathrm{v}\left(\mathrm{k}\right)·\frac{\partial f\left[\mathrm{r},\mathrm{k},t\right]}{\partial \mathrm{r}}+e\mathrm{E}\left(\mathrm{r},t\right)·\frac{\partial f\left[\mathrm{r},\mathrm{k},t\right]}{\partial \hslash \mathrm{k}}=-\frac{\partial f\left[\mathrm{r},\mathrm{k},t\right]}{\partial t}{|}_s$$

(4)

Here the velocity v is determined from the k-dependence of the conduction band energy ε(k):

$$\mathrm{v}\left(\mathrm{k}\right)=\frac{1}{\hslash }\frac{\partial \epsilon \left(\mathrm{k}\right)}{\partial \mathrm{k}}$$

(5)

The second term in Equation (4) represents diffusive effects, the third term represents the drift of the distribution induced by an applied electric field, the RHS represents the change in f due to scattering processes. The scattering integral (actually a sum over k-states) is given by

$$\frac{\partial f\left[\mathrm{r},\mathrm{k},t\right]}{\partial t}{|}_s={\displaystyle {\displaystyle \sum }_{\mathrm{k}{}^{\prime }}}\left\{f\left[\mathrm{r},\mathrm{k}{}^{\prime },t\right]\left(1-f\left[\mathrm{r},\mathrm{k},t\right]\right)R\left[\mathrm{r},\mathrm{k}{}^{\prime },\mathrm{k},t\right]-f\left[\mathrm{r},\mathrm{k},t\right]\left(1-f\left[\mathrm{r},\mathrm{k}{}^{\prime },t\right]\right)R\left[\mathrm{r},\mathrm{k},\mathrm{k}{}^{\prime },t\right]\right\}$$

(6)

We shall be mostly concerned with phonon scattering for which the scattering rate from k to k’ may be written as a sum over the phonon modes b as:

$$\begin{array}{cc}\hfill R\left[\mathrm{r},\mathrm{k},{\mathrm{k}}^{\prime },t\right]=& {\displaystyle \sum _{\mathrm{b}\mathrm{K}}}{\Lambda }_b\left[\mathrm{k},{\mathrm{k}}^{\prime }\right]\left\{{\delta }_{k-q-k^{\prime },K}\left(f_{\mathrm{phonon}b}[\mathrm{r},\mathrm{q},t]+1\right)\delta \left(\epsilon \left({\mathrm{k}}^{\prime }\right)-\epsilon \left(\mathrm{k}\right)+\hslash {\omega }_{\mathrm{bq}}\right)\right.\hfill \\ \hfill & \left.+{\delta }_{k+q-k^{\prime },K}{f}_{\mathrm{phonon}b}[\mathrm{r},\mathrm{q},t]\delta \left(\epsilon \left({\mathrm{k}}^{\prime }\right)-\epsilon \left(\mathrm{k}\right)-\hslash {\omega }_{b\mathrm{q}}\right)\right\}\hfill \end{array}$$

(7)

Here Λ_b is the scattering kernel, K represents a reciprocal lattice vector for Umklapp processes, ħω_bq is the energy of the b mode phonon at wave vector q. The delta-functions in energy and momentum represent the conservation of energy and momentum within a scattering event. The function f_{phonon b}[r,q,t] is the probability distribution function for a phonon in mode b with wave vector q. The scattering rates are local in space and time, and evaluated in the Born approximation and do not depend on the applied electrostatic field. The first term in the brackets {…} of Equation (7) represents electron scattering with the emission of a phonon; the second term describes absorption of a phonon. If the phonon distributions are in internal thermal equilibrium at a well-defined lattice temperature T they assume the form of a Bose-Einstein distribution:

$$f_{\mathrm{phonon} b}\left[\mathrm{r},\mathrm{q},t\right]\to (e^{\frac{\hslash {\omega }_{bq}}{k_{B}T}}-1)$$

(8)

The Boltzmann transport equation is usually solved numerically by using the Monte-Carlo method [28]. Let us now specialise to the steady-state electro-thermal problem for which $\frac{\partial f}{\partial t}=0$.

The local kinetic energy current density is given by:

$${\mathrm{J}}_{kinetic}\left(\mathrm{r}\right)={\displaystyle \sum }_k^{ }\epsilon \left(\mathrm{k}\right)\mathrm{j}\left(\mathrm{r},\mathrm{k}\right)\to {\displaystyle \int }{d}^{3}k {\rho }_k \epsilon \left(\mathrm{k}\right)\mathrm{j}\left(\mathrm{r},\mathrm{k}\right)={\displaystyle \int }{d}^{3}k {\rho }_k \epsilon \left(\mathrm{k}\right)\mathrm{v}\left(\mathrm{k}\right)f\left[\mathrm{r},\mathrm{k}\right]$$

(9)

where the summation over k is replaced by integration by introducing ${\rho }_k$ the density of k-states.

We may obtain an expression for J_kinetic(r) by multiplying the steady-state version of Equation (4) by ε(k) and summing over all k to obtain after an integration by parts and re-arranging the collision integral:

$$\nabla ·{\mathrm{J}}_{kinetic}\left(\mathrm{r}\right)-\mathrm{E}\left(\mathrm{r}\right)·\mathrm{J}\left(\mathrm{r}\right)=-{\displaystyle \sum }_k^{ }\epsilon \left(\mathrm{k}\right)\frac{\partial f\left[\mathrm{r},\mathrm{k}\right]}{\partial t}{|}_s$$

(10)

$$\begin{array}{cc}={\displaystyle \sum _k}{\displaystyle \sum _{k^{\prime }}}\left\{\epsilon \left(\mathrm{k}\right)-\epsilon \left({\mathrm{k}}^{\prime }\right)\right\}f[\mathrm{r},\mathrm{k}]\left(1-f\left[\mathrm{r},{\mathrm{k}}^{\prime }\right]\right)R[\mathrm{r},\mathrm{k},\mathrm{k}]\hfill & =-P\left(\mathrm{r}\right)\hfill \end{array}$$

(11)

$$P\left(\mathrm{r}\right)=-\nabla ·{\mathrm{J}}_{Energy}\left(\mathrm{r}\right)=-\nabla ·{\mathrm{J}}_{kinetic}\left(\mathrm{r}\right)+\mathrm{E}\left(\mathrm{r}\right)·\mathrm{J}\left(\mathrm{r}\right)$$

(12)

where from Equation (7) ε(k’)–ε(k) is the energy transferred to or absorbed from the lattice per unit time when an electron scatters from state k’ to state k by emitting or absorbing a phonon. The explicit expression for the local power density P(r) is

$$P\left(\mathrm{r}\right)={\displaystyle \sum _k}{\displaystyle \sum _{k^{\prime }}}\left\{\epsilon \left({\mathrm{k}}^{\prime }\right)-\epsilon \left(\mathrm{k}\right)\right\}f[\mathrm{r},\mathrm{k}]\left(1-f\left[\mathrm{r},{\mathrm{k}}^{\prime }\right]\right)R\left[\mathrm{r},\mathrm{k},{\mathrm{k}}^{\prime }\right]$$

(13)

showing that for elastic scattering processes ε(k’) = ε(k), the contribution to the power density vanishes: no dissipation. Re-calling our discussion of ballistic transport it is evident from Equation (13) that in semi-classical transport the heat dissipation must occur where the scattering rate R[r,k,k’] is non-zero.

Equation (12) makes the Joule heating explicit, where E(r) is the local electrical field and J(r) = ej(r) is the electrical current density. The Joule power density E.J when integrated throughout the device produces the familiar Joule power I V. The first term on the RHS of Equation (12) tracks the spatial variation of the “kinetic energy” current density of the electron ensemble. Crucially, if the kinetic energy current density J_kinetic(r) does not change with position r the electron energy distribution will follow the bending of the electrostatic potential and the electron mean kinetic energy will not increase relative to the conduction band edge; therefore the energy distribution of the electrons will remain in local equilibrium, there will be no hot electron phenomena. On the other hand, if there is a large change in the spatial dependence of the kinetic energy current density there will be heating or cooling of the electronic sub-system.

Equations (9) and (10) lead directly to the thermodynamically significant detailed balance equations and have much wider generality than the classical picture.

The recognised route to simulating devices operating under non-isothermal conditions is to add a transport model to a standard heat equation. Using the local heat generation rate P(r,t) as the source of excess heat we obtain the heat equation (Fourier’s Equation) for the local lattice temperature T(r,t):

$$C_h\frac{\partial T\left(\mathrm{r},t\right)}{\partial t}=\nabla ·\kappa \nabla T\left(\mathrm{r},t\right)+P\left(\mathrm{r},t\right)$$

(14)

where $C_h$ is the total heat capacity and κ(r,t) is the local thermal conductivity.

For steady-state simulations Equation (8) reduces to:

$$\nabla ·\kappa \nabla T\left(\mathrm{r}\right)=-P\left(\mathrm{r}\right)$$

(15)

There are three main approaches to specify the heat source term P: (i) the Joule-heating model, exemplified by local equilibrium carriers commonly used in drift-diffusion modelling; (ii) assuming non-equilibrium electrons (hot electrons) interacting with equilibrium phonons; (iii) non-equilibrium electrons interacting with local equilibrium phonons.

(i) If we neglect electron-hole recombination and generation processes the Joule heating model [60] assumes:

$$P\left(\mathrm{r},\mathrm{t}\right)=\mathrm{E}\left(\mathrm{r},\mathrm{t}\right)·\mathrm{J}\left(\mathrm{r},\mathrm{t}\right)$$

(16)

where E is the local electric field and J is the local electrical current density.

Evidently, the highest heating rate will occur where the scalar product is at a maximum. In conventional FETs most studies indicate that the maximum heating occurs under the gate where most of the potential drop takes place and where the electrical current density is the largest.

(ii) Here, we require an Equations (12)–(14) which picks up the transfer of energy from the electron system to the phonons in the local lattice. If the phonons are assumed to be in internal thermal equilibrium at temperature T then we can evaluate the net energy loss via the RHS of Equation (14) by using the matrix-elements and phonon distribution functions within the scattering rates R of Equation (7) by their equilibrium values [12,61]. More generally, we may evaluate a local lattice temperature T(r) (from self-consistent solution of the coupled heat balance equation [62,63]).

(iii) To obtain full generality [34,35,36,64,65,66,67,68,69], we must introduce and solve self-consistently the full set of phonon transport equations for each phonon type at each location in the system including internal equilibrium processes such as phonon-phonon scattering and external processes such as boundary scattering. If the phonon distributions are parameterised by a different non-equilibrium value T_b(r) for each mode b then it is still relatively easy to proceed; if not, the problem becomes a major challenge.

Further progress is only possible by resorting to a full space-dependent, energy-resolved quantum statistical approach to which we now turn.

2.2. LDOS and the Transition from Classical to Quantum Modelling

Consideration of the local energy density of states (LDOS) provides a clear example of the need for a transition to a quantum description of electro-thermal transport. Semi-classically, one assigns the particle momentum p to a wave vector k by p = ħk, where the underling implication of a plane wave state exp[ik·r] only has meaning in the evaluation of the scattering rate matrix elements in the collision integral. Its amplitude is normalised to unity. The local density of states is then defined as:

$${\rho }_{LDOS}\left(\mathrm{r},\epsilon \right)={\displaystyle {\displaystyle \sum }_k} \delta \left(\epsilon -\epsilon \left(\mathrm{k}\right)-V\left(\mathrm{r}\right)\right)={\displaystyle \int }{\mathrm{d}}^{3}k {\rho }_k\delta \left(\epsilon -\epsilon \left(\mathrm{k}\right)-V\left(\mathrm{r}\right)\right)$$

(17)

Noting that ρ_k = 2L^N/(2π)^N where N is the number of dimensions and L is the normalisation length, Equation (14) is easily evaluated for 3D, 2D and 1D systems as:

$${\rho }_{LDOS}^{3D}\left(x,y,z,\epsilon \right)=\left(\frac{m}{{\pi }^2{\hslash }^2}\right)\sqrt{\frac{2m\left(\epsilon -V\left(x,y,z\right)-{\epsilon }_c\right)}{{\hslash }^2}}$$

(18)

$${\rho }_{LDOS}^{2D}\left(x,y,\epsilon \right)=\left(\frac{m}{{\pi }^{ }{\hslash }^2}\right)\theta \left(\epsilon -V\left(x,y\right)-{\epsilon }_c\right)$$

(19)

$${\rho }_{LDOS}^{1D}\left(x,\epsilon \right)=\left(\frac{2m}{{\pi }^{ }{\hslash }^2}\right)/\sqrt{\frac{2m\left(\epsilon -V\left(x\right)-{\epsilon }_c\right)}{{\hslash }^2}}$$

(20)

where m is the electron mass for a simple isotropic band structure, ε_c is the conduction band edge; it is assumed that the dimensionality factor L = 1 and $\theta$ is the unit step function. In the Appendix A, Figure A1, Equation (20) is illustrated classically and for a simple quantum generalisation for a typical band edge profile for a 1D NWFET.

The quantum origins of Equation (17) are very simple; consider for example, a general quantum Hamiltonian H that has a complete set of eigenvalues and eigenstates {ε_α,|α >} with completeness condition:

$${\displaystyle {\displaystyle \sum }_{\alpha }}|\alpha ><\alpha |=\mathrm{I}$$

(21)

where I is a unit operator. The density of states in energy per unit volume follows as

$${\rho }_{DOS}\left(\epsilon \right)={\displaystyle {\displaystyle \sum }_{\alpha }}\delta \left(\epsilon -{\epsilon }_{\alpha }\right)={\displaystyle {\displaystyle \sum }_{\alpha }}<\alpha |\delta \left(\epsilon -H\right)|\alpha >$$

(22)

$${\rho }_{DOS}\left(\epsilon \right)=Tr\left[\delta \left(\epsilon -H\right)\right]$$

(23)

$${\rho }_{DOS}\left(\epsilon \right)={\displaystyle {\displaystyle \sum }_{\beta }} < \beta |\delta \left(\epsilon -H\right)|\beta >={\displaystyle {\displaystyle \sum }_{\beta ,\alpha }}|<\beta |\alpha >{|}^2 \delta \left(\epsilon -{\epsilon }_{\alpha }\right)$$

(24)

where Equation (23) shows the representation independence of the trace operation. Equation (24) re-expresses the result in terms of any basis {ε_β,|β >} in particular if we choose the position representation {|r >}, Equation (24) leads to the local density of states

$${\rho }_{LDOS}\left(\mathrm{r},\epsilon \right)={\displaystyle {\displaystyle \sum }_{\alpha }}|<\mathrm{r}|\alpha >{|}^2 \delta \left(\epsilon -{\epsilon }_{\alpha }\right)=<\mathrm{r}|\delta \left(\epsilon -H\right)|\mathrm{r}>$$

(25)

For a simple free electron model the eigenstates $\alpha$are just k-states : <r|k> = exp[ik·r] with eigenvalues ε(k) and we recover Equation (14) for V = 0.

By using the operator identity

$$\delta \left(\epsilon -H\right)={\displaystyle \underset{\eta \to 0^{+}}{\lim }} \frac{1}{2\pi i}\{{(\epsilon -i\eta -H)}^{-1}-{(\epsilon +i\eta -H)}^{-1}\}$$

(26)

Equation (25) may be formulated in terms of the advanced and retarded Green’s operators for the steady-state operator equation

$$\left(\epsilon +i\eta -H\right)G^R=\mathrm{I}; \left(\epsilon -i\eta -H\right)G^A=\mathrm{I}$$

(27)

$$\begin{array}{cc}\hfill {\rho }_{LDOS}(\mathrm{r},\epsilon )=<\mathrm{r}|& \delta (\epsilon -H){\left|\mathrm{r}>=\frac{-1}{2\pi i}\left\{G^R(\mathrm{r},{\mathrm{r}}^{\prime },\epsilon )-G^A(\mathrm{r},{\mathrm{r}}^{\prime },\epsilon )\right\}\right|}_{{\mathrm{r}}^{\prime }\to \mathrm{r}}\equiv A(\mathrm{r},\epsilon )\hfill \\ \hfill =& \left(\frac{-1}{\pi }\right)Im{G}^R{\left(\mathrm{r},{\mathrm{r}}^{\prime },\epsilon \right)}_{\infty }\hfill \end{array}$$

(28)

In the position representation Equations (27) are matrix equations for the Green functions G^R(r,r’,ε), G^A(r,r’,ε); they determine the description of the energy structure of the system.

In thermal equilibrium, the energy distribution of electronic states will follow a Fermi-Dirac distribution

$$f_0\left(\epsilon -\mu \right)={\{1+\exp \left[\frac{\left(\epsilon -\mu \right)}{k_{B}T}\right]\}}^{-1}$$

(29)

so that the equilibrium particle density will be

$$\begin{gathered} n_0\left(\mathrm{r},\epsilon \right)={\displaystyle {\displaystyle \sum }_{\alpha }}|<\mathrm{r}|\alpha >{|}^2 \delta \left(\epsilon -{\epsilon }_{\alpha }\right)f_0\left(\epsilon -{\epsilon }_{\alpha }-\mu \right) \\ =<\mathrm{r}\mid \delta (\epsilon -H)f_0(\epsilon -H-\mu N)\mid \mathrm{r}>\equiv \frac{-i}{2\pi }{G_0^{<}(\mathrm{r},{\mathrm{r}}^{\prime },\epsilon )\mid }_{{\mathrm{r}}^{\prime }\to \mathrm{r}} \end{gathered}$$

(30)

Here we introduce the statistical (superscript “<” denotes “lesser than”) equilibrium Green function. This is easy to evaluate if the eigenstates of H are known. Equations (21)–(30) provide a simple account of how the NEGF formalism arises.

However, if H is split into an interaction free Hamiltonian H₀ (which may include the Hamiltonian for non-interacting phonons) and an interaction term H_int describing the electron-phonon interactions, then the trace operation in Equations (28) and (30) must be widened to include a trace over the phonon states to maintain completeness. A self-consistent perturbation theory therefore requires a simultaneous treatment of the states within thermal equilibrium and the dynamical evolution of the states in non-equilibrium.

The route to doing this was laid out in [41,42,43,44] in the 1960s, but was difficult to apply, partly because of computational resources. Important technical issues and improvements were made much later [40,70,71,72,73,74]

2.3. The NEGF Formalism

The NEGF formalism is ultimately based on statistical quantum field theory; in the following we will draw from references [41,42,43,44,70,71,72,73,74], but we will only review the parts of that theory necessary to consider electro-thermal transport and we avoid the very technical derivations and issues. An excellent practical introduction to the NEGF as applied to devices will be found in works and references therein by Datta [57,58], background many-body perturbation theory is discussed in detail in reference [71]. The Keldysh approach is outlined in 2.3.1 but may skipped to proceed directly to the required steady-state NEGF equations that are described in Section 2.3.2. More technical details will be found in the Appendix, Sections Appendix B, Appendix C, Appendix D, Appendix E, Appendix F and Appendix G.

2.3.1. The Keldysh Picture

Let us start with a generic second-quantised many-body Hamiltonian H = H₀ + H_int describing a system of interacting (via H_int) electrons and phonons. Any property of the system may be represented by an operator $\mathsf{\Xi }$ and its statistical expectation value. < $\mathsf{\Xi }>$ may be computed in thermal equilibrium from the grand canonical ensemble ${\rho }_{eq}$ of statistical mechanics giving:

$$< \mathsf{\Xi }\left(\mathrm{t}\right){>}_{eq}=\frac{Tr\left[\mathsf{\Xi }{\rho }_{eq}\right]}{Tr\left[{\rho }_{eq}\right]}; {\rho }_{eq}=\exp \left[-\frac{H-\mu N}{k_{B}T}\right]{\rho }_{eq}=\exp \left[-\frac{H-\mu N}{k_{B}T}\right]$$

(31)

where T is the absolute temperature and µ is the chemical potential. Here N represent the second-quantised number operator for excitations (electrons and phonons). The trace operation Tr[…] means sum the diagonal matrix elements of the operand in any complete set of many-body states.

An important example of interest is the single electron equilibrium Green function that has a two-point space-time structure:

$$<\mathsf{\Xi }\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right){>}_{eq}\equiv {\mathrm{G}}_{eq}\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right)\equiv <\mathrm{T}\psi \left(\mathrm{r},t){\psi }^{†}({\mathrm{r}}^{\prime },t^{\prime }\right){>}_{eq}$$

(32)

Equation (32) represents the probability amplitude for the insertion of an electron at (r’,t’) and the removal of an electron at (r,t). The operators ψ(r,t), ψ^†(r,t) are anti-commuting electron annihilation and creation field operators in the Heisenberg picture (state vectors constant, operators time-dependent.) T is the Dyson time-ordering operator which for fermions re-orders the field operators so that the term with earliest time is on the right and also multiplies the result by +1 or −1 depending on whether the new order is an even or odd permutation of the original product.

Different time-orderings of Equation (32) give Green functions with different behaviour. The Green functions ${\mathrm{G}}_{eq}^{<}, {\mathrm{G}}_{eq}^{>}$ (cf. Section 2.2) are the specific time-orderings:

$${\mathrm{G}}_{eq}^{<}\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right)={\mathrm{G}}_{eq}^{ }\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right) (t<t^{\prime })$$

(33)

$${\mathrm{G}}_{eq}^{>}\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right)={\mathrm{G}}_{eq}^{ }\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right)(t>t^{\prime })$$

(34)

$${\mathrm{G}}_{eq}^R\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right)=-\mathrm{i}\mathsf{\theta }\left(t-t^{\prime }\right)<\left\{{\psi }^{†}\left(\mathrm{r}{}^{\prime },t\prime ), {\psi }^{ }({\mathrm{r}}^{ },t^{ }\right)\right\}{>}_{eq}$$

(35)

$${\mathrm{G}}_{eq}^A\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right)=\mathrm{i}\mathsf{\theta }\left(t\prime -t\right)<\left\{{\psi }^{ }\left(\mathrm{r},t), {\psi }^{†}({\mathrm{r}}^{\prime },t^{\prime }\right)\right\}{>}_{eq}$$

(36)

The “lesser/greater” Green functions G^</G^> are so named by the conditions t < t’ or t > t’). Here, the brackets { , } represent the anti-commutator.

Using the above Equations (31) and (32) it is relatively straightforward to derive Equations (28) and (30). Equation (31) plus considerable manipulation yields an important boundary condition that relates G^< to G^> in thermal equilibrium:

$${\mathrm{G}}_{eq}^{<}\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right){|}_{t=0}=-\exp \left[\beta \mu \right]{\mathrm{G}}_{eq}^{>}\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right) {|}_{t=-i\hslash \beta }$$

(37)

The latter is the basis for developing finite temperature thermal equilibrium perturbation theory in complex time, see reference [71]

The non-equilibrium problem is to start at some initial time t₀ with a system (such as a device) in thermal equilibrium and then initiate an external perturbation H_ext (applied fields or temperature gradients) that drives the system into non-equilibrium. To establish an infinite order perturbation theory for interacting systems we write the full Hamiltonian as

$$H=H_0+H_{ext}\left(t\right)+H_{int}$$

(38)

where H₀ describes free excitations (electrons, phonons in our case); H_ext (t) describes coupling to external fields switched on at some time > t₀; and H_int represents interaction between excitations (electron-phonon coupling in our case). The system density matrix satisfies the equation of motion:

$$i\hslash \frac{\partial \rho }{\partial t}=\left[H,\rho \right]$$

(39)

Assuming adiabatic switching-on of the interactions, the density matrix at time t is in the interaction picture given by the S-matrix ( evolution operator):

$$\rho \left(t\right)=S\left(t,-\infty \right)\rho \left(-\infty \right)S^{†}\left(t,-\infty \right)$$

(40)

$$S\left(t,t^{\prime }\right)=\mathrm{T}exp\left[-i/\hslash {{\displaystyle \int }}_{t^{\prime }}^{t}dt{}^{″} H_{int}\left(t{}^{″}\right)\right]$$

(41)

$$S^{†}\left(t,t^{\prime }\right)=\mathrm{T}exp\left[-i/\hslash {{\displaystyle \int }}_{t^{ }}^{t\prime }dt{}^{″} H_{int}\left(t{}^{″}\right)\right]$$

(42)

The time ordering is along the real-time integration paths.

$$<\mathsf{\Xi }\left(t\right)>=Tr\left[\mathsf{\Xi }\left(t\right)\rho \left(t\right)\right]=Tr\left[S^{†}\left(t,-\infty \right)\mathsf{\Xi }\left(t\right)S\left(t,\infty \right)\rho \left(-\infty \right)\right]$$

(43)

In a zero-temperature many-body systems with a stable ground state |0> (vacuum state) it is possible to reduce Equation (32) to:

$$<\mathsf{\Xi }\left(t\right)>=\frac{<0|\mathrm{T}S\left(-\infty ,\infty \right)\mathsf{\Xi }\left(\mathrm{t}\right)|0>}{<0|S\left(-\infty ,\infty \right)|0>}$$

(44)

where the time-integrals are along the real time axis from −∞ to ∞. The familiar many-body perturbation theory then follows using Wick’s theorem [71] that provides cancellation of divergent terms in the perturbation expansion of the numerator of Equation (33). This simplification does not hold for non-equilibrium systems. Instead, in the Keldysh formalism Equation (33) is written as

$$<\mathsf{\Xi }\left(t\right){>}_{0eq}=Tr\left[{\mathrm{T}}_{\mathrm{contour}}{S}_{\mathrm{contour}}\mathsf{\Xi }\left(t\right){\rho }_{0eq}\right]$$

(45)

where the S-matrix, S_contour, is evaluated along a time contour C _Keldysh that runs from $-\infty$to t (C⁺) and then back from t to −∞ (C⁻).

Figure 1 displays the Keldysh contour showing the arguments (red) of the contour ordered Green function G^<(t,t’). T_contour is the time ordering operator on the contour C_contour = C_Keldysh. By inserting the identity operator S(t,∞)S(∞,t) the Keldysh contour is extended from (−∞,∞) and back again: (∞,−∞).

In Equation (45) ρ_0eq is the equilibrium density matrix for non-interacting excitations. The subsequent time-dependent many-body infinite perturbation theory then follows a similar Feynman diagram expansion as for zero temperature theory but involves integrations along the Keldysh contour. The perturbation theory is most easily derived using equations of motion for a hierarchy of Green functions.

In the steady-state, neglecting initial state correlations and transients, and choosing a time-independent external perturbation H_ext, the basic single-particle electron Green functions may be expressed in the Heisenberg picture:

$$\mathsf{\Xi }\left(t\right)=\exp \left[\frac{i\left(H_0+H_{ext}+H_{int}\right)\left(t-t_0\right)}{\hslash }\right]\mathsf{\Xi }\left(t_0\right)\exp [\frac{-i\left(H_0+H_{ext}+H_{int}\right)\left(t-t_0\right)}{\hslash }$$

(46)

where now the Hamiltonian H has been augmented to include the perturbing interaction H_ext.

In the case of the single particle electron Green function G(r,r’,t,t’) we find:

$${\mathrm{G}}_{ }\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right)\equiv -i<{\mathrm{T}}_{\mathrm{Keldysh}}\psi \left(\mathrm{r},t){\psi }^{†}({\mathrm{r}}^{\prime },t^{\prime }\right){>}_{0eq}$$

(47)

where as before the terms ψ(r,t), ψ^†(r,t) are annihilation, creation operators in the Heisenberg picture Equation (46) for the electron quantum fields, time-ordered by T_Keldysh (the Dyson time-ordering operator on the Keldysh contour). Equation (47) may be partitioned into four piece-wise defined Green functions chosen according to the location of the time variables t,t’ on the Keldysh contour. A common choice, which is adopted here, is closest to the Kadanoff-Baym definitions (a better choice for constructing equations of motion and self-energies is detailed by Lifshitz and Pitaevskii [70]):

$${\mathrm{G}}_{}^{<}\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right)=\mathrm{G}\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right) (t<t^{\prime })=i<{\psi }^{†}\left(\mathrm{r}\prime ,t\prime ) {\psi }^{ }({\mathrm{r}}^{ },t^{ }\right){>}_{0eq}$$

(48)

$${\mathrm{G}}_{}^{>}\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right)=\mathrm{G}\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right) (t>t^{\prime })$$

(49)

$${\mathrm{G}}_{ }^R\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right)=-\mathrm{i}\mathsf{\theta }\left(t-t^{\prime }\right)<\left\{{\psi }^{†}\left(\mathrm{r}{}^{\prime },t\prime ), {\psi }^{ }({\mathrm{r}}^{ },t^{ }\right)\right\}{>}_{0eq}$$

(50)

$${\mathrm{G}}_{ }^A\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right)=\mathrm{i}\mathsf{\theta }\left(t\prime -t\right)<\left\{{\psi }^{ }\left(\mathrm{r},t), {\psi }^{†}({\mathrm{r}}^{\prime },t^{\prime }\right)\right\}{>}_{0eq}$$

(51)

The Heisenberg equations of motion for these Green functions are coupled differential—integral equations that may be obtained by differentiating with respect to the two variables t and t’ (and r,r’). These equations immediately show up coupling to the two-electron Green functions that in turn are coupled to higher order Green functions to form a well-known hierarchy. Within systematic perturbation theory the hierarchy is curtailed by introducing self-energy terms that close off the hierarchy by simplifying the higher order Green functions. Under these conditions the equations of motion for the piece-wise time sub-domain Green functions are coupled but together form a closed system of equations. Using the Langreth theorem [40], which connects the Keldysh contour ordered Green functions to the real-time Green functions G^<, G^>, G^R, G^A, all the relevant Green functions satisfy a Dyson-like form with a suitable self-energy Σ(r,r’,t,t’):

$$\begin{gathered} {\mathrm{G}}^{\mathsf{\lambda }}\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right)={\mathrm{G}}_0^{\mathsf{\lambda }}\left(\mathrm{r},{\mathrm{r}}^{\prime },t,t^{\prime }\right)+ \\ {{\displaystyle \iint }}^{ }{{\displaystyle \iint }}^{ }{\mathrm{d}}^3\mathrm{r}{}^{″}{\mathrm{d}}^3\mathrm{r}{}^{‴}{\mathrm{d}}^{ }t{}^{″}{\mathrm{d}}^{ }t{}^{‴}\left\{{\mathrm{G}}_0^{\mathsf{\lambda }}\left(\mathrm{r},{\mathrm{r}}^{″},t,t^{″}\right)\mathsf{\Sigma }\left(\mathrm{r}{}^{″},{\mathrm{r}}^{‴},t{}^{″},t^{‴}\right]{\mathrm{G}}^{\mathsf{\lambda }}\left(\mathrm{r}{}^{‴},{\mathrm{r}}^{\prime },t{}^{‴},t^{\prime }\right]\right\} \end{gathered}$$

(52)

Here, ${\mathrm{G}}_{0 }^{\lambda }$is the relevant unperturbed Green function (H_int = 0).

Because the piece-wise Green functions are not independent, the self-energies are functionals of the other self-energies and Green functions.

Considering now a non-equilibrium steady-state, energy-resolved picture let us transform to new variables τ_relative = t–t’; τ_mean = (t + t’)/2; then performing a Fourier transform of the Green functions over the relative time variable τ_relative yields the same forms described in elementary fashion in Section 2.2:

$$G^{<}\left(\mathrm{r},\mathrm{r}{}^{\prime },\epsilon \right)=\left(\frac{1}{\hslash }\right){{\displaystyle \int }}_{-\infty }^{\infty }{\mathrm{d}\mathsf{\tau }}_{relative}\exp \left[\frac{i\epsilon {\tau }_{relative}}{\hslash }\right]G^{<}\left(\mathrm{r},\mathrm{r}{}^{\prime },{\tau }_{relative},{\tau }_{mean}\right)$$

(53)

(where there is no dependence on τ_mean in the steady state) and similarly for G^>, G^A, G^R. The Dyson Equation (52) become matrix equations in the space variables and the self-energies become convolution integrals over intermediate energies.

2.3.2. The Steady State NEGF Equations

For the device Hamiltonian comprising an uncoupled term, an interaction term H_int coupling electrons to phonons, and an applied perturbation H_ext resulting from the applied electric fields and the confinement potential:

$$H=H_0+H_{ext}+H_{int}$$

(54)

The resulting steady-state coupled Green function equations are:

$${\mathrm{G}}^R={\{\epsilon +i\eta -H_0-H_{ext}-{\mathsf{\Sigma }}^R\}}^{-1}={\mathrm{G}}_0^R+{\mathrm{G}}_0^R{\mathsf{\Sigma }}^R{\mathrm{G}}^R$$

(55)

$${\mathrm{G}}^A={\mathrm{G}}^{R†}$$

(56)

$${\mathrm{G}}^{<}={\mathrm{G}}^R{\mathsf{\Sigma }}^{<}{\mathrm{G}}^A; {\mathrm{G}}^{>}={\mathrm{G}}^R{\mathsf{\Sigma }}^{>}{\mathrm{G}}^A$$

(57)

Equation (55) is essentially a Dyson equation involving the self-energy ${\mathsf{\Sigma }}^R$ which is a functional of the other Green functions and self-energies. Equation (57) is deceptively simple but it is the analogue of the classical Boltzmann equation as discussed briefly in Appendix C. Indeed from Equations (55) and (57) we have

$$\{\epsilon +i\eta -H_0-H_{ext}{\}\mathrm{G}}^{<}={\mathsf{\Sigma }}^R{\mathrm{G}}^{<}+{\mathsf{\Sigma }}^{<}{\mathrm{G}}^A$$

(58)

2.3.3. The Projection Theorem

So far we have tacitly treated an open system. A key requirement in NEGF theory applied to device modelling is to project from the open system comprising a device coupled to its environment of contacts and thermal reservoirs into a finite simulation domain (see Appendix B, Figure A3). This non-trivial matter was first solved by Caroli et al in 1971 and their derivation of the theorem is sketched in the Appendix B. The main result is that the finite device domain may be modelled by the device Hamiltonian H_W (W for nanowire) provided that an additional set of non-Hermitian self-energies be added to H_W for each relevant Green function. These contact self-energies ${\mathsf{\Sigma }}_{contacts}^{<},{ \mathsf{\Sigma }}_{contacts}^{>},{\mathsf{\Sigma }}_{contacts}^R,{\mathsf{\Sigma }}_{contacts}^A$are localised interactions on each contact-device boundary surface that describe the injection or extraction of carriers (and consequently the inflow and outflow of charge, current and energy).

2.3.4. Numerical Solution of NEGF Equations

Equations (55)–(57) may be evaluated in the position representation, by inserting complete sets of position eigenstates; if the position variables are discretised the result is a set of matrix equations. The problem is significantly simplified if the self-energies are diagonal. We further require expressions for H_ext and the self-energies ${\mathsf{\Sigma }}_{ }^{\lambda },{\mathsf{\Sigma }}_{contacts}^{\lambda }$.

In solving the NEGF equations, there are three other equations to be solved self-consistently:

(i) The scattering self-energies must be evaluated self-consistently in at least the self-consistent Born approximation (fast algorithms for this and higher order perturbation are described in [75]. Appendix D gives an example of ${\mathsf{\Sigma }}_{}^{<},{\mathsf{\Sigma }}_{ }^R$in the SCBA for a generic electron-phonon interaction.

(ii) Poisson’s equation for the self-consistent applied potential and confinement potential (necessary for image charge induced electrostatic self-energies [10]) must be evaluated.

(iii) If the lattice is considered to be held at constant ambient temperature, there is no further issue; if not, the heat equation(s) driven by the heat dissipation into the phonon system must be considered self-consistently. In the most advanced modelling, this latter step has been replaced by coupling the above equations self-consistently to the phonon Green functions [50,75]. A few groups use local atomic or molecular orbitals as the basis for tight binding style Hamiltonian models, but these techniques are restricted to a few hundred atoms [76]. The effective mass Hamiltonian method [76,77] is derived from tight binding and is discussed further in Section 3.

2.3.5. NEGF Electro-Thermal Modelling

The NEGF equations are very powerful compared to Boltzmann transport theory because they are able to describe quantum phenomena such as tunnelling, resonances, bound states, fully quantised scattering theory, coherent and non-coherent flows. The density of states particularly LDOS can then reveal shift and broadening in the energy levels) due to the scattering self-energies/ contact self-energies) and hot spots in the spatially-resolved energy current densities. In the latter case we may immediately recover the quantum equivalent of the classical Equations (11) and (12):

The NEGF equivalents of the electro-thermally important quantities described in Section 2.1 are as follows.

The energy- and space- resolved charge density involves only the diagonal part of G^<(r,r^′,ε) (following from the limit (r→r^′, t→t^′), of Equations (47) and (48))

$$n\left(\mathrm{r},\epsilon \right)=-i\left(\frac{e}{2\pi }\right){\mathrm{G}}^{<}\left(\mathrm{r},\mathrm{r}{}^{\prime },\epsilon \right){|}_{\mathrm{r}\prime \to \mathrm{r}}$$

(59)

For a general Hamiltonian, H the velocity operator is defined in the Heisenberg picture by $\mathrm{v}\equiv \mathrm{dr}/\mathrm{d}t=\left(\frac{i}{\hslash }\right)\left[H,\mathrm{r}\right]$.

For a simple effective mass Hamiltonian, the (electrical) energy-resolved current density is:

$$\mathrm{J}(\mathrm{r},\epsilon )=e\mathrm{j}(\mathrm{r},\epsilon )={\left.\left(-\mathrm{e}/4\pi m^{*}\right)\left(\nabla -{\nabla }^{\prime }\right)G^{<}\left(\mathrm{r},{\mathrm{r}}^{\prime },\epsilon \right)\right|}_{{\mathrm{r}}^{\prime }\to \mathrm{r}}$$

(60)

$${\mathrm{J}}_{energy}\left(\mathrm{r}\right)={\displaystyle \int }\epsilon \mathrm{j}\left(\mathrm{r},\epsilon \right)\mathrm{d}\epsilon$$

(61)

The local power transfer to the lattice is

$$\mathrm{P}\left(\mathrm{r}\right)=\nabla ·{\mathrm{J}}_{energy}\left(\mathrm{r}\right)$$

(62)

or, separating out the Joule power,

$$\mathrm{P}\left(\mathrm{r}\right)=-\nabla ·{\mathrm{J}}_{kinetic}\left(\mathrm{r}\right)+\mathrm{E}\left(\mathrm{r}\right)·\mathrm{J}\left(\mathrm{r}\right)$$

(63)

where the first term on the RHS of Equation (63) is the divergence of the electron kinetic energy current density, while the second term is the product of the local electric field and the current density and represents the local Joule heat. The LHS is the net energy flow from the electron gas to the phonon assembly and is the source term for a heat equation (s). Equation (63) is derived by writing the kinetic energy ${\epsilon }_{kinetic}\left(\mathrm{r},\epsilon \right)=\epsilon -e\varphi \left(\mathrm{r}\right)$. The kinetic energy current density is then

$${\mathrm{J}}_{kinetic}\left(\mathrm{r}\right)={\displaystyle \int }{\epsilon }_{kinetic} \mathrm{j}\left(\mathrm{r},\epsilon \right)\mathrm{d}\epsilon$$

(64)

and

$${\mathrm{J}}_{Energy}\left(\mathrm{r}\right)={\mathrm{J}}_{kinetic}\left(\mathrm{r}\right)+\varphi \left(\mathrm{r}\right)\mathrm{J}\left(\mathrm{r}\right)$$

(65)

Taking the divergence of Equation (65) and using $\mathrm{E}\left(\mathrm{r}\right)=-\nabla \varphi \left(\mathrm{r}\right)$ plus the continuity equation result $\nabla ·\mathrm{J}\left(\mathrm{r}\right)=0$, recovers Equation (63).

Using Equations (55)–(57), we obtain an expression for the divergence of the total energy current density

$$\begin{gathered} \nabla ·{\mathrm{J}}_{Energy}\left(\mathrm{r}\right)={\displaystyle \int }\mathrm{d}\epsilon \epsilon {\displaystyle \int }{d}^3{r}^{\prime }\left\{\right\{{\mathrm{G}}^R\left(\mathrm{r}, \mathrm{r} {}^{\prime },\epsilon \right)-{\mathrm{G}}^A\left(\mathrm{r}, \mathrm{r} {}^{\prime },\epsilon \right)\}{\mathsf{\Sigma }}^{<}\left( \mathrm{r} {}^{\prime },{\mathrm{r}}^{ },\epsilon \right)+ \\ { \mathrm{G}}^{<}\left(\mathrm{r}, \mathrm{r} {}^{\prime },\epsilon \right){\mathsf{\Sigma }}^A\left(\mathrm{r}{}^{\prime },{\mathrm{r}}^{ },\epsilon \right)-{\mathsf{\Sigma }}^R\left(\mathrm{r}, \mathrm{r} {}^{\prime },\epsilon \right){\mathrm{G}}^{<}\left(\mathrm{r}{}^{\prime },{\mathrm{r}}^{ },\epsilon \right)\} \end{gathered}$$

(66)

$$\nabla ·{\mathrm{J}}_{Energy}\left(\mathrm{r}\right)={\displaystyle \int }\mathrm{d}\epsilon \epsilon <\mathrm{r}|\{{\mathrm{G}}^R\left(\epsilon \right)-{\mathrm{G}}^A\left(\epsilon \right)\}{\mathsf{\Sigma }}^{<}\left(\epsilon \right)+{\mathrm{G}}^{<}\left(\epsilon \right){\mathsf{\Sigma }}^A\left(\epsilon \right)-{\mathsf{\Sigma }}^R\left(\epsilon \right){\mathrm{G}}^{<}\left(\epsilon \right)|\mathrm{r}>$$

(67)

The operator within the position matrix element of Equation (67) is related to the current operator as discussed in the Appendix D, Appendix E and Appendix F.

Since the self-energies here include coupling to the contacts it is possible to use Equation (66) to locate precisely where the energy dissipation occurs (for example, within particular device regions or at the contacts).

The above equations plus explicit forms for the electron-phonon interaction self-energies (see Appendix D) give the quantum equivalent of the classical electro-thermal formulation of Section 2.1, i.e., Equations (11) and (13). It is important to demonstrate that expression $\nabla ·{\mathrm{J}}_{Energy}\left(\mathrm{r}\right)$vanishes for elastic scattering. This is demonstrated in Appendix E, where $\nabla ·{\mathrm{J}}_{Energy}\left(\mathrm{r}\right)$is shown explicitly to vanish everywhere within the device volume using the self-energies for high temperature elastic acoustic phonon scattering in the self-consistent Born approximation. For the elastic case, there is no energy dissipation within the channel but when the contact self-energies are introduced there may still be energy dissipation in the drain contact region.

3. Simulation Methodology for the Silicon Gate-All-Around Nanowire FET

In this section we describe the carrier transport model and our electrothermal model used in the simulation of a gate-all-around nanowire transistor. In a nutshell, our method simultaneously solves the NEGF Equations (55)–(57), the Poisson Equation (3) and heat Equation (15) with the appropriate boundary conditions.

We start by estimating the electrostatic potential and use it to solve the NEGF Equations (55)–(57) for a range of energies; these two equations need to be solved self-consistently as electron-phonon scattering self-energies depend on the electron density. From the calculated Green functions, the electron density is obtained and inserted into the Poisson equation to calculate the new electrostatic potential. At the same time, the current and Green functions are used to calculate the local power dissipation that is the input into to the heat equation in order to find the temperature profile. This process continues until density, current and temperature do not change. In what follows, we will describe in detail the different physical models, boundary conditions and simplifications used in solving the equations.

3.1. Device Structure and the 1D Representation

In this work, the effective mass approximation [78,79,80] has been used to approximate the non-interacting Hamiltonian H₀ of the silicon crystal entering in Equation (54). The effective masses used have been extracted from tight-binding calculations, see reference [80]. In addition, our model for Si considers anisotropic bands and includes the six $\Delta$-valleys, see Appendix A and Figure A2 for the k-space location of the valleys. The axis of the wire is oriented in the [100] direction, for a square cross section, 4 ellipsoids are equivalent as shown in Figure A2.

In devices with cylindrical or rectangular symmetry the spatial dependence of the Green functions may be reduced to 1D forms (in, x,x’) by projecting out the transverse (y,z) dimensional dependence using self-consistent transverse eigenstates (discrete sub-bands n) computed at each point along the x-direction (see Section 2.3 and references therein). The full Green functions are then replaced by effective 1D Green functions: $\mathrm{G}\left(x,x^{\prime },\epsilon \right).$ This simplification is a key result for the efficient modelling of nanowire structures.

The device (typified by Figure 2) is divided in cross-sections perpendicular to the axis of the wire. The standard Schrödinger equation is solved in each cross-section by using the electrostatic potential calculated by the Poisson equation. This results in a set of eigenstates for each valley and cross-section considered. These eigenstates are commonly called sub-bands or modes and are associated with individual conduction band valleys. The longitudinal Hamiltonian is expanded in these sub-band states that are coupled to each other by the relevant overlap integrals. This method is called the coupled mode space method [81,82,83,84] and is very efficient for very narrow nanowires as the number of sub-bands considered can be kept to a minimum. In detail, the full volume of the device simulation domain is discretised in a rectangular mesh; the Schrödinger equation is then solved for each cross-section perpendicular to the wire axis. The potential energy entering in the equation is the electrostatic potential energy calculated from the Poisson equation and the confining potential at the SiO₂ interface. The kinetic energy is given by the anisotropic mass tensor corresponding to the particular valley. From the solution of the Schrödinger equation, a set of eigenvalues and eigenvectors (wavefunctions) are obtained for each cross section. These eigenvectors or modes are used to expand the 3D spatially discretized Green function [84] of size (n_x × n_y × n_z, n_x × n_y × n_z), where n_y and n_z is the number of discretisation points in the cross-section direction and n_x in the transport direction. From these modes and by integration in the transversal spatial coordinates, an equation for the Green function is constructed that couples the modes in one cross section with the modes of the nearest neighbouring cross-sections [83,84]. The size of this Green function matrix is of order (n_m × n_x, n_m × n_x), where n_m is the total number of modes used in the simulations. Note that the 3D spatially resolved Green function is substantially larger in size than the mode space Green function as mentioned before. In the discretised Hamiltonian, the matrix element (or coupling) between mode m in cross section i and mode l in cross section i + 1 is given by:

$$-\frac{{\hslash }^2}{2m_x{\left(\Delta x\right)}^2}{\displaystyle \iint }{\psi }_i^{*m}(y,z){\psi }_{i+1}^l\left(y,z\right)dydz$$

(68)

where ${\psi }_i^m\left(y,z\right)$ is the wavefunction for mode m in cross section i, the y and z are the spatial coordinates representing the cross section. $\Delta x$ is the discretisation step in the longitudinal direction or the transport direction. The asterisk * in the above expression stands for the complex conjugate. The double integral is over the whole cross-section. Equation (68) shows that the probability of an electron going from mode m in cross section i to mode l in cross-section i + 1 depends on the overlap between wavefunctions. It should be noted that if the potential energy in the cross-sections are similar the above integral between two different modes is zero. This means that there is a very small probability (or amplitude) for the electrons during their motion along the wire to switch modes. When the integral between different modes is assumed zero, the simulation method is called the uncoupled mode space approach. In this case the Green function in the mode/longitudinal space reduces to n_m Green functions of size (n_x $\times$ n_x). This reduces the computational load substantially. After the longitudinal Green function is obtained for each valley, mode and energy, the electron density and current is calculated from them and the mode wavefunctions. For details of the coupled mode space as applied to NEGF, we recommend the following references [83,84].

We have considered the acoustic and optical (for f-type and g-type) intervalley scattering mechanisms, as well as the intra-valley elastic acoustic phonons scattering mechanism. The specific parameters and models used can be found in references [85,86]. The electrostatic potential is calculated from the Poisson equation in the mean field approximation. We have used Neumann boundary conditions in the source and drain contacts; and Dirichlet boundary conditions at the gate contact.

Figure 2 shows the schematics of a typical wrap-round gate silicon nanowire with a silicon core comprising source, channel and drain regions. There is a high degree of symmetry. Therefore, we may use the sub-band projection methodology to obtain the 1D density of states throughout the silicon core as a function of energy and position as indicated.

Figure 3 shows the detailed 1D local density of states (LDOS) computed along the nanowire axis for each sub-band energy: it is proportional to (ε−ε_n)^−1/2 where ε_n is the local sub-band energy. Notably, the multiple reflection patterns (a simple example is in Appendix A) that are shown for energies lower that the top of the maximum sub-band energy around the middle of the channel. This is an indication of coherent transport. The wavelength of the electrons increases as the energy increases; this can be confirmed by the shorter separation of the maxima and minima in the figure.

3.2. Non Equilibrium Green Function Formalism: Density and Current Density

Equations (41) and (42), used for the calculation of the Green function, implicitly contain Pauli’s exclusion principle [57,58,87]. There, ${\mathsf{\Sigma }}^{<}$ is a sum of two terms, one is proportional to the electron-phonon scattering rate and the other representing the boundary self-energies that are different from zero only at the contacts. This latter term involves the product of the electron distribution of the contact and the density of states at the contact. This approximation in the contacts is a consequence of the fluctuation-dissipation theorem [88]. In one dimension, for each mode, the energy-resolved electron concentration and current density could be calculated from the mode-space Green functions by [83,84]:

$$n_l\left(x,\epsilon \right)=-\left(\frac{i}{2\pi }\right){\mathrm{G}}_{ll}^{<}\left(x,x,\epsilon \right)$$

(69)

$$J_{nl}\left(x,\epsilon \right)=-\frac{e}{4\pi m^{*}}{\displaystyle \underset{x^{\prime }\to x}{\lim }} \left(\frac{\partial }{\partial x}-\frac{\partial }{\partial x^{\prime }}\right){\mathrm{G}}_{nl}^{<}\left(x,x^{\prime },\epsilon \right)$$

(70)

where x represents the spatial coordinate, e the electron charge snd m* the effective mass. The n and l in the RHS label the modes. Equations (69) and (70) must be modified by a factor of two for spin degeneracy. The total 3D electron density is calculated from the above 1D electron density and the cross-section wave function by:

$$n\left(\mathrm{r},\epsilon \right)={\displaystyle {\displaystyle \sum }_l}{n}_l\left(x,\epsilon \right){\psi }_x^{*l}\left(y,z\right){\psi }_x^l\left(y,z\right)$$

(71)

Figure 4 shows the electron density at low and high gate bias, showing electrical neutrality, i.e., the potential is flat at the contact and the electron concentrations are similar in the source and drain contact. The electron density drops close to the oxide interface showing the volume inversion and the shape of the ground state transversal wave function. For narrow-body nanowires, this concentration of electrons along the wire axis is usually called the volume inversion or electron confinement. At very high gate voltage, the source-drain barrier disappears and, in order to fulfil electrical neutrality, the leads or ends of the source and drain region should be broadened (see Landauer [59]). At high gate bias, the electron potential energy along the wire direction resembles a step function or a linear drop with an energy difference proportional to the applied bias, so the right moving waves from the source reach the drain but the left moving waves from the drain are reflected by the energy drop. This means that the density of electrons in the drain will be always larger than in the source and therefore electrical neutrality would not be fulfilled. However, if simulations are carried out with broadened leads, electrical neutrality in the source/drain are independent of the gate potential and are controlled by the electrostatic of the leads.

Figure 5 shows the energy resolved current spectra calculated by Equations (69) and (70) along the nanowire transistor axis. The source contact is in the left end and the drain is at the right end. The first sub-band is shown as reference. The part of the current with energy lower than the top of the barrier sub-band energy represents tunneling current. As the device has a 15 nm channel length, the tunneling current is not negligible. From the figure it can be seen that electron at the source moving from the source towards the drain comes from two different energy regions:

(i)

energy close to the top of the barrier and

(ii)

energy close to the sub-band energy at the source contact.

Electrons entering the source region (i) will be able to go through the barrier, unless they emit phonons and lose energy. Those in region (ii) need to absorb phonons in order to cross the barrier and in doing so will cool the lattice locally. The temperature profile simulations presented in the next section confirm this statement. At the drain side, electrons injected from the source (at high energy) start to dissipate part of the energy and the current at the drain is more spread downwards in energy compared with the current at the top of the barrier. However, the average energy of the electron at the drain end is approximately equal to 0.25 eV that is substantially larger than the sub-band energy (−0.55 eV).

3.3. Energy Balance and Power

The Equations (61)–(63) represent the energy balance between electrons and phonons and is discussed in several references [89,90,91,92].

These equations represent the total local power dissipated for the electron system into the phonon system. As indicated before, the last term on the RHS of Equation (63) is the local Joule power that integrated through the device produces the well-known Joule power I·V. In Equation (63) the term $\mathrm{E}\left(\mathrm{r}\right)$ represents the local electric field along the nanowire axis. The first term on the RHS is the change in the “kinetic energy” current of the electron ensemble, if this kinetic energy current does not change it means that the electron system is adiabatically following the “bending” of the potential and therefore it is not increasing the electron mean energy relative to the conduction band minima, i.e., no hot electron phenomena occur. In this case, the electron system or more specifically the electrons taking part in the current are in equilibrium with the local potential. On the other hand, a large change in kinetic energy current implies the heating or cooling of the electron system.

Equation (63), for the local power, can be simplified for the case of local electron-phonon self-energies and the resulting equation obtained [56] is essentially a detailed balance equation (see Appendix E, F specificaly Equation (A46)):

$$P\left(\mathrm{r}\right)={\displaystyle \int }\mathrm{d}\epsilon \epsilon \left\{{\mathrm{G}}^{<}\left(\mathrm{r}, \mathrm{r},\epsilon \right){\mathsf{\Sigma }}^{>}\left(\mathrm{r}, \mathrm{r},\epsilon \right)-{\mathrm{G}}^{>}\left(\mathrm{r}, \mathrm{r},\epsilon \right){\mathsf{\Sigma }}^{<}\left(\mathrm{r}, \mathrm{r},\epsilon \right)\right\}$$

(72)

The electron-phonon self-energies used in this paper are also treated in a local approximation (diagonal approximation), so that Equation (72) is the local version of Equation (A46) in the Appendix F. The RHS of Equation (72) expresses the local power as a function of the Green functions and self-energies. It shows that the energy mismatch between the rate of $\epsilon$ -outgoing (leaving the state of energy$\epsilon$) electrons (the first term in the RHS) and the rate of $\epsilon$ -outgoing holes (second term in the RHS) is transferred to the phonon system. Note that for elastic scattering the expression in the curly brackets became zero [91] (Appendix E demonstrates this explicitly). As mentioned before, the first term inside the integral in the RHS represents the outgoing electrons from state$\epsilon$ and the second term is the ingoing electrons into state$\epsilon$. For elastic scattering these rates are equal as the electron does not change its energy after scattering and P(r) vanishes except in the contacts. Equation (72) represents the detailed energy balance between the electron and the phonon system.

Figure 6 shows the kinetic term and the Joule term of Equation (63). Their sum or the total power is also shown. Note that the kinetic term is plotted with reverse sign. The two terms are very similar but with opposite sign therefore in order to obtain an accurate value of the total power, the current needs to be calculated with great accuracy which implies a large number of Born iterations for the self-energies.

The 3D Fourier heat equation used to calculate the local temperature profile is solved in the oxide and in the Si nanowire core. Neumann boundary conditions are used in the source and drain contacts and the Dirichlet boundary condition is used at the metal gate. The heat sources are calculated from the total power or using the Joule heat power, depending of with model we are using. Both sources include the volume inversion effects.

3.4. Results and Discussion

In this section, we report electrothermal simulations of a gate all around nanowire transistors (Figure 2). The cross section of the silicon core is 2.2 × 2.2 nm² and the source, channel and drain have lengths of 20 nm, 15 nm and 20 nm, respectively. The channel is undoped and the doping in the source and the drain are 10²⁰ cm⁻³. In all of this work, we assumed a drain bias of V_d = 0.6 V. We have compared the drain current vs gate voltage characteristic (transfer characteristic) for the case of standard scattering, i.e., assuming the phonon system or the lattice is at equilibrium at 300 K with cases in which the lattice is at higher temperature.

The local phonon temperatures are calculated from Fourier’s law. Two cases have been considered: case (i): in which the electrons are in local equilibrium with the electrostatic potential, so the energy transfer from the electron system to the phonon system is given by the standard Joule law or the second term in Equation (63); and case (ii), in which the electrons relax to equilibrium depending on the energy relaxation length (given by the inelastic electron-phonon scattering length) [92] and therefore the energy transfer are given by the LHS of Equation (63). This latter case considers that the electrons injected from the source into the drain leave the drain at a higher energy than the drain average (equilibrium) energy. These electrons are commonly called “hot electrons”. Consequently, these electrons dissipate less energy inside the device as compared with the electrons in local equilibrium and therefore produced less heating of the lattice in the active part of the device.

In a previous paper [92], we have simulated and examined different sizes of source and drain extensions for a nano-transistor with a similar cross section to the one studied here but self-heating was not considered. In that paper, we found that electrons needed a drain extension length larger than 100 nm in order to guarantee that the total power dissipated inside the device became equal to the Joule power.

However, even for larger source and drain extensions as long as the channel length is around 10 nm, inside the device, the potential energy of the electrons changes hundreds of millielectron volts in just a few nanometers. This length is too short compared with the energy relaxation length. As a consequence, the electron transport around the gate is essentially out of equilibrium or carried out by “hot electrons”. The energy relaxation depends on the density of states and on the electron-phonon scattering rates. However, previous papers [76] for small cross section devices confirm an increase in electron-phonon scattering rates due to an increase in scattering form factors. This fact favours a decrease in the energy relaxation length in small devices but this is not sufficient to allow full electron energy relaxation inside the device.

In general, devices with cross sections larger than 10 × 10 nm² and with channel/source/drain of a few 100 of nanometers are expected to comply with the Joule law but as mentioned previously a more accurate assessment would look at the local spatial variation in the electrostatic potential or electron potential energy, which depends on the device structure and dimensions.

The current-voltage characteristic for these three scattering simulation cases: (i) no self-heating (NSH), (ii) hot electrons self-heating (HSH) and (iii) Joule self-heating (JSH) are shown in Figure 7.

As it is clear from Figure 7 that under self-heating conditions, i.e., in which the lattice has a high temperature, the scattering rate is large and the current at Vg = 0.8 V is 30% smaller for the hot electrons case as compared with the non-self heating case. However, when the Joule self-heating is considered the effect is larger, and the drain current is reduced by 70%, compared to the NSH case. For low gate bias, as the current is small, the impact of self-heating and scattering in general is diminished. In what follows, we will discuss the results at low gate bias and large gate bias.

At low gate bias the current is low and consequently the power dissipated is low; however, the fact that we have a large source-drain bias barrier produce some interesting effects. Figure 8 and Figure 9 show the current spectra at Vg = 0.2 V for the case of the hot electrons. Figure 9 shows an enlargement of the high energy region of Figure 8 where the electrons energy exchange occurs.

The figures show that the electrons at the bottom of the band entering in the source absorb phonons and therefore are able to cross the barrier, this phonon mediated effect tends to transfer energy from the source to the drain of the device, and therefore it has a cooling effect on the source of the device and slight heating at the drain.

As the nanowire cross-section is small the LDOS is quasi-1D and therefore there is a relatively large DOS at the sub-band energy at the source end (and indeed at the drain end); consequently, there is an overwhelming number of electrons at that low energy. This induced a high probability of phonon absorption for the electrons in the source end. This is one of the reasons that the energy-resolved current at the source is split into two components for a large source-drain barrier.

Figure 10 shows the temperature profile across a plane of the device axis including the oxide at Vg = 0.3 V. The figure shows a slight cooling (from 300 K) of the source of the device. This effect is slightly enhanced at Vg = 0.3 V as the height of the barrier has decreased and therefore the electrons are more efficiently transferred lower energy in the source to the top of the barrier. The cooling of the source diminished at higher gate bias as a large fraction of electrons are able to tunnel or directly transit the barrier without the need of phonon assistance. This indicated an optimum length in which the cooling of the source is more efficient. At very low gate bias we have found that the total power dissipated by the electron current in the device is negative; this means that the electrons absorbing energy from the phonon system on the source are not able to release enough energy at the drain to counteract the cooling effect at the source. This net energy absorption of electrons in the source region has been observed in other papers [86,91] through the calculation of the power absorbed along the wire for the electrons using the Equation (9) as will be shown in the Figure 11, here. As mentioned before, it can also be inferred by the plot of the energy current resolved along the nanowire axis (see Figure 8 and Figure 9). However, Figure 10 shows a decrease in temperature (cooling) of the lattice, which depends on the net interchange of heat of the lattice with the surrounding, as it is not sufficient that the electron system absorb energy from the lattice locally in order to cool the lattice. If the drain or other parts are hot enough, heat will be transfer to the source avoiding cooling. The fact that the electron going through the drain does not release enough energy there is one of the facts which allows for the cooling of the source. If we had assumed that the electrons are in equilibrium with the potential or, i.e., assuming Joule heating, the net power dissipation would be always positive, that is a power transfer from the electrons to the lattice. In addition, there would have been a large dissipation in the drain, heating the overall drain and surroundings including the source.

The local power density profile used as a source in the heat source is shown in Figure 10 along the nanowire axis for the case of hot electrons. It is shown that most of the power transfer occurs close to the source-drain barrier. The negative power density at the source end of the barrier indicates that the electron system is taking energy from the lattice or phonon system as described before.

At high gate voltage, the source-drain barrier is small or negligible. Figure 12 and Figure 13 show the profiles of the heat source term for Vg equal to 0.6 V and 0.7 V, respectively. At Vg = 0.6 V, there is a small region in which the heat source has negative values: in this region the electrons are still absorbing energy from the lattice; however, there is a greater heat transfer from the electrons to the lattice at the drain region and as a consequence the overall temperature increases. When the barrier is sufficiently low, as is in the case of Vg = 0.7 V (Figure 13), there is no region of negative-valued heat source; in this case, the potential resembles a step potential and the electrons just release energy when transiting the device. The energy-resolved current density for Vg = 0.6 V and 0.7 V is shown in Figure 14 and Figure 15, for the case when self-heating and hot electron effects are considered. The corresponding temperature profiles are shown in Figure 16 and Figure 17. Note that although the current spectra look very similar, the change in maximum temperature is approximate 60 K. Figure 17 resembles the 1D temperature plots of Figure 7 in reference [93]. We have carried out a more thorough comparison in our previous paper [94]. However, it is worth mention that our nanowire has rectangular cross section and the nanowire in [93] has a circular cross section (not exactly matching the areas). We are using an effective mass approximation, and the reference [93] used full band. The scattering models are also different as they used the full phonon bands. The reference [93] simulated the phonons using a NEGF formalism instead of our simplified heat equation. However, the level of current obtained in our simulations are of similar order of magnitude to those in [93]. Considering all the differences between the models, our calculated temperature profiles qualitatively follow those found in [93]. The temperature profiles there range from 300–500 K and this also matches the range of our results.

The heating of the lattice due to the energy transfer from the electron system into the phonon system produces an increase in the electron-phonon scattering rate and therefore has the effect of reducing the current. This reduction is substantial (at Vg = 0.7 V), compared with the case in which the phonon system is considered at a temperature of 300 K or no self-heating. This is clearly shown in the current voltage characteristics depicted in Figure 7. However, the effect is even large if we neglected hot electron effects and assumed that the electron system is in local equilibrium with the potential, i.e., assuming the Joule law. The Figure 7 shows a very large reduction in current when Joule heating is considered. The corresponding temperature profile and heat source profile for this case are shown in Figure 18 and Figure 19. The heat source in this case (Figure 20) is just proportional to the gradient of the electrostatic potential. Note that the maximum temperature (over 600 K) is substantially larger than the previous case, indicating that the Joule heat grossly overestimates the power dissipation inside small nanowire transistors.

Figure 19 shows the temperature along the device corresponding to the case in which the electron system dissipates the Joule power into the lattice but where we fixed the lattice temperature of 300 K in computing the scattering of electrons. In this case as there is not a reduction in the electron energy current the lattice temperature increases substantially, the maximum peaking around 1000 K. This shows the danger of considering only the energy transfer from the electrons to the lattice, but not the impact of the heated lattice on the electron system for small nano-transistors.

3.5. Perspectives and Challenges

In this section, we describe certain open issues and limitations of the simulation of carrier transport in nano-transistors relevant to our work. Even though substantial progress in nanowire transistor simulation has been achieved, reliable and unified models do not exist yet. The majority of simulation methodology has a simplified description of the contacts. Usually, the contacts have the same cross section as the active part of the device. This is very convenient when using the NEGF formalism as it requires a region of uniformity in order to properly inject and absorb carriers through self-energies. A more realistic description of the contacts would require [59] a broadening of the contacts relative to the device active region. This broadening, of course, would depend on the realistic geometry of the device and needs a large simulation domain as a large part of the extension regions would need to be added. Calculations using different source/drain geometries have been carried out in a ballistic context [95]. In the self-consistent NEGF/Poisson formalism used to simulate devices, the Fermi levels are fixed in the boundary of the source and drain and the potential is allowed to float to achieve charge neutrality. This fixing of the Fermi level would be more accurate as the electrical neutrality in the contact should not be affected by the device active region, as the contacts will ultimately represent reservoirs in equilibrium. However, these reservoirs could be made of or connect 1D, 2D and 3D electron states. In addition, the broadening of the contacts will create new challenges from a computational point of view as scattering needs to be considered in those contacts in order to inject electrons in the narrow regions with the proper equilibrium distribution, the proper quantum states and energies. The transition region between the simulated contacts and the narrow device region will induce unwanted reflections if the transition is abrupt, therefore a smooth change in cross section, from wide to narrow channel is required in conjunction with decoherence and randomization due to electron-phonon scattering. A proper smoothness of the transition will require the use of the finite element method [96] to treat the change of geometry in an intrinsic way. However, this will stretch the use of the recursive algorithm [97], which is the pinnacle of success in NEGF computational versatility. Electron-electron scattering in the interfaces and contacts could play a role but it would be challenging to introduce into reliable, computationally manageable and even theoretical models [51,52]. In addition, electrons in the channel experience remote interactions from interfaces and contacts such as soft optical interface phonon scattering [53,54] and remote Coulomb scattering. The latter references indicate that these mechanisms could have substantial impact on device performance. The modelling and description of these mechanisms are not unique and represent a challenge due to the reduced and changing dimensionality of the nanostructure.

The present study has not considered the influence of atomistic effects on electro-thermal modelling, such as stray discrete impurities, trapped charge at the interfaces or discrete effects of the atomistic interfaces, although we have made extensive studies of these effects in purely electron transport modelling (e.g., [11,79]). However, the high surface to volume ratio for nanowires raises the issue of electronic interface states and their effect on electro-thermal properties which deserves future scrutiny. In addition, the role of localised phonon modes at the interfaces could lead to important thermal channels and local non-equilibrium phonon effects.

However, there is one aspect of simplification that needs caution: it has been usual in most studies to neglect the real part of the self-energies. Unfortunately, this is not justified in general. For example, the retarded self-energy has to satisfy causality requirements and these lead to Kramers-Kronig relations (Hilbert transform) that relate the real part of the self-energy to the imaginary part. As we have shown, elsewhere (see Appendix G and references herein), neglect of the real part may lead to serious violations of the spectral density sum rule and possible violations of charge conservation. Significant differences may be demonstrated in device modelling [77]. In general, the real part of the retarded self-energy leads to energy-dependent level shifts and may lead to bifurcation of the spectral density of states. The drawback is that the procedure for the full self-energy (within the self-consistent Born approximation) is very much more compute-intensive (further details in references Appendix G).

The use of a heat equation in this work and others as an alternative to more complex methodologies and moving towards more computational efficient methods is promising. The similarity in results with other methodologies demonstrated in ref [94] is encouraging in spite of the strong dissimilarities between methodologies. Their use of NEGF description for phonon transport and the full phonon dispersion law, instead of the simple heat equation, made the methodology extremely computationally demanding and therefore limiting its scope. We use a very simplified version of the phonon scattering mechanisms [85,98] that has been proven to produce similar values and trends for low field mobilities [76] and which has been widely used in the Monte Carlo/Greenwood-Kubo methodologies and adapted for nanowire [98]. The calculation of phonon modes in confined nanostructures is also a challenge, as they strongly depend on the boundary conditions used [99]. In our study, the acoustic deformation potential has been parameterised to partially take care of the decrease in mobility for narrow nanowire following [86,98]. Using a heat equation has the advantage of extending thermal simulations in larger regions surrounding the device and interconnects. In addition, interfaces and boundaries between materials of different thermal conductivities can be handled. This opens up the possibility of using more realistic thermal environments [100,101] and compact models as thermal resistance. Thermal conductivities which vary with the temperature could also be incorporated. Simple phonon balance equations can be used instead of the heat law if necessary, alternatively, different heat equations for different phonon types (for example separately for polar and acoustic) could also be used and this could partially address heat transport through interfaces and phonon-phonon interaction. Finally, thermal transport processes in low-dimensional and non-trivial spatial topology nanostructures (so-called “holey” materials) are becoming important [102]for applications in thermo-electrics, photonics and batteries . Reference [102] provides a detailed review of this challenging field and showed that thermal conductivity measurements at room temperature are considerably lower than the Casimir limit for structures smaller than a few hundreds of nanometers. The traditional continuum models of thermal diffusion and conduction no longer apply. The power of the spatially resolved NEGF methodology may prove important in that context.

4. Conclusions

In this paper, we revised and reviewed some of the quasi-classical and quantum methodologies used for the electrothermal simulation of nanoelectronics devices. Details and insights about the physics behind the formalism were also provided and these were complemented with the relevant literature and clarifying appendixes. The rest of this work focused on the application of a blended electrothermal methodology to a gate all around nanowire transistor.

Specifically, we have carried out electrothermal simulations of a narrow gate all-around nanowire transistor using a quantum transport formalism. The non-equilibrium Green function formalism has been used to describe the electron transport and the electron-phonon scattering is described within the self-consistent Born approximation

We have modelled the lattice heating by a Fourier law, and used as a heat source the true loss or gain of energy of the electrons (hot electrons) to the phonons.

These simulations have been compared with the use of the Joule law as a heat source. The drain on-current is reduced around 30% when the hot phonon source is used and 70%, compared to when the Joule law alone is used.

We have found that at low gate bias a slightly cooling of the source region by a few degrees occurs. The effect is large at an intermediate gate bias that depends on the electron-phonon scattering strength and the length of the source. This effect might be used for cooling the device if the source material and the dimension of the device are properly engineered.

At high gate bias, the injected electrons go straight to the drain region without the need of phonon absorption as the source drain barrier height is negligible. The dissipation of energy from the high energy electrons injected from the source induces an overall heating in the active region of the device. This process is more pronounced in the drain side as has been confirmed by the simulations. This effect becomes larger as the current increases. However, as the lattice heats, the average phonon energy or the equivalent temperature of the phonon system increases, and as a consequence, the electron-phonon scattering rate increases, the electron current is decreased from its value when the lattice is at room temperature. This results in less heating of the lattice. Our simulation produces a maximum device temperature of 400 K at large gate bias.

We have also investigated the case that the electrons are in local equilibrium with the lattice and used the Joule power as a source of the heat equation, this approximation is only valid for large devices. This produces a large heating of the lattice and consequently a large decrease in the drain current. The maximum temperature is approximate 600 K.

The methodology presented in this paper is computationally efficient and allows the accurate prediction of the effect of self-heating in small devices and produces very similar results to more sophisticated but computationally intensive methodologies. As indicated, some simple models can reproduce the trends of more sophisticated ones in a more computational efficient way; however, these models require verification or at least justification of the principle used. These models usually provide a simpler and concise explanation of the behaviour of the devices as they used few physical parameters. This work is a small step in this direction of the search for more simplified and phenomenologically-based models allowing the incorporation of more phenomena and larger parts of the device environment in the simulation domain. The solution of the 3D heat equation has a similar computational complexity to the 3D Poisson equations, which are solved quite efficiently by iterative methods. Therefore, quantum transport can be naturally connected with self-heating and thermal transport in a large domain. This provides more realistic boundary conditions, i.e., boundaries at room temperature should be far away from the active region of the device. For a device, designer/engineer accuracy and speed are equally important as substantial numbers of devices and conditions would need to be considered before an optimum decision is taken.

It is anticipated that our methodology may find wide application in the design of power efficient nanodevices.