Robust Simulations of Nanoscale Phase Change Memory: Dynamics and Retention

A robust simulation framework was developed for nanoscale phase change memory (PCM) cells. Starting from the reaction rate theory, the dynamic nucleation was simulated to capture the evolution of the cluster population. To accommodate the non-uniform critical sizes of nuclei due to the non-isothermal conditions during PCM cell programming, an improved crystallization model was proposed that goes beyond the classical nucleation and growth model. With the above, the incubation period in which the cluster distributions reached their equilibrium was captured beyond the capability of simulations with a steady-state nucleation rate. The implications of the developed simulation method are discussed regarding PCM fast SET programming and retention. This work provides the possibility for further improvement of PCM and integration with CMOS technology.


Introduction
Phase change memory (PCM) is one of the emerging non-volatile memories with the most mature technology and the most industrialized prospects at present. It has the advantages of nanoscale cells, high reliability and endurance, high speed, and compatibility with CMOS technology [1]. In particular, integrations at the back-end-of-line of advanced CMOS technology brings exciting new features. The write speed of PCM can possibly be raised further to the sub-nanosecond level through material engineering [2,3] or innovative programming methods [4], comparable to other CMOS-based memories such as static random access memory (SRAM). One of the most widely studied phase change materials thus far is Ge 2 Sb 2 Te 5 (GST). One current concern for PCM is its reliability, e.g., the high-resistance amorphous state shifts to low-resistance states with unintended crystallization [5][6][7]. Nucleation at random sites and subsequent growth leading to the formation of percolation paths might be responsible for retention failure. Further theoretical and experimental studies are useful for maturity of the PCM technology, which enriches and enhances the CMOS technology in various applications such as storage class memory (SCM) and neuromorphic computing.
From a theoretical perspective, nanoscale modeling and simulations play important roles in PCM technology advancements. Empirical formulas and Johnson-Mehl-Avrami (JMA) models [8] are widely used in compact models [9][10][11]. At the numerical simulation level, a classical nucleation and growth (NG) model proposed by Peng et al. [12] is commonly used to reflect the random nucleation and non-uniform distribution of phase states of materials. This model successfully takes the randomness of nucleation/growth and the distribution of phase states into consideration. Nucleation occurs in a small volume if the nucleus formation probability is significant. Growth occurs when there is a crystalline neighbor around the considered volume and the growth probability is high. Together with the NG model, simulations of PCM have been reported [13][14][15][16][17][18][19] by the Monte Carlo method based on geometry discretization from the finite element method. In [13], authors reported the simulation framework as well as the programming operations and the thermal boundary resistance effect. Another numerical simulation algorithm [14] was reported as a basis to propose a phase change module for multi-value storage capability. Further, the thermoelectric effects have been incorporated [15] in the coupling of thermal and electrical simulations. At the same time, simulations are also proved useful in the study of PCM variations and retention failures. By considering the stochastic nucleation, the intrinsic retention as well as its statistics, e.g., the cell-to-cell and cycle-to-cycle variability have been studied [16,17]. The phase field method [18,19] was later included in the simulation framework to address the crystallization from the perspective of total energy (bulk free-energy and interface energy) reduction. It has been shown that the retention time is improved for scaled PCM. However, one basic assumption behind the traditional simulations is a steady-state nucleation rate at the initial stage. The incubation period, i.e., the period in which the nuclei grow to their critical size, is not considered. On the other hand, a critical nuclei size is usually assumed upon nucleation for a discretized volume in the NG model. As such, the preferred volume size is equal to the critical size. Since the critical nuclei size is dependent on temperature [20], the NG model is more suitable for phase changes under a uniform temperature profile [20]. Alternative simulation frameworks are desired for nanoscale PCM.
In this work a simulation framework was developed for the evolution process of the nanoscale PCM embryo distribution from the perspective of dynamic nucleation. The classical nucleation growth model was improved to be suitable for the phase transition under non-uniform temperature conditions. Additionally, the relationship between the transient nucleation rate and the incubation period of the crystal nucleus was analyzed. The accelerating SET operation, as well as the device retention failure, was also discussed.

The Theory and Simulation Algorithm
PCM of the classical mushroom type, as shown in Figure 1a, was used to demonstrate the simulation theory and algorithm. crystalline neighbor around the considered volume and the growth probability is high. Together with the NG model, simulations of PCM have been reported [13][14][15][16][17][18][19] by the Monte Carlo method based on geometry discretization from the finite element method. In [13], authors reported the simulation framework as well as the programming operations and the thermal boundary resistance effect. Another numerical simulation algorithm [14] was reported as a basis to propose a phase change module for multi-value storage capability. Further, the thermoelectric effects have been incorporated [15] in the coupling of thermal and electrical simulations. At the same time, simulations are also proved useful in the study of PCM variations and retention failures. By considering the stochastic nucleation, the intrinsic retention as well as its statistics, e.g., the cell-to-cell and cycle-tocycle variability have been studied [16,17]. The phase field method [18,19] was later included in the simulation framework to address the crystallization from the perspective of total energy (bulk free-energy and interface energy) reduction. It has been shown that the retention time is improved for scaled PCM. However, one basic assumption behind the traditional simulations is a steady-state nucleation rate at the initial stage. The incubation period, i.e., the period in which the nuclei grow to their critical size, is not considered. On the other hand, a critical nuclei size is usually assumed upon nucleation for a discretized volume in the NG model. As such, the preferred volume size is equal to the critical size. Since the critical nuclei size is dependent on temperature [20], the NG model is more suitable for phase changes under a uniform temperature profile [20]. Alternative simulation frameworks are desired for nanoscale PCM.
In this work a simulation framework was developed for the evolution process of the nanoscale PCM embryo distribution from the perspective of dynamic nucleation. The classical nucleation growth model was improved to be suitable for the phase transition under non-uniform temperature conditions. Additionally, the relationship between the transient nucleation rate and the incubation period of the crystal nucleus was analyzed. The accelerating SET operation, as well as the device retention failure, was also discussed.

The Theory and Simulation Algorithm
PCM of the classical mushroom type, as shown in Figure 1a, was used to demonstrate the simulation theory and algorithm.

Electrothermal Module
Joule heat is generated by applying electric pulses of different durations and amplitudes, which is used to realize the transformation of different phases of GST, which is the basic principle of PCM. The coupling of the current continuity equation and heat conduction equation is the key to building an electrothermal module: where J is the current density, σ is the electrical conductivity, ρ is the density of GST, c p is the heat capacity at constant pressure, T is the temperature, and κ is the thermal conductivity.
To reproduce the correct thermal boundary, a large enough domain including the PCM cell and its surrounding dielectric is always simulated. The Dirichlet boundary condition is applied on the electrodes with its temperature of 293.15 K, and the Neumann boundary condition is applied on other surfaces. Thermal boundary resistance (TBR) at the interface between different materials, for example, the interface of GST/BE, is easily incorporated into the simulation. It facilitates a confinement of heating and hence the reduction of the reset current, which may be intentionally introduced in the PCM cell design [21,22]. Figure 1b briefly shows the discretization with cuboidal grids. The electrical conductivity of a mixed phase σ m is calculated with the Wiener upper bound model [23]: where σ c and σ a are the electrical conductivity of the crystalline and amorphous states, respectively, and C f is the crystalline fraction of GST. The thermal conductivity of GST is divided into phonon thermal conduction and electronic thermal conduction. The electrical conductivity of the material is proportional to the thermal conductivity, which can be described by the Wiedemann Franz Law. The phonon thermal conductivity of the mixed state κ phm can be expressed as: where κ phc and κ pha are the phonon thermal conductivity of the crystalline and amorphous states, respectively, and the thermal conductivity of the mixed phase κ m is where L 0 is the Lorenz number.

Dynamic Nucleation
The first step of crystallization requires the formation of new phase nuclei in the parent phase, similar to the classical nucleation and growth theory. The driving force for nucleation is that the volume free energy of the crystalline state is lower than the volume free energy of the amorphous state, and the resistance is due to the formation of a new surface by the crystal embryo, which will cause the surface energy γ to increase. The crystal embryo is approximately a sphere to reduce interface area; it will not become a stable crystal nucleus until the radius r of the crystal embryo exceeds the critical nucleation radius r c . After that, the nucleus can continue to grow. The expression for the number of molecules n c contained in the critical nucleation radius is where V a is the average molecular volume of GST, and ∆G V is the free energy change per unit volume of the liquid-to-solid phase transition, which can be expressed as [24] where ∆H m is the latent heat, and T m is the melting temperature of the GST. A crystal embryo containing n molecules appears in the amorphous phase, and the total free energy difference ∆G n is where r is radius of the embryo and can be expressed as The crystal embryo in the amorphous phase may increase its size with more molecules, or it may decrease its size with less molecules due to dissolution. The expansion and dissolution processes occur in parallel and are described by the following two expressions [14]: where A 1 , A n−1 , A n , and A n+1 are embryos containing 1, n − 1, n, and n + 1 molecules, respectively, and P + n is the probability that the embryo A n will increase by one molecule. P − n is the probability that the embryo A n will decrease by one molecule. The transient expression N n,t of the number of embryos containing n molecules in the system at time t is, then [25], where O n is the number of molecules on the surface of the embryo A n , which is approximately O n ≈ 4n 2/3 ; D is the diffusion coefficient; λ is the atomic jump distance; and ∆g n is the amount of energy change in the system after the embryo A n gains or loses one molecule. The nucleation rate at any time can be obtained by solving the grain flux: With the passage of annealing time, the number of embryos will gradually increase, and the size will continue to increase, and the nucleation rate will continue to rise, but at the same time, more embryos will dissolve, and the number of embryos will develop towards a steady state, while the nucleation rate reaches the steady-state value I s nc . After nucleation, under the driving force of crystallization, the crystal nucleus will grow up. Since the mobility of molecules is affected by temperature, its growth rate is affected by temperature. Kelton et al. developed the GST grain growth rate v g suitable for a wide temperature range [26]:

The Algorithm and Improved Crystallization Model
On the basis of the above model framework, a new algorithm was developed, as shown in Figure 1b. According to the principle of equal distribution of the sample and the population, a dual grid technique was proposed to solve the dynamic nucleation based on the Monte Carlo method. Unlike the classical nucleation growth model, which regards the grid as the element of the crystal nucleus, this model only regards the grid as the carrier of the crystal nucleus, and the change of the grid crystalline fraction is obtained by solving the ratio of the crystallized volume change to the grid volume. The steps are as follows: First, the GST region is divided into many small grids with a volume of V grid (V grid = d 3 grid ). It is assumed that the temperature in the grid is uniformly distributed and considers the grid as a part of a large volume of V bulk whose temperature is uniform. n grid (n grid = V bulk /V grid ) is the number of samples in the population, and the number of nuclei in the population in a time step ∆t can be calculated as N c bulk (N c bulk =I s nc ·∆t·V bulk ); then, the probability P n of crystal nuclei appearing in the small grid is The probability of the appearance of a crystal nucleus is combined with the Monte Carlo algorithm to simulate the nucleation process. The flow chart is shown in Figure 1c. If a crystal nucleus appears in the grid, the size of the crystal nucleus is recorded as the critical nucleation. The crystalline fraction of the updated mesh due to nucleation is If the volume of the crystal nucleus is smaller than the volume of the grid, record the size of the crystal nucleus; if the volume of the crystal nucleus exceeds the volume of the grid, the surrounding grid will evenly divide the grid and the crystalline fraction of the grid needs to be updated. The crystalline fraction of the updated mesh due to growth is where r t−∆t is the radius of the crystal nucleus at time t − ∆t, and N c is the number of nearest surrounding grids with C f = 1. Figure 2 shows PCM R-V characteristics with an amorphous initial state. Programming voltage pulses with different amplitudes and fixed width t p = 170 ns were applied to the PCM cell, in accordance with [15]. The resistance value was obtained by a reading following each programming pulse. With the increasing of voltage amplitudes, the PCM cell went through the SET to low resistance and then the RESET to high resistance. The simulation results agreed reasonably with the experimental data. In particular, the partial RESET of varying degrees was well reproduced. Figure 3 shows the simulated crystallization period as a function of the annealing temperature from the improved crystallization model and the classical nucleation growth model. The results are indeed the average values of multiple measurements. Different grid sizes were set up for the purpose of robustness verification. The purple dotted line represented results of the classical nucleation growth model under the critical nucleation size (which changes according to the temperature). Two other grid sizes, 1 nm and 2 nm, were also used for comparison. With the classical NG model, it could be seen that at the same temperature, the selection of the grid size had a greater impact on the crystallization period, because it is unreasonable to set a uniform grid size for different temperatures corresponding to different critical nucleation sizes. In contrast, the improved crystallization model gave more stable results without significant impacts from grid size, and the results were also close to those of the NG model with a temperature-dependent grid size.

The Incubation Period
Assuming that there are no embryos in the initial state, the function relationship between the number of nucleus J n and time t is shown in Figure 4a. It can be seen that as time goes by, the number of crystal nuclei gradually increases from zero. It gradually stabilizes, indicating that the steady-state nucleus rate is reached. The θ τ in Figure 4a is the nucleus incubation period, and Figure 4b shows the nucleus incubation period at different temperatures, assuming that there is no embryo in the initial state. With the increase in temperature, the incubation period first decreases because of the increase in the mobility of molecules and the probability of their collision. The incubation period then increases slightly, most probably because of the decrease in the degree of subcooling and the increase in the critical nucleation radius, which requires a longer time to reach the critical nucleation molecule number.  Figure 5 shows that dynamic nucleation in the steady-state and the traditional steadystate nucleation rate are in good agreement with a slight deviation at the peak value, and the optimal nucleation temperature corresponding to the peak slightly changes. Nucleation rates in Figure 5 correspond to the slope of linear coordinates in Figure 4a.

Implication on PCM Cell Operations
Based on the above simulation framework, PCM operations were simulated in this section. The default geometry parameters used are summarized in Table 1.  Figure 6 shows the crystallization conditions after annealing at 580 K and 660 K for 30 ns, assuming that there is no embryo in the initial state. In Figure 6a, a large number of small grains were formed, and in Figure 6b, a smaller number of grains was observed with larger sizes. The incubation periods were 5.6 ns and 2.9 ns under 580 K and 600 K, respectively. In the simulation period of 30 ns, the total number of nucleus depended on the slope in Figure 4a, which was higher under 580 K. The growth rate as given by Equation (6) was enhanced under 660 K for a large crystalline fraction.

Symbol
Description Value unit WGST width of GST region 190 nm LGST length of GST region 190 nm HGST height of GST region 90 nm Figure 6 shows the crystallization conditions after annealing at 580 K and 660 K for 30 ns, assuming that there is no embryo in the initial state. In Figure 6a, a large number of small grains were formed, and in Figure 6b, a smaller number of grains was observed with larger sizes. The incubation periods were 5.6 ns and 2.9 ns under 580 K and 600 K, respectively. In the simulation period of 30 ns, the total number of nucleus depended on the slope in Figure 4a, which was higher under 580 K. The growth rate as given by Equation 6 was enhanced under 660 K for a large crystalline fraction.

Cf
LGST Figure 6. Crystallization after annealing at different temperatures, namely (a) 580 K and (b) 660 K, for 30 ns, assuming that there is no embryo in the initial state. Figure 7 shows the application of three different SET pulse schemes and the phase distribution corresponding to points A-F. Applying a low-amplitude pulse will rapidly nucleate the inside of the active area, but due to the low temperature inside the active area, the growth rate of the formed nuclei slows down and the active area cannot be crystallized quickly. The case is shown as the top phase maps A and B in Figure 7b. Applying a higher amplitude pulse, although the growth rate of the crystal nucleus is faster at this temperature, the active area cannot be crystallized quickly from inside due to a lower nucleation rate. The case is shown as the middle phase maps C and D in Figure 7b. Applying a dual-amplitude pulse, a lower-amplitude pulse promotes the nucleation process, and then a higher-amplitude conventional SET pulse promotes the rapid growth of the crystal nucleus and crystallizes almost the entire active region. The case is shown as the bottom phase maps E and F in Figure 7b. The instantaneous temperature profile for A of lowamplitude pulse case and D of high-amplitude pulse case are shown in Figure7c,d.  Figure 7 shows the application of three different SET pulse schemes and the phase distribution corresponding to points A-F. Applying a low-amplitude pulse will rapidly nucleate the inside of the active area, but due to the low temperature inside the active area, the growth rate of the formed nuclei slows down and the active area cannot be crystallized quickly. The case is shown as the top phase maps A and B in Figure 7b. Applying a higher amplitude pulse, although the growth rate of the crystal nucleus is faster at this temperature, the active area cannot be crystallized quickly from inside due to a lower nucleation rate. The case is shown as the middle phase maps C and D in Figure 7b. Applying a dualamplitude pulse, a lower-amplitude pulse promotes the nucleation process, and then a higher-amplitude conventional SET pulse promotes the rapid growth of the crystal nucleus and crystallizes almost the entire active region. The case is shown as the bottom phase maps E and F in Figure 7b. The instantaneous temperature profile for A of low-amplitude pulse case and D of high-amplitude pulse case are shown in Figure 7c,d. Figure 8 plots the annealing of a PCM cell, assuming different amounts of crystal embryos for the initial states. In the simulation, the cell was annealed under a constant temperature of 450 K. Figure 8a shows the phase distribution at the initial moment with a defined active region. Figure 8b is the phase distribution after annealing for 9000 s, without considering the embryos generated in the previous quenching process. The retention process was mainly due to the growth from the outer to inward. This was because when the embryos produced during cooling were not considered, the embryos' incubation periods were longer, and the transient nucleation rate was lower. Figure 8c is the phase distribution after annealing for 6000 s. The failure process of the device was due to both the growth of the outer edge crystal and the nucleation of the amorphous region. The retention failure is accelerated compared with Figure 8b. Figure 8d is the phase distribution after annealing for 3000 s when an even larger number of embryos was considered. The nucleation of the amorphous region was dominant, leading to a percolation path with early failure. In other words, the initial nucleation rate was close to a steady-state rate, and the device would fail prematurely due to the large number of nuclei in the active region.   Figure 8 plots the annealing of a PCM cell, assuming different amounts of crystal embryos for the initial states. In the simulation, the cell was annealed under a constant temperature of 450 K. Figure 8a shows the phase distribution at the initial moment with a defined active region. Figure 8b is the phase distribution after annealing for 9000 s, without considering the embryos generated in the previous quenching process. The retention process was mainly due to the growth from the outer to inward. This was because when the embryos produced during cooling were not considered, the embryos' incubation periods were longer, and the transient nucleation rate was lower. Figure 8c is the phase distribution after annealing for 6000 s. The failure process of the device was due to both the growth of the outer edge crystal and the nucleation of the amorphous region. The retention failure is accelerated compared with Figure 8b. Figure 8d is the phase distribution after annealing for 3000 s when an even larger number of embryos was considered. The nucleation of the amorphous region was dominant, leading to a percolation path with early failure. In other words, the initial nucleation rate was close to a steady-state rate, and the device would fail prematurely due to the large number of nuclei in the active region.   Figure 8 plots the annealing of a PCM cell, assuming different amounts embryos for the initial states. In the simulation, the cell was annealed under a temperature of 450 K. Figure 8a shows the phase distribution at the initial mom defined active region. Figure 8b is the phase distribution after annealing for 900 out considering the embryos generated in the previous quenching process. The process was mainly due to the growth from the outer to inward. This was beca the embryos produced during cooling were not considered, the embryos' incub riods were longer, and the transient nucleation rate was lower. Figure 8c is the p tribution after annealing for 6000 s. The failure process of the device was due to growth of the outer edge crystal and the nucleation of the amorphous region. T tion failure is accelerated compared with Figure 8b. Figure 8d is the phase di after annealing for 3000 s when an even larger number of embryos was consid nucleation of the amorphous region was dominant, leading to a percolation early failure. In other words, the initial nucleation rate was close to a steady-state the device would fail prematurely due to the large number of nuclei in the activ  The simulation implementation and efficiency were summarized briefly. MATLAB and COMSOL were used jointly to implement the simulation flow of Figure 1c. While MATLAB was used for meshing and dynamic nucleation calculation with the in-house code, COMSOL was used for the solutions of electrothermal formulations. On a workbench with a CPU (intel(R) Core (TM) i7-4790 @ 3.60 GHz), it took about 7 h to obtain the data point shown in Figure 2, and about 0.5 h to finish the annealing under a constant temperature, shown in Figure 6.

Conclusions
A simulation framework was developed to study the operation of PCM cell operations. The dynamic nucleation process was implemented with the reaction rate theory. An improved crystallization model was also developed with robustness regarding the choice of grid size. Together with the electrothermal module, the PCM simulation framework was used to guide the PCM operations. The incubation period was captured with the developed method, including descriptions of the embryo distributions. A dual-amplitude pulse scheme was then analyzed to accelerate the SET process, i.e., using a lower temperature to induce enough nucleus, and then using a higher temperature to enhance the growth rate. The different failure mechanisms of the device caused by different embryo distributions were analyzed, which may explain the early retention failure driven by nucleation.