Abstract
Phase change materials are a type of emerging materials whose states will change under certain conditions, which then lead to changes in resistance. To study the characteristics of the phase change materials, a numerical simulation model of the resistive change unit based on the finite element method and the classic nucleation/growth theory is established, while the partial differential equations of electricity and heat conduction and the discrete formula of the finite element are also derived. According to the phase transition process of phase change materials, a crystalline-amorphous simulation model is also proposed in this paper to simulate the electrical and thermal properties and phase transition process of the resistive change unit. Simulations of the resistance change unit under single pulses with different amplitudes and widths as well as the simulations under continuous pulses are conducted in this paper. These results verify the characteristics of resistance change and can provide references for selecting the parameters of the resistance change units.
1. Introduction
Phase change materials (PCM) usually have two states, which are crystalline and amorphous. The structural differences between the two states result in electronic performances. Most of the research studies on phase change materials include structural device optimization, phase change material selection, reversible transformations between the crystalline and amorphous states, and the exploration of the connection between materials structure and performance. Among the current known phase change materials, the chalcogenide-based materials are the most mature, and the most popular are Ge-Sb-Te (GST) due to the high rate of crystallization. In 1968, Ovshinsky [1] presented the memory based on phase change theory for the first time. Due to the different optical and electronic characteristics during the transition between the crystalline and amorphous state, phase change materials can be used to store digital information; in other words, the crystalline and amorphous state can represent “1” and “0,” respectively, which is called the Ovshinsky effect.
GST materials have been deeply studied for their high crystallization rate (<100 ns), high ratio of high vs. low resistance, and good thermal stability [2]. The rapid phase change switching between the crystalline and amorphous states leads to the change of reflectance in PCM, where both states are stable at room temperature resulting in long data storage time (>10 years). This virtue allows PCM to be widely used in nonvolatile random-access memories (PCRAM) [3, 4].
GST materials have a higher melting point (600°C) requiring a higher amplitude operating current/voltage pulse to produce a phase change, which makes the working current/voltage of the resistance change unit using GST material larger. Among the many structures of the resistance change unit, the mushroom structure has a simple shape and is relatively easy to prepare, so it is widely studied. When certain electrical pulses are imposed on the resistance change units, the reversible transformation between a low-resistivity polycrystalline state and a high-resistivity amorphous state can be observed.
Even though PCMs are well studied, research on PCMs is still a hotspot, such as the study on their behavior under high pressure [5] and the study on the structural changes of GST upon rapid cooling by molecular dynamics simulations and atomistic cluster alignment analysis [6], and the research on their applications in other fields as well. For instance, Wainstein et al. [7] investigated RF switches based on phase change memory (PCM) and Xu et al. [8] discussed the use of PCMs to implement artificial neurons and synapses.
In this paper, we employed a methodology using numerical calculation to study the characteristics of phase change materials, where the numerical results should help to guide experimental research.
2. Numerical Simulation Model
The resistive change unit undergoes the reversible conversion phenomenon of the resistance as applied by electrical pulses. The physical mechanism behind this is the state changes between crystalline and amorphous of the phase change materials. The materials show low resistance in the crystalline state and high resistance in the amorphous state. Taking the I-V characteristics of the typical phase change material GST as an example (Figure 1), under lower voltages, the resistivities of the materials in two states are clearly distinguished, and the crystalline GST exhibits linear Ohmic electrical characteristics.

Phase change materials exhibit two ways to transition from crystalline to amorphous. One is ion implantation [9], and the other is melt-quenching. The latter is adopted in this paper, which applies a high but narrow electric pulse on the resistive change unit, where a large amount of Joule heat is generated in a short time, and the local temperature rises rapidly above the melting point of the materials. Due to the narrow pulse, the phase change materials do not have enough time to rearrange into a long-range ordered crystalline state and instead change into a disordered amorphous state. Compared with the narrow and intensive pulse mentioned above, if a relatively low but wide pulse is employed on the phase change materials, Joule heat is generated which induces a temperature increase in a range between the crystallization temperature and the melting point . Owing to the wider pulse, the amorphous phase change materials have sufficient time to nucleate, grow, and transform into a crystalline state. This time-varying resistance conversion is called Ovonic memory switching (OMS). In this paper, we have established a physical model and employed numerical calculations to simulate and analyze this phenomenon, and thus we study the resistance change characteristics of phase change materials.
2.1. The Structural Model of Resistance Change Unit
The model of the resistance change unit used in this paper is shown in Figure 2, which includes the bottom electrode (metal W), the bottom electrode contact (or heater TiN), the phase change layer (GST), the top electrode contact (TiN), and top electrode (W), where the material parameters of each layer are from literature [10, 11].

The dimension parameters of the resistance change unit are shown in Table 1, where , , and are the thickness of the top, bottom electrode contact layers, and the phase change layer, respectively. and denote the area of horizontal section of the phase change layer and the lower electrode contact layer, respectively. Those parameters of the model are variable and will directly affect the current density distribution of the resistive element, as well as the resistance, consequently.
2.2. Computing Theory
The process of phase change between the amorphous and crystalline state in a resistance change unit can be described by differential equations of electricity, heat, and phase change [12] and can be calculated using the finite element method. The overall calculation flow chart is shown in Figure 3.

The finite element method (FEM) is one of the most important numerically analytic engineering methods. The FEM divides the simulation object into finite element units (FEUs) whose vertices and connections points are called finite element nodes (FENs). All the FEUs and FENs together form a finite element mesh.
The PCM layers are divided into small cubes, where , , and are dimensions in the , , and directions, and three parameters depend on the demand of simulation; Figure 4 shows an example of a finite element mesh.

2.2.1. Numerical Calculation of Electrical Simulation
The principle of electrical simulation for the resistance change units is based on basic electrical Laplace equations [13, 14]. The pulse imposed on the resistance change unit will generate Joule heat, and this thermal energy will lead to the phase transition. Since the resistance change unit is composed of several materials with different resistivities, the current density throughout the resistance change unit is uneven. Once the pulse is imposed on the resistance change unit, the potential distribution inside the unit can be described by a Laplace equation.where , and are the internal coordinates of the resistance change unit and and are conductivity and potential, respectively. Because the resistance change unit is divided into a finite number of FEUs and the materials inside each FEU are considered to be homogeneous, equation (1) is simplified as follows:
A pulse with an amplitude is imposed on the top and bottom electrodes of a resistance change unit, and the sides of the GST are connected to the insulation materials; thus, the solution of the steady-state differential equation of equation (2) has the following boundary condition:
Let and denote the potential of the two adjacent material films, where the potential inside the phase change unit is continuous, so we have
The potentials of the vertices of each FEU in the resistance change unit can be obtained by solving the differential equation, and the interpolation is employed to further develop the potentials inside the FEUs. With , the electric field of each FEU can be calculated by the following equation:and the current density can then be derived as follows:
2.2.2. Numerical Calculation of Thermal Simulation
The thermal simulation for the resistance change units is based on basic thermal equations [12, 15]. The current flowing through the resistance change unit generates heat and the heat conducts through the materials, causing temperatures in the resistance change unit to vary from space to space, as well as from time to time, which is denoted by , where the temperature distribution satisfies the transient partial differential equation (PDE).where is the material density, is the specific heat at constant pressure, is the thermal conductivity, is time, is the temperature, and is the internal heat source intensity, which relates to the current density and the resistivity of materials as follows:
To solve the PDE, certain boundary and initial conditions are needed. Assume that the temperature of the bottom boundary is the maximum temperature beyond which the CMOS circuits do not properly work (), and the temperature of the top boundary is set to be room temperature () due to the high conductivity of top electrode materials (i.e., W):
The GST layer is surrounded by insulative material with low thermal conductivity. To simplify the analysis, assume that there is no heat conduction through the side boundaries, which yields
Let and denote the temperature of the two adjacent material films, respectively, and and denote the heat conduction of the two adjacent material films, respectively, which is defined in the following equation:
Assuming that the initial temperature of the resistance change unit is which is generally 25°C, then the initial condition can be shown as follows:
2.2.3. Numerical Calculation for the Simulation of the Phase Transition Process
Numerical calculation of the phase transition process of phase change materials is key to the simulation of the resistance change unit, where the phase transition process can be attributed to the crystallization of the phase change material [10, 11]. The calculation uses classical nucleation/growth theory to simulate the crystallization process of the resistance change unit [16].
The crystallization of the phase change material in an amorphous state requires a certain number of nuclei with sizes larger than the critical value, which is a nucleation process. After which, the disordered atoms and molecules begin to grow around the nuclei. Let be the difference of the amorphous and crystalline regions Gibbs free energy, be the amorphous-crystalline interface energy, be a function of the contact angle , while can distinguish homogeneous or heterogeneous nucleation processes in multiphase materials or on their interfaces.
The energy barrier of forming a critical size nucleus is shown as follows:
The relationship between nucleation rate and temperature can be expressed as follows:where is a frequency, is the number of atoms per unit volume of the phase change material, is the Boltzmann constant, and is the nucleation activation energy.
The Gibbs free energy, , is related to the temperature, and it is determined by latent heat and melting point, which is shown as follows:where is the enthalpy difference per unit volume between the crystalline and the liquid phase at the melting point and is the melting point.
According to the frequency , the atomic distance , and the activation energy of the amorphous to crystalline diffusion, the grain growth rate after nucleation is shown as follows:
2.2.4. Resistance Calculation of the Resistance Change Unit
By imposing a low voltage that does not provoke phase change on the phase change unit and selecting a section S as a reference plane in the resistance change unit in a direction perpendicular to the direction of the current flow, the current density in all FEUs can be obtained (see Figure 5). Then, all the FEUs that cut through are denoted by ; thus, the current from all the FEUs cut by and will add up to the total current.where is the current density of and is the sectional area of cut by . According to Ohm’s law = /, the resistance of the resistance change unit is obtained.

3. Simulation Results Analysis
The phase change of the materials is essentially the structural transformation under the excitation of different energies. The level of phase change can be controlled in a certain range by changing the excitation energy, where one or more intermediate states can be obtained consequently.
In the following section, several simulations will be conducted by imposing different pulses on the resistance change unit to verify the procedure of phase changes of GST, as well as the resistance change. The initial states of the GST materials in all the simulations are either fully crystalline or fully amorphous, where the pulses are all ideal square wave pulses, and all the parameters are listed in Table 2.
3.1. Analysis of the Single Pulse Simulation Experiment
The 1st simulation is conducted by imposing pulses with different amplitudes, and the parameters are in the first row of Table 2. The initial states of GST are all fully crystalline, pulse widths are the same, but the amplitudes of the pulses vary from 0.2 V to 4.4 V.
The relationship between the resistances and the pulse amplitude is shown in Figure 6. From the figure, it can be seen that if the amplitude of the pulse is less than 1.5 V, the resistance remains unchanged. If the amplitude continues to increase, the resistance exponentially increases, which strides over more than three orders of magnitude with the pulse amplitude varying from 1.6 V to 3.0 V. However, the resistance barely changes and remains high if the amplitude is more than 3.0 V.

Figure 7 shows the phase distribution inside the GST layer, where red represents the crystalline part and blue represents amorphous part, while the remaining colors indicate intermediate parts between the two states, which are illustrated by the color bar in the right of 7. It can be seen from the figure that under the pulse with an amplitude of 1.6 V, amorphous regions are generated in the lower part of GST layer, which implies that this part has the highest temperature in the unit. The Joule heat relates to the current density, and the bottom electrode is small, so the area around the bottom electrode is supposed to show the highest temperature. However, the bottom electrode is TiN with high thermal conductivity, and the heat is quickly released near the bottom electrode. So, the highest temperature part has a distance to the bottom electrode, as shown in Figure 7(a).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)
The increment of amplitude of the pulse leads to the increment of the excitation energy, which may result in the expansion of the amorphous area, as well as the gradual decrement of the gap between the amorphous area and bottom electrode. The resistance increases along with the increment of the amplitude of the pulse imposed on the resistance change unit, but as long as the amplitude reaches a certain level (about 3.0 V in this simulation), the expansion of the amorphous area contributes little to the increment of the resistance, and resistance of the unit remains stable.
The 2nd simulation is carried out by imposing single pulses with a fixed amplitude but different widths, and the parameters are in the second row of Table 2. The initial states of GST are all fully amorphous and pulse amplitudes are fixed to 1.0 V, while the widths of the pulses vary from 10 ns to 400 ns.
Figure 8 presents the resistances of the unit by imposing pulses of different widths. The resistance significantly decreases with the increment of the pulse width from 30 ns to 150 ns, and if the pulse width is greater than 150 ns, the downward resistance trend slows.

To explain the phenomenon, the phase distribution is plotted as shown in Figure 9. The Joule heat generated by current pulse accumulates at the center of the amorphous area, so this area stays at a high temperature which causes difficult crystallization. The temperature drops from the center point, as long as the temperature is suitable for crystallization, and the crystalline margin around the central amorphous area appears. With increasing pulse duration, the crystalline margin expands outward and the resistance decreases. Once the crystalline margin contacts the bottom electrode, the resistance begins to decrease.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)
3.2. Analysis of a Continuous Pulse Simulation Experiment
A 3rd simulation was conducted by imposing a continuous pulse on the resistance change unit, where the amplitude, duration, and duty cycle of all pulses were 1.5 V, 3 ns, and 37.5%, respectively, which are shown in Table 2. The initial states of the GST for the 3rd simulation are all fully crystalline. In the 1st simulation, a 1.5 V pulse causes the GST material to partially change phase, and if more pulses with the same amplitude are imposed on the unit, the energy will accumulate in GST causing more of the crystalline portion to transition into the amorphous state.
Figure 10 shows the resistance curve imposing a continuous pulse, which shows that the resistance of the unit increases with increasing numbers of pulses, and when the resistance reaches around 1 M, it remains stable. The increment of the resistance under multiple pulses is illustrated by Figure 11, which shows the distribution of the phase in the GST area. An appropriate single pulse allows the GST layer to be partially amorphous, and the amorphous area is near the bottom electrode due to the high thermal conductivity of the bottom electrode, as shown in Figure 11(a). As more pulses are imposed on the unit, the amorphous area gradually expands and narrows the conductive channel, and the resistance decreases. When the amorphous area is large enough to completely cover the bottom electrode, the resistance remains relatively stable.


(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)
A 4th simulation was conducted by imposing a continuous pulse on the resistance change unit, while the initial states of the GST for this simulation are all fully crystalline. The amplitude, duration, and duty cycle of all pulses are 0.8 V, 10 ns, and 70%, respectively. Comparing the parameters for the 3rd simulation, the pulses for this simulation are wider but lower, which benefits higher crystallization of the GST layer.
Figure 12 shows the resistances of the unit by imposing different numbers of pulses on it. As seen in the figure, the resistance decreases along increasing numbers of pulses, and when the number of the pulses exceeds 12, the rate of decrease begins to slow.

The phase distribution is illustrated in Figure 13, where the process of the change of phase distribution is similar to the 2nd simulation. It means that the Joule heat generated by the single pulse causes the middle area to show the highest temperature above , the temperature descends radially around the area, and the area where the temperature is between and begins to crystallize. Once the pulse is withdrawn, the area with temperature above its melting point remains amorphous. With increasing numbers of pulses, the crystalline margin expands outward and the resistance decreases.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)
The type of pulse for the 5th simulation is also continuous, and the initial states of the GST for this simulation are all also fully crystalline too. The amplitude, duration, and duty cycle of all pulses are 0.6 V, 20 ns, and 50%, respectively. The pulses for this simulation are wider but lower than the pulses employed for the 4th simulation.
Figure 14 shows the resistances of the unit by imposing different numbers of pulses on it, and the result is similar to Figure 12 as the resistance decreases with increasing numbers of pulses.

The phase distribution is present in Figure 15, but it is not exactly the same as the phase distribution of the 4th simulation. From the figure, it can be seen that the middle area of the GST layer is fully crystallized after the first 9, which means that the temperature of this area is between and . As more pulses are imposed, the temperature of this area becomes higher than , and once the pulse is withdrawn, the area becomes amorphous. However, the overall resistance decreases with increasing numbers of the pulse due to the expanded crystalline area.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)
Both Figures 12 and 14 show a floor between high and low resistance, and the potential reason behind this is that along with increasing crystalline pulses, the crystalline area begins to contact the bottom electrode, which is where the first resistance decrease occurs. However, due to the high thermal conductivity of the electrode, the crystalline area cannot cover the bottom electrode, and as more pulses are imposed on the unit, the crystalline area expands until it covers the whole electrode, which leads to the second drop in resistance.
4. Application of Resistance Change Unit
Chalcogenide-based phase change material is a well-acknowledged data storage medium; however, this material also has potential applications in other fields, such as power systems. As the power system becomes more complicated, the probability of short circuits increases. Fault current limiter (FCL) is an effective solution to control the current during faults within low levels. For an ideal FCL, it should meet these key characteristics [17] as follows:(i)Very low impedance during normal operation(ii)Zero power loss in normal operation(iii)Large impedance in fault conditions(iv)Quick activation when fault occurs(v)Fast recovery after fault removal(vi)No DC component
In recent years, many circuit topologies of FCLs have been proposed, such as superconducting FCLs (SFCLs), solid-state FCLs, and resonant FCLs (RFCLs). Although SFCLs show excellent performance in normal mode as well as during and after fault, the SFCLs need cooling systems to maintain a superconducting state, which are expensive and highly power consuming. This defect makes SFCLs impractical. Among FCLs, RFCLs as nonsuperconducting FCLs have drawn more attention owing to some exclusive merits. Rezaee and Harley [17] investigated different types of RFCLs, where a series RFCL was found to be effective [18, 19].
Figure 16(a) shows a fundamental circuit topology of SRFCL. There is an inductor connected with a capacitor in series, and there is a switch parallel with the capacitor in this circuit topology. Additionally, a series resistance must be added to represent the resistance of the inductor. There are many options for the switching device, where the most straightforward is a back-to-back thyristor switch.

(a)

(b)
Before the fault, the switch is open and the series circuit carries the load current. If and satisfy the following equation, in which is the system frequency, then the series connection of and shows zero impedance:
Once a short circuit occurs in the power system, the impedance rapidly decreases, the low impedance leads to overcurrent in the system, and then the switch closes which will short circuit the capacitor . Thus, the impedance increases and limits the fault currents within a low level.
In some application scenarios, the switch is replaced by another element; for instance, a resistance change unit can be used to act like a switch in RFCL (see Figure 16(b)). From the result of the simulation above, the difference between high and low impedance is almost four orders of magnitude, so assuming the low resistance is 200 and the high resistance is about 1 , to satisfy the applications of the power system, the capacities of current and voltage should be improved. This requirement will be met by a series-parallel connection of the small units. Let and , when the system is in normal mode of operation, so the resistance change unit has a high resistance and the series circuit shows zero impedance at a system frequency of (see Figure 17). However, during the fault, the capacitor voltage will rapidly increase, as shown in Figure 8, and the resistance of the resistance change unit in parallel with the capacitor will decrease to low levels due to the transition from the amorphous to crystalline phase. As shown in Figure 17, the overall impedance is about 200 because the capacitor is partially short circuited, and thus, the fault current will be reduced.

5. Conclusions
In this paper, the characteristics of the resistance change unit are studied by numerical simulation. Mathematical models of the resistance change unit were devised to characterize the electrical and thermal features by a serial of Laplace equations. If a pulse is imposed on the electrode of the unit, the current density inside the phase change material can be calculated though the mathematical models. Consequently, the Joule heat distribution, as well as the temperature distribution, can also be obtained. Then, according to the temperature distribution, the phase distribution of the phase change material in the unit can be calculated by the classic crystal nucleation-growth principle and melting-quenching process. Finally, a pulse with low amplitude, which does not provoke a change of state, is imposed on the unit, so the current distribution and total current can be obtained and the resistance can be calculated by Ohm’s law.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was funded by the Science and Technology Project of State Grid Zhejiang Electric Power Co., Ltd. (5211JH170008) and the National Natural Science Foundation of China (51872104).