Formulation of Ground States for 2DEG at Rough Surfaces and Application to Nonlinear Model of Surface Roughness Scattering in nMOSFETs

Electron mobility in extremely-thin-body (ETB) nanosheet channels and at cryogenic temperature is known to be dominated by surface roughness scattering. However, the conventional model of surface roughness scattering lacks accuracy because it requires the use of excessive roughness parameters to represent the experimental results. One of the main difficulties for the surface roughness scattering model is that the higher-order perturbations should be accurately included in the model because the surface roughness scattering is a strongly nonlinear phenomenon. Therefore, in this study, the formulation of ground states of two-dimensional electron gas (2DEG) at rough surfaces is derived by introducing a concept of the space-averaged perturbation Hamiltonian. This revised formulation of 2DEG at rough surfaces is different from the conventional solution for 2DEG at the flat surface. The space-averaged perturbation Hamiltonian is invisible in the linearized perturbation system, while its effect is significant in the system with the nonlinear perturbation energy. We combine the revised 2DEG formulation with a nonlinear model of surface roughness scattering and calculate the 2DEG mobility of the bulk Si and ETB Si-on-insulator (SOI) nMOSFETs. As a result, the experimental mobility of bulk and ETB SOI nMOSFETs is well explained in a wide temperature range of 4.2 to 300 K by using the roughness parameters experimentally obtained by transmission electron microscopy (TEM), which also supports the understanding of mobility at cryogenic temperature. The revised nonlinear model reveals that surface roughness scattering under the present model is 13 times stronger than that predicted by the conventional linear model.


I. INTRODUCTION
An extremely-thin-body (ETB) nanosheet channel is one of the most promising structures for the future advanced technology nodes of complementary metal-oxide-semiconductor (CMOS) devices because of its superior immunity to short channel effects and compatibility with the three-dimensional stacking process [1], [2], [3]. Here, thickness scaling of the nanosheet channel is necessary to improve gate controllability, while mobility degradation is a severe challenge for channel thickness scaling [4], [5], [6], [7], [8]. Therefore, the high-mobility materials robust to surface roughness scattering in ETB channels have recently gained great interest. In parallel, CMOS circuits at cryogenic temperature have recently attracted strong interest as the control circuits of the quantum computer [9]. Since the dominant scattering mechanism of electrons in ETB channels and at cryogenic temperature is surface roughness scattering [10], [11], [12], [13], [14], [15], [16], the quantitative understanding of surface roughness scattering becomes increasingly important. However, the conventional model of surface roughness scattering [10], [11], [12], [13], [14], [15], [16] has the problem that the surface roughness is overestimated in simulations to represent the experiments in comparison with the roughness parameters measured by transmission electron microscopy (TEM) [17], [18], [19]. As a result, it is difficult with the conventional model of surface roughness scattering to predict and assess the mobility in ETB channels for future channel materials alternative to Si.
Particularly, in the conventional Prange-Nee model of surface roughness scattering, the linearized perturbation Hamiltonian is employed to simplify the calculation [10], [11], [12], [13], [14], [15], [16]. However, the actual physical picture of surface roughness scattering is strongly nonlinear and asymmetric phenomena, which cannot be quantitatively modeled by such a linear model with only the first-order perturbation. Here, the nonlinearity means that the matrix element is not proportional to the amount of roughness. This fact leads to the excessively large roughness assumed to represent the experimental mobility. On the other hand, a nonlinear model of surface roughness scattering has been recently proposed [20], [21], [22], [23], which can include the higher order perturbations. This model has explained the experimental results with reasonable roughness parameters. These previous studies [20], [21], [22], [23] are valuable as the first introduction of the nonlinear model. However, their formulations are still insufficient because the nonlinear perturbation effect on wavefunctions of confined electrons, which does not occur for the linear perturbation, is neglected. Therefore, it is more difficult to represent the surface-roughness (SR)-limited mobility accurately for electronic systems with large scattering potentials and high kinetic energies of electrons, such as ETB III-V channels with the light effective mass. For example, the kinetic energy term in Hamiltonian, which is large for ETB III-V channels, is unreasonably neglected in the previous nonlinear model [20], [21], [22], [23]. In fact, under a nonlinear perturbation model, where the space-averaged perturbation Hamiltonian is not zero, the steady state of electrons is also affected by the perturbation Hamiltonian.
Therefore, in this study, we present a new formulation of surface roughness scattering of the two-dimensional electron gas (2DEG) at metal-oxide-semiconductor (MOS) surfaces, where the nonlinear model is applied to 2DEG with revised ground states under the rough surfaces. Compared to a model based on first-principles calculations [24], our model within the framework of effective mass approximation and Fermi's golden rule has advantages in terms of the computational cost and the convenience. The SR-limited mobility of bulk and Si-On-Insulator (SOI) n-channel metaloxide-semiconductor field-effect transistors (nMOSFETs) with Si/SiO 2 MOS interfaces is calculated by using our revised nonlinear model and compared with the experimental mobility of Si nMOSFETs, evaluated in a wide temperature range of 4.2 -300 K. It is found that our model can represent the experimental mobility of the bulk and SOI nMOSFETs by using the same and realistic roughness parameters, which are in good agreements with TEM analysis. Also, by comparing calculated and experimental mobility at cryogenic temperature, we have suggested the possibilities that the current model of the screening effect on surface roughness scattering is insufficient and/or that the conduction of electrons in the 4 valley of Si can be dominant by that in the low-mobility tail states.

A. GROUND STATES OF 2DEG AT ROUGH SURFACE
The schematic view of surface roughness scattering is shown in Fig. 1. Electrons are scattered more strongly by the oxide potential, when the height of roughness is positive. Therefore, surface roughness scattering is regarded as a strongly nonlinear and asymmetric scattering process with respect to . In order to calculate the mobility accurately, such nonlinearity needs to be included in the model, indicating that the higher-order perturbation should be included in the model of surface roughness scattering. Here, the Prange-Nee model has been most widely used among the previous studies of surface roughness scattering [10], [11], [12], [13], [14], [15], [16]. In the Prange-Nee model, however, only the first-order perturbation Hamiltonian is considered into account, and the direct scattering by the oxide potential is neglected. Therefore, the Prange-Nee model is lacking in accuracy.
Recently, the need and first construction of nonlinear modeling in surface roughness scattering has been reported [20], [21], [22], [23]. However, as will be discussed in the next Section II-B, the previous model is still insufficient in its original form. This is mainly attributed to the fact that the nonlinear effect on the ground states of the wavefunction is fully neglected. This effect becomes more serious as the nonlinearity of perturbation becomes stronger, for example in III-V and ETB channels. This problem can be solved by considering the nonlinear effect on the ground states of 2DEG. In fact, under a nonlinear perturbation model, where the space-averaged perturbation Hamiltonian is not zero, the steady ground states of electrons are also affected by the perturbation Hamiltonian. In this section, therefore, we derive the ground states of 2DEG with a rough MOS surface. In the next section, it will be shown that our formulation of the ground states is necessary for the nonlinear model of surface roughness scattering.
Since the mean free path of electrons (few tens nm) is much longer than the correlation length of surface roughness (∼ 1 nm), conducting electrons should feel the effective potential with the averaged perturbation. Therefore, the wavefunction of 2DEG, ξ(z), is reasonably approximated as uniform on the conduction plane. Then, the energy in the z direction, E z,i , and the in-plane averaged energy, E z,i r , in subband i can be expressed by where · · · r is an average operator on the plane as Here, X(z, ) is a roughness-dependent physical quantity, and f 1 ( ) is the probability distribution function of roughness . H is only the r-dependent term in (1), which is the Hamiltonian of a channel with roughness (r) as where m z (z, ) is a confinement effective mass at point z, V(z, ) is a potential, and is the reduced Plank constant. When the wavefunction ξ i (z) is the eigenstate of H r , E z,i r corresponding to the total energy in the system is minimized. Therefore, the energy operator of the ground state of Hamiltonian is H gnd = H r , which is the total energy including the perturbation. It should be noted that H gnd = H r is only a pseudo-steady state for electrons that are conducting and scattered. The concept of the averaged Hamiltonian is also verified by the fact that the mean free path of electrons is much longer than the correlation length of surface roughness.
When the perturbation has a linear relationship to the height of roughness as H = H 0 + H 1 × , where H 1 is the first-order perturbation, the averaged Hamiltonian H r amounts to H 0 . This is implicitly assumed in the Prange-Nee model. However, H r = H 0 for the nonlinear perturbation with respect to surface roughness. Here, H r can be derived by Effective mass and potential abruptly change at a semiconductor/oxide interface. Therefore, all X(z, ) r terms in (5) near the interface are different from X(z, 0), while X(z, ) r far from the interface is almost equal to X(z, 0). In the numerical calculation of · · · r , the physical quantity is smoothly connected at the interface by the sigmoid function S(z) as where z 0 and the mesh size is 0.02 nm, which is small enough for the calculation results to converge because the impact of the shape of the Sigmoid function is negligible on the average quantities · · · r . Also, V q and V ox are the electrostatic potential and the energy barrier between a semiconductor and a gate oxide, respectively, m ox and m sz are the confinement effective mass of oxide and semiconductor, respectively, and T B is the channel thickness. The potential V q should be calculated for each by solving the Poisson equation. In order to calculate the in-plane average · · · r in (3) and (5), the probability distribution function f 1 ( ) should be accurately modeled to fit the real MOS interface roughness. The autocorrelation of the roughness is generally modeled by the Gaussian or exponential formulation [10], [11], [12], [13], [14], [15], [16], [17], [18], [19], [20], [21], [22], [23]. On the other hand, f 1 ( ) is actually complex due to the atomic discreteness. However, f 1 ( ) is also modeled by the Gaussian or exponential formulation to simplify the calculation in this study as where rms is a root-mean-square (RMS) value of surface roughness. In this study, the exponential formulation of The examples of the potential and the wavefunction of the lowest subband are shown in Fig. 2. Here, the black and red curves correspond to them under the completely flat and rough interface, respectively. rms of 0.2 nm is used for all calculations. The apparent potential of V(z) r has a less steep slope at the interface because it is the averaged one. The wavefunction of the 2 valley in the bulk (100) Si/SiO 2 interface in Fig. 2(a) moves away from the interface by taking into account the roughness, to minimize the perturbation energy. On the other hand, the change in the wavefunction of the 4 valley in Fig. 2(b) by roughness is much smaller because the wavefunction spreads more widely due to lighter m sz . For the 5-nm-thick InAs/Al 2 O 3 interface in Fig. 2(c), the penetration of ξ(z) into the oxide is large because of the very light m sz . This large penetration makes it more difficult to calculate the matrix element of surface  roughness scattering for III-V in the Prange-Nee and the previous nonlinear models.

B. MODELING METHOD OF SURFACE ROUGHNESS SCATTERING
In this section, the calculation method of the matrix element and the necessity of the averaged Hamiltonian including the effect of surface roughness on the ground states are described. The matrix element is the most important parameter to determine the mobility limited by surface roughness scattering. The matrix element from subband i to j is expressed by where the ground state of Hamiltonian is The (generalized) Prange-Nee model is commonly used in the previous studies, where the linear relationship of M ij ( ) ∝ is assumed [10], [11], [12], [13], [14], [15], [16]. Therefore, the matrix element in the generalized Prange-Nee model is derived by Although the expression of (13) is very simple and easy to calculate, the Prange-Nee model lacks accuracy for the strongly nonlinear and asymmetric nature of surface roughness scattering, shown in Fig. 1.
In the previously proposed nonlinear model [20], [21], [22], [23], M ij ( ) in (10) was numerically calculated for all values by using the formulation of H gnd by (11). Then, any higher order perturbations can be included in the nonlinear model. On the other hand, in our model, M ij ( ) is calculated for all by using the expression of H gnd by (12), which includes the roughness effect on the ground state of Hamiltonian.
Here, i = i / rms for the Gaussian formulation, while i = i / rms for the exponential formulation in (17). It is worth pointing out that both the Gaussian or exponential formulation can be independently applied to C (r) and f 2 ( 1 , 2 , ρ). In this study, the exponential formulations of f 1 ( ), C (r), and f 2 ( 1 , 2 , ρ) are simply used. The power spectrum density of the matrix element is derived by the two-dimensional Fourier transform of C M (r) as where q is the wavenumber and J 0 (qr) is the zeroorder Bessel function [20]. It should be noted that when C M (r → ∞) = 0, the integral in (19) diverges to infinity and is very oscillatory as |M ij (q)| 2 → ±∞. However, in the previous nonlinear model [20], [21], [22], [23], since C M (r → ∞) = 0, the calculation of (19) might not have been correctly done, and the DC component in (19) Therefore, M ij ( ) r = 0 is a necessary condition for the matrix element model. It should be noted that this condition is not visible in the linear Prange-Nee model because it always holds. While the condition of M ij ( ) r = 0 is only derived from the convergence condition in (19), M ij ( ) r has also the important physical meaning of being the expected value of the perturbation energy. Thus, M ij ( ) r should be zero, when applied to the perturbation theory. On the other hand, one may consider that M ij ( ) in (15) can be replaced by M ij ( ) − M ij ( ) r because, statistically, the autocorrelation function is determined by subtracting the average of the original function. Of course, it is possible to calculate However, we believe that this correction lacks the physical adequacy of the perturbation. The matrix elements calculated by the Prange-Nee model [14], the previous nonlinear model [20] and our model are shown in Fig. 3. Here, M ij ( ) of the lowest subband as i = j = 1 is shown. The values of M 11 ( ) calculated by the Prange-Nee model are so different from those calculated by the nonlinear ones. As clearly seen here, the space-averaged perturbation energy, M 11 ( ) r , is positive in the previous nonlinear model, in spite of the fact that M ij ( ) r = 0 is the necessary condition as discussed above. On the other hand, M ij ( ) = 0 at = 0 in our model because the Hamiltonian at the geometrical center of = 0, H 0 , is different from the perturbative average of the Hamiltonian. As a result, M ij ( ) r in our model is strictly equal to zero, as proven in the following.
A characteristic example of the quite nonlinear M ij ( ) is shown in Fig. 3(c) for a 5-nm-thick InAs/Al 2 O 3 interface, where m sz is very light. The matrix element calculated by the generalized Prange-Nee model in Fig. 3(c) is not realistic because the derivative of M ij ( ) at = 0 is coincidentally very small due to a large contribution of the kinetic energy term. In the previous nonlinear model, the amount of the matrix element in < 0 is also quite large and positive because the kinetic energy term of Hamiltonian is very large for InAs due to light m sz . Therefore, such kinetic energy of the first term in (4) was unreasonably omitted in the calculation of the matrix element in the previous studies [20], [21], [22], [23], which is significant for III-V with light m sz . However, this way undermines the physical understanding of surface roughness scattering. On the other hand, since our model can include higher-order perturbations and the kinetic energy term in Hamiltonian, it can be applied to materials with light m sz such as III-V semiconductors.
The autocorrelation of the matrix element, C M (r), for electrons in the 2 valley at the (100) Si/SiO 2 interface, calculated by the previous nonlinear [20] and our models, is shown in Fig. 4. C M (r → ∞) = 0 for the previous nonlinear model because M ij ( ) 2 r = 0 as shown in Fig. 3(a). This nonzero C M (r → ∞) in the previous nonlinear model causes the divergence of |M ij (q)| 2 in (19). Therefore, the averaged Hamiltonian H r should be used for the ground states at rough surfaces. It is confirmed, on the other hand,  for our model that C M (r → ∞) = 0. The energy difference of the lowest subband between H 0 and H r , (E c − E c0 ), is shown in Fig. 5. Here, (E c − E c0 ) is the additional energy in the system caused by the existence of the nonlinear perturbation due to surface roughness. This additional increase in the quantization energy can appear as an additional threshold voltage shift from the expected value at the flat interface, which was actually observed experimentally in Si nMOSFETs [4].

C. SURFACE ROUGHNESS AT TOP AND BOTTOM INTERFACE
In the previous section, we consider only the roughness of the one-side interface in ETB channels. However, in the real ETB channels, there is surface roughness at the top and bottom interfaces. In the previous studies [16], [20], the matrix elements introduced by the roughness of the top and bottom interfaces are added when the roughness at the top and bottom interfaces are uncorrelated. However, it should be noted that, in our model, the matrix elements cannot be simply added even under the assumption that the roughness at the top and bottom interfaces are uncorrelated. This is because the existence of the roughness at both interfaces itself modulates the wavefunction. As similar to the discussion in the previous section, the matrix element is where t and b are the height of roughness at the top and bottom interfaces, respectively. Here, H gnd is the averaged Hamiltonian H t , b r , where the average operator is Also, ρ tb is the Pearson's product-moment correlation coefficient between the top and bottom interfaces. In this study, ρ tb = 0 is reasonably assumed to simplify the calculation. When ρ tb = 0, f 2 can be simplified as where the integration range is (−∞, ∞) in (26) and (27), and f 4 (· · · ) is the four-variate joint probability distribution function. When the top and bottom interfaces are uncorrelated as ρ tb = 0, f 4 (· · · ) can be simplified by where ρ t (r) and ρ b (r) are the Pearson's product-moment correlation coefficient of the top and bottom interfaces, respectively. Although the expression of C M (r) by (27) is a strict one, it is not realistic to directly calculate (27) because it is extremely time-consuming. Therefore, we approximated the autocorrelation C M (r) by where rms,t and rms,b are RMS of the surface roughness at the top and bottom interfaces, respectively. When RMS at the one side interface is taken to be zero, C M (r) can be calculated by (15), which makes the calculation cost realistic. However, it is worth to note that the approximation of C M (r) by (29) underestimates the real value of C M (r).

D. SCREENING EFFECT
The matrix element M scr ij (q) including the screening effect is delivered by using the scalar dielectric function D (q) as [16], [29] where ε s is the dielectric constant of a semiconductor. The definition of the form factor F i,i (q) and the polarization factor i,i (q) are same as in [29].

E. MOMENTUM RELAXATION TIME
Solving the Boltzmann transport equation (BTE) for anisotropic valleys is a complex problem. Therefore, in this study, we use an approximated solution of the linearized BTE under a weak electric field in isotropic valleys. The momentum relaxation time, τ sr,i (k), in subband i is [14], [16] 1 where g v,j is the degeneracy of subband j. Also, v x,i and k i are the velocity of x direction and the wavenumber of subband i, respectively. In an isotropic valley, the momentum relaxation time at energy E can be obtained by solving the following equations [20], [29] where the matrix A is expressed by Here, k i and k j are the wavenumber of the initial state in subband i and the final state in subband j, respectively. k i = |k i |. m d,i is the density-of-state effective mass of subband i. α i and E ci are the nonparabolicity factor and the ground energy of subband i, respectively.

F. PHONON SCATTERING
The contribution of phonon scattering is pronounced for a thick channel and low carrier density at room temperature [11], [30], [31]. Also, phonon scattering in (100) SOI cannot be ignored even in the 3-nm-thick channels [4], [25], [30]. The momentum relaxation time of intravalley scattering by acoustic phonon for electrons in subband i, τ ac ph,i , is [11], [16], [32] 1 τ ac ph,i (E) where U(E) is the step function, k B is the Boltzmann constant, T is temperature, ρ m is the mass density of the crystal, s l is the longitudinal sound velocity, and D ac is the deformation potential of acoustic phonon. The momentum relaxation time of intravalley scattering by surface optical (SO) phonon for electrons in subband i, τ sop ph,i , is [11], [33], [34] 1 τ where s = 1, 2 denotes one of the two SO phonon modes.
where, D k , E k , and N k are the deformation potential, the energy, and the number of the k-th intervalley phonon, respectively. In this study, scattering between the two

G. MOBILITY
The total momentum relaxation time, τ i , and the mobility, μ total,i , for electrons in subband i are [20] 1 where m x,i and N s,i are the conduction effective mass and the areal density of electrons in subband i, respectively. The total mobility is The parameters used for the mobility calculation are summarized in Table 1 [11], [32], [35]. Here, the tunneling effective mass is used as m ox of SiO 2 [35].

A. MOBILITY OF BULK SI nMOSFETs
The influence of the subband number included in the mobility calculation is shown in Fig. 6. It is found that when more than 3-subbands are included, the mobility becomes saturated. Therefore, 5-subbands are considered in each valley for all the calculations in this study. The E eff dependence of SR-limited mobility, μ sr , is shown in Fig. 7. The experimental mobility at T = 300, 77, and 4.2 K is also shown in Figs. 7(a), (b) and (c), respectively [36]. Here, the experimental mobility at 4.2 K in Fig. 7(c) is measured in this study for the identical samples. Although the E eff universality of μ sr for different acceptor concentrations (N A ) is not simply supported by the nonlinear model because of its complicated formulation, the present simulation results show the good E eff universality of μ sr as shown in Fig. 7(a). The constant mobility of μ ni = 2×10 4 cm 2 /Vs limited by neutral impurity scattering is included to fit with the experimental mobility at 4.2 K in Fig. 7(c). It should be noted that only the 2 valley is considered in Fig. 7(c) to eliminate the effect of the unreasonably large mobility of the 4 valley electrons due to the screening effect as discussed later in Fig. 8. The total mobility is in excellent agreement with the experimental one at 77 and 4.2 K by using the same roughness parameters of rms = 0.2 nm and = 0.8 nm. On the other hand, the Prange-Nee model needs to represent the experimental results by using rms = 0.5 nm and = 0.8 nm. The roughness parameters used in the simulation and evaluated by TEM are summarized in Table 2. Here, the exponential formulation in (9) and (16) is used in this study because the roughness spectra measured by TEM are closer to the exponential ones [18], [25]. The roughness parameters in our model are certainly in good agreement with the TEM analysis [17], [18], [19]. On the other hand, it should be pointed out that rms evaluated by TEM might not be too accurate because of the projection effect [17], [18], [19] and the existence of discrete atomic images [25]. TEM images are obtained as averaged ones along the transmission direction of electrons and, thus, the rms value by TEM is apparent and smaller than the real value, which is called the projection effect. In addition, there is a lower limit on rms evaluated by TEM due to the existence of discrete atomic images. The lower limit of measurable rms is roughly a/4 ∼ 0.14 nm, where a is a lattice constant. These two effects make it difficult to evaluate the real small rms by TEM. However, we would believe that the present rms of 0.2 nm is very realistic because the rms of initial substrates evaluated by atomic force microscopy is typically less than 0.2 nm [37], [38], and the thermal oxidation process might not degrade the interface roughness so much [7], [39].
However, it is still questionable the validity of our mobility model at cryogenic temperature. The calculated SR-limited mobility at 300 and 4.2 K with and without the screening effect is shown in Fig. 8(a). Also, the scalar dielectric   function D (q), which is the screening term in (30), is shown in Fig. 8(b). The screening effect is negligible at 300 K and for the 2 valley electrons at 4.2 K because D (q) ∼ = 1 for the large wavenumber q. On the other hand, the mobility of the 4 valley electrons is extremely large at 4.2 K due to the strong screening effect, which is far from the experimental mobility in Fig. 7(c). There are two possible interpretations for the low experimental mobility of the 4 valley electrons as shown in Fig. 8(c). First, the model of the screening effect on surface roughness scattering is inaccurate. In the current screening model, the potential change is induced by the change of the electron distribution functions due to the perturbations. On the other hand, such a potential change by the screening effect also makes the ground state Hamiltonian H r closer to H 0 , resulting in the enhancement of the amount of the matrix element, as shown in Fig. 3. This is because the unscreened wavefunction, which is the eigenstate of H r , moves away from the interface to decrease the scattering by surface roughness. However, the establishment of more accurate screening modeling of the ground state is still a challenge.
The second interpretation, which is more likely, is the significant decrease in the mobility due to the occupation of the electrons in the tail states with low mobility. The electron occupancy and the Fermi level E f − E c in the 4 valley are shown in Fig. 9. The Fermi energy of the 4 valley is less than 1 meV even for E eff < 0.5 MV/cm, while the energy width of the tail states has been estimated to be ∼ 1.3 meV [40]. Therefore, all the electrons can be occupied in the tail states at 4.2 K, where the mobility should be low because the tail states can originate from the Lifshitz tail [41] and are different from electronics sates in a standard extended band. The energy width of 1.3 meV [40] is estimated for the 2 valley, whereas the similar width can be also expected for the 4 valley.
The SR-limited mobility of electrons in each valley at 77 K is shown in Fig. 10(a). μ sr of the 4 valley electrons is higher than μ sr of the 2 valley electrons, even though the conduction mass m x is 0.19m 0 and 0.315m 0 for the 2 and 4 valleys, respectively, where m 0 is the free electron mass. Here, the temperature is taken to be 77 K to simplify the discussion by suppressing the electron occupancy of higher subbands. The high μ sr in the 4 valley is attributed to the two screening mechanisms. First, the screening effect in (30) is strong for the 4 valley electrons because the Fermi wavenumber is very low. Second, the z-direction electric field is also screened by the 2 valley electrons. The effective electric field of each valley is shown in Fig. 10(b). The effective electric field at each subband i, E eff ,i , and the total effective electric field, E eff , are determined by where E(z) is the electric field. E eff in the 2 valley is much higher than E eff in the 4 valley because the spread width of the wavefunction of 2 valley electrons is narrower due to the heavy m sz and the electron occupancy is higher due to the lower quantization energy. Therefore, μ sr in the 4 valley can be higher than μ sr in the 2 valley at low temperature because of lower E eff for electrons in the 4 valley. Experimentally, the total effective electric field can be given by N s and the coefficient η as follows.
where eN dpl is the charge density of the depletion layer. The coefficient η, determined by the slope of E eff , amounts to 0.59, 0.32, and 0.5 for the 2 , 4 , and total valleys, respectively, whereas the experimental η of the total valley is known to be 0.5 for (100) Si [36]. Also, η for the 4 valley is smaller in the low N s region than 0.32. It should be noted that η of the total valley defined by (46) is always 0.5 in theory [16]. Therefore, in the bulk nMOSFETs with the multi-valley electron occupation, the existence of the valley with heavy m x does not simply result in mobility reduction. On the other hand, in ETB channels, μ sr in valleys with light m sz is severely degraded as discussed in the next Section III-B. The temperature dependence of μ sr of electrons in all the valleys and E eff ,i of electrons in the lowest subband is shown in Fig. 11. Here, the screening effect is temporary ignored because of the too strong influence on electrons in the 4 valley at cryogenic temperature. However, it is not a bad approximation, as shown in Fig. 8(a). While μ sr is independent of temperature in the single-valley and single-subband conditions (data is not shown), μ sr for the multi-valley structures is actually dependent on temperature. As temperature decreases, the electron occupancy in the subbands with more widely-spread wavefunctions decreases because of the larger quantization energy, leading to the reduction in E eff ,i in the lowest subband of the 2 valley. As a result, the total mobility increases with a decrease in temperature. Since the screening effect becomes also stronger at lower temperatures, the tendency for μ sr to increase at low temperature still remains unchanged.

B. MOBILITY OF ETB SOI nMOSFETs
The SOI thickness dependence of the phonon-limited, the SR-limited, and the total mobility under N s = 3 × 10 12 cm −2 is shown in Fig. 12. The experimental mobility [4] is well explained by the calculated values. Here, the validity of our model is also supported by the fact that the same roughness parameters are used in Figs. 7 and 10. The SOI thickness dependence of μ sr under N s = 3 × 10 12 cm −2 is shown for various rms in Fig. 13(a). The general trend is similar among various rms . For the channel thickness of less than 4 nm, which is narrower than the spread width of the wavefunction of electrons in the 2 valley, as shown in Fig. 2(a), μ sr rapidly drops in proportion to T r B with r = 5∼6. This strong T B dependence in ETB channels is well known as thickness fluctuation scattering expressed by [42] In the linear model with the infinite energy barrier, μ sr is proportional to T 6 B . On the other hand, in our model with a finite energy barrier, the factor r = d(ln μ sr )/d(ln T B ) is lower than 6, as shown in Fig. 13(b). The rms and dependence of μ sr for bulk (100) Si/SiO 2 nMOSFETs under N s = 5 × 10 12 cm −2 is shown in Figs. 14(a) and (b), respectively. Here, μ sr is proportional to −2 in our model, as similar to that in the linear model in (45). On the other hand, μ sr is proportional to −2.7 rms and −2 rms in our model and the Prange-Nee model, respectively. As a result, surface roughness scattering has a 13 times stronger effect in our model than that predicted in the conventional Prange-Nee model. Therefore, there is little room for improvement of μ sr in Si nMOSFETs because the rms value of 0.2 nm at Si MOS interfaces is sufficiently small. The E eff dependence of μ sr for the 2-6-nm-thick SOI is shown in Fig. 15. In the lower E eff region, μ sr drops more significantly with channel thickness scaling, attributed to the wider spread width of the wavefunction with lower E eff . μ sr increases very slightly with channel thickness scaling in the high E eff region because of the reduction of the electron occupancy in the 4 valley.
The comparison of μ sr between the single-gate and double-gate structures under N s = 3 × 10 12 cm −2 is shown in Fig. 16. The mobility is slightly lower for the double-gate structure because the electron occupancy in the 2 valley  with light m x is lower. On the other hand, the mobility is almost the same between the single-gate and double-gate structures in the ETB channels less than 3 nm because the electron occupancy in the 2 and 4 valleys is determined by the channel thickness. It should be noted that the  formulation of the screening effect by the scalar dielectric function in (30) is not accurate for the double-gate structure, and the tensorial dielectric function should be used [22], [43]. However, the scalar dielectric function is used in this study to reduce the calculation cost. Fig. 17 shows μ sr of 2 and 4 valley electrons. The channel thickness at which μ sr starts to drop is thicker for the 4 valley than the 2 valley because of the wider wavefunction of 4 valley electrons. Therefore, the utilization of the anisotropic valley with the combination of light m x and heavy m sz is important to suppress the influence of surface roughness scattering in ETB channels, as seen in (48). It should be recalled here that there is little room for improvement of μ sr in Si nMOSFETs because the rms value of 0.2 nm at Si MOS interfaces is sufficiently small. Therefore, the introduction of other materials with the anisotropic valley is very important for the future CMOS technology nodes.

IV. CONCLUSION
Since surface roughness scattering is a strongly nonlinear phenomenon, higher-order perturbations must be taken into account to quantitatively predict mobility. We have proposed the new surface roughness scattering model on a basis of the re-consideration of the ground states for 2DEG with rough MOS surfaces and a combination with the nonlinear surface roughness scattering probability. This revised nonlinear model of surface roughness scattering can represent the experimental mobilities well by using the realistic roughness parameters obtained from the TEM analysis. It has been revealed that surface roughness scattering has a much stronger influence than that predicted by the Prange-Nee model. As a result, there is little room to improve μ sr in Si MOSFETs because the rms value of 0.2 nm at Si MOS interfaces is sufficiently small. Thus, the introduction of high mobility materials with anisotropic valleys is much important. Then, the quantitative nonlinear model in this study allows us to predict mobility in ETB nanosheet channels and to provide an assessment for future channel materials in the advanced CMOS technology.