# Numerical Study on Resistance Change Characteristics of Phase Change Materials

Authors: Yifang Su; Yihang Zhang; Kaifeng Lin; Guanjun Zhao; Qinzheng Yu; Zuoqi Hu; Xiongfei Tao
Journal: Advances in Materials Science and Engineering
 (2021)

---

## 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.


---

## Body

## 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 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.


## 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 (Figure1), under lower voltages, the resistivities of the materials in two states are clearly distinguished, and the crystalline GST exhibits linear Ohmic electrical characteristics.Figure 1
I-V characteristic curve ofGe2Sb2Te5.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 Tm 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 Tg and the melting point Tm. 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 Figure2, 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].Figure 2
Cross-sectional structure of a typical resistance change unit.The dimension parameters of the resistance change unit are shown in Table1, where hTEC, hBEC, and hGST are the thickness of the top, bottom electrode contact layers, and the phase change layer, respectively. SGST and SBEC 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.Table 1
Resistance change unit structural model parameters.
hTEC (nm)hGST (nm)SGST (nm2)hBEC (nm)SBEC (nm2)602001000×1000200250×250


### 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.Figure 3
Resistance change unit simulation calculation flow chart.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 intoNx×Ny×Nz small cubes, where Nx, Ny, and Nz 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.Figure 4
An example of 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.(1)∇σx,y,z∇Φx,y,z=0,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:(2)σx,y,z∇2Φx,y,z=0.A pulse with an amplitudeU 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:(3)Φx,y,z|top boundary=U,Φx,y,z|bottom boundary=0,∂Φx,y,z∂n|side boundary=0.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(4)Φi=Φj.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Φx,y,z inside the FEUs. With Φx,y,z, the electric field of each FEU can be calculated by the following equation:(5)Ex,y,z=−∇Φx,y,z,and the current density can then be derived as follows:(6)Jx,y,z=σx,y,zEx,y,z.


#### 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 Tx,y,z,t, where the temperature distribution satisfies the transient partial differential equation (PDE).(7)ρx,y,zCpx,y,z∂T∂t=kx,y,z∂2T∂x2+∂2T∂y2+∂2T∂z2+qv,where ρ is the material density, Cp is the specific heat at constant pressure, k is the thermal conductivity, t is time, T is the temperature, and qv is the internal heat source intensity, which relates to the current density and the resistivity ρelec of materials as follows:(8)qvx,y,z=J2x,y,zρelecx,y,z.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 (Tb), and the temperature of the top boundary is set to be room temperature (Tt) due to the high conductivity of top electrode materials (i.e., W):(9)Tx,y,z,t|top boundary=Tt,Tx,y,z,t|bottom boundary=Tb.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(10)k∂T∂n|side border=0.LetTi and Tj denote the temperature of the two adjacent material films, respectively, and qi and qj denote the heat conduction of the two adjacent material films, respectively, which is defined in the following equation:(11)Ti=Tj,qi=qj.Assuming that the initial temperature of the resistance change unit isTini which is generally 25°C, then the initial condition can be shown as follows:(12)Tx,y,z,t|t=0=Tini.


#### 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ΔGV be the difference of the amorphous and crystalline regions Gibbs free energy, γ 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.The energy barrier of forming a critical size nucleus is shown as follows:(13)ΔG∗=163γ3ΔGV2fθ.The relationship between nucleation rate and temperature can be expressed as follows:(14)Int=α1⋅N0⋅exp−Eα1KBT⋅exp−ΔG∗KBT,where α1 is a frequency, N0 is the number of atoms per unit volume of the phase change material, kB is the Boltzmann constant, and Eα1 is the nucleation activation energy.The Gibbs free energy,ΔGV, is related to the temperature, and it is determined by latent heat and melting point, which is shown as follows:(15)ΔGV=ΔHfTm−TTm7TTm+6T,where ΔHf is the enthalpy difference per unit volume between the crystalline and the liquid phase at the melting point and Tm 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 Vg after nucleation is shown as follows:(16)Vgt=α2⋅d⋅e−Eα2KBT⋅1−exp−ΔGVKBT.


#### 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 sectionS 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 S are denoted by FEUi; thus, the current from all the FEUs cut by S and will add up to the total current.(17)I=∑IFEUi=∑JFEUi×SFEUi,where JFEUi is the current density of FEUi and SFEUi is the sectional area of FEUi cut by S. According to Ohm’s law ℛ = U/ℐ, the resistance of the resistance change unit is obtained.Figure 5
Schematic diagram of calculating the total resistance of the resistance change unit.


## 2.1. The Structural Model of Resistance Change Unit

The model of the resistance change unit used in this paper is shown in Figure2, 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].Figure 2
Cross-sectional structure of a typical resistance change unit.The dimension parameters of the resistance change unit are shown in Table1, where hTEC, hBEC, and hGST are the thickness of the top, bottom electrode contact layers, and the phase change layer, respectively. SGST and SBEC 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.Table 1
Resistance change unit structural model parameters.
hTEC (nm)hGST (nm)SGST (nm2)hBEC (nm)SBEC (nm2)602001000×1000200250×250


## 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.Figure 3
Resistance change unit simulation calculation flow chart.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 intoNx×Ny×Nz small cubes, where Nx, Ny, and Nz 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.Figure 4
An example of 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.(1)∇σx,y,z∇Φx,y,z=0,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:(2)σx,y,z∇2Φx,y,z=0.A pulse with an amplitudeU 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:(3)Φx,y,z|top boundary=U,Φx,y,z|bottom boundary=0,∂Φx,y,z∂n|side boundary=0.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(4)Φi=Φj.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Φx,y,z inside the FEUs. With Φx,y,z, the electric field of each FEU can be calculated by the following equation:(5)Ex,y,z=−∇Φx,y,z,and the current density can then be derived as follows:(6)Jx,y,z=σx,y,zEx,y,z.


### 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 Tx,y,z,t, where the temperature distribution satisfies the transient partial differential equation (PDE).(7)ρx,y,zCpx,y,z∂T∂t=kx,y,z∂2T∂x2+∂2T∂y2+∂2T∂z2+qv,where ρ is the material density, Cp is the specific heat at constant pressure, k is the thermal conductivity, t is time, T is the temperature, and qv is the internal heat source intensity, which relates to the current density and the resistivity ρelec of materials as follows:(8)qvx,y,z=J2x,y,zρelecx,y,z.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 (Tb), and the temperature of the top boundary is set to be room temperature (Tt) due to the high conductivity of top electrode materials (i.e., W):(9)Tx,y,z,t|top boundary=Tt,Tx,y,z,t|bottom boundary=Tb.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(10)k∂T∂n|side border=0.LetTi and Tj denote the temperature of the two adjacent material films, respectively, and qi and qj denote the heat conduction of the two adjacent material films, respectively, which is defined in the following equation:(11)Ti=Tj,qi=qj.Assuming that the initial temperature of the resistance change unit isTini which is generally 25°C, then the initial condition can be shown as follows:(12)Tx,y,z,t|t=0=Tini.


### 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ΔGV be the difference of the amorphous and crystalline regions Gibbs free energy, γ 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.The energy barrier of forming a critical size nucleus is shown as follows:(13)ΔG∗=163γ3ΔGV2fθ.The relationship between nucleation rate and temperature can be expressed as follows:(14)Int=α1⋅N0⋅exp−Eα1KBT⋅exp−ΔG∗KBT,where α1 is a frequency, N0 is the number of atoms per unit volume of the phase change material, kB is the Boltzmann constant, and Eα1 is the nucleation activation energy.The Gibbs free energy,ΔGV, is related to the temperature, and it is determined by latent heat and melting point, which is shown as follows:(15)ΔGV=ΔHfTm−TTm7TTm+6T,where ΔHf is the enthalpy difference per unit volume between the crystalline and the liquid phase at the melting point and Tm 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 Vg after nucleation is shown as follows:(16)Vgt=α2⋅d⋅e−Eα2KBT⋅1−exp−ΔGVKBT.


### 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 sectionS 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 S are denoted by FEUi; thus, the current from all the FEUs cut by S and will add up to the total current.(17)I=∑IFEUi=∑JFEUi×SFEUi,where JFEUi is the current density of FEUi and SFEUi is the sectional area of FEUi cut by S. According to Ohm’s law ℛ = U/ℐ, the resistance of the resistance change unit is obtained.Figure 5
Schematic diagram of calculating the total resistance of the resistance change unit.


## 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.(1)∇σx,y,z∇Φx,y,z=0,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:(2)σx,y,z∇2Φx,y,z=0.A pulse with an amplitudeU 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:(3)Φx,y,z|top boundary=U,Φx,y,z|bottom boundary=0,∂Φx,y,z∂n|side boundary=0.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(4)Φi=Φj.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Φx,y,z inside the FEUs. With Φx,y,z, the electric field of each FEU can be calculated by the following equation:(5)Ex,y,z=−∇Φx,y,z,and the current density can then be derived as follows:(6)Jx,y,z=σx,y,zEx,y,z.


## 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 Tx,y,z,t, where the temperature distribution satisfies the transient partial differential equation (PDE).(7)ρx,y,zCpx,y,z∂T∂t=kx,y,z∂2T∂x2+∂2T∂y2+∂2T∂z2+qv,where ρ is the material density, Cp is the specific heat at constant pressure, k is the thermal conductivity, t is time, T is the temperature, and qv is the internal heat source intensity, which relates to the current density and the resistivity ρelec of materials as follows:(8)qvx,y,z=J2x,y,zρelecx,y,z.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 (Tb), and the temperature of the top boundary is set to be room temperature (Tt) due to the high conductivity of top electrode materials (i.e., W):(9)Tx,y,z,t|top boundary=Tt,Tx,y,z,t|bottom boundary=Tb.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(10)k∂T∂n|side border=0.LetTi and Tj denote the temperature of the two adjacent material films, respectively, and qi and qj denote the heat conduction of the two adjacent material films, respectively, which is defined in the following equation:(11)Ti=Tj,qi=qj.Assuming that the initial temperature of the resistance change unit isTini which is generally 25°C, then the initial condition can be shown as follows:(12)Tx,y,z,t|t=0=Tini.


## 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ΔGV be the difference of the amorphous and crystalline regions Gibbs free energy, γ 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.The energy barrier of forming a critical size nucleus is shown as follows:(13)ΔG∗=163γ3ΔGV2fθ.The relationship between nucleation rate and temperature can be expressed as follows:(14)Int=α1⋅N0⋅exp−Eα1KBT⋅exp−ΔG∗KBT,where α1 is a frequency, N0 is the number of atoms per unit volume of the phase change material, kB is the Boltzmann constant, and Eα1 is the nucleation activation energy.The Gibbs free energy,ΔGV, is related to the temperature, and it is determined by latent heat and melting point, which is shown as follows:(15)ΔGV=ΔHfTm−TTm7TTm+6T,where ΔHf is the enthalpy difference per unit volume between the crystalline and the liquid phase at the melting point and Tm 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 Vg after nucleation is shown as follows:(16)Vgt=α2⋅d⋅e−Eα2KBT⋅1−exp−ΔGVKBT.


## 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 sectionS 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 S are denoted by FEUi; thus, the current from all the FEUs cut by S and will add up to the total current.(17)I=∑IFEUi=∑JFEUi×SFEUi,where JFEUi is the current density of FEUi and SFEUi is the sectional area of FEUi cut by S. According to Ohm’s law ℛ = U/ℐ, the resistance of the resistance change unit is obtained.Figure 5
Schematic diagram of calculating the total resistance of the resistance change unit.


## 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 Table2.Table 2
Pulse parameters.
SimulationNumber of pulsesAmplitude (V)Width (ns)Duty cycleGST initial state110.2∼4.430.375Crystalline211.010∼4000.7Amorphous31∼131.630.375Crystalline42∼150.8100.7Amorphous52∼200.6200.5Amorphous


### 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 Table2. 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 Figure6. 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 6
Pulse amplitude vs. resistance curve of simulation 1.Figure7 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).Figure 7
Phase distribution of simulation 1. The phase distribution of GST after imposing the pulses with amplitudes of (a) 1.6 V, (b) 1.8 V, (c) 2.0 V, (d) 2.4 V, (e) 2.8 V, (f) 3.2 V, (g) 3.8 V, and (h) 4.6 V on the resistance change unit, respectively.
(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 Table2. 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.Figure8 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.Figure 8
Pulse width vs. resistance curve of simulation 2.To explain the phenomenon, the phase distribution is plotted as shown in Figure9. 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.Figure 9
Phase distribution of simulation 2. The phase distribution of GST after imposing pulses with widths of (a) 10 ns, (b) 20 ns, (c) 30 ns, (d) 50 ns, (e) 70 ns, (f) 90 ns, (g) 110 ns, and (h) 150 ns on the resistance change unit, respectively.
(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 Table2. 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.Figure10 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.Figure 10
Pulse number vs. resistance curve of simulation 3.Figure 11
Phase distribution of simulation 3. (a–h) The phase distribution of GST after imposing 2∼9 (subsequently increasing by one) continuous pulses on the resistance change unit, respectively.
(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.Figure12 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.Figure 12
Pulse number vs. resistance curve of simulation 4.The phase distribution is illustrated in Figure13, 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 Tm, the temperature descends radially around the area, and the area where the temperature is between Tm and Tg 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.Figure 13
Phase distribution of simulation 4. (a–h) The phase distribution of GST after imposing 2, 3, 4, 6, 9, 10, 12, and 15 continuous pulses on the resistance change unit, respectively.
(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.Figure14 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.Figure 14
Pulse number vs. resistance curve of simulation 5.The phase distribution is present in Figure15, 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 Tm and Tg. As more pulses are imposed, the temperature of this area becomes higher than Tm, 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.Figure 15
Phase distribution of simulation 5. (a–h) The phase distribution of GST after imposing 1, 2, 4, 7, 9, 10, 14, and 18 continuous pulses on the resistance change unit, respectively.
(a)(b)(c)(d)(e)(f)(g)(h)Both Figures12 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.


## 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 Table2. 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 Figure6. 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 6
Pulse amplitude vs. resistance curve of simulation 1.Figure7 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).Figure 7
Phase distribution of simulation 1. The phase distribution of GST after imposing the pulses with amplitudes of (a) 1.6 V, (b) 1.8 V, (c) 2.0 V, (d) 2.4 V, (e) 2.8 V, (f) 3.2 V, (g) 3.8 V, and (h) 4.6 V on the resistance change unit, respectively.
(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 Table2. 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.Figure8 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.Figure 8
Pulse width vs. resistance curve of simulation 2.To explain the phenomenon, the phase distribution is plotted as shown in Figure9. 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.Figure 9
Phase distribution of simulation 2. The phase distribution of GST after imposing pulses with widths of (a) 10 ns, (b) 20 ns, (c) 30 ns, (d) 50 ns, (e) 70 ns, (f) 90 ns, (g) 110 ns, and (h) 150 ns on the resistance change unit, respectively.
(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 Table2. 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.Figure10 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.Figure 10
Pulse number vs. resistance curve of simulation 3.Figure 11
Phase distribution of simulation 3. (a–h) The phase distribution of GST after imposing 2∼9 (subsequently increasing by one) continuous pulses on the resistance change unit, respectively.
(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.Figure12 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.Figure 12
Pulse number vs. resistance curve of simulation 4.The phase distribution is illustrated in Figure13, 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 Tm, the temperature descends radially around the area, and the area where the temperature is between Tm and Tg 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.Figure 13
Phase distribution of simulation 4. (a–h) The phase distribution of GST after imposing 2, 3, 4, 6, 9, 10, 12, and 15 continuous pulses on the resistance change unit, respectively.
(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.Figure14 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.Figure 14
Pulse number vs. resistance curve of simulation 5.The phase distribution is present in Figure15, 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 Tm and Tg. As more pulses are imposed, the temperature of this area becomes higher than Tm, 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.Figure 15
Phase distribution of simulation 5. (a–h) The phase distribution of GST after imposing 1, 2, 4, 7, 9, 10, 14, and 18 continuous pulses on the resistance change unit, respectively.
(a)(b)(c)(d)(e)(f)(g)(h)Both Figures12 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 componentIn 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].Figure16(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.Figure 16
Fundamental diagram of the series resonance fault current limiter.
(a)(b)Before the fault, the switch is open and the seriesLC circuit carries the load current. If L and C satisfy the following equation, in which fs is the system frequency, then the series connection of L and C shows zero impedance:(18)2πfsL=12πfsC.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 capacitorC. 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 Figure16(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. This requirement will be met by a series-parallel connection of the small units. Let L=633mH 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 fs=50Hz (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.Figure 17
Curve of impedance vs. frequency.


## 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.



---

*Source: 7934173-2021-07-06.xml*