Full-band quantum transport simulation in the presence of hole-phonon interactions using a mode-space k·p approach

Fabrication techniques at the nanometer scale offer potential opportunities to access single-dopant features in nanoscale transistors. Here, we report full-band quantum transport simulations with hole-phonon interactions through a device consisting of two gates-all-around in series and a p-type Si nanowire channel with a single dopant within each gated region. For this purpose, we have developed and implemented a mode-space-based full-band quantum transport simulator with phonon scattering using the six-band k · p method. Based on the non-equilibrium Green’s function formalism and self-consistent Born’s approximation, an expression for the hole-phonon interaction self-energy within the mode-space representation is introduced.


Introduction
The aggressive down-scaling of the complementary metalsilicon-oxide (CMOS) field-effect transistors (MOSFETs) has allowed major technological innovations, among them, novel materials [1], higher-κ gate dielectrics [2] and new structures such as gate-all-around (GAA) nanowire FETs (NW-FETs) [1,3]. The latter is considered a very promising architecture, and a key component, for the next generation of transistors, both lateral and stacked GAA FETs [1,3,4].
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
According to the International Roadmap for Devices and Systems (IRDS) [4], the implementation of NW-FETs is a high priority for the industry. Hence, understanding their working principles/challenges, e.g. short-channel effects, is an urgent matter.
CMOS technology is also emerging as a promising asset for the realization of future quantum technologies, which combined with current nanometer fabrication techniques, can be used to fabricate multi-gate silicon structures with few dopants. A sketch of these devices is presented in figure 1, consisting of two GAA in series and a p-type Si nanowire channel with a single dopant within each gated region. Devices with similar architecture have shown promising applications in planar silicon-based spin qubits [5] or multi-dopant spectroscopy [6]. In order to provide insight into the working principles of these quantum devices, in this work, we study both coherence and dissipative quantum Figure 1. Simplified schematic of a silicon NW-FET nanowire with two gates. Distance between gates and the thickness of oxide layer is 1 nm. The latter is made from SiO 2 . Gate lengths and silicon body thickness are L g1 = L g2 = 7 nm and T b = 3 nm, respectively. Note that in each gated region there is a dopant, whose energy levels are independently controlled by the applied gate voltages. transport phenomena through the device introduced in figure 1.
For accurate simulation of the electrical properties of NW-FETs, a full quantum treatment of the charge carriers is mandatory. A full quantum transport simulation involves the self-consistent solution of both Schrödinger and Poisson equations. Commonly, Schrödinger and its corresponding quantum kinetic transport equations are solved within the framework of the parabolic effective mass approximation (PEMA) of the band structure. However, the PEMA is most likely to fail to reproduce the valence band structure due to its high anisotropy. Quantum transport simulation of p-type nanodevices then requires a full-band description of the valence band structure.
Formalisms like tight-binding (TB) [7,8], density functional theory (DFT) [9,10], and k · p [11] are capable of computing the proper material band structure. While, the simulation of larger devices with atomistic quantum transport simulators becomes computationally costly, the k · p formalism, usually restricted to the Γ-point, offers a trade-off between accurate description of the band structure and much lower computational burden for device transport simulation. The most significant attribute of the k · p method is that it is based on six-bands for describing the valence band of Si and eightbands for both the conduction and valence bands of direct band-gap semiconductors [12], e.g. III-V materials.
The k · p method has been used for simulating the transport properties of double-gate [13] and nanowire MOSFETs [14][15][16]. Comparison with TB simulations has shown that the two approaches are in agreement [17,18]. However, only a few attempts at including phonon scattering within the k · p framework have been made [13,18]. These studies, though, are based on a real-space representation of the k · p Hamiltonian, which is computationally costly.
In this work, we have developed a full-band quantum transport simulator based on the six-band k · p method that efficiently tackles hole-phonon interactions. This was achieved by combining the models proposed in [15] and [16]. Whereas, in the latter, the k · p Hamiltonian is projected into a much smaller subspace constructed by sampling the Bloch modes significant to the transport, the former allows a modespace Hamiltonian via unitary transformation of the reducedorder Hamiltonian. The transport problem is then solved by means of the non-equilibrium Green's function (NEGF) technique and the hole-phonon interactions are included based on the self-consistent Born's approximation (SCBA) [19,20]. Using the resulting mode-space Hamiltonian, a simplified expression for the hole-phonon interaction selfenergy is introduced. This approach has been successfully implemented in the multi-scale transport simulation tool called Nano-Electronic Simulation Software (NESS) [21], which is currently under development at the University of Glasgow.
The paper is organized as follows. In section 2 and appendix A, the six-band k · p method, reduced-order and mode-space models are presented. Section 2 is also complemented by appendix B, where the elements of the k · p Hamiltonian discretized in K-space are listed. Section 3 briefly summarizes the main expressions of both NEGF and SCBA formalisms. Moreover, the phonon scattering self-energy within the mode-space representation is introduced. In section 4, we discuss our simulation results. Finally, conclusions are drawn in section 5.

Six-band k · p Hamiltonian
The three-band k · p Hamiltonian in { |X⟩, |Y⟩, |Z⟩ } basis can be generally expressed as, where I denotes the identity matrix and the A, B, . . . , G ′ matrices contain information about the crystallographic direction, as shown in appendix A. The Kane parameters, L, M, N, are given in terms of γ Luttinger parameters, where t 0 = −ℏ 2 /2m 0 and m 0 is the free electron mass. Then, by choosing the spin-dependent basis |λS⟩ with λ = X, Y, Z and S = ↑, ↓ for spin-up and spin-down spinors, the six-band Hamiltonian takes the form, .
The last term accounts for the spin-orbit interactions with, where ∆ is the spin-orbit splitting parameter.

Reduced-order and mode-space Hamiltonian
By applying an external potential V(r) the periodicity of nanostructures is broken. The hole wavefunctions must then be obtained from the envelop equation, where the total wavevector k has been replaced by its differential operator. The Hamiltonian H R in equation (5) can be discretized in real space and implemented on NEGF solvers [13,14]. Nevertheless, the real-space representation of H R is computationally costly, particularly for the simulation of larger devices. We have adopted an alternative approach for the discretization of H R by exploiting both the reduced-order [16] model and mode-space [15] representation. This is summarized as follows.
The Hamiltonian H R is first transformed to K-space by the unitary transformation proposed in [15]. After a proper discretization in the transport direction, with N x cross-sectional planes, the K-space Hamiltonian can be written as, The unperturbed Hamiltonian h 0 = H 0 K + W K + W † K , the coupling Hamiltonian W K and potential V K (j) of the j-th crosssectional plane are matrices of size N b N K . N b is the number of bands considered in the H k·p Hamiltonian and N K is the number of mesh points in the K-space. N K depends on the nanowire body thickness and the mesh space in the y-and z-directions, considering the x-direction as the transport direction. N K can be optimized by constructing a reduced basis [16], as small as possible, able to accurately generate the original band structure within the energy range of interest, E n ≲ E ≤ E v , from the top of the valence band edge E v up to a cut-off energy E n .
The reduced basis is constructed by solving the eigenvalue problem, coming from equation (6) when no external potential is applied, at different k x sampled points. Filtering the eigen-energies lying inside the energy window of interest, a matrix with their corresponding eigenmodes Ψ m can be constructed, being then QR factorized [16]. Using Q as a unitary transformation matrix, the reduced-order Hamiltonian discretized in K-space can then be obtained:H whereh 0 = Q † h 0 Q. Note that the aforementioned approach becomes equivalent to the one proposed in [15] when only equation (7) is solved at k x = 0. The K-space reduced-order Hamiltonian might still be large for efficient full-quantum transport simulation. Instead, a mode-space representation of equation (8) can be employed, leading to,H where,H The overall approach has been, therefore, transformed into finding the lowest N M ≪ N b N K eigenmodes relevant for transport simulation. For the latter, equations (9) and (10) have been implemented into a modespace NEGF solver, as follows.

Non-equilibrium Green's function approach with hole-phonon interactions
Within the mode-space framework presented in section 2, the retarded and lesser Green's functions for the active device region are written, in matrix notation, as, respectively. The HamiltonianH M and retarded/lesser self- is the self-energies for the semi-infinite leads, whereas, Σ R/< S accounts for all scattering mechanisms. In the case of both acoustic and optical phonons, Σ Holephonon interactions are tackled perturbatively within the SCBA [19,20].
In practice, the retarded/lesser Green's functions are computed by exploiting a recursive algorithm that allows us to access their most relevant elements for the calculation of the physical quantities, such as the hole density, where ∆x j = x j+1 − x j is the distance between adjacent planes. In equation (13), the following transformations, must be performed to recover the hole density from modespace to real-space representation. U K (r nm ) is the unitary matrix [15] to transform from K-space to real-space and r nm = (y m , z n ) was used for brevity.
In the case of the current, on the other hand, no transformation is needed. It can be straightforwardly computed from, ] .
The Green's function elements G < j,j ′ are N M × N M block matrices.

Hole-phonon interaction self-energy within the mode-space representation
Assuming bulk hole-phonon interactions, the retarded/lesser phonon self-energy in mode-space can be expressed as, where n (n ′ ) and p (p ′ ) are mode indices. For acoustic phonons, whereas, are the lesser and retarded Green's functions for the holeoptical phonon interactions, respectively. The term N ph is the Bose-Einstein distribution for the phonon energy ℏω ph . During numerical evaluation of equations (17) to (20), the principal value of G R and the real part of the self-energies are generally neglected. It was shown that this approximation introduced no significant error in the transport properties [20]. The charge and current criteria of convergence were set to 1%.

Evaluation of coefficients F n,p n ′ ,p ′
In the hole-phonon interaction self-energy introduced in equation (17), the coefficients F carry information about the type of hole-phonon coupling and form factors due to the basis transformations, expressed as, where the x j was omitted for brevity and d = x, y ,z. Being a vector of size N b × 1, ζ † n is the real-space wave function for the mode n at x j and . It is retrieved from the following basis transformation, where Λ n = U M (n) is the mode-space wave function for the mode n. Finally, the integral in equation (21) is carried out over the cross/confinement section/direction of the nanostructure. Furthermore, in this work, D d describes the hole-(acoustic/optical)-phonon interactions within the six-band k · p framework, being a 6 × 6 matrix. Herein, they are as given in [13], following the valence band deformation potential theory [22].

Simulations
The aforementioned mode-space full-band quantum transport approach has been implemented into the multi-scale transport tool called NESS [21]. We have efficiently simulated coherent and dissipative transport through a device consisting of two GAA in series and a p-type Si nanowire channel with a single dopant within each gated region. Both gates' work functions are 4.55 eV. The external source/drain bias is set to -0.1 V during this whole study. Figure 1 provides a detailed description of the device. The Si nanowire is 3 nm thick covered by 1 nm layer of SiO 2 . The gate lengths are L g1 = L g2 = 7 nm, being 1 nm distance from each other. The source/drain is highly doped with N A = 10 20 cm −3 and the region under both gates is considered intrinsic. The Coulomb potential of the dopants is described by adding a negative electron charge on the impurity site in the Poisson's equation. This approach correctly describes the long-range Coulombic tail, which is known to have the largest impact on the transport [23]. Electron-electron interaction is then only treated within the Hartree approximation by solving the NEGF transport equations self-consistently with Poisson's equation.
The source and drain Fermi energies are fixed at E FS = 0 and E FD = 0.1 eV, respectively. Then, the energy window considered for quantum transport simulations is chosen as E FS − 1.0 eV < E < E FD + 0.2 eV in order to ensure a good convergence. Luttinger parameters used for the six-band k · p Hamiltonian were taken from [15], where they were calibrated to sp 3 s * TB electronic calculations performed with the TCAD tool OMEN [19].
Finally, it is worth mentioning that the performance of the mode-space k · p approach presented in previous sections is limited by the number of modes N M , size of the nanowire, and by how dense the energy grid is for the NEGF recursive algorithm. For instance, thicker nanowires would require us to increase N M and to consider a denser energy grid. However, in comparison to real-space k · p, the mode-space approach would still offer higher performance.

Band structure calculations
The K-space Hamiltonian has been discretized using a uniform mesh in all directions with mesh step a 0 = 0.2 nm, leading toÑ K = 225. By applying the QR factorization approach, as explained in section 2, the size of the original problem is reduced, for instance, by approximately 60% in the case of a window energy of E v − E n = 0.6 eV, finding that only N K = 92 modes are needed to fully reproduce the band structure with very good accuracy. The latter can be observed in figure 2, where the comparison of the full band structures of a 3 nm body thickness Si nanowire in the ⟨111⟩ crystallographic direction was computed by means of the K-space and mode-space Hamiltonians. Note that for energies below E ≈ −1.7 eV the mode-space approach becomes less accurate, particularly around k x = 0. Regarding the transport problem in section 3, we have considered N M = 12 (≪ N K ) modes after directly diagonalizing equation (8). Although including additional modes could lead to more accurate results, no substantial improvement of the accuracy was obtained. In particular, the latter number of modes was a good trade-off between accuracy and computational cost.   to the effect of the external source-to-drain bias. In equilibrium conditions, i.e. when V DS = 0 V (not shown here), DELs from both dopants perfectly align to each other in terms of energy.

Ballistic limit
In figure 4(b), where V GS1 = −0.6 V and V GS2 = −1.0 V, the highest DELs of each dopant, are separated by approximately −0.18 eV despite the absolute value of V GS2 being increased by 0.4 V. This counter-intuitive phenomenon can be explained by the different screening mechanisms of each dopant. The dopant on the left is closer to the source, getting screened faster than the dopant on the right, which is mostly populated by the electrons injected from the drain. First, the decrease of the gate voltage increases the initial DEL, which is progressively populated by source/drain holes. This screening mechanism produces an opposite effect of level lowering, and the DEL slowly saturates towards the average potential in the source/drain, as observed in figure 5. When the gate voltage is strong enough, the DEL can go above the average potential in the source/drain. This is the case of V GS1 = −0.9 V and V GS2 = −1.0 V, in figure 5, where one can see how the DEL of the first dopant abruptly goes above the average potential in the source. Meanwhile, the DEL of the second dopant slightly changes over the V GS1 range for all values of V GS2 considered in this work.

Impact of hole-phonon interactions
This section reports on the impact of hole-phonon interactions on the ballistic simulation results presented above. The quantum transport simulations, in the presence of phonon scattering, were carried out by means of the NEGF and SCBA formalisms within a mode-space representation of the sixband k · k Hamiltonian introduced in section 3. Simulations are presented below.   Figure 6 shows the LDOS in the presence of hole-phonon interactions for the same gate voltage configuration as its ballistic counterpart in figure 4(b) at V GS1 = −0.6 V and V GS2 = −1.0 V. Due to phonon scattering, the LDOS is globally broadened, as observed around the DELs. Note that in the ballistic regime ( figure 4(b)), three sharp DELs exist. However, by switching the phonon scattering on, the highest two DELs are merged in one broader energy level. One can also observe that there is an overlapping of energy levels of the second dopant.
A stronger effect of phonon scattering can be seen in the current spectrum, as reported in figure 7, which compares spectra corresponding to the ballistic regime ( figure 7(a)) and the one in the presence of hole-phonon interactions (figure 7(b)) for V GS1 = −0.6 V and V GS2 = −1.0 V. In the case of the ballistic, we can observe three resonant current branches within the window energy of −80 and −150 meV. Each branch corresponds to the energy levels of the first dopant, as shown in the LDOS plotted on figure 4(b). The fourth and broadest branch in figure 7(a) below the potential barrier carries the current contribution of thermionic holes injected from the source and drain regions. In the presence of phonon scattering, holes injected from the source can relax into the resonant energy levels by emitting a single or multiple phonons. The latter is observed in figure 7(b), where the resonant sites help holes to relax to higher energies, starting from the first dopant site up to the drain region. Interestingly, holes that are backscattered can also relax to higher energies in the direction of the source region. Figure 8 shows the impact of hole-phonon interactions on the I D − V GS1 characteristics for two different configurations: V GS2 = −0.6 V (blue lines) and V GS2 = −1.0 V (red lines). When compared to the corresponding ballistic case, the current is reduced over the entire range of V GS1 shown here. Their losses are plotted in the inset figure with the same color code and V GS1 range. The percentage loss is defined as |I ball − I scatt |/I ball × 100%. Similar to current p-type MOSFETs, the lower the V GS1 , the stronger the impact of phonon scattering causing a reduction in the current and therefore enhancing the losses. For the transistor considered in this work, this trend only happens for high voltages, i.e V GS1 ⩾ −0.7 V, within the sub-threshold part of the I D − V GS1 characteristics. The loss  has a slight decrement at V GS1 = −0.9 V and V GS1 = −0.7 V in the case of V GS2 = −0.6 V and V GS2 = −1.0 V, respectively. The latter could be explained by a substantial contribution of the tunneling current flowing through the resonant site of the first dopant when its energy level gets closer to the Fermi level. For lower V GS1 , optical-phonon scattering is expected to become dominant, contributing to an increase in the loss, which reaches up to 60% in the case of V GS1 = −1.0 V and V GS2 = −1.0 V. Figure 9 shows the current spectrum for the previous configuration where the strong effect of the phonon scattering can be seen. In this case, the DEL is below the source band edge, creating a resonant site where holes injected from the source can tunnel to the drain. A second current branch around the energy level of the first dopant is then generated. Due to the strong electric field within the gate regions, it quickly fades due to the large number of holes relaxing to higher energies until they reach the top of the drain band edge.

Conclusion
A mode-space-based full-band quantum transport simulator with phonon scattering using the six-band k · p method has been developed and implemented in Nano-Electronic Simulation Software (NESS) from the University of Glasgow. Based on the non-equilibrium Green's function formalism and selfconsistent Born's approximation, an expression for the holephonon interaction self-energy within the mode-space representation was introduced. We efficiently simulated a transistorlike device, which may be relevant for future quantum technologies, formed by two GAA in series and a p-type Si nanowire channel with a single dopant within each gated region. The dopant's Coulomb potentials were mimicked by adding a negative electron charge on the impurity site in the 3D Poisson's equation.
Our findings suggest that hole-phonon interactions have a significant detrimental impact on the transistor performance when compared to ballistic simulations. Here, we have reported that the current could be reduced by up to 60%, depending on the gate voltage configuration. U KL (r ⊥ ) = U(r ⊥ , k L ) = 2 L y L z sin(k p y) sin(k q z).
(32) k L = π(p/L y , q/L z ) is evaluated in the k-space grid point L = (p, q). L y and L z are the cross-section lengths in y-and z-directions, with real-space coordinates r ⊥ = (y, z), respectively. Carrying out the integral over the cross-sectional plane in equation (31), one straightforwardly obtains, After a proper discretization in the transport x-direction, with N x cross-sectional planes, the K-space Hamiltonian can be written as, which briefly can be re-expressed as follows, for the j-th cross-sectional plane. In the latter, the spin-orbit coupling interactions have been included in its first term on the right-hand side. Ultimately, in equation (36), and the potential V ′ K in K-space can be numerically evaluated, for instance, as, where r i = (y m , z n ) is the i-th grid point in real space.