Reduced k-points based computationally efﬁcient band structure approach for UTB device simulation

—The accurate calculation of channel electrostatics parameters, such as charge density and potential, in ultrathin body (UTB) devices requires self-consistent solution of the Poisson’s equation and the full band structure, which is channel material and thickness dependent. For cubic crystals like silicon, the semi-empirical sp 3 d 5 s ∗ tight-binding (TB) model is preferred in device simulations, over the density functional theory, to obtain the full band structure because of being computationally less intensive and equally accurate. However, the computational time of the TB model scales non-linearly with the channel thickness and becomes cumbersome for silicon, beyond 5 nm, primarily because of the increasing size of the TB hamiltonian that needs to be solved over the entire k-space, in the irreducible Brillouin zone. In this work, we precisely identify those k-points corresponding to the energies close to the band minima, where the Fermi-Dirac probability signiﬁcantly affects electrostatics parameters. This enables us to demonstrate a computationally efﬁcient approach based on solving the hamiltonian only on those reduced number of k-points. The rigorous benchmarking of the channel electrostatics parameters obtained from this approach is performed with results from accurate full band structure simulations showing excellent agreement over a wide range of channel thicknesses, oxide thicknesses, device temperatures and different channel orientations. By showing that the approach presented in this work is computationally efﬁcient, besides being accurate, regardless of the number of atomic layers, we demonstrate its applicability for simulating UTB devices.


I. INTRODUCTION
U LTRA-thin body (UTB) devices offer better immunity from short channel effects and their planar geometry makes them a good candidate to succeed bulk MOSFETs at nanometer technology nodes [1]. Additionally, the fully depleted intrinsic channel provides immunity to UTB devices from random dopant fluctuations [2]. Hence, these devices show excellent prospects for low power digital [3] as well as analog applications [4].
The channel of the UTB device consists of a stack of atomic layers in the thickness direction. The number of atomic layers affects the energy dispersion relation of the channel material and therefore the electronic properties of the UTB device depends considerably on the channel thickness [5] [6]. Apart from this, the calculation of the channel electrostatics parameters of the UTB device requires consideration of various other effects such as: (i) quantum confinement effect [7] channel orientations and surface passivations, (iii) projections of dispersion valley from bulk to surface brillioun zone [9], and (iv) effective mass variation as a function of channel thickness and orientation [10], that can be inherently incorporated using atomistic or full band structure simulations. The accurate full band structure simulations are done using density functional theory (DFT), however, the computational burden it costs limits its usage for device simulations [11]. The semi-empirical sp 3 d 5 s * tight-binding (TB) model is an alternative, which generates fairly accurate full band structure for materials with cubic crystal symmetry, with reducing computational time compared to DFT. Given these advantages, the TB model is widely used to obtain channel electrostatics parameters for UTB and nanowire devices [12] [13].
To obtain the channel charge density and potential, the full band structure obtained from the TB model needs to be solved self-consistently with the Poisson's equation, using an iterative scheme [13] [14]. This scheme becomes time consuming and computationally cumbersome with increase in number of atomic layers along the channel thickness, thus limiting its applicability for larger channel thicknesses. A simplification in the scheme was proposed by multiplying the Boltzmann factor exp(qφ/kT ) in the charge density expression, in order to consider the effect of potential (φ) on the band structure, thus, avoiding solving the TB hamiltonian in the self-consistent loop [14]. This approach works well below the threshold voltage, however, it starts predicting semi-classical behaviour at higher gate voltages and higher channel thicknesses. Besides this approach, by utilizing parallelism, Luisier et. al [15] have shown a multi-scale TB simulation for channels thicknesses upto 12 nm, which requires high computational resources.
The main reason for high computational time of the TB model, with large number of atomic layers, is the increasing dimension of the TB hamiltonian that needs to be diagonalized for the entire k-space in the irreducible Brillouin zone (BZ). Furthermore, this problem will become more serious for some materials of recent interest that have lower lattice constants, such as Diamond (3.567Å [16]), cubic-BN (3.615Å [16]), BAs (4.777Å [17]), SiC (4.359Å [18]), AlN (4.394Å [19]), and GaN (4.512Å [19]), for which even small channel thicknesses will have a large number of atomic layers. Therefore, given the well documented advantage of the sp 3 d 5 s * TB model, we try to extend its usability for larger number of atomic layers, while being computationally efficient compared to the full band structure approach, without sacrificing its accuracy and versatility.
In this work, we propose an approach based on the systematic reduction of k-points that results in the reduction of the number of times the TB Hamiltonian needs to be diagonalized. This is made possible by precisely identifying the k-points corresponding to the range of energies where the Fermi-Dirac probability is significant, which in-turn, impacts the channel charge density. This significantly reduces the computational time in the proposed approach.
As a way of demonstrating the validity of the proposed reduced k-points based band structure approach, we apply it to a double gated (DG) device with a UTB silicon channel for various channel thicknesses ranging from 2 nm to 10 nm over (100) and (111) channel orientations, oxide thicknesses ranging from 1 nm to 5 nm and over different device temperatures from 15 K to 300 K. The channel electrostatics parameters calculated are rigorously benchmarked with the accurate full band structure simulation results and are found to be in excellent agreement, while also showing 90% reduction in the computational time.
Additionally, the proposed band structure based simulation results are also compared with the results obtained from other approaches in literature, such as the density gradient model (DGM) and the effective mass approximation based Schrodinger-Poisson's model (EMA), clearly showing good qualitative agreement over an entire range of channel thicknesses and gate voltages.
The rest of the paper is organized as follows: In section II, the details of getting the self-consistent solution of the full band structure and the Poisson's equation is discussed and the reduced k-points based band structure approach is proposed. In section III, the channel electrostatics parameters obtained using the proposed approach benchmarked with the full band structure is provided first for the (100) as well as (111) channel orientation. The comparison of the proposed approach with other approaches in literature is shown in section IV. Finally, the conclusion is presented in section V.

II. APPROACH
In this section, firstly a procedure is discussed to obtain the self-consistent solution of the full band structure and the Poisson's equation for a UTB DG MOS device with an intrinsic Si (100) surface, which is shown in Fig. 1.

A. Full band structure approach
The sp 3 d 5 s * tight-binding (TB) model considering 10orbitals has been used to obtain the full band structure of the thin Si channel [14] [20]. In this model, the two centre overlap integrals considering four nearest neighbours per atom are taken into account, while the spin-orbital interactions have been neglected [21] [22]. Infinite crystal periodicity and thus Bloch's theorem is assumed to hold good along the channel length and width directions. However, given the ultra thin nature of the channel, the crystal cannot be assumed to be periodic along the thickness direction, z. The TB hamiltonian has been constructed using the well known procedure outlined in [14][20] [23]. The TB fitting parameters for silicon given by Boykin et al. [24] are used at room temperature, while those from Jancu et al. are used at lower temperatures [25]. The top  and bottom surfaces are assumed to be hydrogen passivated [26].
Given the direction of confinement being along the z direction, the wave-vector k z in the bulk Brillioun zone (BZ) is considered to be zero, thus, reducing the 3D BZ to a 2D one. The 2D BZ is discretized along the x and y directions using the step size ∆k x and ∆k y , which are taken as 0.05 × 2π a , where a is a lattice constant of Si. The TB hamiltonian is assembled on each k-point in the irreducible BZ and eigenvalues of the hamiltonian plotted along k gives the band structure of the channel material.
The algorithm of the full band structure approach is shown in Fig. 2, where the full band structure is solved selfconsistently with the Poisson's equation, which is given as where, φ is the electrostatic potential, is the channel's dielectric constant and ρ(z) is the electron density along the direction of the channel thickness (T si ). It may be noted that, the channel thickness is a function of number of atomic layers (N ), which, for the case of (100) surface, is given as (N − 1) × a/4. The boundary conditions at the top and bottom oxide-channel interface are the same as used by Majumdar et al. [14].
The electron density ρ(z) used in Eqn. 1 obtained from [27] is given below (2) where, n is the number of energy levels in the conduction band, the factor F is the Fermi-Dirac probability, g k is the degeneracy factor at each k-point, E F is the gate Fermi level and ψ k n (z) is the normalized wavefunction along z. Due to two possible spin states at each k-points, the factor of 2 is multiplied in the above equation. The double summation takes into account the contribution of various energy levels in the conduction band (n) at all k-points in the BZ for electron density calculation, thus, termed as the full band structure approach. The potential energy (-qφ) obtained by solving the Poisson's equation is added with the onsite energies at the diagonal of the TB Hamiltonian to consider the effect of potential on the bandstructure. The resulting TB Hamiltonian is solved self-consistently with the Poisson's equation to obtain the channel charge density and potential at a given gate voltage, as shown in Fig. 2.

B. Proposed reduced k-points based bandstructure approach
In this section, a computationally efficient approach is proposed where k-points of interest are identified and utilized in the band structure simulations to accurately represent the electrostatics of the UTB MOS device. By restricting the simulation to relevant k-points, instead of the entire k-space used in the full band structure approach, higher computational efficiency is achieved, which is demonstrated over a wide range of channel thicknesses, oxide thicknesses, device temperatures and surface orientations.
The conduction band of the full band structure obtained using the algorithm given in Fig. 2, for (100) Si channel with T si = 2 nm at V g = 0 V, is shown in Fig. 3a. From this figure, we notice two valleys (Γ and X), out of which the band minima can be found at the Γ point, for each energy level n, across all k-points. The energies in the figure are specified with respect to the valence band maxima, which is also located at the Γ point for the case of Si. This indicates that, for T si = 2 nm, Si behaves as a direct bandgap semiconductor. For an intrinsic channel, at zero gate voltage (V g ), the Fermi level lies in the middle of the band gap, while, with increase in V g , the conduction band moves towards the Fermi level, which beyond the threshold voltage shows insignificant movement, as can be seen in Fig. 3b. It may also be seen that the band minima remains at Γ point irrespective of V g . Similarly, for T si = 10 nm, where the confinement effect substantially recedes, both the X and Γ valleys have significant contributions at V g = 0 V, as seen in Fig. 4a. However, it is interesting to notice from Fig. 4b, that the separation in energies between the Γ point and X point goes on increasing with V g , thus, direct bandgap transition can be prominently seen for T si = 10 nm at higher gate voltages. The variation in Fermi-Dirac (FD) probability across all kpoints is shown for n = 1 and 2, for T si = 2 nm and 10 nm, in Fig. 5, when V g = 0 V. From Fig. 5a, it can be seen that the maximum value of FD probability is seen at the Γ point for T si = 2 nm, while apart from the Γ point for T si = 10 nm, FD probability is also significant near the X point, as shown in Fig. 5b. Therefore, from Figs. 3, 4 and 5, it may be seen that in order to accurately simulate the electrostatics of the UTB MOS device, k points near the band minima play an extremely important role. The identification of those k-points at V g = 0 V, T = 300 K, and by solving the band structure selfconsistently, at only those k-points of interest, as opposed to the entire k space, will enable a bandstructure approach to be computationally efficient, while also being accurate. Thus, it becomes important to identify those k-points whose energy is located close to the band minima and whose FD probabilities are likely to be significant.
Given that the maximum value of the FD probability (F max ) is obtained for the energy at the band minima (E min ), the F max can be given by the following expression, It may be noted that the FD probability reduces exponentially across all k-points for those energy levels which are far away from the Fermi level. By specifying the range of energies above the band minima, based on the magnitude of the FD probability, we propose an approach to identify the k-points of importance. Suppose the FD probability of any energy (E k n ) is η multiplied by F max , then the energy range (∆E = E k n − E min ) for which the FD probabilities are between η×F max and F max , can be expressed as Equation 4 can be simplified by considering the values of F max obtained from the full band structure approach as shown in Fig. 5, which is of the order of 10 −10 , at V g = 0 V. Therefore, by substituting this value in Eqn. 4, we can write ∆E ≈ kT ln 1 η . It may be noted here that η is a semi-empirical parameter, whose choice is likely to impact the accuracy as well as computational time of the proposed approach. A larger value of η will result in reduction of the number of significant k-points, which will in-turn reduce the computational time, but may also lead to substantial error in comparison with the full band structure approach. Therefore, an optimal value of η is required, such that it remains valid over a wide range of device conditions. Based on this choice of η, in Fig. 6, we propose a modified algorithm, which we term as the reduced band structure approach (RBA), as this is based on using reduced number of k-points. One of the key features of the proposed approach is that the full band structure needs to be solved only once to identify the relevant k-points based on the value of η, which can then be utilized in the self-consistent loop to speed up the iterative scheme.
In the next section, through detailed benchmarking of the RBA with the accurate full band structure approach (FBA), the optimal value of η is determined for different device conditions.

APPROACH
The benchmarking of the channel electrostatics parameters for different channel and oxide thicknesses is carried out first for Si-channel with (100) surface orientation. The benchmarking and validity of the approach for lower device temperatures is performed on the (100) surface and then over different channel orientations i.e. (111) surface.

A. Channel with (100) surface orientation
The comparison between the electron density profiles obtained from the full band structure approach (FBA) and the reduced band structure approach (RBA), for different channel thicknesses and gate voltages, are shown in Fig. 7. In case of RBA, three values of η = 0.1, 0.05 and 0.01 corresponding to ∆E of 2.19kT , 3kT and 4.6kT , respectively, are considered, to determine the optimal value of η, which will enable the best agreement with FBA, for V g = 0 V and 2 V, respectively. From these figures, it can be clearly seen that for η = 0.01 excellent agreement can be seen between the electron density profiles obtained using the FBA and RBA for both V g = 0 V and 2 V, and T si = 2 nm and 10 nm at T = 300 K (maximum percentage error less than 0.5%).
To show the reduction of computational time using the RBA, the time taken by both the approaches to simulate electron density, for various T si , at V g = 0 V, on a machine with 16  GB RAM, Intel i7 CPU, is given in Table. I. It can be seen that beyond 37 atomic layers (ALs), which corresponds to 5 nm T si for Si (100) surface, the FBA starts to become computationally cumbersome. On the other hand, precisely chosen but fewer number of k-points, enable the RBA to have substantially lower simulation time compared to FBA. The huge reduction in computational time makes the proposed RBA viable over large number of ALs, while ensuring that it remains as accurate as the FBA.
To establish the validity of RBA over different oxide thicknesses (T ox = 1 nm and 5 nm), the integrated channel electron density is benchmarked with results from FBA in Fig. 8a and 8b, for T si = 2 nm and 10 nm, respectively. It may be noted that for η = 0.01, RBA shows excellent agreement with FBA over the entire range of channel and oxide thicknesses at room temperature.  Now, we establish the validity of the proposed RBA for simulating the electrostatics of UTB DG MOS device at low temperatures. At low temperatures, the TB hamiltonian is constructed using the Si TB parameters given by Jancu et al. [25]. The variation of bandstructure at 50 K and 300 K for T si = 2 nm at V g = 0 V and 2 V, is shown in Fig. 9. It may be noted that apart from the change of band gap, which can be seen at Γ point, there is an insignificant variation in the curvature of the band.
The full band structure approach is then implemented using the algorithm shown in Fig. 2 and channel electrostatics parameters are obtained. It may be noted from Eqn 4 that η ≈ exp − ∆E kT . Therefore, for a particular value of ∆E, η must sharply decrease with decreasing temperature. In Fig. 10, we compare the integrated electron density obtained from FBA and RBA (with η appropriately determined) showing excellent agreement, over the entire range of gate voltages and channel thicknesses, at T = 15 K and 50 K, respectively. The η required to obtain good match with the FBA for device temperature is determined and tabulated in Table II. In this section, we show the validity of the RBA for simulating the electrostatics of UTB DG MOS device, with Si (111) surface, whose structure and atomic arrangement is shown in Fig. 11. Before implementing the FBA, the coordinates of the nearest neighbours for this atomic arrangement are determined and the TB hamiltonian is constructed using the Si TB parameters [24]. The surface is passivated with hydrogen atoms, similar to the case of (100) surface. The total number of k-points in the irreducible BZ for (111) surface, is 201, as compared to 121 for (100) surface. As the average inter-plane spacing is more in the case of (111), the channel thickness of 2 nm and 10 nm, corresponds to 14 and 64 ALs, which were 17 and 75, respectively, in case of Si (100) surface. The FBA is then implemented using the algorithm shown in Fig. 2 and channel electrostatics parameters are obtained. The band structure for 2 nm and 10 nm T si are shown in Fig.  12, where the band minima is seen near X point for all T si . Then, the information about significant k-points is extracted, using Eqn. 4, while following the algorithm shown in Fig. 6. The electron density profiles obtained using RBA for different values of η, for channel with T si = 2 nm and 10 nm, are shown in Fig. 13, where it can be seen that electron density profiles obtained from RBA shows excellent match with FBA for η = 0.01 over a wide range of gate voltages and channel thicknesses. Finally, in order to show the applicability of RBA (with η = 0.01) for different oxide thicknesses, we compare the integrated electron density obtained from RBA with the FBA, showing very good agreement, over the entire range of device parameters considered, at T = 300 K. In the case of (111) surface orientation as well, the computational time using RBA is reduced by 90%, compared to FBA for the simulation of electron density profile of T si = 10 nm, at V g = 0 V. Thus, through rigorous benchmarking, we established that the proposed reduced band structure approach is as accurate as the full band structure approach, while being far more computationally efficient over a wide range of atomic layers and other device conditions. In order to demonstrate the greater effectiveness of the RBA for accurately simulating the effects of quantum confinement on the electrostatics of UTB devices, in the next section, we compare the results from the RBA with other methodologies in literature, for a typical case of Si (100) surface.

IV. COMPARISON WITH OTHER METHODOLOGIES
The channel electron density and integrated electron density obtained from the RBA is compared with the density gradient model (DGM) used to simulate the quantum confinement effect in Sentaurus TCAD [28] and the effective mass Schrodinger equation approach (EMA) used in the SCHRED tool [29]. The DGM is employed to include the quantum confinement effects by adding potential like function Λ n in the classical charge density [30] [31], where as, EMA is based on the parabolic E −k relationship, which is assumed to be an accurate description near the band minima. This assumption of parabolic E − k relationship is violated in a devices having strong confinement and therefore, not valid for T si < 3 nm [10] [32].  The channel electron density profiles for T si = 2 nm, where strong structural confinement is seen, and for T si = 10 nm, where weaker confinement is seen, is shown in Fig. 15a and 15b, respectively, for Si (100) surface at V g = 2 V. It may be noted, that the RBA agrees much better with DGM compared to EMA for T si = 2 nm, while showing much better agreement with EMA for T si = 10 nm. In Fig. 16, the percentage error in integrated electron density from RBA, for T si = 2 nm and 10 nm is compared with results from EMA and DGM for various gate voltages. The considerable error between various approaches (EMA and DGM) and the results from the RBA (already benchmarked with results from the full band structure approach), particularly for V g < V th , shown in Fig.  16, highlights the importance of the bandstructure approach for accurately simulating the electrostatics of UTB devices, particularly in the strong volume inversion region.

V. CONCLUSION
In this work, the reduced k-points based band structure approach was presented for simulating the electrostatics of UTB devices. Based on the Fermi-Dirac probabilities, we modeled a semi-empirical parameter, which was used to obtain the condition to reduce the k-points, in the first irreducible Brillioun Zone. Using this, the band structure was only calculated on those significant k-points, which resulted in substantial reduction in the computational time of this approach over the full band structure approach, while being equally accurate. We show that the value of this semi-empirical parameter remains unchanged over different device conditions at a particular temperature, while changing only with device temperatures. The versatility of the proposed approach was established through comprehensively benchmarking it against the full band structure approach over a wide range of channel thicknesses, oxide thicknesses, device temperatures and surface orientations, showing excellent agreement throughout. With this, we show that the proposed reduced k-points based band structure approach is extremely beneficial for simulating the electrostatics of UTB devices independent of the number of atomic layers.