A numerical approach to quasi-ballistic transport and plasma oscillations in junctionless nanowire transistors

A numerical framework for DC and RF small-signal simulations of nanowire transistors is presented, which is based on the self-consistent solution of the Poisson, Schrödinger, and Boltzmann transport equations and is stable for the entire range from weak to strong particle scattering. The proposed approach does not suffer from the deficiencies due to the transformation of the Boltzmann transport equation into the energy space and can handle the quasi-ballistic case. This is a key requirement for the investigation of plasma resonances and other high-mobility phenomena. The in-house solver is validated with results of a previously developed simulator based on the H-transformation for a conventional $$\hbox{N}^+\hbox{NN}^+$$


N
+


NN
+


 silicon transistor with strong scattering. Then, its results are compared with those of moments-based models and it is shown that these do not provide a satisfactory description of the electron dynamics in the quasi-ballistic transport regime. Furthermore, the internal boundary conditions of the transport models at the contacts are found to have a significant impact on plasma resonances and the physics-based thermal-bath boundary condition strongly suppresses them.


Introduction
Recently, scientific interest in solid-state terahertz (THz) devices (operating at frequencies between 300 GHz and 3 THz) has increased significantly [1][2][3][4][5]. The main idea, based on the shallow-water analogy put forward by Dyakonov and Shur [6][7][8], is to use field-effect transistors with high-mobility channels as nonlinear elements for the resonant detection of incident THz radiation. Moreover, the possibility of plasma instabilities leading to the emission of THz radiation by planar field-effect transistors, as pointed out by Dyakonov and Shur, has resulted in investigations by many groups (e.g., [9][10][11][12][13]). The question arises whether similar plasma instabilities could occur in nanowire or nanotube transistors. Although the methods developed previously for planar transistors could be applied to nanowire transistors, it is not clear whether the assumptions on which these models are based hold in the case of quasi-ballistic transport, which is a prerequisite for plasma instabilities [14]. Therefore, we discuss in addition to the drift-diffusion (DD) model, which is similar to the transport model used by Dyakonov and Shur, also the Boltzmann transport equation (BE) from which the DD model can be derived under rather stringent approximations [15].
Since the stabilization scheme based on the H-transformation for the BE fails for quasi-ballistic transport, a scheme based on the phase-space trajectories of the electrons is presented [16]. Furthermore, in the case of quasi-ballistic transport the stationary drain current of a MOSFET is limited by the finite injection velocity of the electrons from the source into the channel (e.g., [17,18]). This leads to a much more complex behavior of the source contact at high frequencies that can be described by the simple approximation made by Dyakonov and Shur and that the small-signal electron density must satisfy a Dirichlet boundary condition (D-BC) at the source side of the channel [6]. In order to describe the quasi-ballistic transport in such devices, it is not only necessary to go beyond the drift-diffusion approximation but to also include a more realistic contact model.
In the semiclassical framework, the BE provides the most precise information about the distribution of the charge carriers in the multidimensional phase space [19]. The most frequently applied method to solve the BE is the Monte-Carlo method [19], but there are no Monte-Carlo methods known to the authors that can perform small-signal simulations at complex frequencies required by the Dyakonov-Shur approach to plasma instabilities. On the other hand, deterministic solvers for the BE have been successfully used in the case of nanowire transistors (e.g., [20]), even including small-signal calculations [21]. In these cases, a transformation from the wavenumber to energy is used, which causes trouble in the case of a 1D k-space, because the density of states diverges at zero energy. This exacerbates the problems of small-signal calculations already encountered in the case of a 2D k-space [22]. Furthermore, in the ballistic limit the distribution function can be strongly asymmetric and discontinuous w.r.t. the wavenumber. This leads to further problems with the usual stabilization schemes based on the maximum entropy dissipation scheme [23] and the H-transformation [24]. These deficiencies are fatal in the numerical analysis of plasma resonances, where a numerically robust solver for the BE is required for handling the ballistic limit. In order to solve the above-mentioned problems, the BE is discretized directly in the phase space [16] and the quasiballistic limit of a nanowire transistor is investigated. The particular focus is on demonstrating a stabilization method, which is based on the method of characteristic curves and matrix exponentials. Thermal bath BCs, with the assumption of thermal equilibrium only for the electrons entering the device, give a more realistic description of the contacts [25]. The results are compared to the results of the DD equations.
The paper is organized as follows. Section 2 provides a description of the numerical framework used to evaluate the DC and small-signal behavior extending to the ballistic limit. In Sect. 3, the device structure is presented and the simulation results concentrating, in particular, on the small-signal drain self-admittance are discussed. Finally, in Sect. 4 the conclusions are drawn.

Theory
The numerical models are mostly constructed with the usual methods [16,[26][27][28][29][30], and only the new aspects are discussed. Since it is difficult to include the Pauli exclusion principle properly in moments-based transport models and it was neglected by Dyakonov and Shur, we will also neglect it in this investigation.
The device model for the silicon nanowire transistor consists of three partial differential equations, the Poisson, Schrödinger, and Boltzmann equations [20]. The Poisson equation is formulated for the 3D real space where (⃗ r) is the permittivity at position ⃗ r in real space, (⃗ r, t) the quasi-stationary potential at time t, q the positive electron charge, n(⃗ r, t) the electron density and N D (⃗ r) the donor concentration. The quasi-stationary approach is sufficient in this case, because even at 1 THz the wave length of light in silicon is about 100 μm and thus much larger than the device dimensions.
The device is sliced into N sec sections perpendicular to the transport direction (z-direction), and the quasi-stationary Schrödinger equation is solved for 2D cross sections parallel to the xy-plane at fixed positions for z under the assumption that the wave functions vanish outside of the silicon region. ℏ is the reduced Planck constant, m xx , m yy the electron masses in x and y direction assuming a ⟨100⟩ orientation of the silicon nanowire and parabolic bands, (⃗ r, t) the envelope function and sub (z, t) the subband energy. The superscript = (v, s) contains the valley index v of the minima of the conduction band and the subband index s. Again, the quasi-stationary approximation is used, because the largest part of the confining potential, which results in the boundary conditions, does not depend on time and due to the high symmetry of the device, which suppresses acoustic-wave-like plasma oscillations perpendicular to the transport direction.
The multi-subband 1D Boltzmann equation along the transport direction (z-axis) is given by where f (z, k, t) is the electron distribution function at wavenumber k, v (k) the group velocity and S{f } the scattering integral for single electron interactions [20] L sys is the 1D system volume and W , � (z, k, k � , t) the rate of a transition from state (k � , � ) into (k, ) calculated by Fermi's golden rule [31]. Unless otherwise stated, we use the macroscopic relaxation time approximation (RTA) to improve the consistency of the BE and the moments-based models where RTA is the macroscopic relaxation time and p eq (k) the one-particle Maxwell-Boltzmann distribution. The electron density is evaluated by The Poisson and Schrödinger equations are discretized with the standard methods (finite volume method, dimensional splitting, finite differences etc. [26]). The resultant linear systems of equations are solved with PARDISO [32] and the eigenvalue problems with FEAST [33]. The BE is directly discretized in the phase space, which is the Cartesian product of the position along the channel axis z and the axial wavevector k. An equidistant spatial grid z i = iΔz with i = 1, … , N z is introduced. In the k-space, a nonuniform grid with constant Δk near k = 0 and a Δk at higher energies corresponding to a constant step in energy is chosen. The number of grid nodes for negative wavenumbers is N k and the same for positive ones. For the sake of legibility, the time arguments are dropped from now on.
We start our formulation by considering the ballistic limit ( W , � (z, k, k � ) = 0 ), because this is the most critical case w.r.t. stability. Assuming a linear approximation for the subband energy sub (z) between the adjacent grid nodes z i and z i+1 , the force F i = −( sub,i+1 − sub,i )∕(z i+1 − z i ) is constant and the characteristic curves (ballistic electron path in the phase space) are parabolas. If the force is positive, the electrons are accelerated into the positive z-direction, and for k > 0 , an electron can move from z i to z i+1 (blue trajectories in Fig. 1). The wavenumber of an electron in state starting at z i with k z i is at position z i+1 2 m zz is the electron mass in transport direction. Thus, k z i+1 ≥ K i+ 1 ∕2 holds. For electrons starting at z i+1 and moving into the negative direction, we get The electrons can pass the barrier only, if their kinetic energy is larger than sub,i − sub,i+1 and thus k z i+1 ≤ −K i+ 1 ∕2 (red trajectories in Fig. 1). If this is not the case , the electron returns to position z i+1 , where its wavenumber has reversed sign (black trajectories in Fig. 1). The change in the wavenumber depends only on the difference in the subband energy and not on the exact force profile between both grid nodes. We can therefore assume that the subband energy is piecewise constant and changes abruptly at position z i+ 1 ∕2 = (z i+1 + z i )∕2 . Due to the Liouville theorem, the electron flux at position z i+ 1 ∕2 on the lefthand side of the step must equal the one on the right-hand side. For k L > 0 , where the superscript L indicates quantities on the LHS and R the ones on the RHS, we get with Thus, for the piecewise constant subband energy profile the characteristic curves are no longer continuous at z i+ 1 ∕2 and the values of the distribution function on the left and right sides of the step have to be distinguished. For k R < −K i+ 1 ∕2 , the electrons can pass the barrier from right to left resulting in the negative root for the wavenumber together with (10). For |k R | < K i+ 1 ∕2 , we get reflection Discrete equations are obtained by box integration in the k-space. The k-grid has 2N k grid nodes, which are symmetrically distributed around zero and the value zero itself is excluded from the grid. The jth box is given by k − j < k < k + j with its boundaries k ± j = (k j + k j±1 )∕2 and j = 1, … , 2N k . The interpolation function for the jth box is one inside the jth box and zero outside where Θ(x) is the Heaviside function. The distribution functions are assumed to be piecewise constant with the value f ,j (z) for the jth box On the LHS/RHS of the step in subband energy, we have Particle flux conservation is mapped onto the discrete k-grid by box integration. This allows to discretize (10) together with (13) for j > N k (positive wavenumbers) The first term on the RHS describes the flux of electrons from the left to the right and the second term electrons reflected at the step. The integration limits are defined as If we assemble the N k values of the distribution function for positive wavenumbers into a vector on the LHS of the step ⃗ f Lp,i+ 1 ∕2 , on the RHS ⃗ f Rp,i+ 1 ∕2 and for negative wavenumbers p is invertible and does not depend on the position in real space for j, j � = 1, … , N k . The other two matrices are Flux conservation for particles moving from right to left results with (10) and (12) in The second term on the RHS is due to the assumption of a maximal positive wavenumber k + 2N k , and all particles, which were to exceed this value while moving from left to right √ , are instead reflected to ensure flux conservation (dashed black line in Fig. 1). Furthermore, a similar set of equations can be formulated for a positive step in energy. Under stationary conditions, the ballistic distribution function does not depend on position for a constant subband energy and within a box of the real space grid Due to the TB-BCs, the distribution function of the incoming particles is given by the equilibrium one ⃗ f eq [25] The Fermi energies of the contacts are chosen such that the space charge densities are zero within the contacts at equilibrium. This closes the system of equations, and in most cases, a unique solution can be obtained for ⃗ f p,i , ⃗ f n,i . 1 Since the Pauli exclusion principle occurs in the BE only in the scattering integral, it can be included in the ballistic case by using Fermi functions instead of Boltzmann distributions in the boundary conditions.
The solution becomes unique for arbitrary subband profiles, if scattering is taken into account, because scattering couples all states for a given position [20]. Since scattering couples different subbands and valleys, the index is dropped and the vector of the distribution function contains from now on the distribution functions for all values of . Due to the assumption of a piecewise constant subband energy, the force is zero within a box of the real space. Integration of the BE over a box in k-space yields with the scattering integral Due to scattering, the solution depends on the z-coordinate within a box of the real space, even if we assume that the transition rate itself does not depend on position within the box. With ⃗ f i = ⃗ f (z i ) , Ŝ i =Ŝ(z i ) and the matrix exponential [35], the solution is within the ith box of the real space grid With the abbreviation Ĉ i (z) = (z i − z)D −1 vŜi , the solution at the boundaries of the ith box in real space is These equations replace (25), (26) of the ballistic case. Thus, it is possible to calculate the distribution function on the leftand right-hand side of the steps based on the distribution function on the ith grid node and a complete set of equations can again be assembled. The matrix exponential can be calculated with the methods described in Ref. [35], but for large matrices the numerical accuracy might not be sufficient. In this case, Expokit together with an iterative solver, which requires only multiplications of the matrix exponential with a vector, yields results with a higher precision and is more CPU efficient [36]. If the Pauli exclusion principle were included in the scattering integral, (29) would become nonlinear and could be integrated by ODE-solvers.
For small-signal analysis, we assume a time dependence for the solution quantities given by x(t) = x + ℜ{ x exp(st)} , where x is the large-signal stationary solution, x the complex phasor of the small-signal solution and s = + i the complex frequency, where is the real part describing a damped ( < 0 ) or growing ( > 0 ) solution and = 2 f the imaginary part corresponding to a harmonic oscillation with the frequency f. The time derivative in the BE can therefore be replaced by s, and within a box of the spatial grid, the linearized BE has two additional terms compared to the stationary case (29) D k is a diagonal matrix, which is due to box integration in k-space and contains the k-space box volumes Since the scattering integral depends on the subband energies and wave functions, linearization of the scattering integral yields the term on the RHS of (34). The corresponding derivatives are calculated by the chain rule and stationary perturbation theory for the SE [30]. With the abbreviation

the solution of the small-signal BE is obtained
With this result, the values of the small-signal distribution function at the boundaries of the real-space box can be calculated by numerical means. Linearization of (20) yields and (24) similarly. These equations together with the smallsignal PE and the small-signal version of the boundary conditions (27), (28) can be used to assemble a linear system of equations for the complete system. For the calculation of the small-signal terminal currents, a formulation of the Ramo-Shockley theorem is used that is consistent with the presented numerical framework [21,37]. Thus, admittance parameters can be calculated as a function of the complex frequency s, as required by the analysis of plasma instabilities.

Results and discussion
The test structure with a gate all-around architecture is shown in Fig. 2. The rectangular channel cross section is 5 nm × 4 nm with a gate oxide thickness of 1 nm (further details are given in Table 1). First, simulation of the ballistic case is demonstrated, because this is the numerically most challenging case. Since scattering is neglected, particles cannot transfer between valleys and subbands, and only a single parabolic subband with m xx = m zz = 0.19m e , m yy = 0.98m e is simulated ( m e is the electron rest mass). The k-space grid has a constant spacing of Δk = 0.0225 nm −1 at low energies and at higher ones a spacing corresponding to an energy of 2.5 meV. The maximal energy is 1.01 eV. The structure described in Table 1 is simulated with a grid spacing in real space of 0.5 nm, and the doping profile is shown in Fig. 3. To demonstrate the stability of the approach, the electron density in the device is shown for different gate and drain voltages under ballistic conditions (Fig. 3). The electron density can vary by more than 10 orders of magnitude without any problems. In Fig. 4, distribution functions in the k-space are shown near the drain and the source/channel barrier (top of the subband energy). The distribution function is positive in the whole phase space and shows no spurious oscillations. Due to the larger number of electrons moving from the left to the right than into the other direction, the distribution function is strongly asymmetric. Since the electrons do not loose energy by scattering when moving from source to drain, they gain qV DS = 100 meV in energy and a second maximum occurs in the distribution function near the drain at positive wavenumbers. Although the cutoff by the source/channel barrier is abrupt, the left flank of the second maximum near the drain is not as abrupt as at the top of the barrier. This is due to the finite k-space grid, which leads to artificial diffusion in the k-space, because a k-space box on the LHS of a step of the subband energy in real space is connected to multiple boxes on the RHS and vice versa (see (17)).   This is a fundamental problem of tensor-product grids in phase space and occurs regardless of the treatment of the derivative in k-space, because any formulation of the derivative in k-space involves multiple grid nodes. In principle, a finer grid in k-space could reduce the diffusion, but a finer k-space grid requires a finer real space grid to prevent numerical problems and a certain level of diffusion cannot be avoided. On the other hand, the impact of the artificial diffusion is weak, the phase-space solver based on the tensor-product grid is stable in the ballistic case, and the problems due to the H-transformation in the case of the small-signal parameters are avoided. In silicon nanowires, the mobility is rather low due to strong scattering resulting in diffusive transport [34], and it is not clear whether the numerical methods developed for the ballistic case will work. To investigate this, a silicon nanowire is simulated. The electron bandstructure in silicon is modeled in the conventional way accounting for six ellipsoidal and parabolic valleys, and five subbands are considered in each valley. As for the scattering mechanisms, various inter-valley and intra-valley electron-phonon interactions are included according to [38] and the deformation potential of elastic phonon scattering is adjusted in highdoped regions to reproduce the results of [20]. The simulations are performed at room temperature, assuming a ⟨100⟩ orientation in transport direction. The calculated I DS − V DS curves are presented in Fig. 5 (red curves), and the results are smooth.This shows that also strongly diffusive transport can be simulated. In order to validate the model, the results are compared to consistent simulations by other means. The BE is solved in addition based on the H-transformation [21] and by projection onto the first 10 Hermitian polynomials (TM) [14], and the results are almost identical. At V DS = 0.2 V , the calculations by the H-transformed BE and moments equations differ by ≈ 1.5% for V GS = 0.6 V and ≈ 0.8% for V GS = 0.5 V , where an energy grid with a uniform spacing of approximately 2.6 meV has been used for the H-transformed BE. Figure 6 shows the I DS − V GS characteristics of the device for V DS = 0.1 V and different discretization parameters. The artificial diffusion in phase space translates to artificial carrier heating and increases the subthreshold slope of the device. This effect, however, is not very significant, and our simulations show that even for a coarse k-grid ( Δk = 0.055 nm −1 ), the subthreshold slope changes from SS = 60.5 mV∕dec to SS = 66.5 mV∕dec as we move from a very coarse spatial grid ( Δz = 2 nm ) to a very fine one ( Δz = 0.2 nm ). Refining the k-grid alleviates this problem and we have SS = 62.5 mV∕dec for Δk = 0.0225 nm −1 . Although Δk cannot be reduced below a certain value for a given Δz , the spurious increase in the subthreshold slope can be reduced to such a level that it is negligible (at least at room temperature). Keeping this in mind, our results confirm the validity of all three approaches for the simulation of conventional devices with relatively low mobilities.
Next, we investigate the quasi-ballistic limit of the nanowire transistor. In the rest of this paper, single-subband and single-valley transport together with the RTA for the scattering integral is assumed for the sake of CPU efficiency and consistency of the BE and the momentsbased models. In the case of high mobilities, the impact of scattering is weak anyway and the error due to the RTA is small for small V DS . In order to check the accuracy of our BE solver, the ballistic small-signal drain self-admittance Y DD is calculated for different grid refinement factors and the relative error w.r.t. the finest grid is shown in Fig. 7. It is important to note that the refinements of the z-and k-grids have contrary impacts on the numerical stability of the discretized equations, which leads to a Courant-Friedrich-Lewy-like condition for Δz and Δk . While the k-grid should be fine enough to capture discontinuities, the z-grid must be refined accordingly to prevent numerical instabilities. Hence, Fig. 7 presents the relative error with simultaneous refinement of both z-and k-grids. The error decreases for finer grids and the method converges.
We investigate the drain self-admittance as a function of the frequency, because it is an important quantity to determine the impact of plasma waves on the device behavior. For example, a negative real part of the drain self-admittance corresponds to an instability and could enable the generation of THz waves [7,39]. Results of the BE are compared to the ones calculated with projections onto different numbers of Hermitian polynomials and boundary conditions, where the first two Hermitian polynomials result in the DD model. In order to obtain a plasma instability, specific bias conditions have to be applied to the contacts of the device. The gate/source port should be short-circuited and the drain/source port open (i.e., the drain self-admittance should be zero). These bias conditions should not be confused with the boundary conditions of the transport models at the contacts inside the device, which are determined by the contact model.  Fig. 8 (strong damping of plasma waves) and = 10 6 cm 2 ∕Vs (quasi-ballistic case) in Fig. 9 for TB-BCs and 10 for D-BCs. For lower mobilities, good agreement between the two methods is observed over the entire frequency range and the details of the source and drain boundary conditions do not play an important role, i.e., even the drift-diffusion (DD) model with D-BCs can give reasonably accurate predictions in the diffusive regime for frequencies f < 3 THz . Although the DD model deviates from the results of the BE solver at higher frequencies (e.g., it has an error of ≈ 8% at f = 8 THz ), it performs well in predicting the overall behavior of the real part. Truncation of the Hermite polynomial expansion at higher orders and TB-BCs improve the accuracy of the results. However, when the mobility is increased to = 10 6 cm 2 ∕Vs , moments of the BE fail to even provide a qualitative description of the device's small-signal behavior (Fig. 9). Although the model with ten moments (TM) does a better job than the DD model, the overall error below 5 THz is rather large. The DD model not only shows an artificial peak at about 5 THz, it also yields a far too large admittance at zero frequency. This is a well-known deficiency of the DD model [40,41]. The situation becomes worse, if we change the boundary conditions of the moments-based models from TB-BCs to D-BC, which were used by Dyakonov and Shur (Fig. 10). The unrealistically large values of ℜ{Y} at zero frequency increase by another 2-3 orders of magnitude. In addition, the peaks of the real part of the admittance become much sharper and higher and the minima smaller. The latter result is important, because the Dyakonov-Shur instability corresponds to a zero of the drain self-admittance (pole of the drain self-impedance for a short-circuited input). In Fig. 11, the absolute value of the drain self-admittance is shown for the DD model with different BCs, where the maxima and minima correspond to the poles and zeros, respectively. For D-BCs, the poles and zeros are lined up at = −1∕2 RTA [6], which is very small for = 10 6 cm 2 ∕Vs and results in sharp peaks on the imaginary axis. The three poles with an imaginary part between 4 and 10 THz correspond to the three peaks of the DD model in Fig. 10 in that frequency range. On the other hand, TB-BCs shift these poles (and zeros) to the left and the plasma resonances are strongly damped. This reduces their impact on the self-admittance at = 0 which is why the results in Fig. 9 are much smoother than in Fig. 10. For a plasma instability poles or zeros with > 0 are required. Unfortunately, the more realistic TB-BCs move the poles and zeros far to the left making plasma instabilities even more improbable. Since the DD model fails for quasi-ballistic transport anyway, the absolute values of the drain self-admittance calculated by the BE are shown in Fig. 12. The poles and zeros of the BE also occur at large negative real parts of s, and even a nonzero drain/source bias does not result in a significant shift to the right.
In Fig. 13, the nonequilibrium behavior of ℜ{Y DD (i2 f )} is shown for various drain/source voltages as a function of frequency. The BE results are compared to the tenth-order model. As it is evident at zero frequency, the transistor operates in the saturation regime for V DS > 0.2 V ( Y DD (0) → 0 ). Increasing the drain bias shifts the peak at f = 8 TH z to higher frequencies (see Fig. 12). The momentsbased model follows this behavior up to V DS = 0.1 V , whereas for higher biases it becomes unstable and cannot produce meaningful results. The BE, on the other hand, can be solved without problems for larger voltages and its results do not show any active behavior (i.e., ℜ{Y DD } < 0 ) for drain voltages larger than zero. The above results clearly show that the internal contact BCs of the transport model play an important role and that the more realistic TB-BCs strongly damp plasma resonances.
Since the damping of the plasma resonances depends on the type of the BCs, it would be interesting to investigate the impact of the D-BCs onto the BE results. The D-BCs are applied in the case of the moments-based models to the even moments (density, energy density (temperature), etc.). In the case of the BE, this would correspond to the application of D-BCs to the even parts of the distribution functions. Without scattering, such a boundary condition violates the Liouville theorem, because the ratio of the even part of the distribution function at the source and drain at energies above the top of the barrier is fixed by transport for an inversion symmetric bandstructure and cannot be imposed by the BCs. Thus, D-BCs for the even part of the distribution function and the ballistic BE are incompatible.

Conclusion
In this paper, a numerically stable approach for phase space discretization of the BE was proposed, which is applicable to the diffusive as well as ballistic regime. The presented BE solver was used to calculate the small-signal self-admittance of an N + NN + nanowire transistor up to THz frequencies, and the results were compared to the ones by moments-based models. D-BCs were shown to result in spurious plasma resonances for weak scattering, and moments-based models fail at describing the smallsignal behavior of quasi-ballistic devices. This calls into question the credibility of DD or higher-order models for the investigation of high-mobility phenomena in electronic devices.