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

The aggressive down-scaling of the complementary metal-silicon-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]. 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 transport phenomena through the device introduced in figure 1.

Zoom In Zoom Out Reset image size

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 nano-devices 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 eight-bands 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–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 mode-space Hamiltonian via unitary transformation of the reduced-order 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 self-energy 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.

The three-band k · p Hamiltonian in $\big\{| X \rangle, | Y \rangle, | Z \rangle \big\}$ basis can be generally expressed as,

$$\begin{align} H’_{k \cdot p} = & Mk^2 \mathbf{I} + (L-M) \left( A k^2_x + B k^2_y + C k^2_z + D k_x k_y \right.\nonumber \\ &\left. + F k_y k_z + G k_x k_z \right) -N \left(A’ k^2_x + B’ k^2_y + C’ k^2_z + D’ k_x k_y \right.\nonumber \\ & \left.+ F’ k_y k_z + G’ k_x k_z \right), \end{align} \tag{ 1 }$$

where $\mathbf{I}$ denotes the identity matrix and the $A, B, \dots, 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,

$$\begin{equation} L = t_0 \big( \gamma_1 + 4 \gamma_2 \big), \,\, M = t_0 \big( \gamma_1 - 2 \gamma_2 \big), \,\, N = 6t_0\gamma_3, \end{equation} \tag{ 2 }$$

where $t_0 = -\hbar^2/2m_0$ and m₀ is the free electron mass. Then, by choosing the spin-dependent basis $| \lambda S \rangle$ with λ = X, Y, Z and $S = \, \uparrow,\downarrow$ for spin-up and spin-down spinors, the six-band Hamiltonian takes the form,

$$\begin{equation} H_{k \cdot p} = \left(\begin{array}{cc} H’_{k \cdot p} & 0\\ 0 & H’_{k \cdot p} \end{array}\right) + \left(\begin{array}{cc} \mathbf{M} & \mathbf{N}\\ -\mathbf{N}^* & \mathbf{M}^* \end{array}\right). \end{equation} \tag{ 3 }$$

The last term accounts for the spin–orbit interactions with,

$$\begin{equation} \mathbf{M} = \frac{\Delta}{3} \left(\begin{array}{ccc} 0 & -i & 0\\ i & 0 & 0 \\ 0 & 0 & 0\end{array}\right), \quad \mathbf{N} = \frac{\Delta}{3}\left(\begin{array}{ccc} 0 & 0 & 1\\ 0 & 0 & -i \\ -1 & i & 0\end{array}\right), \end{equation} \tag{ 4 }$$

where Δ is the spin–orbit splitting parameter.

Reduced-order and mode-space Hamiltonian

By applying an external potential $V(\mathbf{r})$ the periodicity of nanostructures is broken. The hole wavefunctions must then be obtained from the envelop equation,

$$\begin{align} H_R(-i\nabla, \mathbf{r}) \Psi(\mathbf{r}) & = E\Psi(\mathbf{r}) \nonumber \\ & = \left[ H_{k \cdot p} (-i\nabla) + V(\mathbf{r}) \mathbf{I} \right]\Psi(\mathbf{r}) , \end{align} \tag{ 5 }$$

where the total wavevector $\mathbf{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,

$$\begin{equation} H_K(j) = h_0 + V_K(j). \end{equation} \tag{ 6 }$$

The unperturbed Hamiltonian $h_0 = H_{0~{\rm K}} + W_{K} + W^{\dagger}_{K}$, the coupling Hamiltonian W_K and potential V_K (j) of the j-th cross-sectional plane are matrices of size $N_b N_K$. N_b is the number of bands considered in the $H_{k \cdot 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} \lesssim E \leq E_{\mathrm{v}}$, from the top of the valence band edge $E_{\mathrm{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,

$$\begin{equation} \left( H_{0~{\rm K}} + W_{K} \, \mathrm{e}^{i k_x \Delta_x} + W^{\dagger}_{K} \, \mathrm{e}^{-i k_x \Delta_x} \right) | \Psi_m \rangle = E | \Psi_m \rangle, \end{equation} \tag{ 7 }$$

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:

$$\begin{align} \tilde{H}_K(j) & = \tilde{h}_0 + \mathbf{Q}^{\dagger} V_K(j) \mathbf{Q} \nonumber \\ & = \tilde{h}_0 + \tilde{V}_K(j), \end{align} \tag{ 8 }$$

where $\tilde{h}_0 = \mathbf{Q}^{\dagger} h_0 \mathbf{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,

$$\begin{eqnarray} &&\tilde{H}_{M}(j) = \tilde{H}_{0~{\rm M}}(j) + \tilde{W}_M(j) + \tilde{W}^{\dagger}_M(j) ,\end{eqnarray} \tag{ 9 }$$

where,

$$\begin{align} \tilde{H}_{0~{\rm M}}(j) & = U^{\dagger}_M(j) \left( \tilde{H}_{0~{\rm K}} + \tilde{V}_K(j) \right) U_M(j), \nonumber \\ \tilde{W}_M (j) & = U^{\dagger}_M(j) \tilde{W}_K U_M(j+1) \end{align} \tag{ 10 }$$

are $N_M \times N_M$ matrices and the unitary transformation matrix U_M (j) consists of N_M mode wave functions $| \Psi_m(j) \rangle$ of size $N_bN_K$ obtained from the reduced eigenvalue problem $\tilde{H}_K(j) | \Psi_m(j)\rangle = E |\Psi_m(j) \rangle.$ The overall approach has been, therefore, transformed into finding the lowest $N_M \ll N_bN_K$ eigenmodes relevant for transport simulation. For the latter, equations (9) and (10) have been implemented into a mode-space NEGF solver, as follows.

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,

$$\begin{eqnarray} G^R(E) & = \left[ E\mathbf{I} - \tilde{H}_M - \Sigma^R_C(E) - \Sigma^R_S(E) \right]^{-1}, \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} G^ \lt (E) & = G^R(E) \left[ \Sigma^ \lt _C(E) + \Sigma^ \lt _S(E) \right] G^{\dagger R}(E), \end{eqnarray} \tag{ 12 }$$

respectively. The Hamiltonian $\tilde{H}_M$ and retarded/lesser self-energies $\Sigma^{R/ \lt }$ are $N_xN_M \times N_xN_M$ matrices. $\Sigma^{R/ \lt }_{C}$ is the self-energies for the semi-infinite leads, whereas, $\Sigma^{R/ \lt }_{S}$ accounts for all scattering mechanisms. In the case of both acoustic and optical phonons, $\Sigma^{R/ \lt }_{S} = \Sigma^{R/ \lt }_{ac} + \Sigma^{R/ \lt }_{opt}$. Hole-phonon 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,

$$\begin{eqnarray} p(x_j,y_m,z_n) & = \frac{-i}{\Delta x_j} \int \frac{\mathrm{d} E}{2\pi} \mathrm{{Tr}}\left[ p_R(x_j,y_m,z_n; E)\right], \end{eqnarray} \tag{ 13 }$$

where $\Delta x_j = x_{j+1} - x_j$ is the distance between adjacent planes. In equation (13), the following transformations,

$$\begin{eqnarray} p_R(x_j,\mathbf{r}_{nm}; E) & = U_K(\mathbf{r}_{nm}) p_K(x_j; E) U^{\dagger}_K(\mathbf{r}_{nm}), \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} p_K(x_j; E) & = U_M(j) G^ \lt _{j,j}(E) U^{\dagger}_M(j), \end{eqnarray} \tag{ 15 }$$

must be performed to recover the hole density from mode-space to real-space representation. $U_K(\mathbf{r}_{nm})$ is the unitary matrix [15] to transform from K-space to real-space and $\mathbf{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,

$$\begin{eqnarray} &&I(x_j) = \frac{e}{\hbar} \int \frac{\mathrm{d} E}{2\pi} \mathrm{{Tr}} \left[ \tilde{W}_M(j) G^ \lt _{j+1,j}(E) - G^ \lt _{j,j+1}(E)\tilde{W}^{\dagger}_M(j) \right]. \end{eqnarray} \tag{ 16 }$$

The Green's function elements $G^ \lt _{j,j’}$ are $N_M \times N_M$ block matrices.

3.1. 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,

$$\begin{equation} \Sigma^{R/ \lt }_{p,p’}(x_j,E) = i \, \mathrm{{Im}} \bigg\{ \sum_{n,n’} F^{n,p}_{n’,p’}(x_j) \, \overline{G}^{R/ \lt }_{n,n}(x_j,E) \bigg\}, \end{equation} \tag{ 17 }$$

where n ($n’$) and p ($p’$) are mode indices. For acoustic phonons,

$$\begin{equation} \overline{G}^{R/ \lt }(x_j,E) = G^{R/ \lt }_{j,j} (E), \end{equation} \tag{ 18 }$$

whereas,

$$\begin{align} \overline{G}^{ \lt }(x_j,E) = & N_{ph} G^{ \lt }_{j,j} (E-\hbar\omega_{ph}) + \left( N_{ph} + 1\right) G^{ \lt }_{j,j} (E+\hbar\omega_{ph}), \nonumber\\ \end{align} \tag{ 19 }$$

$$\begin{eqnarray} \overline{G}^{R}(x_j,E) = &N_{ph} G^{R}_{j,j} (E-\hbar\omega_{ph}) + \left( N_{ph} + 1\right) G^{R}_{j,j} (E+\hbar\omega_{ph}) \nonumber \\ & + \frac{1}{2}\left[ G^{ \lt }_{j,j} (E-\hbar\omega_{ph}) - G^{ \lt }_{j,j} (E+\hbar\omega_{ph}) \right], \end{eqnarray} \tag{ 20 }$$

are the lesser and retarded Green's functions for the hole-optical 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%.

3.2. 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,

$$\begin{equation} F^{n,p}_{n’,p’} = \sum_{d} \int \mathrm{d} \mathbf{r}_{\perp} \zeta^{\dagger}_n(\mathbf{r}_{\perp}) D_d \zeta_p(\mathbf{r}_{\perp}) \cdot \zeta_{p’}(\mathbf{r}_{\perp}) D^{\dagger}_d \zeta_{n’}(\mathbf{r}_{\perp}), \end{equation} \tag{ 21 }$$

where the x_j was omitted for brevity and d = x, y ,z. Being a vector of size N_b × 1, $\zeta^{\dagger}_n$ is the real-space wave function for the mode n at x_j and . It is retrieved from the following basis transformation,

$$\begin{equation} \zeta_{n}(\mathbf{r}_{\perp}) = U_K(\mathbf{r}_{\perp}) \Lambda_n ,\end{equation} \tag{ 22 }$$

where $\Lambda_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].

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₂. The gate lengths are $L_{\mathrm{g1}} = L_{\mathrm{g2}} = 7$ nm, being 1 nm distance from each other. The source/drain is highly doped with $N_{\mathrm A} = 10^{20}$ cm⁻³ 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_{\mathrm{FS}} =$ 0 and $E_{\mathrm{FD}} =$ 0.1 eV, respectively. Then, the energy window considered for quantum transport simulations is chosen as $E_{\mathrm{FS}}-1.0$ eV $\lt E \lt$ $E_{\mathrm{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~{\rm 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.

4.1. Band structure calculations

The K-space Hamiltonian has been discretized using a uniform mesh in all directions with mesh step a₀ = 0.2 nm, leading to $\tilde{N}_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_{\mathrm{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 $\langle 111 \rangle$ 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 \, (\ll 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.

Zoom In Zoom Out Reset image size

4.2. Ballistic limit

Figure 3 shows $I_{\mathrm{D}}-V_{\mathrm{GS1}}$ characteristics for different applied voltage at gate two, i.e. $V_{\mathrm{GS2}} = -$0.5, −0.6, −0.7and −1.0 V. The lower the $V_{\mathrm{GS2}}$, the higher the current in the whole range of $V_{\mathrm{GS1}}$ shown here. One can also observe that all curves start saturating at $V_{\mathrm{GS1}}\lesssim-0.6$ V, whereas in the ON-state, the current increases approximately one-order of magnitude at each value of $V_{\mathrm{GS2}}$. By applying a voltage at each gate, the effective potential barrier in the two-gated region is lowered and the dopant energy levels (DELs) are pushed up to higher energies. The latter are observed in the local density of states (LDOS) plotted in figure 4 for different states of the device in figure 1. One can observe how the DELs can be controlled by the applied gate voltages. The white line is the potential averaged with the hole concentration: $V_{\mathrm{av}}(x) = \frac{1}{p_{K}(x)} \int \mathrm{d} \mathbf{r}_{\perp} p (x,\mathbf{r}_{\perp}) V (x,\mathbf{r}_{\perp}).$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4(a) shows the LDOS when the applied voltage is the same on both gates, $V_{\mathrm{GS1}} = V_{\mathrm{GS2}} = -0.6$ V. One can observe that the DELs of the second dopant are slightly shifted upwards by approximately 17 meV. The latter is due to the effect of the external source-to-drain bias. In equilibrium conditions, i.e. when $V_{\mathrm{DS}} = 0$ V (not shown here), DELs from both dopants perfectly align to each other in terms of energy.

In figure 4(b), where $V_{\mathrm{GS1}} = -0.6$ V and $V_{\mathrm{GS2}} = -1.0$ V, the highest DELs of each dopant, are separated by approximately −0.18 eV despite the absolute value of $V_{\mathrm{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_{\mathrm{GS1}} = -0.9$ V and $V_{\mathrm{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_{\mathrm{GS1}}$ range for all values of $V_{\mathrm{GS2}}$ considered in this work.

Zoom In Zoom Out Reset image size

4.3. 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 six-band 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_{\mathrm{GS1}} = -0.6$ V and $V_{\mathrm{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.

Zoom In Zoom Out Reset image size

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_{\mathrm{GS1}} = -0.6$ V and $V_{\mathrm{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 back-scattered can also relax to higher energies in the direction of the source region.

Zoom In Zoom Out Reset image size

Figure 8 shows the impact of hole-phonon interactions on the $I_{\mathrm{D}}-V_{\mathrm{GS1}}$ characteristics for two different configurations: $V_{\mathrm{GS2}} = -0.6$ V (blue lines) and $V_{\mathrm{GS2}} = -1.0$ V (red lines). When compared to the corresponding ballistic case, the current is reduced over the entire range of $V_{\mathrm{GS1}}$ shown here. Their losses are plotted in the inset figure with the same color code and $V_{\mathrm{GS1}}$ range. The percentage loss is defined as $|I_{ball} - I_{scatt}| / I_{ball} \times 100\%$. Similar to current p-type MOSFETs, the lower the $V_{\mathrm{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_{\mathrm{GS1}} \geqslant -0.7$ V, within the sub-threshold part of the $I_{\mathrm{D}}-V_{\mathrm{GS1}}$ characteristics. The loss has a slight decrement at $V_{\mathrm{GS1}} = -0.9$ V and $V_{\mathrm{GS1}} = -0.7$ V in the case of $V_{\mathrm{GS2}} = -0.6$ V and $V_{\mathrm{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_{\mathrm{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_{\mathrm{GS1}} = -1.0$ V and $V_{\mathrm{GS2}} = -1.0$ V.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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 self-consistent Born's approximation, an expression for the hole-phonon interaction self-energy within the mode-space representation was introduced. We efficiently simulated a transistor-like 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.

The authors are most grateful to Mathieu Luisier for granting us access to the IIS computer infrastructure to perform the simulations leading to the results presented in this work. This research was funded by the Engineering and Physical Sciences Research Council (EPSRC), throughGrant Nos. EP/S001131/1 and EP/P009972/1. This project has also received funding from the EPSRC Impact Acceleration Account scheme under Grant Agreement No. EP/R511705/1 (Nano-Electronic Simulation Software (NESS)—creating the first open source TCAD platform in the world and Fast Track - Development boost for the Device Modelling group open-source NESS computational framework).

Conceptualization, HCN; methodology, HCN; software writing and implementation, HCN and CMB; original draft preparation, HCN; writing review and editing, HCN; CMB - writing and editing; supervision, VPG and AA; funding acquisition, CMB, VPG and AA. All authors have read and agreed to the published version of the manuscript.

In section 2, the six-band k · p Hamiltonian hasbeen given in terms of the three-band k · p Hamiltonian in $\big \{| X \rangle, | Y \rangle, | Z \rangle \big\}$ basis. As expressed in equation (1), information about the crystallographic direction is contained in the following group of matrices:

$$\begin{equation} A = \left(\begin{array}{ccc} c_t^2c_p^2 & 0 & 0\\ 0 & c_t^2s_p^2 & 0 \\ 0 & 0 & s_t^2\end{array}\right), \: B = \left(\begin{array}{ccc} s_p^2 & 0 & 0\\ 0 & c_p^2 & 0 \\ 0 & 0 & 0\end{array}\right), \end{equation} \tag{ 23 }$$

$$\begin{equation} C = \left(\begin{array}{ccc} s_t^2c_p^2 & 0 & 0\\ 0 & s_t^2s_p^2 & 0 \\ 0 & 0 & c_t^2\end{array}\right), \: D = 2c_t c_p s_p \left(\begin{array}{ccc} 1 & 0 & 0\\ 0 & -1 & 0 \\ 0 & 0 & 0\end{array}\right), \end{equation} \tag{ 24 }$$

$$\begin{equation} F = 2 s_p c_p s_t \left(\begin{array}{ccc} -1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 0\end{array}\right), \: G = -2 c_t s_t \left(\begin{array}{ccc} c_p^2 & 0 & 0\\ 0 & s_p^2 & 0 \\ 0 & 0 & -1\end{array}\right), \end{equation} \tag{ 25 }$$

$$\begin{equation} A’ = \left(\begin{array}{ccc} 0 & c_t^2 c_p s_p & -c_t c_p s_t\\ `\text{'} & 0 & c_t s_t s_p \\ \text{`}\text{'} & `\text{'} & 0\end{array}\right), \: B’ = \left(\begin{array}{ccc} 0 & -c_p s_p & 0\\ `\text{'} & 0 & 0 \\ 0 & 0 & 0\end{array}\right), \end{equation} \tag{ 26 }$$

$$\begin{align} C’ = \left(\begin{array}{ccc} 0 & s_t^2 c_p s_p & c_t s_t c_p\\ `\text{'} & 0 & -c_t s_t s_p\\ `\text{'} & `\text{'} & 0\end{array}\right),\nonumber\\ D’ = \left(\begin{array}{ccc} 0 & c_t(s_p^2-c_p^2) & -s_t s_p\\ `\text{'} & 0 & -c_p s_t \\ `\text{'} & `\text{'} & 0\end{array}\right), \end{align} \tag{ 27 }$$

$$\begin{eqnarray} F’ = \left(\begin{array}{ccc} 0 & s_t (c_p^2-s_p^2) & -s_p c_t\\ `\text{'} & 0 & -c_t c_p \\ `\text{'} & `\text{'} & 0\end{array}\right),\nonumber \\ G’ =\left(\begin{array}{ccc} 0 & -2 c_t c_p s_t s_p & c_p(s_t^2-c_t^2)\\ `\text{'} & 0 & s_p(c_t^2-s_t^2) \\ `\text{'} & `\text{'} & 0\end{array}\right), \end{eqnarray} \tag{ 28 }$$

where $s_\alpha$ and $c_\alpha$ represent sinα and cosα, respectively. For the sake of brevity, symmetric matrix elements have been omitted (''). In the spirit of a spherical coordinate system (k, θ, φ), φ (α = p) is the azimuthal angle between k_x -axis and the projection of the k − vector onto the $k_x-k_y$ plane. Whereas, θ (α = t) is the angle formed by the latter and the direction to which the k − vector points. Both angles can be then extracted from,

$$\begin{equation} \cos\phi = \frac{1}{l_x} \frac{1}{d_{l_x l_y \infty}}, \quad \cos\theta = \frac{d_{l_x l_y \infty} }{d_{l_x l_y l_z}} ,\end{equation} \tag{ 29 }$$

where l_i is the corresponding Miller index for i = x, y, z and the reciprocal distance $d_{l_x l_y l_z}$ has been defined as,

$$\begin{equation} d_{l_x l_y l_z} = \sqrt{\sum_{i = x,y,z} l^{-2}_i }. \end{equation} \tag{ 30 }$$

For the particular case in this work, l_x  = 1, l_y  = 1 and l_z  = 1, one can find that $c_p = 1/\sqrt{2}$ and $c_t = \sqrt{2/3}$. Finally, by means of the set of equations (23)–(30), other relevant crystallographic directions can be then easily incorporated in equation (1) for quantum transport simulations.

The real-space Hamiltonian given by equations (1)–(5) is first unitary transformed to K-space via [15],

$$\begin{equation} H_{K,L,L’}(x) = \int \mathrm{d} \Omega \, U_{K_L}(\mathbf{r}_{\perp}) H_R(-i\nabla, \mathbf{r}) U_{K_{L’}}(\mathbf{r}_{\perp}) ,\end{equation} \tag{ 31 }$$

where $\mathrm{d} \Omega = \mathrm{d} y \mathrm{d} z$ and

$$\begin{equation} U_{K_L}(\mathbf{r}_{\perp}) = U(\mathbf{r}_{\perp},\mathbf{k}_L) = \frac{2}{L_yL_z} \sin(k_p y) \sin(k_q z). \end{equation} \tag{ 32 }$$

$\mathbf{k}_L = \pi (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 $\mathbf{r}_{\perp} = (y,z)$, respectively. Carrying out the integral over the cross-sectional plane in equation (31), one straightforwardly obtains,

$$\begin{eqnarray} K^2_{y(z)}& = \int \mathrm{d} \Omega \, U_{K_L} k^2_{y(z)} U_{K_{L’}} = k^2_{p(q)} \delta_{pp’(qq’)} ,\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray} K_y & = \int \mathrm{d} \Omega \, U_{K_L} k_y U_{K_{L’}} = -\frac{4i}{\pi} k_{p’} \frac{p}{p^2-p’^2}\delta_{qq’} ,\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray} K_z & = \int \mathrm{d} \Omega \, U_{K_L} k_z U_{K_{L’}} = -\frac{4i}{\pi} k_{q’} \frac{q}{q^2-q’^2} \delta_{pp’}. \end{eqnarray} \tag{ 35 }$$

After a proper discretization in the transport x-direction, with N_x cross-sectional planes, the K-space Hamiltonian can be written as,

$$\begin{eqnarray} H_K = & \left(\begin{array}{cc} H’_{0~{\rm K}} & 0\\ 0 & H’_{0~{\rm K}} \end{array}\right) + \left(\begin{array}{cc} W’_K + W’^{\dagger}_K & 0\\ 0 & W’_K + W’^{\dagger}_K \end{array}\right) \nonumber \\ &+ \left(\begin{array}{cc} \mathbf{M} & \mathbf{N}\\ -\mathbf{N}^* & \mathbf{M}^* \end{array}\right) \delta_{pp’}\delta_{qq’} + \left(\begin{array}{cc} V’_K & 0\\ 0 & V’_K \end{array}\right), \end{eqnarray} \tag{ 36 }$$

which briefly can be re-expressed as follows,

$$\begin{equation} H_K(j) = H_{0~{\rm K}}(j) + W_K(j) + W_K^{\dagger}(j) + V_K(j), \end{equation} \tag{ 37 }$$

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),

$$\begin{align} H’_{0~{\rm K}} =& M \left(\frac{2}{\Delta x^2} + k^2_p + k^2_q \right) \delta_{pp’}\delta_{qq’} \mathbf{I_{3\times 3}} \nonumber \\ & + (L-M) \bigg( A \frac{2}{\Delta x^2} \delta_{pp’}\delta_{qq’} + B K^2_y + C K^2_z + F K_yK_z \bigg)\nonumber \\ & - N \bigg(A’ \frac{2}{\Delta x^2} \delta_{pp’}\delta_{qq’} + B’ K^2_y + C’ K^2_z + F’ K_y K_z \bigg), \end{align} \tag{ 38 }$$

$$\begin{align} W’_{K} =& -\frac{M}{\Delta x^2} \delta_{pp’}\delta_{qq’} \mathbf{I_{3\times 3}} \nonumber \\ & + (L-M) \bigg( A\frac{-1}{\Delta x^2} \delta_{pp’}\delta_{qq’} + D \frac{-i}{2\Delta x} K_y + G \frac{-i}{2\Delta x} K_z \bigg) \nonumber \\ & - N \bigg( A’\frac{-1}{\Delta x^2} \delta_{pp’}\delta_{qq’}+ D’ \frac{-i}{2\Delta x} K_y + G’ \frac{-i}{2\Delta x} K_z \bigg), \end{align} \tag{ 39 }$$

and the potential $V’_K$ in K-space can be numerically evaluated, for instance, as,

$$\begin{eqnarray} V_{K,L,L’}(j) = \sum^{N_yN_z}_{i = 1} \Delta \Omega_i U_{K_{L}}(x_j,\mathbf{r}_i) V(x_j,\mathbf{r}_i) U_{K_{L’}}(x_j,\mathbf{r}_i), \end{eqnarray} \tag{ 40 }$$

where $\mathbf{r}_i = (y_m, z_n)$ is the i-th grid point in real space.