Numerical Study on Resistance Change Characteristics of Phase Change Materials

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.*ese results verify the characteristics of resistance change and can provide references for selecting the parameters of the resistance change units.


Introduction
Phase change materials (PCM) usually have two states, which are crystalline and amorphous. e 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 chalcogenidebased 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]. e 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).
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 abinitio 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.

Numerical Simulation Model
e resistive change unit undergoes the reversible conversion phenomenon of the resistance as applied by electrical pulses.
e physical mechanism behind this is the state changes between crystalline and amorphous of the phase change materials. e 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. e 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 T m 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 T g and the melting point T m . Owing to the wider pulse, the amorphous phase change materials have sufficient time to nucleate, grow, and transform into a crystalline state. is 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.

e Structural Model of Resistance Change Unit.
e 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]. e dimension parameters of the resistance change unit are shown in Table 1, where h TEC , h BEC , and h GST are the thickness of the top, bottom electrode contact layers, and the phase change layer, respectively. S GST and S BEC denote the area of horizontal section of the phase change layer and the lower electrode contact layer, respectively. ose parameters of the model are variable and will directly affect the current density distribution of the resistive element, as well as the resistance, consequently.

Computing eory.
e 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. e overall calculation flow chart is shown in Figure 3. e finite element method (FEM) is one of the most important numerically analytic engineering methods. e 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.   Advances in Materials Science and Engineering e PCM layers are divided into N x × N y × N z small cubes, where N x , N y , and N z are dimensions in the x, y, and z directions, and three parameters depend on the demand of simulation; Figure 4 shows an example of a 6 × 4 × 3 finite element mesh.

Numerical Calculation of Electrical Simulation.
e principle of electrical simulation for the resistance change units is based on basic electrical Laplace equations [13,14]. e 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 x, y, and z 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 U 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 Φ i and Φ j denote the potential of the two adjacent material films, where the potential inside the phase change unit is continuous, so we have e 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 Φ(x, y, z) inside the FEUs. With Φ(x, y, z), the electric field of each FEU can be calculated by the following equation: and the current density can then be derived as follows:

Numerical Calculation of ermal
Simulation. e thermal simulation for the resistance change units is based on basic thermal equations [12,15]. e 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 T(x, y, z, t), where the temperature distribution satisfies the transient partial differential equation (PDE).
where ρ is the material density, C p is the specific heat at constant pressure, k is the thermal conductivity, t is time, T is the temperature, and q v is the internal heat source intensity, which relates to the current density and the resistivity ρ elec 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 (T b ), and the temperature of the top boundary is set to be room temperature (T t ) due to the high conductivity of top electrode materials (i.e., W): e 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 T i and T j denote the temperature of the two adjacent material films, respectively, and q i and q j 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 T ini which is generally 25°C, then the initial condition can be shown as follows:

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]. e calculation uses classical nucleation/growth theory to simulate the crystallization process of the resistance change unit [16]. e 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 ΔG V be the difference of the amorphous and crystalline regions Gibbs free energy, c be the amorphous-crystalline interface energy, f be a function of the contact angle θ, while f can distinguish homogeneous or heterogeneous nucleation processes in multiphase materials or on their interfaces. e energy barrier of forming a critical size nucleus is shown as follows: e relationship between nucleation rate and temperature can be expressed as follows: where α 1 is a frequency, N 0 is the number of atoms per unit volume of the phase change material, k B is the Boltzmann constant, and E α1 is the nucleation activation energy. e Gibbs free energy, ΔG V , is related to the temperature, and it is determined by latent heat and melting point, which is shown as follows: where ΔH f is the enthalpy difference per unit volume between the crystalline and the liquid phase at the melting point and T m is the melting point. According to the frequency α 2 , the atomic distance d, and the activation energy E α2 of the amorphous to crystalline diffusion, the grain growth rate V g after nucleation is shown as follows:

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). en, all the FEUs that cut through S are denoted by FEU i ; thus, the current from all the FEUs cut by S and will add up to the total current.
where J FEU i is the current density of FEU i and S FEU i is the sectional area of FEU i cut by S. According to Ohm's law R � U/I, the resistance of the resistance change unit is obtained.

Simulation Results Analysis
e phase change of the materials is essentially the structural transformation under the excitation of different energies. e 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. e 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.

Analysis of the Single Pulse Simulation Experiment.
e 1st simulation is conducted by imposing pulses with different amplitudes, and the parameters are in the first row of Table 2. e 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. e 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. e 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). e 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. e 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. e 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. e 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. e 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. e Joule heat generated by  Figure 5: Schematic diagram of calculating the total resistance of the resistance change unit. Advances in Materials Science and Engineering current pulse accumulates at the center of the amorphous area, so this area stays at a high temperature which causes difficult crystallization. e 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.

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. e 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. e 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). 6 Advances in Materials Science and Engineering 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 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. e 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. e 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 T m , the temperature descends radially around the area, and the area where the temperature is between T m and T g 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. e type of pulse for the 5th simulation is also continuous, and the initial states of the GST for this simulation  Advances in Materials Science and Engineering are all also fully crystalline too. e amplitude, duration, and duty cycle of all pulses are 0.6 V, 20 ns, and 50%, respectively. e 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. e 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 T m and T g . As more pulses are imposed, the temperature of this area becomes higher than T m , 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.
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.

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    [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. ere 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. ere are many options for the switching device, where the most straightforward is a back-to-back thyristor switch. Before the fault, the switch is open and the series LC circuit carries the load current. If L and C satisfy the following equation, in which f s is the system frequency, then the series connection of L and C 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 C. us, 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 MΩ, to satisfy the applications of the power system, the capacities of current and voltage should be improved. is requirement will be met by a series-parallel connection of the small units. Let L � 633 mH and C � 16 μF, when the system is in normal mode of operation, so the resistance change unit has a high resistance and the series LC circuit shows zero impedance at a system frequency of f s � 50 Hz (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.

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. en, 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
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare there are no conflicts of interest regarding the publication of this paper. Advances in Materials Science and Engineering 11