Multi-physics modeling of phase change memory operations in Ge-rich $Ge 2 Sb 2 Te 5$ alloys

Phase change memories (PCM) relying on the electrical contrast between amorphous and crystalline phases of chalcogenide materials have been identified as a promising solution for embedded non-volatile memory technologies.¹ In such devices, the state of individual memory cells is controlled using short but intense pulses of electric current. The Joule heating generated by these pulses can locally melt the material, and the shape of the pulse controls the cooling rate. The latter determines whether the material ends up in a crystalline or amorphous state.

Materials for PCM must, therefore, fulfill several requirements: they must have a strong contrast of resistivity between the amorphous and the crystalline state, but at the same time, they must have a rapid crystallization kinetics to allow for fast switching of memory states. The stoichiometric compound $Ge 2 Sb 2 Te 5$ (GST), originally developed for rewritable optical disks,² is the most widely used active material. However, for some applications (namely, for the automotive market), memories must withstand high temperatures (⁠ $150 ° C$⁠) for several years. GST is not compatible with such requirements because its crystallization temperature, the temperature at which the amorphous phase spontaneously crystallizes, is relatively low (around $150 ° C$⁠). In the case of a spontaneous and unintentional crystallization, the stored information is lost.

One possibility to increase the crystallization temperature is to enrich GST with additional germanium.³ While this Ge-rich GST (GGST) allies excellent crystallization temperature with a reasonably high switching speed, it exhibits germanium segregation,⁴ which leads to the nucleation of a new crystalline phase, almost pure germanium. The active material becomes a two-phase polycrystal, and this complex structure can have a large impact on device performance. Predictive modeling of phase change and segregation would, thus, be of great benefit to aid technological development and to improve the understanding of the underlying mechanisms.

In a previous work, a multi-phase field model (MPFM) has been developed to simulate the evolution of the GGST microstructure.⁵ In this formalism, the local presence of each phase is indicated by a scalar phase field, and the grain structure is described by a field that indicates the local crystal orientation. The evolution of all the fields is coupled to the diffusion of chemical elements, which is driven by gradients of appropriate chemical potentials. The latter depend on temperature, such that the time-dependent temperature field is required for a realistic simulation of memory operations. In our previous work,⁵ the temperature field for a given electric pulse shape was obtained by a separate electrothermal solver. This weak (one-way) coupling neglects the influence of the phase distribution on the electric current, which induces inhomogeneous heating. This precludes a quantitative modeling of memory operations.

In this article, we present a more complete model for phase change memories, obtained by a full coupling of the MPFM with an electrothermal model. More precisely, the equations for electrical and heat conduction are added to the MPFM model, with phase-dependent transport coefficients. This model can account for the multi-physics nature of memory operations and the influence of the material structure on the performances of the memory cell. Since the coupled system of equations is numerically stiff, numerical methods with multiple grids and time steps are required for efficient simulations.

The model contains a large number of material parameters, which are gathered from the literature or taken from in-house measurements. Comparison of the simulations with measurements performed on memory devices allows us to assess the quality of the model and to improve several parametrizations. Notably, it turns out that the thermal boundary resistances (TBR) at the interfaces between different materials have to be taken into account. TBR lead to discontinuities in temperature between the two sides of an interface when the interface is crossed by a heat flux. In macroscopic systems, TBR are usually negligible, but in PCM devices with characteristic size of a few tens of nanometers and high heat fluxes, considerable temperature jumps can occur.⁶ The experimental device characteristics can only be reproduced by our model if TBR are included both on the interfaces between the different parts of the device and on the internal interfaces between different phases of the active material. This is one of the main new results of our work.

Some preliminary aspects of this work have already been published;⁷ here, the full details of the model as well as new results are presented. The paper starts by the description of the simulation setup and the model equations in Sec. II. The equations contain numerous parameters, and Sec. III discusses our choices, taking into account the existing literature. Some details concerning the numerical implementation are discussed in Sec. IV, and the main results are presented in Sec. V.

The geometry of PCM devices has been optimized to limit power consumption while operating the memory cell. The so-called wall architecture, schematically shown in Fig. 1, was found to maximize the current density inside the active volume and, thus, to reduce the programming current. While the exact geometric details of the structure are not important for our purpose, the main features should be reproduced to make a comparison with experimental measurement meaningful. Therefore, we set up a two-dimensional simulation cell that represents a cross section through the wall structure, as shown in Fig. 2.

The electric currents that are used to probe or to change the state of the memory cell flow between a high and narrow electrode shown at the bottom of the figure, the so-called heater, and an extended counterelectrode on the top. The Joule heating used to melt the material mainly occurs in the heater but also inside the active material. The generated heat is evacuated to the surroundings of the memory cell by heat conduction. A complete modeling of memory operations must, therefore, include equations for the electric current, for heat generation and dissipation, and for the structural changes in the active material. Whereas the latter need to be solved only within the active material, electric current and heat flow need also to be taken into account outside the GGST layer. In order to optimize the performance of our code, we solve the different equations on different domains with different simulation grids, as will be described in detail below in Sec. IV.

The domain must be large enough so that the temperature field is not significantly disturbed by boundary conditions. A total domain of $300 nm$ width $×$ $240 nm$ height is used as a good trade-off between resolution speed and realistic thermal profile.

All the details of the crystallization model have been presented in a previous publication⁵ and need not be repeated here; for reference, the model equations and parameters are summarized in  Appendix A.

The crystallization model has been formulated as a multi-phase field model (MPFM). In this mesoscopic approach, the simulation domain can be occupied by different phases that are separated by diffuse interfaces. To capture the segregation effect, three phases are considered: crystalline germanium (Ge), crystalline $Ge 2 Sb 2 Te 5$ (GST), and a disordered phase corresponding either to amorphous or to liquid GGST. Each phase is associated with a scalar phase field $p i$ (⁠ $i=1$ for Ge, $i=2$ for GST and $i=3$ for the disordered phase) that ranges between 0 and 1 (1 when the phase is present and 0 when it is not).

In addition, a concentration field $c$ is introduced, which corresponds to the local excess of germanium with respect to GST composition, meaning that $c=0$ is pure GST and $c=1$ is pure germanium. This amounts to a pseudobinary approximation that restricts the possible compositions of the ternary mixture to the straight segment linking $Ge 2 Sb 2 Te 5$ and pure Ge in the ternary phase diagram.

Interface motion and chemical diffusion are driven, respectively, by grand potential differences between the phases and by gradients of the chemical potentials. Once the free energy functions for all the phases are specified, the equations of motion can be derived from the thermodynamics of irreversible processes. Finally, an orientation-field model is also included to account for the polycrystalline aspect of the GGST material and the formation of grain boundaries. This adds two scalar orientation fields to the crystallization model, but they do not couple to the thermal and electrical models presented below.

The additional dependence of $σ A / L$ on the modulus of the electric field $E$ is important to reproduce the ovonic threshold switching effect¹¹ that is key to PCM operations: at high electric field, the conductivity of the amorphous phase strongly increases to reach a value similar to the one of the crystalline material. This makes it possible to use programming pulses of a constant intensity, regardless of the previous memory state.

In addition to the latent heat and molar volume that were already needed for the crystallization model⁵ (see Tables VI and VII), the thermal and the electrical models require the knowledge of multiple physical parameters: thermal conductivities, specific heats, thermal boundary resistances and electrical conductivities. In the following, we review the data that can be gathered from the literature and describe our choices for the parametrization of our model. For some parameters, measurement done internally on GST and GGST are also exploited.

The thermal conductivities of the three phases (Ge, GST, and amorphous/liquid) in Eq. (2) range over two orders of magnitude. For germanium, the temperature dependence from Ref. 12 is approximated by several linear segments over the 300–1200 K range. They can be constructed from the values listed in Table I. Those values serve as a starting point, but they are reduced by a factor of 8 to match the lower thermal conductivity in germanium thin films¹³ (see discussion in Sec. V F 1).

For the specific heat, the same constant value is used in the three phases, $C m= 26.7 J mo l − 1 K − 1$⁠. It is close to the Dulong–Petit law and agrees with literature data available for germanium^21,22 and GST^23,24 (less than $20 %$ error). This value should hold well for amorphous GGST, but it is probably underestimated for liquid GGST.²³ Nevertheless, as a first approach, we have chosen to also use it for the liquid.

The surroundings of the GGST layer are made of other materials. For simplicity, constant values of $k t h$ and $C m$ are used in these regions; the values used are summarized in Table II. For the top electrode, $k t h$ is provided in Ref. 25. Due to the high Debye temperature of $TiN$⁠, the Dulong–Petit law does not hold and $C m$ increases with temperature.²⁶ However, during simulations, the material stays in the 300–400 K range, and therefore, we use $18.5 J mol − 1 K − 1$ [ $37 J mol − 1 K − 1$ read in Ref. 26 on Fig. 4(c), divided by two because there are two atoms per formula unit of $TiN$].

The bottom electrode temperature is almost constant, too, because of the TBR and its high thermal conductivity; $k t h$ and $C m$ come from Ref. 27.

Contrary to the two electrodes, the temperature in the oxide layers varies a lot. $k t h$ of $Si 3 N 4$ has been measured at CEA-Leti up to 650 K, and it increases linearly. To keep the model simple, and since we do not have any information on the evolution at higher temperatures, we use the same constant value as other authors.^17,28 It corresponds to measurements made at 560 K. $C m$ is also reported to increase with temperature.²⁹ For consistency, we use the value at 560 K, $125.4 J mol − 1 K − 1$⁠, and divide it by the number of atoms.

The heater is made of a material close to $TiN$³⁰ with a higher electrical resistivity. No data are available for this material; however, an article reports a decrease in the thermal conductivity of $TiAlSiN$ as the amount of $AlSi$ increases.³¹ To take this into account, we divide the $TiN$ thermal conductivity by 2. For the specific heat, we could not find any reference in the literature, so we rely on data for $TiN$⁠. Since the temperature of the heater stays quite high for most of the simulation, we consider a higher temperature range than for $TiN$⁠, 600–2000 K.

As for $k t h$ and $C m$ in additional materials, TBR are taken constant, but they vary depending on the materials on both sides of the interface. The values used are summarized in Table III. For the phase change material, we distinguish the amorphous, the crystalline, and the liquid state. The crystalline state correspond to both GST and germanium as both phases are not distinguished. Indeed, most data on TBR at interfaces with $TiN$⁠, $Si 3 N 4$⁠, or the heater correspond to amorphous and crystalline GST (and not GGST). Also, since no data were found for the interface with the heater material, values corresponding to $TiN$ are used instead.

For interfaces between $TiN$ and GGST, the value for the amorphous phase comes from Ref. 32 and for the crystalline phase from Refs. 32 and 33. For $Si 3 N 4$⁠, an author measures the same value for interfaces with amorphous and with crystalline GST.¹⁷ This is inconsistent with what is reported for $TiN$ and $SiO 2$⁠;³⁴ the resistance at an interface with an amorphous phase should be higher. As a consequence, for crystalline GGST, we use the lower resistivity value from Ref. 34. TBR at interfaces with a liquid being impossible to measure, we use a value roughly divided by 20 with respect to the amorphous phase for both $TiN$ and $Si 3 N 4$⁠. The reasoning behind this choice is that we see an increase of the same order of magnitude in the thermal conductivity of GST when the material melts.

For $Si 3 N 4$/ $TiN$ and $Si 3 N 4$/ $W$ interfaces, an article finds a correlation between dielectric/metals TBR and the difference of the Debye temperatures of the two materials.³⁵ Debye temperature of tungsten, $TiN$⁠, and $Si 3 N 4$ are 400,³⁶ 900,²⁶ and 1150 K,³⁷ respectively. This yields the two values reported in Table III. Finally, for the $TiN$/ $W$ interface, data corresponding to $TiN$/ $Al$³³ are used.

The large uncertainties associated with the reported values (see Ref. 32, for instance) and the general lack of data on those parameters make them more prone to significant errors. Therefore, even though some articles report on temperature dependencies of the TBR, we chose to consider only constant values.

Contrary to thermal parameters, electrical conductivities were mostly provided by unpublished in-house material characterization. Since the electrical conductivities of Ge and GST are very sensitive to the presence of impurities, for a meaningful comparison between simulations and experiments, it is necessary to use the values measured on the material that is actually used in the devices. In the three phases present in the active material, the main dependence considered for electrical conductivities is the temperature dependence. In addition, in the amorphous/liquid phase, several additional effects are included, as described below.

The electrical conductivity of the heater is taken as constant and has been calibrated to $5 × 10 4 S m − 1$ (see discussion in Sec. V F 2 below).

Finally, above the melting temperature, the conductivity of the liquid phase is equal to $5 × 10 5 S m − 1$⁠, as reported in Ref. 16. As for the thermal conductivity, amorphous and liquid values are connected linearly. Figure 6 summarizes this information.

While the thermal model is solved on the full domain, the two other models are solved on specific sub-parts. The MPFM is solved in the “active GGST” region (see Fig. 2), the only area where the temperature rise is sufficient to trigger a phase change. Its width has been chosen according to experimental observations and early simulation results. Similarly, the electrical model is solved on a reduced domain comprised of this active region and of the heater. Oxides are approximated as perfect insulators and the two electrodes as perfect conductors.

Boundary conditions also vary from one equation to another. For the crystallization model, the continuity of the microstructure between active and non-active parts of the GGST region is maintained at the left and right interfaces. For the temperature, Dirichlet boundary conditions are applied all around the domain and fixed at 323 K. Finally, for the electrical model, two potentials are applied on the top and the bottom electrodes (if a fixed current is desired, the potentials can be varied with time during the simulation), and zero-flux conditions are imposed at oxide interfaces and other domain boundaries. A summary of the simulation domains and boundary conditions associated with each model is presented in Fig. 7.

The complete coupled model is a numerically stiff multi-scale problem. The characteristic size of the grains and domains in the active material is a few nanometers. Therefore, the thickness of the diffuse interfaces of the phase-field model must be substantially smaller, that is, a fraction of a nanometer. This, in turn, sets the maximal grid spacing that can be used, since the interfaces must be resolved by several grid points. With this constraint, the thermal diffusion inside the bottom electrode (made of tungsten) limits the time step to $dt= 10 − 16 s$ (using the parameters in Table II). For the simulation of a memory operation, which takes place on a microsecond time scale, the required number of time steps would be prohibitively large.

To overcome this limitation, we first tried to switch to an implicit (backward Euler) time scheme. However, improvement made on the time step was offset by the numerical cost of inverting a large matrix. A much better efficiency gain was obtained by two improvements that introduce different grid spacings and time steps for the different model equations.

The first improvement is to solve the thermal model on a coarser grid with $d x ′= r Gdx$⁠, since $dx$ is only limited for the MPFM. The gain on the timestep is significant (⁠ $dt→dt× r G 2$⁠) even for reasonably small values of $r G$⁠. The position of the two grids, one relative to the other, can be seen in Fig. 8. The use of this method requires the transfer of data from one grid to the other. For instance, the temperature is calculated on the coarse grid, but it is also needed on the fine grid to evaluate physical parameters in the MPFM, such as the free energies. Passing a field from the fine to the coarse grid is done by averaging the $r G 2$ fine nodes that are closest to each coarse node. The transfer from the coarse to the fine grid is done with a double linear interpolation (to have every fine node surrounded by four coarse nodes, even at domain borders, we duplicate coarse nodes in accordance with the respective boundary conditions). The two operations are shown in Fig. 9.

With this two-grid method, the actual improvement in computational efficiency is close to $r G 2$⁠, the expected gain from increasing the time step. During the evaluation of the method with $r G$ values from 5 to 11, the simulation times were improved by $88 %$ to $100 %$ of $r G 2$⁠. Therefore, the additional operations required to manage the two grids (see above) do not affect performances significantly.

The second improvement is to solve the thermal model with a reduced time step $d t t h$ compared to the crystallization model (⁠ $d t c r= N d td t t h$⁠). The time step of the thermal model has to satisfy the most stringent stability criterion, while the time step of the crystallization model can be significantly higher. This method implies to solve the thermal model multiple times between each resolution of the crystallization model. Since the time scale of thermal conduction is much lower than the one of crystallization, it is a good approximation to consider that the fields of the MPFM (⁠ $p i$⁠, $c$⁠, etc.) do not evolve between each step of the thermal model. Good performance improvements are obtained when it is much faster to solve the thermal model than the crystallization model (which is the case in our equation system). Indeed, compared to the reference case where both models are solved in each time step, the former is solved the same number of times while the latter is solved $N d t$ times less often.

Several $N d t$ values were evaluated, and the results are shown in Table V. For high $N d t$ values, the speed-up reaches a maximum. This is because the method is intrinsically limited by the ratio of simulation times for the thermal model and the crystallization model.

In addition to the improvements described above, we implemented parallelization with OpenMP (shared memory parallelization). We tuned $r G$ and $N d t$ (⁠ $r G=5$ and $N d t=30$⁠) to ensure an acceptable error on the results. Figure 10 compares temperature profiles with a reference case where all the fields are solved on the fine grid with the same time step. The differences between the two are negligible. After all of these optimizations, most simulations run in only a few hours of wallclock time.

Equation (5) being stationary, the electrical model is not plagued by stability issues. However, it requires the solution of a matrix system, which can be computationally expensive. To reduce the size of the matrix and to improve performance, the coarse grid presented above for the thermal model is also used for the electrical model. In addition, early simulation results showed that the current density and the potential gradient are constant along the heater. Consequently, the domain corresponding to the heater is reduced to one single horizontal line, instead of its entire length. This approximation was tested against simulation on the fully discretized domain, and the error made on the electrical potential is negligible (⁠ $<1 %$⁠).

To solve the matrix, the conjugate gradient method provided by the Eigen C++ library³⁹ is used. At each step, electrical conductivities are determined with the relevant fields (temperature and electric field) calculated in the previous time step, and then the matrix is updated and the system is solved. In Sec. IV A, it was already mentioned that the boundary conditions could be adjusted to force a fixed current instead of a fixed voltage. To achieve this, after the first solution, the current is computed, and if it does not match the target value, one of the two voltage boundary conditions is changed and the system is solved again.

For memory operation simulations, the initial conditions for the thermal and electrical models are simple. The temperature is constant in the domain at 323 K, and no electric current flows through the memory (the potentials at the two electrodes are the same).

The initial microstructure, on the other hand, is a fully crystallized layer of GGST (see Fig. 11), with crystalline germanium and GST. This microstructure is generated by an annealing simulation that aims at reproducing the crystallization of the material (and subsequent segregation) occurring during the fabrication process of the device. As described in more detail in our previous work,⁵ this crystallized layer is obtained by annealing an amorphous domain (here, at 723 K) in which germanium and GST grains nucleate and then grow using the crystallization model. A noteworthy aspect is that during the nucleation process, crystalline seeds are placed randomly, such that different runs lead to different grain distributions. The microstructure presented in Fig. 11 is used for most of the simulations presented in Sec. V.

The model contains a large number of parameters. Some of those (as, for example, the TBR) are not known with great precision, and others (as, for example, the electrical conductivities) can vary with the presence of impurities. Therefore, not surprisingly, simulations with all the values of the parameters as gathered from the literature did not match the experimental data. However, guided by the comparison between simulations and experiments, we were able to refine some of our choices and, thus, to calibrate our model. In order to properly discuss these adaptations, it is useful to have the main simulation results in mind. Therefore, we will first present our main simulation results, obtained with the optimized parameter set, and then discuss below (in Sec. V F) all the changes that have been made to obtain these results and the reasoning behind these modifications.

The programming pulse used during a RESET simulation (similar to the industry standard) is presented in Fig. 12. After the beginning of the pulse, the material heats up, and if the current density is high enough, it melts. After $50 ns$⁠, the current is set to zero, and the memory cools down within $30 ns$⁠.

The evolution of the microstructure during this simulation is displayed in Fig. 13. From the initial fully crystallized layer, a dome-shaped domain (in black) is first melted and then quenched by the fast decrease of the temperature. Between steps 3 (at $t=50 ns$⁠) and 4 (end of simulation), a slight recrystallization begins but quickly stops as the temperature drops. This is visible through the small reduction in the amorphous domain. The dependence of the interface mobility on temperature is crucial to obtain this behavior. As described in our previous paper,⁵ the mobility is obtained from measurements of the crystallization velocity of pure GST as a function of temperature.⁴⁰ As the temperature decreases from the melting temperature, the velocity first increases because of the growing thermodynamic driving force but then decreases again because all atomic rearrangement processes slow down. Therefore, when, as in Fig. 13, the temperature decreases very rapidly, the interface motion cannot keep up with the temperature field, and the liquid is quenched into the solid amorphous state. The final structure is consistent with experimental measurements.⁸ 

The programming pulse used to simulate a SET operation is presented in Fig. 14. Due to the progressive decrease in the imposed current, the evolution of the microstructure is different (see Fig. 15). A similar melting occurs, but the temperature of the crystal/liquid interface stays high enough during cooldown to enable full recrystallization of the active material.

Comparing the snapshots 1 and 4 of Fig. 15, we can see that all the germanium grains that were located at the center of the melted zone have disappeared, whereas some germanium grains that were at the edge of the melted dome end up being larger than before the operation. Thus, there has been a net germanium redistribution from the center to the border of the melted dome. This is also observed experimentally.⁸ To characterize this effect more precisely, we consider three annuli with increasing radius (see Fig. 16) and we track the evolution of the average excess germanium concentration during the SET operation, which is displayed in Fig. 17. The time axis is divided into three parts corresponding to the three phases of the operation. During the initial melting of the material, the outer annulus (only partially melted) sees its concentration slightly increase, while the two other areas (fully melted) homogenize. After that, during the regrowth of the material, as the crystal/liquid front moves toward the center of the memory cell, first, the concentration of the outer annulus increases, and then one of the intermediate annulus. At the same time, the concentration of the liquid (green and orange curves first, then only orange) steadily decreases. This is due to two factors: to grow, germanium grains must absorb germanium from the liquid, and, moreover, the concentration of regrown GST is higher than at the beginning [see the lighter shade of gray on Fig. 16(b)]. The last phase of regrowth starts slightly before the third snapshot picture of Fig. 15, when the growth of the germanium grains stops. Since, in the further evolution, only GST recrystallizes, germanium segregates into the liquid, whose concentration increases in consequence. This explains why, in Fig. 16(b), the closer to the heater (where the last liquid area recrystallized), the higher the concentration. Figure 17 confirms that after the operation, a germanium redistribution toward the edge of the melted dome occurred. The localized accumulation of germanium in the liquid at the end of the recrystallization process is of second order to this global redistribution.

The model is able to reproduce the trend observed experimentally: the germanium concentration increases at the edge of the melted dome while the concentration of the core is lower than initially. However, the magnitude of the simulated concentration redistribution is lower than in experiments. Moreover, the large increase of the liquid concentration at the end of recrystallization is not seen in the experiments. More detailed comparisons would be necessary to determine the exact origins of these discrepancies.

The coupled models can be exploited to visualize the electrothermal fields during memory operations in the GGST segregated microstructure. First, we can look at the temperature profile in the memory. Figure 18 shows the temperature in the full domain at $50 ns$⁠, in two situations: at the top, the programming current is too low to melt the material (Fig. 11 is the corresponding microstructure); at the bottom, the current is higher (same simulation as Fig. 13). In Fig. 18(a), the large differences in thermal properties of germanium and GST are visible: germanium grains clearly stand out. In Fig. 18(b), the temperature is more homogeneous in the melted dome. Another important difference is that if the material stays crystalline, the temperature is almost equal at the top of the heater and at the bottom of GGST, whereas if it melts, the temperature drops in the phase change material. These two figures also highlight the impact of TBR: a sharp drop of several hundred kelvins that occurs at the interfaces of the memory element, in line with values obtained in the literature.^6,41

We can also look at spatial maps of current density and Joule heating. In Fig. 19(a), it can be seen that the electric current avoids germanium grains because of their higher electrical resistivity. This leads to a reduced Joule heating in these grains [see Fig. 19(b)]. In Fig. 20, after the melting of the material, the Joule heating has drastically decreased due to the high electrical conductivity of the liquid phase.

To avoid any confusion, one must understand that unlike the two temperature maps of Fig. 18 that correspond to different programming current (but are both obtained at the end of a $50 ns$ RESET pulse), Figs. 19(b) and 20 correspond to the same programming current. The former is obtained close to the beginning of the pulse (after less than $1 ns$⁠), whereas the latter is obtained after $50 ns$⁠.

We have performed simulations to reproduce two figures of merit that are commonly used to characterize the memory devices: a $R(I)$ plot in Fig. 21 and a $I(V)$ plot in Fig. 22. The $R(I)$ curve is obtained by performing several RESET operations at different programming current intensities $I$⁠, followed by the reading of $R$⁠, the total resistance of the memory cell. Each data point corresponds to one simulation. Similarly, the $I(V)$ curve is obtained by performing several simulations with different $I$⁠. Those simulations are longer RESET pulses (⁠ $100 ns$ instead of $50 ns$ as shown in Fig. 12), which lead to an almost stabilized microstructure. The voltage $V$ across the cell is read before the end of the pulse. The result of this procedure is $V$ as a function of $I$⁠, which can then be inverted to yield the more customary current–voltage characteristic $I(V)$ for comparison with data from experiments. The two pulses are illustrated in Fig. 23.

In Fig. 21, the model reproduces experimental measurements well. At a certain current, the resistance starts to increase because the Joule heating is high enough to melt (and, after the cooling, amorphize) a part of the active material. The effect saturates at higher currents since the dome radius reaches the GGST layer thickness. The resistances at low and high currents correspond, respectively, to the values of the fully crystallized layer and of the amorphous phase. The agreement between simulated and experimental $I(V)$ curves in Fig. 22 is also excellent. At high current, when the active material melts and becomes significantly more conductive, the drop of the applied voltage essentially occurs in the heater. This can be used to determine the electrical properties of this material (see Sec. V F 2).

The impact of the initial microstructure has also been studied. As explained in Sec. IV D, all the memory operation simulations presented so far have used Fig. 11 as a starting point. However, variations in the microstructure, in particular, the initial position of the germanium grains, can lead to different electrical behavior of the memory. A total of nine additional microstructures were generated randomly by multiple annealing simulations (see Fig. 24), and both $R(I)$ and $I(V)$ simulations were performed with each of them. In Fig. 25, the curves associated with the different microstructures are averaged and compared to experimental data. Individual curves are also displayed with thinner lines to visualize the variability of the results. On the $R(I)$ plot, individual curves are scattered; in particular, the onset current leading to the formation of the melted dome varies significantly. However, when taking the average of the ten curves, a good agreement with experimental data emerges.

Now, we will discuss the process that has allowed us to calibrate the model parameters and to obtain the good agreement with experiments that was displayed above.

The thermal conductivity of crystalline GGST, that is, the two-phase polycrystalline material, has been measured up to 700 K.¹⁷ With the values of the thermal conductivities gathered from the literature, simulations with our model, using the microstructure of Fig. 11, yielded an effective thermal conductivity (that is, the ratio of the total heat flux through the entire sample, divided by the temperature difference between top and bottom) that was far too large. This problem was cured in two steps.

First, it was realized that the thermal conductivity of germanium thin films¹³ is significantly lower than the one initially used in our model, corresponding to bulk germanium samples.¹² Unfortunately, the literature does not provide thermal conductivity values for germanium thin films over the whole 300–1200 K temperature range. Instead, it was decided to divide the values reported for bulk germanium samples by a factor of 8; this division is already included in the data shown in Fig. 4. Even with this modification, the model still yields an effective thermal conductivity that is more than $60 %$ higher than the value measured at 300 K. To obtain the experimental value, a reduction in the germanium conductivity by a factor of 40 would have been necessary, which does not seem reasonable. Instead, TBR are added to the model at germanium/GST interfaces. Indeed, just as the “device” interfaces between dissimilar materials, the “internal” interfaces inside the active material should add thermal resistance.

The value of $R 12$ has been calibrated to yield the correct value of the effective thermal conductivity of the two-phase polycrystalline material at 300 K (⁠ $0.8 W m − 1 K − 1$⁠, as in Ref. 17). As shown in Fig. 26, simulations with this value of TBR match well with experimental measurement over an extended temperature range. Furthermore, the value $R 12= 8 × 10 − 9 K m 2 W − 1$ is of the right order of magnitude when compared to the values for the other interfaces provided in Table III. In summary, including TBR in the GGST microstructure seems to be important to accurately reproduce the thermal behavior of the material.

As mentioned along Fig. 22, the electrical properties of the heater can be determined from the $I(V)$ curve. Initially, the electrical conductivity of the material was set at $6.5 × 10 4 S m − 1$⁠, according to in-house measurements done at low temperatures (up to a few hundred degrees Celsius). However, with this value, the simulated current was systematically too high. Since the temperature is much higher during memory operations (especially for high currents), an increase in the electrical resistivity can be expected, as in most metals. Good agreement with the experimental $I(V)$ curve was obtained with the value of $5.0 × 10 4 S m − 1$ already given in Sec. III D.

With the initial parameters, we observed that most germanium grains located inside the liquid dome did not completely melt during a RESET pulse. Indeed, since the melting temperature of GST is lower than the one of germanium, the melting process was faster in the GST matrix, and isolated germanium grains remained within the liquid. After cooling, thus, crystalline germanium inclusions were present in the amorphized dome. This was never observed in the experiments. We had initially thought that the inclusion of the phase-dependent electrical conductivity would resolve this issue, since the higher resistivity of germanium should lead to locally stronger Joule heating. However, this is offset by the inhomogeneous current distribution shown in Fig. 19(a). Therefore, the melting of the germanium grains must be accomplished by heat generated in the heater and transported through the liquid. To promote this melting process, we have lowered the thermal conductivity of the liquid by a factor of 5; this leads to higher temperatures close to the heater inside the liquid, since the heat is more strongly confined. The data shown in Fig. 4 do not include this factor.

The last adjustments concern the $R(I)$ curve (Fig. 21). Initially, the transition between small and large values of $R(I)$ was very abrupt, even if at the same time the size of the amorphous dome was increasing gradually with the current. To obtain the more progressive behavior shown in the figure, it was necessary to take into account the Poole–Frenkel effect as given in Sec. III D. Indeed, this effect reduces the resistivity of small amorphous domes, which leads to a more progressive increase in resistance with the size of the amorphous dome.

On this aspect, the electrical conductivity of the liquid is also key. Indeed, thanks to its high value, there is almost no Joule heating in the melted dome as shown in Fig. 20. This enables more progressive dome sizes by removing most of the Joule heating in the GGST as the melted area expends.

In order to model the evolution of the GGST microstructure during memory operations, a fully self-consistent coupling of the MPFM with a thermal model and an electrical model has been presented. The thermal model accounts for heat generation and its propagation in the memory, while the electrical model provides the electric current density (used to compute the Joule heating) and is used to determine the electrical resistance of the memory. The numerous additional parameters needed have been obtained from literature data and in-house material characterization. The complete model has then been calibrated to match the experimental data.

SET and RESET operation simulations demonstrate the ability of this model to qualitatively reproduce the GGST amorphization and crystallization depending on the programming conditions. It is also possible to exploit the model to study the evolution of electrothermal fields (temperature, current density, etc.) during operations. For instance, we showed that the melting of germanium grains cannot be explained by a higher Joule heating (compared to GST) due to their high electrical resistivity. Beyond that, our simulations of electrical figures of merit are in good agreement with experimental data, even when considering changes in the microstructure of the material. Finally, we have also evidenced that thermal boundary resistances were needed at germanium/GST interfaces to properly reproduce the thermal conductivity measured on crystalline GGST.

Despite these good results, some caveats are in order, and some lines for future improvement can be outlined. Most importantly, our simulations are two-dimensional. The real three-dimensional grain structure is very different from the columnar structure that would be obtained if our two-dimensional heat maps were extended along the third dimension. Furthermore, the geometry of the conduction paths inside the GST phase would be quite different in two and three dimensions. However, it is known that grain coarsening exhibits the same scaling laws in two and three dimensions, only with different prefactors. Since the low-conductivity germanium phase forms isolated islands and does never percolate, there should be a quantitative, but no qualitative difference concerning conduction. This explains why, although two-dimensional, our model still captures many relevant features. For simulations in three dimensions, the orientation-field model needs to be adapted. Indeed, in three dimensions, the orientation of a grain needs to be specified by three scalar parameters, which can be the classic Euler angles or the elements of a unit quaternion.⁴² Furthermore, keeping acceptable simulation times despite the large number of discretization points would require more sophisticated parallelization methods, such as MPI or GPU methods. This is an interesting subject for future work. As an additional benefit, such three-dimensional models could potentially capture the full scope of three-dimensional complexities inherent to many real-world device architectures with arbitrary boundary shapes.

Another important point is the lack of knowledge about parameters. We had to determine many of them approximately due to the lack of more precise information in the literature; this mainly applies to the liquid phase but not exclusively. A systematic sensitivity analysis could reveal which of these parameters are most influential for the device performance. However, in view of the large number of parameters, this represents an enormous amount of work, which is outside the scope of the present article. As an alternative, advanced simulation methods such as molecular dynamics with neural network-trained interatomic potentials⁴³ could be exploited to obtain more precise values at least for some parameters.

Another area of improvement is the crystallization model. In this paper, we mainly focused on the thermal and electrical models and on their calibration; however, some aspects of the crystallization model, such as its kinetics, play a key role in the physics of the device. To this end, simulation results could, for instance, be compared to experimental measurement of the SET resistance as a function of the quenching time. Also, the pseudobinary approximation limits the composition of the system (see Sec. II B). In particular, it is not possible to include pure Sb lamellas as reported recently in the literature.⁴⁴ However, extending the current multi-phase field model to the full ternary phase diagram is very challenging. Set aside the complex mathematical formulation, it would require a lot of additional physical parameters, notably to model the thermodynamic behavior of all the different phases.

Currently, this model already provides an excellent framework to investigate diverse programming conditions and material change. In particular, it is able to simulate the forming operation, the first electrical activation needed by GGST alloy in order to exhibit good switching behavior. This forming operation is known to deeply influence the microstructure.⁴⁵ While the results presented in this paper have not considered post-forming microstructure as starting point for SET and RESET operations, future exploitation of the models could reveal how these operations are impacted by the forming step. Moreover, the capability of the models to explore various programming pulses extends the model's utility to address the multistate capability of the PCM. This capacity of multistate is of major importance for in-memory computing⁴⁶ and neuromorphic applications, where the gradual change of resistance is instrumental in mimicking the analog nature of biological synapses.

This work was financially supported by the Association Nationale Recherche Technologie (ANRT), France, and STMicroelectronics through the CIFRE Contract No. 2020/1070.

The authors have no conflicts to disclose.

Robin Miquel: Investigation (lead); Software (lead); Visualization (lead); Writing – original draft (lead). Thomas Cabout: Conceptualization (equal); Funding acquisition (lead); Investigation (supporting); Methodology (equal); Project administration (equal); Software (supporting); Supervision (equal); Writing – review & editing (supporting). Olga Cueto: Conceptualization (equal); Funding acquisition (supporting); Investigation (supporting); Methodology (equal); Project administration (equal); Software (supporting); Supervision (equal); Writing – review & editing (supporting). Benoit Sklénard: Investigation (supporting); Methodology (equal); Supervision (supporting); Writing – review & editing (supporting). Mathis Plapp: Conceptualization (equal); Investigation (supporting); Methodology (equal); Project administration (equal); Supervision (equal); Writing – review & editing (lead).

The data that support the findings of this study, and in particular the simulation code, are available from the corresponding author upon reasonable request.

In this appendix, all the equations of the crystallization model that are used in the simulations are summarized for completeness; full details and further explanations can be found in our previous paper.⁵ 

For the definition of the free-energy and grand potential functions for each phase, we use a regular solution model. We use the same functions as in Ref. 5, with the same nomenclature but with slightly different values for some of the coefficients; all the parameters used to generate the phase diagram are listed in Table VI, and the phase diagram is displayed in Fig. 27.