Phase change memory^[1-3] has been one important non-volatile memory among several emerging technologies such as RRAM, MRAM and FeRAM, etc. It relies on the switching between amorphous state (high resistance) and crystalline state (low resistance) for memory applications. In general, phase change memory (PCM) is one type of memristor. The need to handle massive amount of data in the era of big data drives the extensive research and developments of these memory technologies. Among various applications, the storage-class-memory (SCM)^[4] calls for candidates with high capacity and high band-width to fill the gap between DRAM and NAND Flash in the state-of-the-art memory hierarchy. PCM is one of the strong candidates for SCM applications^{[1, 3]}.

As of today, PCM has entered the market of standalone memory, for example the 3D XPoint memory by Intel^{[5, 6]} and embedded non-volatile memory (eNVM), for example, that of STMicroelectronics^[7-9]. This comes from years of research and development all around the world. The original concept of phase change memory was proposed early in 1970s^[10]. With the success of phase change material Ge₂Sb₂Te₅ in the optical storage, utilizing the phase change material for electrical storage revived in later 1990s^[11]. Several milestones in the process include the mushroom-type PCM realized with 180 nm CMOS technology^[12], the line type PCM with 90 nm technology node^[13], the embedded 90 nm 1 transistor 1 resistor PCM technology^[14], the dash-type confined PCM with sub-20 nm technology^[15], and the stackable cross point PCM with the Ovinic threshold switching (OTS) selector^[5], and the two-deck cross point PCM for 128 Gb SCM^[16]. With the advancement of technology scaling and material innovations, the state-of-the-art PCM features 4F² cell area, tens of nanoseconds set time, 100 μA reset current, 10¹² endurance cycles, 10 years’ retention at temperatures higher than 120 °C. Fig. 1 shows the schematics of two typical PCM cell structures: the mushroom type and the confinement type. Alternative applications of PCM are also emerging, such as the multi-level cell (MLC) and analog memories^[17-20]. One key factor is to program the PCM into intermediate states between the highest resistance and lowest resistance states. Both SET and RESET processes could be utilized for this purpose. As it takes some time to complete the phase change process, intermediate states are then possible by partial phase changes^[20]. The analog memory properties are being used to represent an electrical synapse with plasticity or an electrical neuron for neuromorphic computing^[21-25] systems. In most cases, the PCM cells are connected together with CMOS devices, especially for realizations of the spiking neuron network.

Along with advancements of the PCM technology, their compact models had been playing important roles. Memory device models are needed for the same purpose of all other device models: to be used for circuit simulations. Although the memory cells are generally repeated in an array or a chip and the reading operation is more or less regarded as a single cell operation, a memory model is still necessary for the memory technology development^[26]. A PCM cell works with a selector device to form the basic bit unit, and a simulation with models facilitates the device matching. The neighboring effects from adjacent PCM cells, as well as the parasitic effects, are accounted for in the model-based simulations. Further, the circuitry besides of the memory array including the sense amplifiers, are designed taking into consideration with the cell properties and the word/bit line voltage drops, etc. Memory reliability^[27] is another concern that could be addressed in circuit simulations with memory models. In the era of neuromorphic computing^[28] with the analog phase change memories, the need for a compact memory model becomes even stronger. The model needs to reproduce not only the high/low resistance states, but also the intermediate states, the state transition time, their voltage dependence, and so on. Reproductions of these properties with the model allow evaluations of the analog cell characteristics, as well as the co-evaluations with the CMOS transistors for analog-type circuits. For advanced technologies, device-technology co-optimizations (DTCO) of the PCM-based computing applications are also essential. All these require a dedicated compact model supporting SPICE circuit simulations.

In this work we review compact models of the phase change memory. Device physics and principle operations as the basis of compact modeling are explained in Section 2. Then reported compact models of different categories from behavior macro models, to physics-based models with different approximations are reviewed in Section 3. Finally, further challenges on PCM modeling and our perspective are explained in Section 4.

There are three distinct states in the operations of PCM as shown in Fig. 2(a): the crystalline state (C), the amorphous state (A) and the melted liquid state (M), although only two of them (C and A) are possible at room temperature. The transition from state A to C is called crystallization. The transition from state C to A is not a direct process but via the state M, and is usually called amorphization. Fig. 2(b) shows the programing temperature for crystallization (also called SET) and amorphization (also called RESET). The SET process is triggered by increasing the temperature above the crystallization temperature (T_cryst or simplyT_c). There are two possibilities in the process, a nucleation kinetics and a growth kinetics, that re-arrange atoms to form an ordered structure. The RESET process is triggered by increasing the temperature above the melting temperature (T_melt or simplyT_m) and then quenching, i.e., decreasing the temperature within a short period. In the process the crystalline structure is melted with the chemical bonds broken, and the atoms re-arrange forming a disordered structure. For memory applications, the energy in the SET and RESET process comes from the Joule heat of the electrical pulse. Some further details in addition to the above essential ones of PCM operations are reported in excellent reviews^{[5, 29-31]}.

Another key property for the SET process is the threshold switching. It refers to an increase in the conductance after the voltage across the amorphous region reaches a critical value. Only with the threshold switching property, Joule heating for increasing the temperature above T_c happens with a reasonable voltage, which otherwise will be too large for applications. In the current–voltage curves, a snap-back is defined as the region with decreasing voltage and increasing current, i.e., a negative differential resistance. The threshold switching effect provides a dynamic resistance which is comparable to the crystalline resistance to enable the SET process. While the physics of threshold switching is still yet to be fully investigated, there are a few proposals such as a thermal model^[33], an impact ionization model^[34], and a thermal runaway model^[35].

A popular theory for the crystallization is the classical nucleation and growth (CNG) theory^{[36, 37]}. In this theory framework, the crystallization is composed of two processes: a nucleation process and a growth process, which have different temperature dependence. The nucleation pictures the formation of a nucleus with a critical size r_c: after the nucleus is larger than this size, it becomes stable and further grows without dissolution. As a result, the nucleation is indeed a stochastic process. With a surface energy γ, and a bulk free energy gain dgdue to the supercooling, the energy of a nucleus with the critical size r_c is:

$ \Delta G = \frac{{16}}{3}\frac{{\pi {\gamma ^3}}}{{d{g^2}}} . $  (1)

The nucleation probability P_n^{[36, 38, 39]} during a time interval Δt is then expressed with the activation energy at the surface E_a1:

$ {P_{\rm n}} = \alpha \Delta t \, \exp \left[ { - \beta \left( {{E_{\rm{a}1}} + \Delta G} \right)} \right] , $  (2)

in which α is a frequency factor,β = 1/(kT) with k the Boltzmann constant and T the temperature. Detailed derivations have been reported in the literature^[39]. Fig. 3(a) plots a schematic energy diagram for the nucleation process.

A deterministic model, which can be regarded as a variant of the above probabilistic model, has a similar formulation of the phase transition rate^[40] based on the maximum transition frequency .

$ {c_N} = r_0^N \; \exp \left[ { - \beta \left( {{E_{{\text{a}}1}} + \Delta G} \right)} \right] . $  (3)

At the same time, the stabilized nucleus grows at any interface between the crystalline and amorphous. The growth is regarded as a combination of an atom diffusion Arrhenius process (with the activation energy E_a2) and a bi-directional process of liquid-solid transition. As a result, the growth velocity is given by^{[36, 39]}:

$ {v_{\rm g}} = f{a_0}\alpha \left[ {1 - \exp \left( { - \beta \left| {dg} \right|} \right)} \right]\exp \left( { - \beta {E_{\rm{a}2}}} \right), $  (4)

in which is the atomic jump distance, andf is a temperature-dependent parameter. In the probabilistic picture, the growth probability during a time interval Δt is:

$ {P_{\rm g}} = pa\Delta t\left[ {1 - \exp \left( { - \beta \left| {d{{g}}} \right|} \right)} \right]\exp \left( { - \beta {E_{\rm{a}2}}} \right), $  (5)

which is transformed from Eq. (4) with f replaced by an exponential factor, and p is a fitting parameter. Fig. 3(b) plots a schematic energy diagram for the growth process. Alternatively, the growth rate similar to the above nucleation rate is derived:

$ {c_G} = r_0^G\left[ {1 - \exp ( - \beta \left| {dg} \right|)} \right]\exp \left( { - \beta {E_{\rm{a}2}}} \right), $  (6)

in which is the frequency factor.

The above Eqs. (2) and (5) have been used in a probabilistic framework^[38], and Eqs. (3) and (6) have been used in a deterministic and local framework^[40]. Fig. 4 plots the nucleation and growth rate which show peaks with non-monotonic temperature dependences.

Based on the CNG model, a phase space framework has been reported^[40] covering all the three phases as shown in Figs. 5(a) and 5(b). At the same time, it is noted that the melting, either from crystalline state or amorphous state, is usually treated as a simultaneous process. When the temperature is raised above the melting temperature T_m, the state is transformed to molten state without a noticeable delay. This is reflected in the phase space framework as shown in Fig. 5(b) with sudden increases of the transition rates of A-to-M and C-to-M.

Another theory to describe the crystallization is the Johnson–Mehl–Avrami–Kolmogorov (JMAK) theory^[41-43]. It assumes isothermal and isokinetic (same temperature dependence of nucleation and growth) conditions, and describes the transformed fraction X(t):

$ X\left( t \right) = 1 - \exp \left[ { - {{(At)}^n}} \right] , $  (7)

$ A\left( T \right) = w\,\exp \left( { - \frac{{{E_{\rm a}}}}{{kT}}} \right) , $  (8)

in which w is a frequency factor and n is the Avrami coefficient, E_a is the activation energy. The JMAK theory may be regarded as an approximation of the CNG theory in a specific temperature range.

The PCM state, either SET, partial SET or RESET, can be achieved by designing the melting and cooling (usually termed as quenching) process. The RESET to high resistance can be done by a “melting and fast cooling” process after melting the memory active region, as shown in Fig. 2. Noticing the quenching time-dependent chalcogenide atoms re-organizations, the SET to low resistance can be done by a “melting and slow cooling” (MSC) method^[44]. It means by controlling the quenching time for atoms reorganizations partial or full SET states can be achieved. Fig. 6 also repeats the traditional “solid phase crystallization” (SPC) method of Fig. 2, in comparison with the MSC method.

In this section, models with different approaches and different levels of abstractions are described, from the macro models, to more physics-oriented models.

One key feature of the PCM is the switching from the amorphous state to crystalline state, or from the amorphous state to a dynamic on state. In the most common current–voltage characteristics, the above switching shows a snap-back. One essential starting point for a macro model is to find an appropriate representation for the device, either a pure mathematic one or a physics-based one for other well-studied devices. In the latter case, a lumped SPICE model may be generated by borrowing the formulations of devices with similar I–V and the switching.

A macro model of PCM was reported^[45]. The core block of the model is shown in Fig. 7(a), and extensions to this block to include the ‘conductive on’ and transient simulation capability were also reported. There are also variations to the core block in latter models, but the main idea is clearly reflected in Fig. 7(a). Including in the ‘word circuit’, interlinked switches are included in the Zin part, determining its outputs (Vset, and Vreset) to the RSff which is a flip-flop. When the input to the Zin part is a lower voltage, Vset is high and Vreset is low. The RSff part generates a high voltage, which is used to control the ‘Zmodel’. For a high Vffout, the switches S1 and S2 respond correspondingly: S1 closed and S2 open. As a result, the resistance sensed from the Isense terminal is different depending on the state of switches. In the S1 branch, the amorphous resistance is included, together with another compensation resistance, to form the high-resistance PCM state. In the S2 branch, the crystalline resistance is included, and another diode is also included, to form the low resistance PCM state which shows a non-linear dependence. The memory state changes when the input to the Zin part goes high, as the flip-flop RSff outputs a low Vffout and the switches S1 and S2 changes their status. As is seen, the macro model used the flip-flop logic for the core memory function, and switches to control the transitions between PCM resistive states.

To model the snap-back, a silicon-controlled rectifier (SCR) type formulation was used for a softer transition from the high-resistance to low-resistance state. Some considerations on the transition smoothness are incorporated in the model with additional circuit components. Fig. 7(b) shows the I–V curves from the macro model in comparisons with the experimental data. For transient simulation, it is necessary to account for the duration of the input pulse, especially for the SET process. An extension to the Word circuit with an integrator is included to detect the input amplitude and duration, and only activates the switching when the crystallization/amorphization is expected to finish after a period.

A similar macro model was reported^[46] with a similar approach but extended features. Fig. 8 shows the flow chart of the macro model. A thermal circuit and a crystallization module were introduced on top of the bi-stable circuit. The two modules were both implemented with circuit macros like amplifiers based on the physics of Joule heating and crystallization kinetics. The Joule heat is calculated for the temperature, which is fed into the crystallization module for a trace of the crystal fraction. As a result, multi-level storage is supported in addition to the binary storage for the previous model. The trade-off between pulse height and width for programming was demonstrated. Similar macro models were reported in the literature^{[29, 47-49]}. It is interesting to note that the macro model^[48] is used to demonstrate a spiking neuron circuit based on PCM.

Macro models are relatively straightforward to be developed without knowing the device physics. As seen above, existing circuit blocks are connected to reproduce the observed PCM properties. The PCM macro models have been assisting the technology development at the early stage. One general drawback of macro models is a high complexity and long computational time. Thus physics-based models are usually pursued when the technology becomes more mature.

A single piece model was reported^[50] for the first time. The model description of the conduction in the amorphous and crystalline phases are both based on the avalanche of chalcogenide materials, Eq. (9). n_status is a multiplication factor in the unit of V^–1, and R₀ describes the low field resistance. Both are regarded as functions of the phase change material status, i.e., the crystalline fraction. For lower voltages, Eq. (9) is simplified as I =V₁/R₀, hence linear. For higher voltages, the current depends on voltage exponentially. The current equation (9) takes the temperature dependence into consideration, which satisfies the Arrhenius law. T_ref is reference temperature. E_a,status is the activation energy which can be extracted from the temperature measurement. Fig. 9 plots the I–V characteristics with two states.

$\begin{array}{l} {I_\alpha }\left( {{V_1},{V_{\rm{status}}}} \right) = \dfrac{1}{{{n_{\rm{status}}} {R_0}\left( {{V_{\rm{status}}}} \right) \cdot \exp \left[ {\dfrac{{{E_{\rm{a,status}}}}}{k}\left( {\dfrac{1}{T} - \dfrac{1}{{{T_{\rm{ref}}}}}} \right)} \right]}} \\ \qquad\qquad\qquad\times\left[ {\exp \left( {{n_{\rm{status}}} {V_1}} \right) - 1} \right] . \end{array}$  (9)

A smooth function F(x, x_t) is proposed for the threshold switching behavior. It can be any function which switches from 0 to 1 around a threshold value x = x_t. In practice, a Fermi-like function is a good option. For the PCM modeling, a parameter of threshold current (I_th) is used to represent the critical condition of threshold switching, and the current–voltage relationship is described by Eq. (10). Starting from the high resistance (e.g. fully amorphous) state, the current is given by Eq. (9) before it reaches I_th. Once the current reaches I_th, the smoothing F(I₁, I_th) goes from zero to one, and the PCM switches to the crystalline state. The conduction is then described by the on state currentI_on, which is part of the I_αSET. The ‘versatile’ function provides a flexible description of the threshold switching, which is also termed as a voltage snap-back. Anyway, the current always increases, unlike the voltage drop across the PCM.

$ {I_1} = {I_\alpha }\left( {{V_1},{C_x}} \right) + F\left( {{I_1};{I_{\rm{th}}}} \right) \cdot \left[ { - {I_\alpha }\left( {{V_1},{C_x}} \right) + F\left( {{V_1};{V_{\rm{hold}}}} \right) \cdot {I_{\rm{on}}}\left( {{V_1}} \right)} \right] . $  (10)

The cross point of the snap-back branch and theI_αSET branch is (V_hold, I_hold), which satisfies the crystalline stateI_αSET of Eq. (9) as well as the load line below:

$ I = - \frac{V}{{{R_{\rm{load}}}}} + \frac{{{V_{\rm{th}}}}}{{{R_{\rm{load}}}}} , $  (11)

in whichR_load is the PCM load resistance (the resistance in series, e.g. the heater resistance).V_th is the threshold voltage. In the implementation, the PCM together with a given load resistance R is treated as a voltage controlled current source (VCCS). At the threshold voltage, the passing current changes from the I_αRESET to I_on.

The model^[50] implemented a mechanism to monitor the quenching time after melting temperature has been reached for the MSC method in Fig. 6. Auxiliary circuits including several blocks are proposed for this purpose. A table lookup approach is used to determine the SET resistance according to the quenching time. Advantages include the ability to model the multi-level cell (MLC) operation. Note that for partial SETs the parameter of V_status between that of the RESET and full SET goes into Eq. (9) to determine the high-resistance current conduction. The threshold swing is chosen to be reorganized by a constant threshold current I_th, which means a reduced threshold voltage for partial SET PCM.

An improved model without involving the table-based approach for crystallization was reported^[51]. A more physics meaningful crystal fraction C_x is used instead of V_status to represent the degree of crystalline state. The time evolution of C_x, i.e. the crystallization kinetics, is developed as follows with a temperature-dependent time constantτ:

$ {C_x} = 1 - \exp \left( { - \frac{t}{\tau }} \right) . $  (12)

It is indeed identical to the JMAK theory with Avrami coefficient n = 1. The phase change material resistance is determined with C_x:

$ {R_x} = {R_{0\rm{c}}} + \left( {{R_{0\rm{a}}} - {R_{0\rm{c}}}} \right) \left( {1 - {C_x}} \right), $  (13)

in which R_0c and R_0a the resistance of state which are completely crystallized and amorphized, respectively.

Fig. 10 shows that Eq. (12) agrees well with the experimental data of resistance evolution in the time domain under a certain temperature. With Eqs. (12) and (13), the crystallization process can be described step-by-step with the temperature at each step. For the implementation, auxiliary circuits in Fig. 11 to integrate for the crystal fraction C_x is introduced. A temperature circuit calculates the heat and obtains the temperature. Taking the temperature as inputs, the timer circuit accumulates time when the temperature is between the crystallization temperature T_c and melting temperature T_m. The crystal fraction circuit implements the Eq. (12) with inputs from the timer circuit. At the same time, C_x goes to zero once the temperature reaches T_m. The crystal fraction C_x is fed back into the current–voltage equations, and a self-consistency loop is implemented in the model.

With the model^[51] pioneering the physics-based PCM modeling, both SET programming themes of MSC and SPC are in principle supported. Compared to macro models, the model calculation efficiency is highly enhanced due to the reduced complexity. There are variations in this model category, for example, with or without switches in the implementations.

A model with rate equations of crystallization and amorphization was reported^[52]. For the dynamic phase transitions, the threshold swing was modeled not only in the voltage/current domain but also in the time domain. By observing a switching time in the range of τ_f = 0.15 ns, a dynamic version of the “versatile” function, Eq. (14) was developed. It is playing a role of smoothing in the transient simulations of the PCM, i.e. the setting process. When I_gst crosses the threshold current I_th, the variable F asymptotically approaches zero or one in the time domain. Fig. 12 plots the waveforms of the PCM voltage, current as well as the variable F.

$ \frac{{\rm{d}F}}{{\rm{d}t}} = - \frac{{F - \theta \left( {{I_{\rm{gst}}} - {I_{\rm{th}}}} \right)}}{{{\tau _f}}} , $  (14)

in which θ(x) is the step function, which is zero or one when x < 0 or x > 0.

The amorphous volume V_a is formulated with the classical nucleation theory instead of the crystallization kinetics Eq. (12) of the JMAK theory used before. In a nucleation-driven case, the crystallization rate is given by P_nV_nV_a/V_m, where P_n is the nucleation probability per unit time, V_m is the volume of a monomer, and V_n is the volume of a nucleus. In a growth-driven case the crystallization rate is given by S_av_g, where S_a is the area of amorphous/crystalline interface, v_g is the growth velocity. For the purpose of model simplicity, the nucleation and growth processes are not treated as successive but concurrent processes. By summing these two crystallization processes, a kinetic equation for C_a is developed in the model with:

$ {\left( {\frac{{\rm{d}{V_{\rm{a}}}}}{{\rm{d}t}}} \right)_{\rm{c}}} = - \left({P_{\rm{n}}}{V_{\rm{n}}}\frac{{{V_{\rm{a}}}}}{{{V_{\rm{m}}}}} + {S_{\rm{a}}}{v_{\rm{g}}} \right) . $  (15)

The model is also self-limited, and the RHS approaches zero when V_a comes from 1 to 0. If the growth term is ignored, the above Eq. (15) resembles the JMAK-type formulation of Eq. (12), where the term τ = V_m/(P_nV_n) characterizes the time constant of the nucleation process (e.g. the time needed for the amorphous volume to change the amount of a monomer V_m).

$ {\left( {\frac{{\rm{d}{V_{\rm{a}}}}}{{\rm{d}t}}} \right)_{\rm{c}}} = - \frac{{{V_{\rm{a}}}}}{\tau } . $  (16)

In principle, the total growth frontier is the summation of that of the scattered nucleus and that of the amorphous/crystalline regional interface. For simplicity, only the later area enters into the growth term of Eq. (15). The advantage of using the rate equation is to capture the essential temperature dependence in the crystallization process.

At the same time, the amorphization process is also modeled with a rate equation. Considering that the temperature difference between the PCM local temperature and the melting temperature is the driving force for the amorphization (melting then quenching), the rate equation is given by:

$ {\left( {\frac{{\rm{d}{V_{\rm{a}}}}}{{\rm{d}t}}} \right)_{\rm{a}}} = \frac{{{T_{\rm{a}}} - {T_{\rm{m}}}}}{{{R_{\rm{t}}}\Delta {h_{\rm{1}}}}} , $  (17)

in which Δh₁ is the latent heat of solid-to-liquid transition and T_a is the temperature at the amorphous/crystalline interface.

Assuming an initial state of full crystallization (V_a= 0), the reset current raises the BE/PCM interface temperature T_b, which is also T_a at the moment, above T_m. Then the RHS drives the melting region interface towards the TE. With a larger V_a hence a new melting/crystalline interface, T_a evaluated at the interface decreases (roughly T_a = (T_b– T_ref) ×(1 – V_a/V_amax) + T_ref ) with a linear temperature profile), also the melting rate decreases. Assuming the peak temperature T_b is constant, the melting region frontier moves until T_a is equal to T_m, which eventually determines the size of the amorphous region. Hence, in the rate formulation the term P_a = (T_a–T_m)/R_t can be regarded as the frontier local melting power, with Δh₁ the heat of fusion (unit of J/cm³). As a result, the formation time of melting volume is Δh₁/P_a. The rate equation is also made self-limiting by a theta function, as the maximum V_a is fixed for a technology. The two rate equations, Eqs. (15) and (17), are combined as a unified rate equation Eq. (18) for the phase dynamics. Eq. (18) is evaluated in a time evolution manner. Together with the other equations, the phase change memory operations including read, set and reset, are modeled.

$ \frac{{\rm{d}{V_{\rm{a}}}}}{{\rm{d}t}} = {\left( {\frac{{\rm{d}{V_{\rm{a}}}}}{{\rm{d}t}}} \right)_{\rm{c}}}\theta \left( {{T_{\rm{m}}} - {T_{\rm{a}}}} \right) + {\left( {\frac{{\rm{d}{V_{\rm{a}}}}}{{\rm{d}t}}} \right)_{\rm{a}}}\theta \left( {{T_{\rm{a}}} - {T_{\rm{m}}}} \right) . $  (18)

Compared to previous models, Sonoda’s model formulates the threshold switching and memory switching in the form of rate-equations. It also considers the different temperature dependences in the nucleation and growth process. With these, the model marks a reference for later rate-equation based models.

Later, another rate-equation based model was reported^[53]. Three variables, the crystalline fraction F_c, the amorphization fraction F_a, and the melting fraction F_m, are used to represent the state of the PCM cell, with F_m the new variable compared to the previous models. The melting fraction increases and may come closer to 1 in a dynamic process depending on the temperature above the melting T_m. When the temperature goes below T_m, the rate equation also describes the process of quenching, i.e. F_m decreases and comes closer to zero. It is noted that the fractions above are not corresponding to the PCM cell physics state. For the three-dimensional active region of the mushroom type, the fractions are in fact calculated assuming a one-dimensional profile. The resistance calculations are also done in a one-dimensional way without going to the details of current distributions in the 3D structure.

The rate equation for F_m is given by:

$ {\tau _{\rm{m}}}\frac{{\partial {F_{\rm{m}}}}}{{\partial t}} + {F_{\rm{m}}} = {\left[ {1 + \exp \left( {\frac{{{T_{\rm{m}}} - T}}{{{\sigma _{\rm{m}}}}}} \right)} \right]^{ - 1}} , $  (19)

in which τ_m is melting time constant. While σ_m is treated as a fitting parameter. When the temperature T goes beyond T_m, Eq. (19) gives the transition of F_m from zero to non-zero depending on σ_m. Thus, the parameter σ_m is also used to control the abruptness of the reset resistance of the R–I curve. A larger σ_m decreases the melting fraction F_m for a given temperature T, which then becomes the amorphous fraction F_a after the reset.

The rate equation for F_c follows the theory of JMAK. As the melting fraction is introduced, the differential version of Eq. (20) is revised a bit and goes as follows:

$ {\tau _{\rm c}}\frac{{\partial {F_{\rm c}}}}{{\partial t}} + {F_{\rm c}} = 1 - {F_{\rm m}} , $  (20)

in which τ_c is the characteristic time of crystallization.

Lastly, the fraction of amorphization F_a is the remaining part besides of the crystalline and melting fractions.

$ {F_{\rm a}} = 1 - {F_{\rm m}} - {F_{\rm c}} . $  (21)

After considering the non-Arrhenius crystallization, Eq. (20) is revised to the following with a fitting parameter b:

$ {\tau _{\rm{set}}}\frac{{\partial {F_{\rm c}}}}{{\partial t}} = {F_{\rm a}} b \, \exp \left( {1 - b {F_{\rm a}}} \right) , $  (22)

in whichτ_set is the temperature-dependent crystallization time. It is used to compensate the simplification of using the highest temperature for the crystallization rate. As the F_a is decreasing, RHS of Eq. (22) first increases then decreases with a peak at a certain F_a. Further, to account for the non-Arrhenius crystallization the set time τ_set is revised by including to activation energy:

$ {\tau _{\rm{set}}} = {\tau _{\rm{HT}}} + {\tau _{\rm{LT}}} = {\tau _{0\rm{HT}}}\exp \left( {\frac{{{E_{\rm{AHT}}}}}{{kT}}} \right) + {\tau _{0\rm{LT}}}\exp \left( {\frac{{{E_{\rm{ALT}}}}}{{kT}}} \right), $  (23)

in whichτ_LT and τ_HT are crystallization time for low temperature and high temperature, τ_0LT and τ_0HT are crystallization time prefactor for low temperature and high temperature, E_ALT and E_AHT are activation energy for low temperature and high temperature.

At the same time, caution is needed to handle the three-variable system as none of F_c, F_a and F_m can be negative, and reduction of F_c in the reset process is not reflected with Eq. (20) or (22). In Ref. [53], the current is modeled with a Poole–Frenkel conduction by assuming a complete voltage drop across the amorphous part.

$ {I_{\rm{PF}}} = {A_{\rm{kPF}}} F \left( { - \frac{{{\Phi _{\rm{PF}}} - {\beta _{\rm{PF}}}\sqrt F }}{{kT}}} \right) , $  (24)

in whichA_kPF,Φ_kPF and β_kPF are fitting parameters.

And the barrier dependence on temperature is given by:

$ {\Phi _{\rm{PF}}} = {E_{\rm{a}0}} - \frac{{{a_{\rm{va}}}{T^2}}}{{{b_{\rm{va}}} + T}} , $  (25)

in which E_a0 is the activation energy at 0 K, considered as a fitting parameter, and b_va are the thermal parameters.

By increasing the voltage, the current first increases with only the barrier lowering term as the heat is not significant. After the current goes above a certain level, the temperature rising comes into the picture and the transport barrier is reduced upon the heat. As a result, the current increases sharply with the voltage, or the voltage required to maintain the current is reduced for a snapback behavior. Fig. 13(a) plots the sharp switching from the amorphous state to crystalline state. With the core model equations Eqs. (19)–(25), the model captures theR–I characteristics as shown in Fig. 13(b).

Compared to previous models, Pigot’s model implements a physics-based threshold switching module, and a simplified crystallization module with rate-equations. The model incorporates related parameters for its flexibility, and has been verified by the industry’s PCM data in STMicroelectronics.

A model incorporating the phase transitions between three possible states was reported^[54]. The key rate equation used for the phase transitions is Eq. (26), with c the capture rate and e the emission rate. Similar to Figs. 5(a) and 5(b)^[40], the normalized ratios of three phase states are given first following the Arrhenius law. The transition rates including those from amorphous to crystalline, from crystallization to melting, and from melting to amorphous are explicitly given. The rates in the opposite directions are calculated from the equilibrium cases. Based on these rates, the dynamic transition equations form the governing equation set for the phase transition:

$ \frac{{\rm{d}\Omega_i}}{{\rm{d}t}} = \sum\limits_{i \ne j} {\sum\limits_{e,c \in P} {\left( {{c_{j \to i}}{\Omega _j} - {e_{i \to j}}{\Omega _i}} \right)} } , $  (26)

in which Ω_i is the normalized ratio of phase state i.

The resistance module incorporates an energy gain theory^[35] based on the Poole–Frenkel transport. Two trap states, a shallow trap state n_T2 and a deep trap state n_T1, are considered with possible re-populations between them. With a low field, the two states are generally conductive separately. With increasing the field, electrons in the deep trap state are possibly transmitted to the shallow trap states with tunneling, forming the non-equilibrium carriers n_T2’. The key equations of the current are summarized:

$\begin{array}{l} I = 2eA \dfrac{{\Delta z}}{{{\tau _0}}}\left[ {n_{\rm{T}1}}\exp \left( { - \dfrac{{{E_{\rm c}} - {E_{\rm{T}1}}}}{{kT}}} \right) + \left( {{n_{\rm{T}2}} + {n_{\rm{T}2}}'} \right)\exp \left( { - \dfrac{{{E_{\rm c}} - {E_{\rm{T}2}}}}{{kT}}} \right) \right] \\ \qquad\times\, {\rm{sinh}} \left( {\dfrac{{dU}}{{kT}}} \right) , \end{array}$  (27)

$ {H_{\rm{a,OFF}}} = \frac{{{E_{\rm{T}2}} - {E_{\rm{T}1}}}}{{kT}} , $  (28)

$ {H_{\rm{a,ON}}} = {H_{\rm{gst,a}}} - {H_{\rm{a,OFF}}} , $  (29)

$ V = {F_{\rm{OFF}}} {H_{\rm{a,OFF}}} + {F_{\rm{ON}}} {H_{\rm{a,ON}}} , $  (30)

in which E_T1 is deep trap energy and E_T2 is shallow trap energy, and Δz is trap distance.

The increase of the electron population in the shallow trap states increase the conduction current significantly. In the process, positive feedback is indeed happening. An increase of the current in the non-equilibrium region increases the electric field in the equilibrium region, which then increases the electron tunneling rate into the shallow trap states, and leads to a further increase of the current. This positive feedback is observed as a switching from the high resistance amorphous state to the dynamic low resistance state. In the model implementations, a current-controlled voltage source is used with simulation convergence considered. For the Joule-heating, the thermal resistance of different regions including those at the interfaces are all considered. Fig. 14(a) shows the schematic of the developed model including all the modules: the transport module, the thermal module and phase transition module, and Fig. 14(b) plots the comparisons between the models and experimental data.

Compared to previous models, Xu’s model implements a physics-based threshold switching, as well as a dynamic phase-transition model of Fig. 5. The model has been used for cross-point memory array simulations in Samsung.

Towards a fully-functional PCM model for circuit designs especially neuromorphic circuits, there are still some challenges ahead. From our perspective, these challenges include at least three parts: modeling for scaled PCM, modeling for variations^[55], and modeling for reliability.

The models in Section 3 are mainly developed for the classical mushroom PCM cells. Along with the technology scaling, alternative memory cell structures are developed such as the confined one^{[5, 15, 16]} with high integration density. The thermal boundary has changed significantly, for example, the temperature profile peaks at completely different locations in the cross-point PCM^[56]. The crystallization speed also changes due to absence of crystal seeds when a full reset of the confined region is triggered. These changes call for extensions of the current models or new generations of PCM models. It is in some sense similar to the scenarios of CMOS modeling which evolves with device geometry changes from planar to FinFET and gate-all-around structures. A PCM model that covers different cell geometries is highly desirable.

Relaxation oscillations in PCMs have been observed^[57-60]. The voltage across a PCM cell shows an oscillation, which features a current dependent frequency. Furthermore, the peak voltages are being relaxed in the process. A conductive filament theory^[61] was proposed to describe the oscillation in chalcogenide OTS devices, which are nowadays used as selectors^[62-65]. Although PCM is based on totally different chalcogenide glass which could undergo rapid crystallization, the relaxed oscillation may be explained with a similar reason. Later, the voltage/current dependent frequency was explained with a simple RC circuit assuming a holding current^[57]. A field-induced nucleation theory was proposed^[59] to describe the threshold voltage and holding voltages, as functions of material and geometry parameters. To utilize the oscillation properties of PCM in circuits for example the neuromorphic circuits, a further extension of the PCM model^[66] is to be developed by incorporating the related theories.

For a complete compact model, its parameter extractions are also needed to represent the technology with a process design kit. Different from the MOSFET model in which the parameters are extracted based on the DC and frequency domain properties, the PCM model requires an extraction also from the time domain properties such as the crystallization dynamics. Thermal-related parameters like the thermal resistance and thermal boundary resistance in a possible thermal network, are also to be extracted. An extraction strategy together with algorithms should be figured out.

A statistical model in an evolutionary manner for PCM conductance in the SET process has been developed^[67]. The model is intended to capture the statistical characteristics of PCM crystallization. Statistical properties are experimentally characterized: a pulse train is applied to a group of PCM cells with their initial conductance set as a high resistance state 0.1 μS. The first observation is that the increment of conductance after a programming pulse is decreased, as seen in Fig. 15(a). The same trend has been captured well in a physics-based model, as shown in Fig. 15(b)^[68]. It is understood from the JMAK equation, suggesting that the crystallization rate decreases with the accumulation of crystal fraction. At the same time, experimental characterizations also reveal a significant variation of the conductance among the PCM cells. Fig. 15(a) shows, for example, after the first partial SET pulse of some PCM cells does not further crystallize though the average conductance seems following the physics-based model prediction. It is speculated that this originates from the randomness of the crystallization process, e.g. the nucleation randomness as explained in Section 2. The statistical nature of PCM programming^[69], sometimes termed as a cycle-to-cycle variation, is not covered in previous compact models.

The model developed by Nandakumar et al.^[67] is summarized as Eqs. (31)–(34). The subscript N denotes the Nth programming pulse. The incremental conductance is a sampling of a Gaussian distribution, with the distribution parameters, the average and standard deviation, given by Eqs. (31) and (32), which are also evolving with the conductance itself. While the other parameters are for the purpose of fitting, χ represents a Gaussian random number of mean zero and variance 1, the variable P_mem is used to capture the programming history effect. The model reproduces very well the variations of measured PCM conductance as shown in Fig. 15(a).

$ {\mu _{\Delta {G_N}}} = {m_1}{G_{N - 1}}\left( {{t_0}} \right) + \left( {{c_1} + {A_1}{P_{\rm{mem}}}} \right) , $  (31)

$ {\sigma _{\Delta {G_N}}} = {m_2}{G_{N - 1}}\left( {{t_0}} \right) + \left( {{c_2} + {A_2}{P_{\rm{mem}}}} \right) , $  (32)

$ \Delta {G_N} = {\mu _{\Delta {G_N}}} + {\sigma _{\Delta {G_N}}}\chi , $  (33)

$ {G_N}\left( {{t_0}} \right) = {G_{N - 1}}\left( {{t_0}} \right) + \Delta {G_N} . $  (34)

While this model has been used for a two-PCM based spiking neural network, it is not directly applicable to CMOS-PCM integrated circuit simulations. Examples of the CMOS-PCM integrated neuromorphic computing circuits are reported^[19] for unsupervised SNN learning with the spiking time-dependent plasticity (STDP). In the example, the spiking time difference between the pre-synaptic and post-synaptic neurons are translated into the voltage pulse amplitude which induces the different amount of conductance change. A statistical transistor model, together with a statistical PCM model, will favor the design and simulations of neuromorphic circuits in a similar way as CMOS circuits. The transistor model has been in use for long but such a statistical PCM model is not yet available. At the same time, it is interesting to note that circuit level innovations are also being proposed to mitigate the variations, for example, by the multi-memristive synapse^[70]. It has been shown that the SNN classification performance is largely improved with the developed approach. However, variations cannot be fully eliminated. A further extension to incorporate the variation and stochastic crystallization is essential for the circuit evaluations, based on a physics-based device-level model of current–voltage characteristics such as those in Section 3. That will be one possible next step for the PCM compact modeling.

PCM conductance drifts after programming, especially after the RESET. Fig. 16(a) shows that the high resistance increases with time, while the low resistance does not change much. Physics-based models^[71] and statistical models^[72] have been developed to describe the resistance drift behaviors, in the single cell and PCM array levels. It is attributed to the structural relaxation after quenching, where the further rearrangement of local atomic structure leads to change of local states with the long-range disorder. A power law dependence of the resistance is empirically used with t₁ as reference time t of PCM programming^[71]:

$ R\left( t \right) = {R_1}{\left( {\frac{t}{{{t_1}}}} \right)^v} , $  (35)

in which R₁ is the resistance value at t = t₁ and v is the power-law exponent. The physics meaning of the exponential factor v is derived^[71], which is related to the increase of action energies for structure relaxation and conduction with time. The theory also explains the resistance dependent v factor observed in Fig. 16(a). Experimental characterizations in Ref. [72] further verifies the empirical expression, and a read noise is also included in the statistical model.

In additional to the resistance increase due to structural relaxation, resistance decrease^[74] due to spontaneous crystallization is another concern of PCM reliability. For memory array application, a long-term data retention at 85 °C is required. Fig. 16(b) shows the high-resistance retention under accelerated test conditions. For a binary or multi-level cell, it is usually a no-or-yes issue, with the retention failure or not. For the neuromorphic circuits, whether and to how much extent the resistance shifts changes circuit performances require further considerations. Other observed reliability issues include cycle endurance, thermal disturbance, and switching threshold voltage drifts. Cycle endurance describes the capacity to re-write in a PCM cell without failure^[75]. The PCM cell may not be programmed into the high resistance from the SET state, or may remain a high resistance without successful SET. Possible mechanisms are the change of intrinsic chalcogenide material properties after undergoing high temperature^[76]. The thermal disturbance is spontaneous crystallization due to heat from neighbor cells’ operations^[77]. Finally, threshold voltage shift also happens^[78], which is expected to be correlated to the resistance shifts. As stated earlier, analog memory property is used to represent the synapse connectivity in neuromorphic circuits. The resistance/threshold voltage shifts in the time domain, i.e. the change in synapse weight, in principle shifts the circuit characteristics like accuracy.

Innovative device level designs such as a projected PCM^[79] are also being excitingly explored to suppress the resistance shifts. The active phase change region is projected to non-phase-change materials for the high resistance state. With it, the reliability of the PCM-based cross-bar network is greatly enhanced. Based on the working principle, the high-resistance state eventually degrades when the resistance shifts in the amorphous region increase to an extent altering the current bypass. Moreover, projected PCM has been implemented in a single-layer neural network^[80], which could keep a high classification accuracy at elevated temperatures without temperature compensated way. Generally, for reliability-aware neuromorphic circuit designs reliability modeling of PCM cells^[81] deserves further effort.

Phase change memory has been on the stage of emerging memory technologies for some time, and recently promises applications in neuromorphic computing. In this work, we review the developments of PCM compact models, and try to provide an analysis of the challenges from our perspective. Generally, developing a compact model based on the device physics is a trend, especially when the device variations and reliability are important considerations. We start from the basic operation principles of PCM, and describe the known theoretical framework of crystallization (including nucleation and growth) and amorphization. Then we explain the evolution of the PCM models from macro models to more physics-oriented models, with different approximations in the model development. While the state-of-the-art PCM models provide a basis for circuit simulations, future model developments are facing challenges from the statistical modeling, to design-for-reliability models to facilitate the PCM applications in neuromorphic circuits.

This work is supported in part by the National Natural Science Foundation of China (62074006, 91964204), in part by the Major Scientific Instruments and Equipment Development (61927901), the Shenzhen Science and Technology Project (GXWD20201231165807007-20200827114656001), Strategic Priority Research Program of the Chinese Academy of Sciences (XDB44010200), Science and Technology Council of Shanghai (19JC1416801), the Shanghai Research and Innovation Functional Program (17DZ2260900), and in part by the 111 Project (B18001).