Impact of plasma induced damage on the fabrication of 3D NAND flash memory

A physical process model for dry plasma etching is presented and applied to simulate vertical channel hole etching, a critical fabrication step in modern three-dimensional (3D) NAND flash memory. The presence of physical etching with high energy ions is shown to induce damage in the underlying silicon, which results in the formation of voids during the subsequent selective epitaxial growth (SEG) step. In this manuscript, we present a model for ion induced damage by storing it as a surface property during the plasma etching simulation. A specialized advection algorithm is subsequently applied to simulate silicon SEG on the bottom source line. The model clearly shows the damage caused by the high energy particles, on the crystal nature of silicon, resulting in poor coverage during the SEG step. The removal of this damaged layer using lower energy plasmas results in highly crystalline epitaxially grown silicon. The simulation results show excellent agreement with experiments in the formation of undesired voids without the low-energy pre-epitaxial plasma treatment.


Introduction
Modern three-dimensional (3D) flash memory (NAND) [1,2] employs a large number of stacked control gates and insulating layerscomprised of alternating metallic (usually tungsten) and silicon dioxide (SiO 2 ) thin films, leading to a recent increase in storage size up to 1.33 TB [3] on a single die [4]. The fabrication of these stacks relies on using silicon nitride (Si 3 N 4 ) as a sacrificial layer, which is etched away and replaced by tungsten in a later processing step [5]. Therefore, when the channel is being etched, the etching proceeds through a stack of alternating SiO 2 and Si 3 N 4 thin layers. Until now, etching through 128 layers was possible, while anything beyond that becomes a problem for presently-available plasma etching techniques [6]. This rapid scaling drives new challenges in the etching of these stacks in order to meet the demands of further increasing memory capacity by introducing even more stacked layers, thereby dramatically increasing aspect-ratios of the feature sizes. The high aspect ratios (HARs) and small feature sizes result in the increase in the plasma densities, which also makes the process more prone to fabrication-induced damage. Due to the different framework requirements for modeling topography motion and plasma induced damage (PID) during a sequence of fabrication steps, a combined study of these two effects was not previously available. In this work we present, for the first time, a process simulation methodology which is able to store PID effects after a topography plasma etching simulation, which is then applied in subsequent simulation steps, such as selective epitaxy. The developed method shows qualitative agreements with experimental values.

Plasma Induced Damage
Fabrication using plasma processes, such as ion-enhanced etching, play a significant role in today's advanced technologies. They have been developed to fabricate very fine patterns with anisotropic features in metal oxide semiconductor field effect transistors (MOSFETs). As the transistors and features scaled down over time, the plasma density required increasing. While advancements have been made to improve plasma processing, there is nevertheless material degradation associated with it, usually referred to as plasma-induced damage (PID) [7,8].
PID has been widely studied in silicon-based electronic devices [9]. The degradation induced by plasmas is generally categorized into three types, namely, physical, charging, and radiation damage [10]. In MOSFETs, the physical damage is induced by high-energy ion bombardment on silicon substrates or other surfaces [11]. Charging damage is induced by the current from a plasma flowing into the insulator or dielectric material [12], while radiation damage corresponds to the breaking of bonds in a material due to high-energy photon interactions [13]. The most relevant source of damage for the devices of interest in this manuscript, namely 3D NAND, is physical PID, which is visually represented in Fig. 1. From the figure, it can be observed that one of the damaging aspects of the process is that the crystalline nature of silicon is disrupted near the region of ion impact.

Physical PID in 3D NAND
The specific challenge in NAND fabrication is the dry etching of highaspect-ratio vertical channel holes through the stacked SiO 2 /Si 3 N 4 layers using high energy plasmas. As previously mentioned, these types of plasma etch chemistries introduce damage in the silicon layers by implanting ions in the crystalline silicon substrate which forms the source line [11]. The physical PID leads to the introduction of a damaged layer, several nm in thickness, with a severe disorder in the crystal structure, as shown in Fig. 2, resulting in poly-crystalline selective epitaxial silicon growth (SEG). This leads to undesirable electrical properties in the source line due to the poor crystal quality, as well as the inclusion of voids in the grown silicon. Recent studies have shown that removing the damaged silicon layer using a post etch plasma treatment can effectively remove these impurities, leading to highly crystalline SEG of silicon at the bottom [14,15].
We propose a physical model which is capable of simulating the channel hole dry etch process, including the resulting source material damage through ion implantation in a high energy plasma, leading to voids in the source contact after SEG. We show how storing the damage information on surface points and including an exponential attenuation parameter during propagation of the surface, leads to results which show excellent qualitative agreement with experimental studies. The model uses an in-house developed level set powered topography simulator ViennaLS [16] in combination with top-down Monte Carlo (MC) ray tracing provided by ViennaRay [16] to accurately describe the surface kinetics during the etching process. For the simulation of SEG, a specialized advection algorithm developed by Toifl et al. [17] is applied. The simulation tools account for multiple material regions in 3D geometries and are combined and implemented in the ViennaTools software ecosystem [16].
This manuscript is organized as follows. In Section 2 the physical process model is introduced, starting with a general explanation of level set-based surface descriptions and how the evolution of these surfaces is modeled. Subsequently, flux calculations using Monte Carlo ray tracing are described as applied for the modeling of surface kinetics and the ioninduced damage of the substrate. Finally, the model and specialized advection algorithm for the description of SEG is discussed. Section 3 shortly introduces the plasma chemistry which determines our various model parameters, after which the simulation setup is described in Section 4. The obtained results are presented in Section 5 and a conclusion is provided in Section 6.

The Level Set Method
To accurately represent the substrate surface and its time evolution during the etching process, a level set based description is used [18]. With this method, the surface is described implicitly by a level set function ϕ( x → ) which is defined at every point x → in space. This function is obtained using signed distance transforms, describing the surface S as the zero level set: In order to propagate the surface during an advection step, the level set equation is solved in time t, given the scalar velocity field v( x → ) which describes the surface normal velocity at each point. This is achieved by discretizing the level set function on a regular grid and applying a finite difference scheme to solve Eq. (2). The velocity field is generated by extending the calculated surface rates to all points in the simulation space. Theses surface rates are calculated using surface kinetics models, described in the next section.

Surface Kinetics
To analyze the surface kinetics which describe the plasma etch process, a model based on the theory of active surface sites is applied [19]. The surface propagation is a result of concurrent etching and deposition phenomena. For this reason, three different types of particle species must be considered:   3. Energetic ions which physically sputter the film and enhance the dissociation of the volatile etch products. It should be noted that these high energy ions are the primary source of PID in silicon.
The rates of particle types impinging on the surface can be summed to give the coverages θ x of different particle types at all surface sites, where x represents the etchant (e), polymer (p), or etchant on polymer (pe). Since the etching time is usually much longer than the surface reaction time scales, we can assume that the coverages reach a steady state on the surface. This means that the amount of the different species found at all times can be expressed in terms of constant fluxes J x . The coverages can therefore be expressed with the following site balance equations [19]: • Polymer sites balance • Etchant on polymer surface balance Here, J x and S x represent the different particle fluxes and sticking probabilities, respectively. Y ie is the ion-enhanced etching yield for etchant particles, Y p is the ion-enhanced etching yield on polymer, Y sp gives the physical ion sputtering yield, and k ie and k ev are the stoichiometric factors for ion-enhanced etching and evaporation, respectively, which are determined by the chemical etching reaction. By solving these steady state equations for the coverages, one can determine etch or deposition rates on the surface. If deposition of polymer dominates, the surface normal velocity is positive and is given by where ρ p is the atomic polymer density. The first term J p S p gives the rate of polymer particles reaching and adsorbing on the surface, while the second term Y p J i θ pe describes the removal of polymer by ion-enhanced etching. Together, these terms describe the deposition of polymer material on the surface, which acts as passivation layer for the chemical etching process. If, on the other hand, etching of the substrate dominates, the negative surface velocity of the substrate is given by where ρ m is the atomic density of the etched material and depends on which layer in the stack is being etched. In Eq. (7) each term accounts for a different type of surface reaction. The first term, J ev θ e , describes the chemical etching process, where etchants bind chemically with the substrate to form volatile etch products which dissolve thermally from the surface. Thus, the evaporation flux J ev is a parameter proportional to the etchant flux J e and depends on the chemical gas and surface composition and temperature of the etching plasma. It is given by where K is a process parameter describing the volatility of the chemical etching process, E a is the activation energy for thermal etching, k B is the Boltzmann constant, and T is the temperature. The second term in Eq. (7), J i Y ie θ e , describes the contribution of ion-enhanced etching. In this surface reaction, volatile etch products which do not dissolve from the surface thermally, absorb energy from impinging ions and consequently dissolve from the surface. Finally, the last term, J i Y sp (1 − θ e ), describes physical sputtering of the substrate by highly energetic ions. Since both chemical and ion-enhanced etching involve etchants, they are proportional to the etchant coverage θ e , while physical ion sputtering takes place directly on the substrate and is thus proportional to the fraction of the surface not covered by the etchant. These different types of processes, leading to a deformation of the surface, are illustrated in Fig. 3. If polymer forms on the top layer which is being etched, only the removal by ion-enhanced etching is considered and the negative surface velocity is again given by Eq. (6), since ion-enhanced etching proceeds on all surfaces, including those covered by the polymer. Together, Eq. (6) and Eq. (7) describe the temporal evolution of the surface, given the particle fluxes J x at each location on the etched substrate, which are calculated using MC ray tracing methods.

Flux Calculation
In order to calculate the particle fluxes, the particle transport moving towards the sample surface inside the etching reactor must be simulated. Since the particle's traversal inside the reactor depends on the specific geometry of the reactor and the processing conditions prevailing inside it, the reactor is divided into reactor-scale and feature-scale regions during simulations, separated by the so-called source plane [20]. The feature-scale region describes the space directly above the wafer, being comparable in size to the geometric features generated by processing the wafer.
The reactor-scale region, on the other hand, is much larger than the feature-scale region and particle transport inside the reactor-scale region is governed by Maxwell-Boltzmann statistics. This is the case, since the reactor-scale region is much larger in comparison to the mean free path of the particles, which is in the range of 100 μm to 10 mm. The feature scale region, however, is small compared to the mean free path, on the order of a few hundred nm, meaning that collisions with the surface are much more frequent than those with other particles in the gas phase. Therefore, ballistic transport can be used to describe the propagation of particles through the feature scale region, which can be simulated in a straightforward manner using ray tracing methods. In particular, a top-down Monte Carlo ray tracing approach is used, where a large number of particles are launched from a source plane and are traced towards the substrate surface in order to determine the point of impact. Each particle will usually represent several hundred molecules or ions of a particular species which correspond to the particle species described in the previous section, i.e., etchant, polymer, or energetic ions. Since the particles in the feature-scale region are independent of each other, the surface intersection calculations can be parallelized straight-forwardly, which renders this method computationally efficient [21].
For the calculation of the precise location of the surface impact, the implicit surface, given by the level set representation, is converted to an intermediate explicit representation by placing a tangential disc at each surface point [22]. Each disc is angled so that its normal is equal to the surface normal and its radius d r is determined by the grid spacing d g . To ensure, the discs form a closed surface with no holes, the disc radius is given by However, this implementation leads to partially overlapping disks and thus intersections can occur on multiple disks. A schematic example of a surface represented by discs, including intersections, is given in Fig. 4. Particles are first initialized at random positions on the source plane with random directions according to particle-specific distributions. Both polymer and etchant are neutral particles and thus are not affected by the electric sheath potential inside the plasma reactor. Therefore, their initial directions can be approximated by a simple cosine distribution. Ions, on the other hand, are accelerated towards the surface by the sheath potential. Therefore, their initial distribution is represented with a power cosine distribution, meaning that their movement is more directional towards the surface [23].
The ion yield efficiencies on the point of impact are dependent on the ion energy E, as well as the incident angle α between the ion particle and the surface normal. The initial energy for each ion particle is assigned based on the process used. In general, the ion yield efficiency for physical sputtering and ion-enhanced etching upon impact on the surface can be expressed as where A is a process-dependent constant, describing the particle yield per unit of energy, E th is the minimum energy ions must attain to etch the substrate, referred to as threshold energy, and f(α) is a function of the incident angle. For physical ion sputtering, this function can be approximated using the expression while for ion-enhanced etching, the expression is applied, as described in [19]. If the energy of an ion is below the threshold, it is reflected specularly, until it reaches a point with a larger incident angle or until it leaves the simulation domain. The ion energy is however adjusted at each surface reflection according to the angle-dependent yield efficiencies.
Polymer and etchant particles, on the other hand, reflect diffusely from the surface. A weight factor which describes the probability of the particle taking the modeled trajectory is reduced after each reflection with particle dependent sticking probabilities. The particles are tracked until the weight factor drops below a certain limit or the particle leaves the simulation domain.

Plasma Impact Damage
In HAR etching processes, the use of high energy ions in a plasma is necessary to achieve the desired channel aspect ratios. However, this introduces the problem of high energy ions being able to penetrate into the bottom substrate, consequently leading to a disordered silicon layer on the bottom of the channel hole, thereby destroying the crystal purity. Experimental studies of ion induced etching damage show that the damaged layer is composed of impurities related to the etchant chemistry [15,24]. To model this damaged layer in the crystalline silicon substrate, the energy of traced ions is recorded on the substrate surface at each flux calculation step. Since the highest ion energies occur at normal incidence, the recorded energies can be applied to model the ion damage in the material directly below the surface. At each surface point, the impinging ion energy E is assumed to decrease through the substrate due to scattering off the silicon atoms, following an exponential attenuation given by where E i is the initial ion energy upon surface impact and d is the normal distance to the surface inside the material, equivalent to the penetration depth. Given the observed thickness of the damaged layer d th = 8 nm, as noted in Fig. 2 [14], and the threshold energy for ion-enhanced etching E th , the attenuation length λ is determined as Subsequently, an ion damage coefficient D(d) is defined and stored for Fig. 4. Schematic representation of a 2D surface using disks. If a particle trajectory intersects any tangential disk, it contributes to the rates of the corresponding grid point. Reflected particles are then traced, starting from the first disk intersection point. Here, reflection of a particle is illustrated with an incidence angle α to the surface normal n → . We also note that the incident particle ray will intersect with two disks.
each surface point of the geometry. The coefficient is proportional to the ion energy at a depth d in the substrate: In order for silicon to grow epitaxially, the material must not be damaged and hence, the ion damage coefficient at the surface must fulfill D⩽0 for growth to take place, which equivalently describes E(d)⩽E th . To remove the damaged layer in a post etch treatment process, the surface is etched using low energy ions [14]. Assuming these low energy ions do not cause any additional implantations in the substrate, the surface is etched until the damage coefficient D at the bottom drops below 0, thereby forming a suitable highly crystalline interface for the subsequent selective epitaxial growth. The surface damage coefficient D is stored as a local parameter at each surface point on the surface level set. During the surface propagation in the post etch treatment, the damage coefficient is reduced according to the exponential attenuation given in Eq. (15), where d is the advected distance in the direction of the surface normal. In order to accurately transport the damage coefficient to the new surface after each advection step, we use the Manhattan normalization. First, we assume that along a surface normal all local surface properties remain unaltered. From the level set equation, we know that information only propagates along the surface normal and thus every new grid point acquires the properties from the grid point from which it arises. Similar to level set advection, this method is very accurate if the surface propagates along its grid lines, but less accurate if it moves diagonal to its grid lines.

Selective Epitaxial Growth
For the simulation of SEG, the specialized numerical method presented by Toifl et al. [17] is applied. In this approach the numerical integration of the level set equation is enhanced by introducing numerical dissipation coefficients based on the Stencil Lax-Friedrichs flux [25], which ensure a monotonic and thus stable evolution of the surface. The rates for the velocity function, which describe the surface evolution are extracted from experimental studies, as described in Section 3.2. The scheme does not rely on numerical calibration parameters and can therefore be employed for generic velocity functions which are purely surface normal-dependent. Additionally, the level set advection is carried out on an additional deposition top layer to robustly handle multiple material regions.
In the SEG step, the algorithm is applied at all surface points describing the bottom Si layer, where the surface damage coefficient fulfills D⩽0. This way it is ensured that SEG commences only on a crystalline Si surface.

Plasma Etching
Etching is modeled in a fluorocarbon based etch process, including Ar and O 2 as inhibitor, as is commonly used for the etching of SiO 2 /Si 3 N 4 layers [26]. As described in Section 2.2, three different physical processes lead to the etching of a substrate: Chemical etching, ion-enhanced etching, and physical sputtering. Fluoride reacting chemically with the silicon surfaces leads to etch products which evaporate back in the gas phase [27]. The fluorocarbon radicals CF + x together with Ar + provide ion bombardment which enables ion-enhanced etching and physical sputtering, where etch products evaporate from the surface due to their now smaller binding energies, or as a result of receiving momentum from the collision cascade induced by the incident ion [28]. The threshold energy E th for physical sputtering is related to the binding energy of the substrate. The addition of O 2 enables the formation of a passivation layer on the sidewall of the substrate, which inhibits lateral etching. The parameters describing these processes in the etching model are extracted from [19,29] and summarized in Table 1. For simplicity, the ion energies are randomly initialized based on a Gaussian distribution around an initial mean energy E ion corresponding to the sheath potential, approximating the results for inductively coupled plasma chemistries from [30].
For the post etch treatment, Luo et al. proposed a chemistry based on an inductively coupled plasma using low energy Cl 2 and NF 3 /CH 2 F 2 gases [14]. By lowering the plasma energy, chemical etching dominates the etching process and ions do not penetrate deep into the surface as is the case with the first etch step. The basic etch mechanisms, however, are the same as in fluorocarbon plasmas. Thus, we can reuse the surface kinetic model with adjusted parameters. In order to account for Cl 2 plasma, the sticking probability for etchant particles is set to the corresponding values for Cl 2 . The initial ion energy is lowered, while the initial particle fluxes for ions and neutral particles are kept the same. The parameters which are different to the fluorocarbon plasma are summarized in Table 2.

Selective Epitaxial Growth
During selective epitaxial growth, silicon is grown on a clean crystalline silicon surface, through the adsorption from Si atoms from the gas phase. On the microscopic scale epitaxial growth is governed by the movement of individual atoms mainly through adsorption, surface diffusion, nucleation, and film coalescence [31]. For a continuum model of epitaxy, the atomistic nature of matter is neglected and microscopic Table 1 Parameters used for the simulation of the vertical channel hole etching process. Values are extracted from [19,29].  Table 2 Parameters used for the simulation of the post etch treatment process, which differ from those in the vertical etching process and parameters to model the ion damage.

Parameter Value Description
Se 0.1 Etchant particle sticking probability Eion 10 eV Mean initial ion energy E th 4 eV Threshold energy for ion damage to occur d th 8 nm Depth of the damaged layer as observed in [14] processes are abstracted with distributed quantities and material properties. Since the incoming atom flux is determined by a single neutral particle species, the flux is assumed to follow Knudsen diffusion in the hole and is thus independent of the lateral position on the wafer. For this reason, the incoming atom flux and incorporation can be represented by growth rates which depend only on the crystal orientation and growth conditions. Growth rates of crystal facets can be experimentally characterized through techniques such as atomic force microscopy (AFM) or cross-sectional transmission electron microscopy (TEM). Therefore, an empirical growth rate function which provides the growth rate of all crystal orientations is constructed. The corresponding rates are extracted from experimental studies by Dutartre et al. [32].

Simulation Setup
The structure of the simulation domain is based on the multimaterial representation presented by Ertl and Selberherr [33]. In this approach, the stacked material layers are represented using individual level sets which are connected by sequential union operations. This multi-material representation enables the resolution of thin layer regions (i.e., thickness below a single grid spacing) and thus an additional layer, placed on top of all layers, is used to capture any deposited polymer.
The geometry of the 60 nm diameter hole is resolved on a regular grid with lateral extensions 64 × 64 and a grid spacing of 1.5 nm. Reflective boundary conditions are assumed in both lateral directions. The etching simulation and post etch treatment process each consist of multiple steps of alternating ray tracing and surface advection. In the ray tracing step, the required fluxes at each surface point are calculated and used to find the resulting surface velocity which is applied for the subsequent level set advection step. The advection step solves the level set equations and moves the surface, effectively identifying the new location of the surface. In the advection algorithm, the time-discretization Δt is adjusted to the previously calculated rates. The algorithm determines the maximum time step possible, according to the Courant-Friedrichs-Lewy (CFL) condition [34] for a sparse level set [22] max|ϕ( x This ensures a stable time integration during level set advection. The number of rays traced is also adjusted during the evolution of the surface. For each grid point a minimum of 3000 rays is simulated to ensure a statistically converged result. This results in anywhere between 3 and 11 million trajectories of each particle species being calculated during each time step. For this reason the main plasma etching simulation was performed on a 40-core Intel Xeon Gold node of a cluster with a computation time of about 2 h. The simulation of the post-etch treatment and the subsequent SEG was carried out on a single workstation with an Intel i7-3770 CPU, requiring approximately 15 min.

Vertical Channel Hole Etching
The geometric profile of the stacked sheets after the etching simulation is shown in Fig. 5. The etch process is performed until the hole is slightly over-etched in the bottom silicon substrate. On the sidewall one can observe the thin passivation layer in red, formed during the etching process.
The damage coefficient on a 3D clip of the surface after the etching process is shown in Fig. 6. Due to the high directional trajectories of ions and the angle-dependent etching yield, the resulting ion damage is predominately confined to the bottom regions of the trench, which is also confirmed experimentally in [14].
An additional effect caused by energetic ions is the charging of dielectric materials, as reported by Memos et al. [35]. The charging of SiO 2 /Si 3 N 4 layers could lead to small non-spectral deflections of the ion trajectories and could affect the etching profile [36,37]. However, most of the damage recorded on the bottom surface is caused by energetic ions which directly hit the surface and are not reflected from the sidewalls. The particles which do reflect off the sidewalls and still have enough energy to damage the crystalline substrate must reflect at very high angles to the surface normal. This results in highly specular reflections with only a minimal loss of energy, as is assumed in our reflection model. Therefore, the surface charging effect on the sidewall is not expected to have a large impact on the ion damage recorded at the  bottom. However, some charging can also take place at the bottom of the hole, ultimately resulting in a reduction in the impacting ion energy and in the surface etch rate, an effect which was excluded in this study. Nevertheless, the implemented model provides for the flexibility of in including any additional effects, specifically because the surface charge can be stored as a parameter within the framework, similar to the way that the damage parameter is stored. For this, only Eq. (10) would need to be adapted to include the effects of surface charging.
During the post-etch treatment, the surface is etched with low energy ions, until the damaged layer has been removed and the crystalline silicon material underneath is revealed, in order for the SEG step to commence on a pure crystal surface. During the surface propagation in the post etch treatment, the damage coefficient is reduced according to the exponential attenuation given in Eq. (15), where d is the advected distance in the direction of the surface normal.

Selective Epitaxial Growth of Silicon
Next, selective epitaxial growth of silicon at the bottom of the channel hole is carried out on the final profile after the post etch treatment. When the damaged layer is present on the silicon substrate a large void is observed, as shown in Fig. 7. Since the disordered layer covers most of the bottom, as observed in Fig. 6, there are only certain points on the sidewall, where SEG can initiate. Such voids lead to an illformed contact to the silicon source line, which will heavily impact the bottom select gate characteristics and reduce the overall quality of the memory stack.
However, when the silicon substrate is cleaned prior to SEG, as described in the previous section, the desired SEG growth is obtained with the grown layer providing full contact with the silicon source line, as can be seen in Fig. 8. The simulation results can be qualitatively compared to the experimental results observed by Lung et al. [15], showing the TEM images of the final structure after SEG with and without a post-etch treatment to clean the silicon substrate, respectively.

Conclusion
A physical process model is applied to simulate vertical channel hole  etching in sequential SiO 2 /Si 3 N 4 3D NAND flash memory stacks. The model, based on a level set surface representation and Monte Carlo ray tracing for flux calculation, is able to simulate the damage induced in a crystalline silicon substrate by ion-enhanced etching. The damage caused on the bottom silicon substrate is modeled using a surface damage coefficient, relating the etch damage to the relative position on the surface using an ion penetration depth. Subsequent silicon SEG is shown to form undesired voids with the damaged layer present, similar to experimental findings [15]. In a pre-epitaxial post etch treatment process, the damaged substrate is removed by low energy ions until the crystalline silicon at the bottom is exposed. The subsequent SEG of silicon provides full contact with the source line, leaving no undesirable voids. This is modeled using the surface damage coefficient which is decreased according to an exponential attenuation, during the etching of the bottom substrate. Due to the physical nature of the presented model, relevant physical effects and underlying etch mechanisms are represented appropriately. Therefore, the proposed model allows to analyze the physical behavior of the etch process in order to optimize the fabrication conditions during channel hole etching and the post etch treatment fabrication steps. Hence, the model serves as a basis for an enhanced understanding of source contact formation in 3D NAND memory cells.
The broad applicability of the developed tools to deal with topography and volume problems provides a sound basis for future studies in damage induced effects during etching processes and other detrimental phenomena taking place during fabrication. The ability to store and process surface information during the simulation of a sequence of fabrication steps allows for the analysis of the impact of damage in various regions of the modeled device.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.