Analytical Modeling of 3D NAND Flash Cell With a Gaussian Doping Profile

The incessantly increasing demand for highly dense storage medium in this era of big-data has led to the development of 3D NAND Flash memories. 3D NAND Flash based SSDs have revolutionized edge storage and become an integral part of the data warehouses and the cloud storage systems. The conventional 3D NAND flash memory with greater than hundred stacked word-line (WL) layers suffer from channel tapering and non-uniformity in the threshold voltage of cells in different WL layers which necessitates the use of different programming voltages for different WL layers and a complex error-correction circuitry. A non-uniform vertical (Gaussian) channel doping profile alleviates the threshold voltage non-uniformity and leads to a significant reduction in the complexity of the error-correction circuitry and a simple (uniform) programming scheme for all WLs. Although behavioral models have been proposed to utilize 3D NAND flash memory for circuit and system-level applications, an analytical model is important to understand the intricate details of the device physics and propose design guidelines for efficient cell design. To this end, in this paper, for the first time, we have proposed an analytical model for the characteristic length, surface potential and inner potential of the Macaroni-body 3D NAND flash cell with vertical Gaussian doping profile. A strong agreement between the analytical model and the TCAD simulations for different gate lengths (down to 25 nm), inner and outer radius, channel and oxide thickness of the 3D NAND flash cell validates the efficacy of the developed model.


I. INTRODUCTION
The dramatic shift from the compute-centric paradigm to the data-centric paradigm in this era of big data has led to a surge in the flow of asynchronous data [1], [2]. Ultra-dense storage technologies are urgently required to cope up with this data explosion [3]. Therefore, to achieve the bit density levels required for efficient storage of big data, conventional planar 2D NAND flash cells were incessantly scaled [3], [4], [5]. However, the performance of planar flash cells deteriorates significantly with shrinking of the feature size owing to the considerably large neighboring cell coupling disturbances and crosstalk [6], [7], [8], [9]. To overcome these challenges, several architectures exploiting vertical (3D) stacking The associate editor coordinating the review of this manuscript and approving it for publication was Cristian Zambelli .
Moreover, to achieve an ultra-high bit density, several (≥176) word line (WL) layers are stacked in commercial 3D NAND flash arrays [20]. This leads to a significant increase in the string height (mold height). Etching a narrow hole in such stacks with a high aspect ratio (typically greater than 40:1) without substantial tapering of the channel region from the source select line (SSL) to the bit select line (BSL) is a technological challenge [17], [21]. This undesirable tapering of the channel region results in cells with different channel thickness at different location along the string [18]. Therefore, the fabrication process for 3D NAND flash induces an inherent non-uniformity in the intrinsic threshold voltage of the cells at different WL layers and necessitates a complex error correction circuitry [17]. To mitigate this non-uniformity in the threshold voltage across the charge-trap flash cells in a string, several techniques such as a non-uniform vertical (Gaussian) channel doping profile, compensated blocking oxide thickness, positiondependent programming/erase voltage or duration, etc. were proposed [17].
Furthermore, to facilitate the utilization of Macaroni body 3D NAND flash array for circuit and system-level applications, behavioral compact models of 3D NAND flash memory cells with uniform channel doping were proposed [19], [22], [23], [24]. Although the behavioral compact model [19] accurately captures the static string current characteristics and the parasitic capacitance components applicable to the Macaroni body 3D NAND flash memory cells with uniform channel doping, it fails to describe the intricate device physics. Therefore, an analytical model is required to understand the electrostatic characteristics and propose design guidelines for optimum performance of Macaroni body 3D NAND flash cells. Analytical models are more accurate since they are derived from the first principles utilizing the fundamental laws of physics such as Maxwell equations and Poisson's equation along with appropriate boundary conditions emanating from Gauss law. These models provide additional information regarding the electrostatics and microscopic properties which are not evident from the experimental results. The behavioral compact models treat the device under test (DUT) as a black box and utilize curve fitting techniques to mimic the macroscopic experimental characteristics by tweaking the empirical parameters. Although even a behavioral compact model may also provide some physical insight, its accuracy while performing any parametric analysis such as investigating impact of change in the structural parameters on the output/transfer characteristics is dubious since the empirical parameters are fitted only to a limited experimental data. Therefore, analytical models are more important for device physicists and technologists to understand the microscopic properties, whereas behavioral compact models are more useful for circuit designers [25].
To this end, an analytical model for the surface potential, inner potential, and the characteristic length (which provides a measure of the electrostatic integrity) of a GAA Macaroni body MOSFET with a uniform doping profile in the channel region were also formulated [12], [13], [26]. However, to mitigate the impact of process-induced substantial tapering of the channel region and the non-uniform threshold voltage in 3D NAND flash cells without introducing significant fabrication overheads or complexity, a vertical Gaussian doping profile in the channel region appears promising [17]. In this technique, to realize a uniform threshold voltage distribution across the string, a high dopant dose is applied at the source end of the channel region to increase the threshold voltage of the cells and a low dopant dose is applied at the drain end of the channel region to reduce the threshold voltage of the cells [17], [27], [28], [29]. Although the analytical model developed in [12] predicts the electrostatic behavior of the conventional uniformly doped 3D NAND flash, the model cannot be extended to provide the design guidelines for the 3D NAND flash cell with non-uniform doping profile presented in [17]. Moreover, to the best of our knowledge, an analytical model for Macaroni body 3D NAND flash cell with a vertical Gaussian doping profile in the channel region is still elusive.
Therefore, in this work, for the first time, we develop an analytical model for the 3D NAND flash cell with uniform threshold voltage distribution and non-uniform (vertical Gaussian) doping profile which may help in providing the necessary design guidelines for further optimizing its performance. We solve the 3D Poisson's equation in cylindrical coordinates with appropriate boundary conditions to derive an analytical model for the characteristic length, surface potential, and inner potential of a Macaroni body 3D NAND flash cell with a vertical Gaussian doping profile in the channel region. Furthermore, we analyze the accuracy of the developed model by comparing results obtained from the analytical model against the data obtained from 3D TCAD simulations. A strong agreement between our model and the TCAD simulations for different structural parameters such as doping concentration of the channel region, gate length, oxide thickness, inner and outer radius validates the efficacy of the developed model. Although the insights provided by the developed analytical model can be also be obtained with the help of TCAD simulations, considering the computational complexity and the resource (processing power and memory)-intensive requirement of the TCAD simulations, it is better to perform the design exploration across the structural and design parameters utilizing the proposed analytical model which provides accurate results while consuming only a fraction of time and resources of the TCAD simulations. Moreover, as the BSIM-CMG model does not allow the user to define a non-uniform doping profile, the compact model developed for conventional 3D NAND flash cell with uniform doping profile utilizing BSIM-CMG model [19] may not be used for design exploration of 3D NAND flash cell with the non-uniform (vertical Gaussian) doping profile.

II. DEVICE STRUCTURE
The 3D schematic view and the cross-sectional view of a Macaroni body 3D NAND flash cell with a vertical Gaussian doping in the channel region are shown in Fig. 1(a) and Fig. 1(b), respectively. As can be observed from Fig. 1, in a Macaroni body cell, the polysilicon channel is etched at the center and filled with a core dielectric filler to realize an ultrathin channel. This ensures a minimal threshold voltage variation due to the inherent grain boundaries and associated traps owing to a reduction in the volume of grain boundaries and traps present within the channel region. Since the thickness of the polysilicon layer is less than the depletion region width contributed by the gate electrode [10], the Macaroni body cell exhibits an enhanced gate control and a high electrostatic integrity.
The structural parameters used for the Macaroni body 3D NAND flash cell with a vertical Gaussian doping profile in this work are summarized in Table 1. The parameters of the Gaussian doping profile ( 2σ 2 )) were tuned to reproduce the non-uniform doping profile reported in [17] for minimizing the threshold voltage non-uniformity as shown in Fig. 2. As shown in Fig. 2 and Table 1, the peak doping concentration at the source end is N D ∼ 10 18 cm −3 and the channel doping concentration at the drain end is N D ∼ 10 15 cm −3 [17]. Moreover, such a doping profile can be realized experimentally utilizing the multi-implant process and tuning the implant energy and ion dose according to the desired doping levels as done in [17].

III. MODEL DESCRIPTION
To understand the impact of vertical Gaussian doping profile in the channel region of a Macaroni body 3D NAND flash cell, starting from the first principles, we have modelled the key parameters dictating the device performance such as the surface potential, inner potential and the characteristic length (which is a measure of scalability and electrostatic integrity [30], [31]). Since the structure of the Macaroni body 3D NAND flash cell is cylindrical, we use the cylindrical coordinate system instead of Cartesian coordinate system to formulate the analytical model. To obtain the characteristic length and the potential profiles, we solved the three-dimensional Poisson's equation.

A. CHARACTERISTIC LENGTH
The 3D Poisson's equation in cylindrical coordinate system can be expressed as: where r, θ , and z are the cylindrical coordinate axes along the radial, angular, and axial (along the channel length) directions, ψ is the channel potential, q is the electron charge, si is the dielectric constant of silicon and N D (z) is the donor concentration in the channel region of the Macaroni body cell. As discussed in Section II, the donor concentration N D (z) varies from the drain end to the source end along the vertical channel in the proposed 3D NAND flash cell as: where k and σ are constants tuned to reproduce the optimal doping profile in [17]. Due to the symmetric structure along the axial direction, the channel potential does not vary with the angle in the radial plane. Therefore, ψ is independent of θ i.e. ∂ψ(r,θ,z) ∂θ = 0. Equation (1) can then be simplified as a 2D Poisson's equation given by: Now, Young's parabolic potential approximation [32] can be utilized to simplify the 2D electrostatic potential distribution along the radial direction inside the channel region as: Since the silicon body is located between r 1 and r 2 (r 1 ≤ r ≤ r 2 ), we shift the origin to r = r 1 and express equation (4) as: where C 1 (z), C 2 (z) and C 3 (z) are arbitrary coefficients.
To find out the unknown coefficients in equation (5), we have assumed that the electric field lines inside the channel terminate at the silicon body-core filler dielectric interface i.e. the electric field at r = r 1 is zero: ∂ψ(r 1 ,z) Furthermore, we have considered the boundary condition i.e. the electric flux at the interface between the silicon film and the gate oxide must be continuous at r = r 2 : si Si tsi V gs − V fb − ψ s where, t Si is the channel thickness, V fb is the flat-band voltage, ψ s is the surface potential (ψ s = ψ (r = r 2 , z)) and C ox is the effective gate oxide capacitance per unit area of cylindrical Macaroni body cell given by [12]: where ox is the dielectric constant of oxide and t ox is the gate oxide thickness. Now, we define the potential at the channel-core filler interface (r = r 1 ) as the inner potential, ψ o = ψ (r = r 1 , z).
Utilizing these assumptions and boundary conditions, equation (5) can be simplified as: To find an expression for the surface potential, we use r = r 2 in equation (7) to get: Now, substituting the value of ψ s obtained from equation (8) in equation (7), we get: Substituting ψ (r, z) from equation (9) in equation (3) and using r = r 1 , we get a second-order differential equation for the inner potential ψ 0 (z) as: where λ is characteristic length of the Macaroni body NAND flash cell expressed as: To obtain an expression for the inner potential, we solved the differential equation (10) using the superposition principle where the solution consists of a complimentary function and a particular integral (PI ). The inner potential (ψ 0 (z)) can then be expressed as: where C 1 and C 2 are constants and PI represents the particular integral. The particular integral (PI ) of equation (12) can be obtained from equation (10) as: where D is the derivative with respect to z (along the axial direction). Although calculation of PI including the Gaussian term is mathematically challenging [29], considering the fact that the channel doping has a constant value at any particular position in the channel along the z-axis, the PI can be approximated as: The details regarding the calculation of PI are discussed in the Appendix. Now, to determine the constants C 1 and C 2 , we have used two boundary conditions: (a) the inner potential at the source terminal ψ 0 (z = 0) = V R , where V R is the difference between the donor Fermi level and the intrinsic Fermi level, expressed as V R = φ t ln N D n i and (b) the inner potential at the drain end ψ 0 z = L g = V R + V ds [12], where L g is the gate length. Applying these boundary conditions in equation (12), C 1 and C 2 are obtained as: 15) VOLUME 10, 2022 where K 1 and K 2 are constants that can be expressed as: Now, the exact expression of inner potential ψ 0 (z) can be obtained by substituting the value of C 1 , C 2 , and PI in equation (12) as: As can be observed from equation (19), the inner potential has a strong dependence on the doping density, gate voltage (V gs ), oxide thickness, and the difference of outer and inner radius (r 2 − r 1 ) (channel thickness).

C. SURFACE POTENTIAL, ψ s (z)
Now, an analytical model for the surface potential (ψ s = ψ (r = r 2 )) distribution of the Macaroni body 3D NAND flash cell can be obtained from equation (9) as: Si 8λ 2 (20) where ψ 0 (z) is given as: and C 1 and C 2 are constants. To estimate these constants, we have again used two boundary conditions: (a) the surface potential at the source terminal ψ s (z = 0) = V R and (b) the surface potential at the drain end ψ s z = L g = V R + V ds [12]. PI can be obtained following the same methodology (equation (14)) as discussed in the Appendix. Applying the aforementioned boundary conditions, C 1 and C 2 can be obtained as: where K 3 and K 4 are constants expressed as: Now, substituting the values of C 1 , C 2 and PI in equation (21) ψ 0 (z) can be expressed as: where K 5 and K 6 are constants expressed as: Si 8λ 2 (28) The exact expression for the surface potential ψ s (z) can then be obtained by substituting the value of ψ 0 (z) from equation (26) in equation (20) as: where K 7 is a constant given by:

IV. RESULTS AND DISCUSSION
To validate the developed analytical model for the inner potential and the surface potential of the Macaroni body 3D NAND flash cell, we have compared the results obtained from the analytical model against the TCAD simulations for different device parameters such as outer radius (r 2 ), inner radius (r 1 ), gate length L g , oxide thickness (t ox ) and gate bias V gs . Sentaurus TCAD [33] was utilized for performing device simulations. Carrier transport was modelled using the drift-diffusion formalism. The degradation in the carrier mobility due to dopant scattering and high electric field was considered using the doping-dependent mobility model and the normal field-dependent mobility model. The model parameters of the TCAD simulation setup were calibrated by reproducing the experimental static characteristics of the 3D NAND flash string with 10 WLs [34] as done in [19]. As can be observed from Fig. 3, the simulation results match with the experimental data of [34] validating the accuracy of the TCAD simulation setup used in this work.
To analyze the variation in the inner potential along the channel length (axial direction) of the Macaroni body 3D NAND flash cell for different channel thickness, we varied the outer radius, r 2 (inner radius, r 1 ) without changing the inner radius (outer radius) and other device parameters such as oxide thickness, channel length and doping profile. First, we fixed the inner radius at 13.5 nm and varied the outer radius to 17.5 nm t si 2 = 4 nm , 19.5 nm t si 2 = 6 nm , 21.5 nm t si 2 = 8 nm , and 23.5 nm t si 2 = 10 nm as shown in Fig. 4(a). Then, we varied the inner radius to 13.5nm t si 2 = 10 nm , 15.5nm  Fig. 4(b). As can be observed from Fig. 4, a reduction in the channel thickness leads to a larger dynamic range of the inner potential and a better gate control [35]. Although the characteristic length which represents the electrostatic integrity (equation (11)) is not a direct function of r 1 or r 2 individually, it depends significantly on the channel thickness (r 2 − r 1 ). A smaller value of (r 2 − r 1 ) implies a lower characteristic length representing an efficient gate control, a higher scalability and a lower leakage current [12]. Therefore, an ultra-thin Macaroni body should be utilized while designing the 3D NAND flash cell with vertical Gaussian doping profile with significantly enhanced electrostatic integrity. Moreover, an ultra-thin Macaroni body also ensures a minimal spatial (cell-to-cell) threshold voltage fluctuation owing to the reduced number of polysilicon grains in the channel region [10]. Also, r 2 dictates the footprint of the vertical Macaroni body 3D NAND flash cell and should be kept as minimum as possible. Furthermore, the results obtained from the analytical model match exactly with the TCAD simulations validating the proposed model.
We have also analyzed the impact of oxide thickness (t ox ) on the distribution of inner potential along the channel length (axial direction) as shown in Fig. 5(a). As t ox increases, the electrostatic integrity reduces due to an increase in the characteristic length (equation (11)). Therefore, the inner potential at the center of the channel region increases leading to a reduction in the dynamic range indicating a poor gate control. This leads to a complex design guideline for the 3D NAND flash cell with vertical Gaussian doping profile with respect to the oxide stack thickness. While an ultra-thin gate oxide stack (yielding a smaller characteristic length) should be used to increase the immunity against the short channel effects, incessant scaling of the gate oxide may degrade the retention characteristics owing to the stress-induced leakage current (SILC) [10], [36].
Furthermore, to understand the scaling prospects, we have also analyzed the impact of gate length L g on the VOLUME 10, 2022 distribution of inner potential along the channel length (axial direction) as shown in Fig. 5(b). As the gate length reduces, the influence of drain electric field on the channel region increases resulting in a reduction in the electrostatic integrity [35]. Therefore, the dynamic range of inner potential distribution is larger for the Macaroni body 3D NAND flash cell with higher gate lengths indicating a better gate control of the channel region. Since 3D NAND flash cell is essentially a vertical structure, we may utilize cells with a larger gate length and improved electrostatic integrity without increasing the area overhead. However, using 3D NAND flash cell with large gate length may restrict the number of word line layers which may be stacked owing to the fabrication constraints while etching deep trenches with high aspect ratio.
To understand the influence of inner and outer radius on the channel electrostatics, we investigated the surface potential distribution by selecting different combinations of the inner radius (r 1 ) and the outer radius (r 2 ) while keeping a constant channel thickness r 2 − r 1 = t si 2 in an attempt to decouple the impact of channel thickness and the inner and outer radius. We observe from Figs. 6 (a) and (b) (also   the channel thickness is constant. This clearly indicates that the surface potential is independent of the value of inner and outer radius and depends only on the channel thickness as also suggested by our analytical model (equation (29)).
While the channel thickness [2(r 2 − r 1 )] needs to be small, the outer radius (r 2 ) should also be kept as minimum as possible since r 2 dictates the footprint of the vertical Macaroni body 3D NAND flash cell. Moreover, as the gate voltage increases, the surface potential reduces gradually approaching its minimum value at the threshold voltage [13].

V. EXTENDING THE MODEL FOR CONSTANT DOPING PROFILE
The developed model can be extended to obtain an analytical model for 3D NAND flash cell with a uniform doping profile [12] by simply using z = 0 in equation (2). The final expressions obtained for the characteristic length, the inner potential and the surface potential corroborates well with the results presented in [12]. The results of the analytical model and TCAD simulations for the inner potential distribution for different parameters such as channel thickness, oxide thickness, and gate length are shown in Fig. 7. It is evident from Fig. 7 that the dynamic range of the inner potential in 3D NAND flash cell with uniform doping profile is significantly lower as compared to the dynamic range of the inner potential of 3D NAND flash cells with vertical Gaussian doping profile (Figs. 4 and 5) for same silicon thickness, oxide thickness, and gate length. This clearly indicates the enhanced electrostatic integrity in 3D NAND flash cell with a non-uniform (vertical Gaussian) doping profile. Therefore, the proposed analytical model helps in providing necessary design guidelines for realizing 3D NAND flash cell with optimal performance.

VI. CONCLUSION
In this work, for the first time, we derived an analytical model for the characteristic length, surface potential, and inner potential distribution of Macaroni body 3D NAND flash cell with vertical Gaussian doping profile. To analyze the efficacy of the developed model, we compared the results obtained from the analytical model against the TCAD simulations. A strong agreement between the analytical model and the TCAD simulations for different device parameters such as oxide thickness, channel length, inner and outer radius, channel thickness and doping profile validates the accuracy of the developed model. The design guidelines obtained from the developed model may provide incentive for efficient design of 3D NAND flash cells.

APPENDIX
The particular integral (PI) of equation (10) can be obtained as: is non-trivial and complex [29], considering the fact that the channel doping has a constant value at any particular position in the channel along the z-axis, the exponential term in equation (A4) yields a constant value. Therefore, the higher order terms of the differential operator reduce to zero and the final expression of particular integral (PI ) can be obtained as: