I. INTRODUCTION

The rapid growth of the big data era has driven an unprecedented increase in the demand for NAND flash memory. To achieve higher storage capacities, the number of stacked layers in 3D-NAND flash memory has steadily increased. Accurate prediction of current-voltage characteristics for next-generation NAND devices with such complex structures is essential. However, the unique architecture of 3D-NAND flash memory, distinct from conventional CMOS transistors, limits the applicability of standard BSIM-based SPICE models. As a result, technology computer-aided design (TCAD) simulations have become a widely adopted and reliable tool for analyzing the electrical and physical behavior of 3D-NAND devices ^[1-3]. However, TCAD simulation is time-intensive and lack the capability to evaluate memory operations in integration with peripheral circuits. Therefore, fast and accurate NAND modeling that can be used in circuit simulations has been actively researched. SPICE-based NAND modeling studies on conventional 2D planar structures have been sufficiently conducted ^[4,5], and research on 3D-NAND modeling is also being actively conducted ^[6,7]. In most approaches, the modeling process begins with a unit-cell using commercial transistor models and is subsequently extended to a NAND string.

Various models have been proposed for modeling transistors. Among them, the Berkeley short-channel IG-FET model (BSIM) ^[8] has been predominantly used. The BSIM model is a physics-based transistor model that accurately reflects the electrical characteristics of short-channel transistors. It provides a standard parameter extraction procedure, allowing for the adjustment of parameters to model transistors with desired characteristics. The cylindrical GAA structure of 3D-NAND can be perfectly implemented using the BSIM-CMG transistor model. However, as the BSIM model has evolved, the number of parameters has increased (${\sim}$1000 ea.), requiring skilled professionals and considerable time to fit them.

To address this, various techniques have been developed to automate parameter extraction ^[9-11]. For example, genetic algorithms (GA) and machine learning (ML)-based optimization techniques can solve this problem. The Genetic Algorithm (GA) is an optimization technique inspired by natural selection, which iteratively evolves solutions by selecting, crossing, and mutating high-fitness candidates from a population of random parameter sets. Using GA, BSIM-CMG parameters matching the desired data can be determined within a few hours. However, while GA-based parameter extraction for single devices has been extensively studied ^[9-11], its application to 3D-NAND strings remains unexplored.

This study presents a method for modeling the current-voltage (I-V) curve of 3D-NAND string using BSIM-CMG and optimizing its parameters through GA. To validate the proposed methodology, SPICE simulation results were compared with those obtained from Sentaurus TCAD ^[12] simulations. The method accurately represents the I-V characteristics of 3D-NAND structures with up to 500 layers including tapered channel holes.

II. CURRENT-VOLTAGE MODELING USING GENETIC ALGORITHM

Fig. 1 illustrates the modeling of the 3D-NAND string using BSIM-CMG employed in this study. As depicted in Fig. 1(a), the 3D-NAND string has a gate-all-around (GAA) structure. To describe the GAA structure of 3D-NAND, each unit cell was implemented through BSIM-CMG (GEOMOD=3) which describes the cylindrical gate transistor. String can be modeled by connecting unit cells in series and biasing bit line (BL) and source line (SL) voltages to each end node as depicted in Fig. 1(b).

Fig. 1. (a) The structure of the 3D-NAND string and (b) the modeling of the 3D-NAND string for circuit simulation using BSIM-CMG.

Fig. 2 presents the structure of the 3D-NAND string, and the parameters used in this study ^[6,13]. Synopsys Sentaurus TCAD was employed to simulate the reference I-V characteristics. The simulation reflected the actual 3D-NAND string geometry, using a tunneling oxide/trap layer/blocking oxide (SiO${}_{2}$/Si${}_{3}$N${}_{4}$/SiO${}_{2}$) gate stack within a macaroni channel structure filled with SiO${}_{2}$. For simplicity, the channel consisted of single-crystalline silicon, and a metal with a work function of 4.8 eV was used as the word-line (WL). Both the WL and spacer were set to 50 nm in length, as were the string select line (SSL), and ground select line (GSL). The simulation employed mobility degradation model, high-field saturation mobility model, and drift-diffusion model in cylindrical coordinates system.

Fig. 2. The diagram of the 3D-NAND string and parameters used in simulation.

The parameters of BSIM-CMG used can be optimized through GA. By selecting parameters that significantly affect the current-voltage characteristics (parameters related to gate capacitance, mobility, channel resistance, etc.), a total of 24 parameters were targeted for modeling. Table I lists some key BSIM-CMG parameters optimized through GA. By employing GA, a parameter set that best matches the reference I-V data can be obtained. The operational procedure of GA is shown in Fig. 3, and each step is described as follows.

Fig. 3. (a) The flow chart of GA and (b) description of the crossover and the mutation process in GA.

1. Initialization: The algorithm creates an initial population for GA by selecting the BSIM-CMG parameters to be optimized and generating 500 individuals with random parameter sets.

2. Fitness Evaluation: The algorithm evaluates the similarity of each individual to the reference data. A fitness score is assigned to each individual based on the following formula, where ${I}_{\rm BL,SPICE}$ represents the bit-line current obtained using BSIM-CMG in SPICE simulation, and ${I}_{\rm BL,reference}$ represents the bit-line current extracted from TCAD simulation. The calculation is performed across the ${V}_{\rm GS}$ range of the I-V graph.

3. Selection: The algorithm selects individuals for the next generation based on their calculated fitness scores. First, it ranks all individuals according to their fitness scores. Then, only the top 50% of individuals are retained and passed on to the next generation.

4. Crossover: The algorithm randomly selects two individuals from the chosen group and generates offspring through the crossover process. For randomly selected parameters, the parameter values of the two individuals are swapped to create two offspring. This process is repeated until the population size is restored. The procedure is illustrated in Fig. 3(b).

5. Mutation: The algorithm induces genetic mutations with a probability of 0.1%, as shown in Fig. 3(b). By introducing random changes to the genetic information of individuals, the algorithm helps maintain genetic diversity and prevents convergence to local minima ^[14].

6. Replacement: The algorithm replaces the previous population with the newly generated individuals. To enhance the performance of the GA, the top 4% of individuals from the previous population are retained and carried over to the next generation.

7. Termination Check: If the termination criteria (fitness score or iteration) are met, the process ends. Otherwise, the above steps are repeated.

Fig. 4 shows the procedure and results of parameter optimization using GA. To simulate a NAND flash memory array, it is more effective to use reference I-V data from a string structure with multiple WLs, rather than a single-cell structure, due to structural differences that necessitate additional calibration. However, simulating strings with a large number of WLs significantly increases the number of circuit instances during optimization. Therefore, appropriate string structures should be selected to ensure a balance between accuracy and computational cost. In this study, as shown in Fig. 4(a), reference I-V data for NAND strings with 1, 3, and 5 WLs were extracted through TCAD device simulation, and GA optimization was performed based on this data. Circuit simulation using the extracted BSIM parameters was conducted with Cadence Spectre 21.1 ^[15]. The average fitness score, defined as the total fitness score divided by the number of I-V curves used, improves over generations, as shown in Fig. 4(b). The I-V characteristic modeling was validated using the parameters of the highest-scoring individual (average fitness score: 0.124) from the final generation with. Fig. 4(c) shows the circuit simulation results using the GA-optimized BSIM parameters, which exhibit excellent agreement with the reference data obtained from the TCAD device simulation.

Fig. 4. The procedure and the results of the parameter optimization through GA. (a) The reference structure used in GA. The N indicates the number of WL$_{\rm s}$. (b) Evaluation of the fitness scores of individuals. Each symbol represents an individual. (c) The reference data extracted through TCAD (square symbol) and the calibrated data with optimized parameters through GA (red line). Unread WLs including GSL and SSL are biased with $V_{\rm PASS}$.

III. SIMULATION RESULTS

We conducted modeling for specific cases (1, 3, and 5 WLs under fixed operating voltage conditions), but we examined whether these results could be extended to strings with different numbers of WLs or other operating voltage conditions.

First, the accuracy of the I-V characteristics of NAND strings with a larger number of WLs was verified. As shown in Fig. 5, the I-V characteristics were successfully modeled for WL counts of 32, 64, 128, 300, and 500. This result demonstrates that the BSIM parameters were carefully and comprehensively selected to accurately reflect the physical behavior of the NAND string. As the number of WLs (stacked layers) increases, the string resistance increases, leading to a corresponding decrease in the current level. The inset clearly illustrates that the on-current decreases inversely with the number of layers. To evaluate the accuracy, fitness scores were calculated for all I-V curves. The highest score was 0.1362, which deviates only slightly from the initial G.A. optimization result of 0.124, indicating the consistency of the model. Furthermore, the on-current levels showed a maximum error of 1.5%, demonstrating that the model effectively captures structural variations induced by different numbers of WLs. Here, the on-current was extracted as the current value at a WL voltage of 6 V. The BL voltage was fixed at 0.7 V unless explicitly stated otherwise.

Next, simulations were performed for different BL voltages, as shown in Fig. 6. The results confirm that the I-V characteristics of the NAND string were well modeled across various BL voltage conditions.

Finally, the variation in I-V characteristics at different selected WL positions was observed. The I-V characteristics were extracted at the top (WL31), middle (WL16), and bottom (WL0) of a 32-layer string. Figs. 7(a) and 7(b) show that the modeling of I-V characteristic variation according to the WL position was well implemented in both linear and log scales. Fig. 7(c) illustrates the variation in ${G} _{\rm m}$ according to the WL position. It can be observed that the peak ${G}_{\rm m}$ increases as the WL position moves closer to the bottom of the string. This is a result of the decrease in the effective source resistance (${R}_{\rm S}$ in Fig. 7(d)) of the cell as the channel length below the cell becomes shorter ^[16]. Fitness scores were also evaluated for the case where the BL voltage and WL positions were modified. The highest fitness score among the curves was 0.1530, which does not deviate significantly from the initial G.A. optimization result of 0.124, indicating the robustness of the model.

When modeling the I-V characteristics of strings in commercially available 3-D NAND flash memory, an additional factor to consider is the shape of the channel hole. Due to the limitations of the etching process, the actual 3-D NAND string exhibits a tapered structure where the diameter of the channel hole decreases toward the lower layers shown in Figs. 8(a) and 8(b) ^[17,18]. The degree of taper in the string is represented by $\theta$, the angle between the vertical axis and the channel hole. When the channel hole is ideally etched, $\theta$ is 0$^\circ$; however, in actual 3D NAND, it can have a value of approximately 0.1$^\circ$ to 0.2$^\circ$. Due to the tapered shape of the string, the channel diameter (D) decreases toward the bottom of the string, resulting in an increase in channel resistance. The current characteristics of each cell according to the taper angle ($\theta$) can be implemented by additionally introducing the diameter (D) parameter in the BSIM-CMG transistor model. As shown in Fig. 8(c), accurate I-V modeling is achieved by representing the diameter value of each corresponding cell layer using the D parameter. The modeling results are presented in Fig. 9. Fig. 9(a) demonstrates that the variation in I-V characteristics of the NAND string due to the taper angle is well captured. Additionally, Fig. 9(b) confirms that the I-V characteristics are accurately modeled based on the cell position within the string at a specific taper angle. When the tapered structure was considered, the highest fitness score was 0.2307, and the on-current error was within 1%, indicating that the model effectively captures the impact of tapering

Fig. 5. I-V Characteristics of NAND Strings with different numbers of WLs. The inset shows only the on current.

Fig. 6. I-V Characteristics of NAND Strings for different bit line voltage ($N = 32$).

Fig. 7. I-V curve for different WL position (a) in linear scale and (b) in log scale. ($N = 32$) (c), (d) Variation of G$_{\rm m}$ due to different effective source resistance with the cell position.

Fig. 8. (a) Tapered structure of 3D-NAND flash memory. (b) Different channel area profile with its location. (c) Implementation of tapered structure in BSIM-CMG.

Fig. 9. (a) I-V curves for different tapered angle. (b) I-V curves for different WL position ($N = 32$).

IV. CONCLUSIONS

In this study, the I-V characteristics of 3D-NAND string were modelled using the BSIM-CMG transistor model combined with a GA-based parameter optimization technique. Reference data for 3-D NAND string with one, three, and five WL layers were obtained using TCAD device simulations, and GA optimization was applied to derive optimal values for a total of 24 parameters. The proposed modelling approach was verified to be applicable to strings with various WL counts, up to 500 layers, and accurately predicted the I-V characteristics for different selected WL positions and voltage conditions. Notably, the method successfully captured the impact of tapered channel structures caused by etching process limitations. The developed modelling technique can be effectively applied to full-chip circuit simulations integrated with peripheral circuits and can further be utilized to predict the characteristics of next-generation 3D-NAND flash memory as the number of layers continues to increase.