Optimization of Bilayer Resistive Random Access Memory Based on Ti/HfO2/ZrO2/Pt

In this paper, the electrothermal coupling model of metal oxide resistive random access memory (RRAM) is analyzed by using a 2D axisymmetrical structure in COMSOL Multiphysics simulation software. The RRAM structure is a Ti/HfO2/ZrO2/Pt bilayer structure, and the SET and RESET processes of Ti/HfO2/ZrO2/Pt are verified and analyzed. It is found that the width and thickness of CF1 (the conductive filament of the HfO2 layer), CF2 (the conductive filament of the ZrO2 layer), and resistive dielectric layers affect the electrical performance of the device. Under the condition of the width ratio of conductive filament to transition layer (6:14) and the thickness ratio of HfO2 to ZrO2 (7.5:7.5), Ti/HfO2/ZrO2/Pt has stable high and low resistance states. On this basis, the comparison of three commonly used RRAM metal top electrode materials (Ti, Pt, and Al) shows that the resistance switching ratio of the Ti electrode is the highest at about 11.67. Finally, combining the optimal conductive filament size and the optimal top electrode material, the I-V hysteresis loop was obtained, and the switching ratio Roff/Ron = 10.46 was calculated. Therefore, in this paper, a perfect RRAM model is established, the resistance mechanism is explained and analyzed, and the optimal geometrical size and electrode material for the hysteresis characteristics of the Ti/HfO2/ZrO2/Pt structure are found.


Introduction
Resistive random access memory (RRAM) is a new type of memory device, which uses the characteristics of transition metal oxide (TMO) film material of a resistive layer to change its internal electronic structure under the action of electric field and converts back and forth between high resistance and low resistance to store information [1,2].Many research groups have extensively studied ZnO, HfO 2 , ZrO 2 , TiO 2 , CeO 2, Ta 2 O 5 , Nb 2 O 5 , and other TMOs for RRAM applications [3][4][5][6][7][8][9][10][11][12][13][14].Among them, ZrO 2 and HfO 2 are considered to be the two most important materials due to their compatibility with complementary metal oxide semiconductors, high specific heat capacity, thermal conductivity, and dielectric constant [15].Arun et al. [16] systematically studied the characteristics of resistance switches based on monolayer HfO 2 using four different metal base electrode materials: Au, Al, Pt, and Cu.However, it was found that the resistance value of monolayer structure is not controllable, and there is no multi-value regulation effect.Moreover, Lee et al. [17] reported that compared with HfO x monolayer devices, the ZrO x /HfO x bilayer has a lower reset current, tight distribution of switching parameters, good switching durability up to 10 5 cycles, and good data retention at 85 • C.After this report, Ismail et al. [18] researched an RRAM with a TaN/HfO 2 /ZrO 2 /Pt structure, where a bilayer HfO 2 /ZrO 2 thin film structure was deposited by radio frequency sputtering at room temperature (RT) to investigate the resistive switching characteristics and the switching ratio > 10 for 10 3 DC cycles was calculated.Therefore, in this paper, HfO 2 and ZrO 2 are used as the transition layers of the bilayer RRAM structure to study the resistance mechanism behind it.
In order to study the resistance mechanism of RRAM devices, many previous RRAM models with different characteristics and accuracies have been proposed [19][20][21][22][23][24][25].Menenzel et al. [25] focused on the effect of the electric field on ion migration, considering the thermal properties of the upper and lower electrode materials: Ti/Nb:STO (Nb-doped SrTiO 3 ) of RRAM devices.They estimated the temperature characteristics of conductive filaments based on experimental data and established an electrothermal 2D finite element model.In the model, the effect of temperature and electric field on the change in SET time with the change in voltage is obtained by simulating the change in the resistivity in the filament region, but the bipolar I-V conversion characteristics of the device were not obtained.Kim et al. [26] modeled and simulated a bilayer Pt/Ta 2 O 5 /TaO x /W RRAM device.The model based on temperature and electric field accelerated oxygen vacancy migration, considering the effect of oxygen vacancy generation rate on the resistance characteristics of the device, obtained the DC bipolar I-V characteristics of the device.However, the recombination of oxygen vacancies and oxygen ions is neglected.Moreover, the transient properties of the Joule thermal model are ignored, and the temperature of the resistive lattice does not change with time.Huang et al. [27] proposed a physical analysis model of metal oxide RRAM under both DC and pulse operating modes.The model covers the transport of oxygen vacancies and oxygen ions, the temperature distribution under heat conduction, and also studies the influence of the width of conductive filaments and the gap between conductive filaments and electrodes, but the thickness of the transition layer is not considered.
Therefore, based on the model of Huang et al., this paper refers to a bipolar RRAM numerical model based on temperature field accelerated ion migration proposed by Larentis et al. [23] and takes into account the generation and recombination of local oxygen vacancy in the medium layer, and then establishes a perfect Ti/HfO 2 /ZrO 2 /Pt electricthermal coupling model.On the Ti/HfO 2 /ZrO 2 /Pt model, the scanning voltage rate is set to 100 mV/s by applying a cycle of forward and reverse scanning voltage.The distribution of oxygen vacancy concentration (n D ), temperature (T), and electric potential (Ψ) in Ti/HfO 2 /ZrO 2 /Pt RRAM with an applied voltage device is discussed in detail.The two-dimensional distribution diagram and one-dimensional curve diagram are used to show the distribution details of the three physical quantities.Thus, the RESET and SET process and mechanism of RRAM under the Ti/HfO 2 /ZrO 2 /Pt double structure are verified.In this paper, not only are the width of the bilayer conductive filament studied but also the thickness of the transition layer and the influence of the commonly used metal top electrode materials (Ti, Pt, and Al) on the hysteresis characteristics of the device.By comparison, it is found that under the conditions of the width ratio of the conductive filament to the transition layer (6:14) and the thickness ratio of HfO 2 to ZrO 2 (7.5:7.5),Ti is selected as the top electrode material.In this case, the RRAM has a stable high-low resistance state, that is, the maximum switching ratio and obvious hysteresis characteristics.Thus, the optimization of RRAM under the Ti/HfO 2 /ZrO 2 /Pt structure is completed by changing the electrode structure of RRAM and the size between the switching layer and the conductive filament in the simulation work.

Simulation Methods
In this paper, the three differential equations of oxygen vacancy concentration (n D ), temperature (T), and electric potential (Ψ) are simulated by using custom coefficient-type partial differential equations.The following are the physical fields and partial differential equations involved in the simulation model of RRAM devices.

Electrical-Thermal Coupling Physical Model
In order to study the change in oxygen ion concentration during the change in resistance value, oxygen ion migration is represented by the migration equivalent of oxygen vacancy.The movement of oxygen vacancies mainly includes three aspects: the diffusion movement of oxygen vacancies, the drift movement under the action of voltage, and the formation and recombination of local oxygen vacancies in the dielectric layer.The migration of oxygen vacancy is expressed by the flow rate of oxygen vacancy per unit volume, including drift and diffusion, so the diffusion and drift flux can be used to represent the migration flux of oxygen vacancy.When the voltage is 0 V, the drift motion does not need to be considered, and the oxygen vacancy will diffuse from the place with high concentration to the place with low concentration, and the diffusion motion will dominate at this time.Once the voltage is applied, the atomic transition barrier will decrease so that the oxygen vacancy can obtain enough energy to jump out of the original position.The reduction in the potential barrier is related to the applied potential, and the proportional relationship is qbE, where q is the amount of charge, b is the minimum size of the partition lattice, which is also the size of the grid in the simulation, and E is the electric field strength.The atomic transition probability is exponentially related to the barrier level, and the oxygen vacancy will jump rapidly after the voltage is applied, thus increasing the drift current, and the drift motion is dominant at this time.Based on the above analysis, the expression of oxygen vacancy migration is as follows [23]: where j diff and j dift represent diffusion and drift flux, respectively, the unit is cm 2 s −1 , n D is the oxygen vacancy concentration, the unit is cm −3 , D is the oxygen vacancy diffusion coefficient, v is the oxygen vacancy drift rate, and the unit is cm 2 s −1 .Therefore, the rate of change in oxygen vacancy per unit volume with time, namely the drift-diffusion equation of oxygen vacancy, can be expressed as follows [23]: and because the diffusion coefficient of ions is related to temperature, it satisfies the Arrhenius equation [23]: D 0 is the exponential factor of the diffusion coefficient, E a is the oxygen vacancy diffusion activation energy, k is the Boltzmann constant, T is the temperature, a is the ion transition distance in nm, and f is the escape frequency.
The drift rate of oxygen vacancy also obeys the Temperature-dependent Arrhenius equation [28]: The reason for the generation and recombination of local oxygen vacancies is that the concentration of oxygen vacancies in the medium layer is uneven, resulting in excessive electric field intensity at the local position after the voltage is applied, which breaks the covalent bond of HfO 2 molecules and forms HfO x and free oxygen ions [24].Oxygen ions can combine with each other to form oxygen molecules and free to the trap at the interface between the oxide and the metal, and HfO x will introduce new oxygen vacancies so that the concentration of oxygen vacancies in the local area becomes higher, that is, the effect.On the contrary, the recombination of oxygen vacancies will occur [29].The expression of oxygen vacancy generation and recombination is as follows [27]: where A is the constant related to vibration frequency, E b is the average activation energy of oxygen vacancy, and E c is the relaxation energy in the recombination process.Finally, as mentioned above, Equations ( 3)-( 7) are brought into the drift-diffusion Equation (2) for calculation, and the values of the parameters involved in COMSOL are shown in Table 1.Relaxation energy during recombination 1 The values of E a and a refer to ref. [29]. 2 The values of f refer to ref. [30]. 3The values of A, b, and E c refer to ref. [31].

Electrical Conduction Model
Current continuity equation [31]: where σ is the conductivity of the conducting filament, expressed in Ω −1 cm −1 , and Ψ is the potential.In order to solve (2)-( 8), we need to give a model of the conductivity σ, which can be expressed as follows: where σ 0 is the pre-exponential factor, the unit is Ω −1 cm −1 ; E AC is the conductive activation energy.Suppose that σ 0 and E AC both depend on the oxygen vacancy concentration.σ 0 is the equation for oxygen vacancy concentration; the specific setting is explained in the following text.

Joule Heat Model
The heat conduction equation is as follows [31]: where ρ is the density, C p is the heat capacity, k th is the thermal conductivity of conductive filament, and its unit is Wm −1 K −1 .It is also assumed that the change in k th depends on the oxygen vacancy concentration.The steady-state mode is exclusively considered, disregarding ρC P ∂T ∂t .By applying Fourier's law, the heat flux formula can be derived as follows: J T − k th • ∇T.The left side of the equation is to find its divergence.Because there is a passive field in the medium layer, the steady-state heat conduction equation is obtained as follows: The COMSOL simulation calculation process is the solution process of the partial differential Equations ( 2), ( 8) and (11), so as to calculate the distribution of oxygen vacancy concentration n D , electric potential Ψ, and temperature T during the device conversion process.

Model Establishment in COMSOL Geometric Model
In this paper, the 2D axisymmetric structure of COMSOL is used to create a geometric model of RRAM, as shown in Figure 1, where the left boundary (r = 0) is the axis of rotational symmetry.The device width, thickness, and top electrode to adjust material size parameters are shown in Table 2.

(
) The COMSOL simulation calculation process is the solution process of the partial differential Equations ( 2), ( 8) and (11), so as to calculate the distribution of oxygen vacancy concentration nD, electric potential Ψ, and temperature T during the device conversion process.

Geometric Model
In this paper, the 2D axisymmetric structure of COMSOL is used to create a geometric model of RRAM, as shown in Figure 1, where the left boundary (r = 0) is the axis of rotational symmetry.The device width, thickness, and top electrode to adjust material size parameters are shown in Table 2.In the model, the same width and thickness of the top and bottom electrodes were set according to the nanoscale dimensions in ref. [32] where the width is twice greater than described in the literature, and the thickness is half that.The width and thickness of CF1 and HfO2, CF2, and ZrO2 were initially set, but these parameters will be changed in our simulation work.The sum of the thickness of the Y2O3 switching layer and Y metal layer in ref. 32 is taken as the thickness of the HfO2 layer and ZrO2.

Definition of Material Properties
According to the parameters required in the analysis of the electro-thermal coupling model, the material properties of each layer entity used in the simulation of  In the model, the same width and thickness of the top and bottom electrodes were set according to the nanoscale dimensions in ref. [32] where the width is twice greater than described in the literature, and the thickness is half that.The width and thickness of CF1 and HfO 2 , CF2, and ZrO 2 were initially set, but these parameters will be changed in our simulation work.The sum of the thickness of the Y 2 O 3 switching layer and Y metal layer in ref. [32] is taken as the thickness of the HfO 2 layer and ZrO 2 .

Definition of Material Properties
According to the parameters required in the analysis of the electro-thermal coupling model, the material properties of each layer entity used in the simulation of Ti/HfO 2 /ZrO 2 /Pt RRAM devices are given in Table 3.In the simulation, it is assumed that the heat capacity C p , density ρ, and relative permittivity of CF1 and CF2 are equal to those of HfO 2 and ZrO 2 , and their conductivity σ and thermal conductivity k th are defined as a function of oxygen vacancy concentration n D .
According to the calculation formula of atomic density, The molar mass of HfO 2 was 210.49g/mol, ZrO 2 was 123.218 g/mol, and N A was Avogadro's constant 6.02 × 10 23 /mol.The atomic density P 1 of HfO 2 = 2.77 × 10 22 cm −3 and P 2 of ZrO 2 = 2.88 × 10 22 cm −3 can be obtained by Equation (12).Since the maximum doping concentration can reach 4.3% of the atomic density [23], the maximum oxygen vacancy concentration of HfO 2 in the simulation is 1.2 × 10 21 cm −3 , and the maximum oxygen vacancy concentration of ZrO 2 is 1 × 10 21 cm −3 .According to Equation ( 9), the conductivity σ can be determined by defining σ 0 in relation to E AC and oxygen vacancy concentration n D .The σ 0 of HfO 2 is estimated from 10 Ω −1 cm −1 to 3300 Ω −1 cm −1 , respectively, by high and low resistance resistivity.The σ 0 of ZrO 2 varies from 10 Ω −1 cm −1 to 3000 Ω −1 cm −1 .Figure 2a shows the relationship between σ0 and oxygen vacancy concentration.The E AC of HfO 2 activation energy in high and low resistance states are 0.087 eV and 0.018 eV, respectively.And the E AC of ZrO 2 is 0.087 eV and 0.018 eV, respectively.Figure 2b shows the relationship between E AC and oxygen vacancy concentration.When the oxygen vacancy concentration reaches half of the maximum concentration, conductive filaments are formed, and the device is converted to a low resistance state [27].In the simulation, the thermal conductivity k th at the filament changes from the thermal conductivity of HfO 2 and ZrO 2 to that of Hf and Zr with the increase in oxygen vacancy concentration.Because the oxygen vacancy concentration n D is assumed to enhance thermal conductivity due to the free-carrier contribution to heat conduction.As shown in Figure 2c, based on the Wiedemann-Franz law, we assumed linear dependence of k th on n D [22,23].In addition, the above processes are considered to be linear changes.

Results and Discussion
RESET is a process where the device shifts from a low resistance state to a high resistance state.When the COMSOL simulation model is established, the initial conditions of the device are that the conductive filament is fully connected; the boundary temperature and initial temperature are room temperature 300 K; the top electrode has a scanning voltage of 0 V→1 V→0 V, where the scanning second rate is 100 mV/s; the bottom

Results and Discussion
RESET is a process where the device shifts from a low resistance state to a high resistance state.When the COMSOL simulation model is established, the initial conditions of the device are that the conductive filament is fully connected; the boundary temperature and initial temperature are room temperature 300 K; the top electrode has a scanning voltage of 0 V→1 V→0 V, where the scanning second rate is 100 mV/s; the bottom electrode is grounded.In order to reflect the resistance process, CF1 and CF2 regions are considered, and the distribution diagram of oxygen vacancy concentration, potential, temperature, and other parameters of the device is obtained under the forward bias voltage of 0.4 V-0.8 V (step by 0.1 V).
When the bias voltage is 0.4 V, the conductive wire morphology in Figure 3a shows no significant change, and CF1 is still in a state of conduction.When the voltage is increased to 0.5 V, the oxygen vacancy concentration in CF1 changes, but no fracture occurs at this time.It can be seen from Figure 3b that the lowest oxygen vacancy concentration in the CF1 layer is about 1.05 × 10 21 cm −3 (>0.6 × 10 21 cm −3 ).Due to the change in oxygen vacancy concentration, CF1 begins to gradually break down.The voltage continues to increase to 0.6 V, and it can be intuitively seen from Figure 3b that CF1 has a clear fracture zone (the region with a blue concentration of 0 in CF1).When CF1 continues to increase the voltage to 0.8 V after fracture, Figure 3d,e show that the fracture region further expands.We simulated the distribution of oxygen vacancy concentration in the conductive filament region of RRAM at a forward bias voltage of 0.4 V-0.8 V (step by 0.1 V) in the z direction, as shown in Figure S1.Combined with the one-dimensional distribution curve in Figure S1, it can be found that the conductive wire morphology clearly breaks at z = 28 nm, the maximum fracture width is about 2 nm, and the V RESET voltage is 0.6 V.
In order to further study the temperature distribution, it can be seen from Figure 4a,b that when the applied voltage is 0.4 V and 0.5 V, the temperature distribution is uniform and mainly distributed within CF1. Figure 4c-e show that as the voltage increases, the temperature is no longer uniformly distributed.We further calculate the distribution of temperature and electric potential in the z direction when the forward bias is 0.4 V-0.8 V (step by 0.1 V), as shown in Figure S2.The specific distribution law can be clearly seen from the one-dimensional distribution curve in Figure S2a: When the temperature increases caused by Joule heat, the one-dimensional distribution curve of temperature is different from the parabola in the initial stage, and it is no longer smooth and distorted.Specifically, the temperature near the ZrO 2 interface is low, and the temperature near z = 28 nm is high.In the region where CF1 is broken, the temperature can reach up to 950 K.The results show that in this region, the migration rate of oxygen vacancy is fast, and CF1 will break at the highest temperature.In order to further study the temperature distribution, it can be seen from Figure 4a,b that when the applied voltage is 0.4 V and 0.5 V, the temperature distribution is uniform and mainly distributed within CF1. Figure 4c-e show that as the voltage increases, the temperature is no longer uniformly distributed.We further calculate the distribution of temperature and electric potential in the z direction when the forward bias is 0.4 V-0.8 V (step by 0.1 V), as shown in Figure S2.The specific distribution law can be clearly seen from the one-dimensional distribution curve in Figure S2a: When the temperature increases caused by Joule heat, the one-dimensional distribution curve of temperature is different from the parabola in the initial stage, and it is no longer smooth and distorted.Specifically, the temperature near the ZrO2 interface is low, and the temperature near z = 28 nm is high.In the region where CF1 is broken, the temperature can reach up to 950 K.The results show that in this region, the migration rate of oxygen vacancy is fast, and CF1 will break at the highest temperature.Then, the distribution of potential in the RESET process is studied.Figure 5 shows  In order to further study the temperature distribution, it can be seen from Figure 4a,b that when the applied voltage is 0.4 V and 0.5 V, the temperature distribution is uniform and mainly distributed within CF1. Figure 4c-e show that as the voltage increases, the temperature is no longer uniformly distributed.We further calculate the distribution of temperature and electric potential in the z direction when the forward bias is 0.4 V-0.8 V (step by 0.1 V), as shown in Figure S2.The specific distribution law can be clearly seen from the one-dimensional distribution curve in Figure S2a: When the temperature increases caused by Joule heat, the one-dimensional distribution curve of temperature is different from the parabola in the initial stage, and it is no longer smooth and distorted.Specifically, the temperature near the ZrO2 interface is low, and the temperature near z = 28 nm is high.In the region where CF1 is broken, the temperature can reach up to 950 K.The results show that in this region, the migration rate of oxygen vacancy is fast, and CF1 will break at the highest temperature.Then, the distribution of potential in the RESET process is studied.Figure 5 shows the distribution of potential in the dielectric layer under different forward bias voltages.As can be seen from Figure 5a,b, when the forward voltage is 0.4 V and 0.5 V and CF1 is Then, the distribution of potential in the RESET process is studied.Figure 5 shows the distribution of potential in the dielectric layer under different forward bias voltages.As can be seen from Figure 5a,b, when the forward voltage is 0.4 V and 0.5 V and CF1 is not broken, the oxygen vacancy concentration in the medium layer is the same, so the potential distribution is uniform, and the maximum potential is at the top of the upper half of the region.Figure 5c,e show that after CF1 breaks, the inner potential of the upper half of the region increases, indicating that the resistance here increases and the partial voltage increases.Figure S2b shows that the potential is uniformly distributed when CF1 is switched on.When the fracture is complete, the electric potential in the fracture area increases sharply, and the resistance value becomes large, showing a high resistance state.
In summary, in the RESET process, combined with the distribution of oxygen vacancy concentration, temperature, and potential, the whole process of CF1 from conduction to fracture can directly or indirectly indicate the fracture partitions the top electrode and the bottom electrode regions, forming a high resistance state.The device completes the RESET process from the low resistance state to the high resistance state.
not broken, the oxygen vacancy concentration in the medium layer is the same, so the potential distribution is uniform, and the maximum potential is at the top of the upper half of the region.Figure 5c,e show that after CF1 breaks, the inner potential of the upper half of the region increases, indicating that the resistance here increases and the partial voltage increases.Figure S2b shows that the potential is uniformly distributed when CF1 is switched on.When the fracture is complete, the electric potential in the fracture area increases sharply, and the resistance value becomes large, showing a high resistance state.In summary, in the RESET process, combined with the distribution of oxygen vacancy concentration, temperature, and potential, the whole process of CF1 from conduction to fracture can directly or indirectly indicate the fracture partitions the top electrode and the bottom electrode regions, forming a high resistance state.The device completes the RESET process from the low resistance state to the high resistance state.
We further verified that not only the RESET process but also the SET process exists in the Ti/HfO2/ZrO2/Pt structure.In this paper, the bias voltage is set to the reverse value, and the bias voltage of the top electrode is set to 0 V→−1 V→0 V, where the sweep second rate is still 100 mV/s, and the oxygen vacancy concentration distribution under the reverse bias voltage is obtained, as shown in Figure 6. Figure 6a,b show that the device is still in the high resistance state of the previous process, CF1 is in the broken state, and the oxygen vacancy concentration in the broken region is low.With the increase in voltage, the concentration of oxygen vacancy in the fracture gap of the filaments increases, and then CF1 recombines in the fracture zone and re-conducts.As shown in Figure 6c, the conductive wire breaks again when the voltage is increased to −0.7 V. We further verified that not only the RESET process but also the SET process exists in the Ti/HfO 2 /ZrO 2 /Pt structure.In this paper, the bias voltage is set to the reverse value, and the bias voltage of the top electrode is set to 0 V→−1 V→0 V, where the sweep second rate is still 100 mV/s, and the oxygen vacancy concentration distribution under the reverse bias voltage is obtained, as shown in Figure 6. Figure 6a,b show that the device is still in the high resistance state of the previous process, CF1 is in the broken state, and the oxygen vacancy concentration in the broken region is low.With the increase in voltage, the concentration of oxygen vacancy in the fracture gap of the filaments increases, and then CF1 recombines in the fracture zone and re-conducts.As shown in Figure 6c, the conductive wire breaks again when the voltage is increased to −0.7 V.The temperature distribution in the SET process of RRAM is shown in Figure 7.In the initial stage, the temperature distribution during CF1 fracture is mainly concentrated at the fracture of the red high-temperature region in Figure 7a,b.The current at the frac--0.85V -0.8 V -0.7 V -0.6 V -0.5 V -0.7 V -0.6 V -0.5 V The temperature distribution in the SET process of RRAM is shown in Figure 7.In the initial stage, the temperature distribution during CF1 fracture is mainly concentrated at the fracture of the red high-temperature region in Figure 7a,b.The current at the fracture is small, so less Joule heat is generated.With the increase in voltage, the current through z = 28 nm increases, and the oxygen vacancy accumulates toward the fracture, so Joule heat is also generated.Then, the high-temperature area at CF1 moves toward the ZrO 2 layer, and the temperature is gradually evenly distributed throughout the device.As shown in Figure S3, after calculating the 2D distribution of oxygen vacancy concentration, we simulated the distribution of temperature and potential in the z direction when applying the reverse bias voltage of −0.85 V, −0.8 V, −0.7 V, −0.6 V, and −0.5 V. Figure S3a is the one-dimensional temperature distribution curve of the device during the SET process.Combined with Figure S3a, it can be determined that the specific location of the fracture region is in CF1, where the temperature is low near z = 28 nm in the z direction, while in CF1, the temperature is high at z = 24 nm, and the migration rate of oxygen vacancy is accelerated.The temperature distribution in the SET process of RRAM is shown in Figure 7.In the initial stage, the temperature distribution during CF1 fracture is mainly concentrated at the fracture of the red high-temperature region in Figure 7a,b.The current at the fracture is small, so less Joule heat is generated.With the increase in voltage, the current through z = 28 nm increases, and the oxygen vacancy accumulates toward the fracture, so Joule heat is also generated.Then, the high-temperature area at CF1 moves toward the ZrO2 layer, and the temperature is gradually evenly distributed throughout the device.As shown in Figure S3, after calculating the 2D distribution of oxygen vacancy concentration, we simulated the distribution of temperature and potential in the z direction when applying the reverse bias voltage of −0.85 V, −0.8 V, −0.7 V, −0.6 V, and −0.5 V. Figure S3a is the one-dimensional temperature distribution curve of the device during the SET process.Combined with Figure S3a, it can be determined that the specific location of the fracture region is in CF1, where the temperature is low near z = 28 nm in the z direction, while in CF1, the temperature is high at z = 24 nm, and the migration rate of oxygen vacancy is accelerated.After the above discussion of the temperature distribution of RRAM devices, further attention should be paid to the distribution of potential.Figures 8a,b and S3b are two dimensional and one-dimensional distribution curves of potential distribution.Figure 8a,b show that the electric potential is mainly concentrated at the red high temperature near the top electrode Ti, showing a high resistance state.Thereafter, Figure 8a-c show that with the increase in voltage, the potential in the lower half region of CF1 is no longer evenly distributed but increases in the inner potential, indicating that the resistance value of the resistance in the lower half region becomes larger.Figure S3b can directly reflect the position change in the CF1 fault region: the abrupt increase in the fault region changes from the previous steep increase region at z = 28 nm to z = 24 nm.The uniform distribution of potential in the previous fracture region indicates that CF1 has a process from fracture to recombination fracture.
In summary, the distribution of oxygen vacancy concentration, temperature, and potential under reverse bias voltage can all indicate that CF1 breaks during RESET to re-conduction during SET, and the device converts from the above high resistance state to a low resistance state, completing the SET process.
value of the resistance in the lower half region becomes larger.Figure S3b can directly reflect the position change in the CF1 fault region: the abrupt increase in the fault region changes from the previous steep increase region at z = 28 nm to z = 24 nm.The uniform distribution of in the previous fracture region indicates that CF1 has a process from fracture to recombination fracture.In summary, the distribution of oxygen vacancy concentration, temperature, and potential under reverse bias voltage can all indicate that CF1 breaks during RESET to reconduction during SET, and the device converts from the above high resistance state to a low resistance state, completing the SET process.
After verifying the existence of REST and SET resistive processes, this paper further studies the effect of conductive filament size on RRAM.CF1 and CF2 in the 2D geometric model of the device are set to different widths and thicknesses; that is, the radius of the conductive filaments is set to 5 nm, 6 nm, 7 nm, and 8 nm, and the thickness takes the thickness of the HfO2 layer as an example (the total thickness of HfO2 layer and ZrO2 layer is 15 nm) and is set to 6.5 nm, 7 nm, 7.5 nm, and 8 nm, respectively, for forward and reverse scanning of the device.The results are shown in Figure 9a,b.The radius of CF1 and CF2 reflects the limiting current of the SET process and the resistance of the low resistance state.Figure 9a shows that with the increase in the radius of CF1 and CF2, the RESET conversion current increases, and the limiting current and low-resistance resistance of the SET process also increase.Therefore, it can be concluded that the larger the RESET conversion current, the larger the SET limiting current.However, in Figure 9a, the hysteresis curves under the reverse bias voltage almost coincide at 7 nm and 8 nm, which also indicates that the influence of excessive conversion current is the instability of the high-low resistance state and the reduction in the switching ratio.After verifying the existence of REST and SET resistive processes, this paper further studies the effect of conductive filament size on RRAM.CF1 and CF2 in the 2D geometric model of the device are set to different widths and thicknesses; that is, the radius of the conductive filaments is set to 5 nm, 6 nm, 7 nm, and 8 nm, and the thickness takes the thickness of the HfO 2 layer as an example (the total thickness of HfO 2 layer and ZrO 2 layer is 15 nm) and is set to 6.5 nm, 7 nm, 7.5 nm, and 8 nm, respectively, for forward and reverse scanning of the device.The results are shown in Figure 9a,b.The radius of CF1 and CF2 reflects the limiting current of the SET process and the resistance of the low resistance state.Figure 9a shows that with the increase in the radius of CF1 and CF2, the RESET conversion current increases, and the limiting current and low-resistance resistance of the SET process also increase.Therefore, it can be concluded that the larger the RESET conversion current, the larger the SET limiting current.However, in Figure 9a, the hysteresis curves under the reverse bias voltage almost coincide at 7 nm and 8 nm, which also indicates that the influence of excessive conversion current is the instability of the high-low resistance state and the reduction in the switching ratio.Figure 9b shows that with the increase in HfO 2 thickness, the RESET conversion current gradually decreases, but the hysteretic characteristics have no specific trend, and the most obvious high-low resistance state and the maximum switching ratio exist at 7.5 nm.Therefore, 7.5 nm is considered to be the best thickness for HfO 2 and ZrO 2 and has the most stable high and low resistance state at this thickness.shows that with the increase in HfO2 thickness, the RESET conversion current gradually decreases, but the hysteretic characteristics have no specific trend, and the most obvious high-low resistance state and the maximum switching ratio exist at 7.5 nm.Therefore, 7.5 nm is considered to be the best thickness for HfO2 and ZrO2 and has the most stable high and low resistance state at this thickness.In the model of RRAM devices, in addition to the influence of the size of CF1 and CF2 on its performance, the material selection of the top electrode also affects the I-V characteristics of the device.In this paper, the influence of commonly used metal top electrode materials on the resistance transition characteristics of the device is studied.As shown in Figure 10, the electrode combination needs to be changed into three structures

(b) (a)
Figure 9. Influence of size on RRAM device performance: (a) the width of the CF1 layer at 5 nm, 6 nm, 7 nm, and 8 nm; (b) the thickness of the CF1 layer at 6.5 nm, 7 nm, 7.5 nm, and 8 nm.
In the model of RRAM devices, in addition to the influence of the size of CF1 and CF2 on its performance, the material selection of the top electrode also affects the I-V characteristics of the device.In this paper, the influence of commonly used metal top electrode materials on the resistance transition characteristics of the device is studied.As shown in Figure 10, the electrode combination needs to be changed into three structures of Pt-Pt, Ti-Pt, and Al-Pt in the analysis process to draw I-V curves for comparison.The results show that the peak current of the Ti/HfO 2 /ZrO 2 /Pt device is about 0.25 mA in the SET process and 0.14 mA in the RESET process with Ti-Pt as an The switching ratio R off /R on of devices with electrode combinations of Pt-Pt, Ti-Pt, and Al-Pt are 1.97, 11.67, and 2.35, respectively.Ti-Pt has the largest switching ratio, so the Ti/HfO 2 /ZrO 2 /Pt structure has a more stable high-low resistance state in comparison.In the model of RRAM devices, in addition to the influence of the size of CF1 and CF2 on its performance, the material selection of the top electrode also affects the I-V characteristics of the device.In this paper, the influence of commonly used metal top electrode materials on the resistance transition characteristics of the device is studied.As shown in Figure 10, the electrode combination needs to be changed into three structures of Pt-Pt, Ti-Pt, and Al-Pt in the analysis process to draw I-V curves for comparison.The results show that the peak current of the Ti/HfO2/ZrO2/Pt device is about 0.25 mA in the SET process and 0.14 mA in the RESET process with Ti-Pt as an electrode.The switching ratio Roff/Ron of devices with electrode combinations of Pt-Pt, Ti-Pt, and Al-Pt are 1.97, 11.67, and 2.35, respectively.Ti-Pt has the largest switching ratio, so the Ti/HfO2/ZrO2/Pt structure has a more stable high-low resistance state in comparison.Based on the above discussion of the size of CF1 and the top electrode material, a bilayer RRAM with a Ti/HfO2 (thickness 7.5 nm)/ZrO2 (thickness 7.5 nm)/Pt structure is adopted in this paper, where the radius width of CF1 is 6 nm.The hysteresis curve, as shown in Figure 11, can be obtained.The forward voltage of 0 V→1 V→0 V is applied to Based on the above discussion of the size of CF1 and the top electrode material, a bilayer RRAM with a Ti/HfO 2 (thickness 7.5 nm)/ZrO 2 (thickness 7.5 nm)/Pt structure is adopted in this paper, where the radius width of CF1 is 6 nm.The hysteresis curve, as shown in Figure 11, can be obtained.The forward voltage of 0 V→1 V→0 V is applied to the device, and the V RESET voltage of RRAM is 0.6 V from conduction to disconnection.When the voltage is increased to 1 V, the RESET process of the device ends, completing the transition from a low resistance state to a high resistance state and then maintaining a high resistance state.Then, the reverse voltage is increased from 0 V to −1 V, and when the voltage is increased, the transformation process from the high resistance state to the low resistance state begins until the CF1 fracture region moves down to the ZrO 2 interface when the voltage increases to −0.7 V.The SET process from high resistance state to low resistance state is completed when the voltage is −0.7 V, that is, the V SET is −0.7 V. Combined with the curve in Figure 11, it can be calculated that the average resistance value of the high resistance state (R off ) of the device is 43.10 kΩ, the average resistance value of the low resistance (R on ) state is 4.12 kΩ, and the switching ratio is R off /R on = 10.46.Compared to the experimental data of Ismail et al.'s TaO 2 /ZrO 2 /HfO 2 /Pt bilayer structure [18] in Table 4, the simulated switching ratio, SET voltage, and RESET voltage for a similar structure are shown in the table.By comparison, it can be observed that in spite of choosing a smaller top electrode size (15 nm in this study compared to Ismail et al.'s 70 nm), the switch ratio differs by approximately 3.Not only that but our model also has lower SET and RESET voltages.

Conclusions
In this paper, firstly, a bilayer RRAM model of Ti/HfO 2 /ZrO 2 /Pt was established using COMSOL Multiphysics, and the RESET and SET resistance processes of RRAM, as well as the distribution of oxygen vacancy concentration (inside CF1 and CF2), temperature, and potential of the two processes were verified.Then, the effects of the size of the conductive filament and the top electrode material on the performance of the device were investigated.Under the conditions of the width ratio of the conductive filament to the transition layer (6:14) and the thickness ratio of HfO 2 to ZrO 2 (7.5:7.5), the RRAM has a stable high-low resistance state.In addition, the results also show that the RRAM's larger CF1 and CF2 size is not better.Moreover, a larger conversion current will appear at the same time that the switching ratio is lower, and the high and low resistance states are not obvious.Next, we set Pt-Pt, Ti-Pt, and Al-Pt as the electrodes and calculated that the switching ratios were 1.97, 11.67, and 2.35, respectively.So, the comparison showed that the switching ratio of the Ti electrode was the largest, which was 11.67.Finally, we found that the switching ratio R off /R on of Ti/HfO 2 /ZrO 2 /Pt RRAM is 10.46 by taking the optimal size and selecting the best top electrode material.And comparing our simulation results with the previous three HfO 2 /ZrO 2 bilayer RRAM [7,18,38], we can see that the smaller top electrode can also achieve an approximately equal switching ratio (about 3 less than Ismail et al.'s RRAM structure), which is more advantageous and valuable in low-power applications and integrated circuits.

Figure 1 .
Figure 1.Schematic diagram of the 2D axisymmetric geometric model of the device.In COMSOL, the Z-axis is used as the rotation axis, and the 2D model is rotated 360° to obtain the solid model.

Table 2 .
List of model size parameters in COMSOL[32].

Z r 1 Figure 1 .
Figure 1.Schematic diagram of the 2D axisymmetric geometric model of the device.In COMSOL, the Z-axis is used as the rotation axis, and the 2D model is rotated 360 • to obtain the solid model.

Table 2 .
List of model size parameters in COMSOL[32].

Figure 3 .
Figure 3. Two−dimensional distribution of oxygen vacancy concentration inside CF1 and CF2 of Ti/HfO2/ZrO2/Pt at the bias voltage of 0.4 V-0.8 V (step by 0.1 V) in the SET process.

Figure 4 .
Figure 4. Two−dimensional temperature distribution in the SET process of Ti/HfO2/ZrO2/Pt at bias voltage 0.4 V-0.8 V (step by 0.1 V).

Figure 3 .
Figure 3. Two−dimensional distribution of oxygen vacancy concentration inside CF1 and CF2 of Ti/HfO 2 /ZrO 2 /Pt at the bias voltage of 0.4 V-0.8 V (step by 0.1 V) in the SET process.

Figure 3 .
Figure 3. Two−dimensional distribution of oxygen vacancy concentration inside CF1 and CF2 of Ti/HfO2/ZrO2/Pt at the bias voltage of 0.4 V-0.8 V (step by 0.1 V) in the SET process.

Figure 4 .
Figure 4. Two−dimensional temperature distribution in the SET process of Ti/HfO2/ZrO2/Pt at bias voltage 0.4 V-0.8 V (step by 0.1 V).

Figure 4 .
Figure 4. Two−dimensional temperature distribution in the SET process of Ti/HfO 2 /ZrO 2 /Pt at bias voltage 0.4 V-0.8 V (step by 0.1 V).

Figure 5 .
Figure 5. Two−dimensional distribution of potential in the SET process of Ti/HfO2/ZrO2/Pt at bias voltage 0.4 V-0.8 V (step by 0.1 V).

Figure 5 .
Figure 5. Two−dimensional distribution of potential in the SET process of Ti/HfO 2 /ZrO 2 /Pt at bias voltage 0.4 V-0.8 V (step by 0.1 V).

Figure 9 .
Figure 9.Influence of size on RRAM device performance: (a) the width of the CF1 layer at 5 nm, 6 nm, 7 nm, and 8 nm; (b) the thickness of the CF1 layer at 6.5 nm, 7 nm, 7.5 nm, and 8 nm.

Figure 9 .
Figure9.Influence of size on RRAM device performance: (a) the width of the CF1 layer at 5 nm, 6 nm, 7 nm, and 8 nm; (b) the thickness of the CF1 layer at 6.5 nm, 7 nm, 7.5 nm, and 8 nm.

Table 1 .
Model parameters in COMSOL.

Table 3 .
Lists the material properties of entities at each layer of the model in COMSOL.(Both HfO 2 and ZrO 2 used in the model are polycrystalline structures.)The parameter values of the cited references are referred to the corresponding reference, and the values of unlabeled parameters are provided by the COMSOL material library.