Simulation method for growth morphology evolution of cbram conductive filament based on thermoelectric coupling

Abstract

This invention provides a simulation method for the growth morphology evolution of conductive wires in CBRAM based on thermoelectric coupling. It employs a two-dimensional matrix to simulate the CBRAM device, and uses a sparse matrix solver to process the nodal voltage equations of the grid points based on a resistive network model to obtain the potential distribution inside the CBRAM device. The finite difference method is used to calculate the temperature field during the growth of the conductive wires inside the CBRAM device. Based on the potential distribution and temperature field, the rate of various reactions occurring for each particle in the CBRAM device is calculated, and the reaction occurring for each particle is determined based on the KMC method. After each particle performs a reaction, the two-dimensional matrix, local resistive network, potential, and temperature field of the CBRAM device are updated. A resistive switching event is triggered when a connection between the top and bottom electrodes of the CBRAM device is detected through a conductive wire. This invention improves the accuracy of conductive wire growth morphology evolution.

Description

Technical Field

[0001] This invention relates to the field of programmable resistive memory technology, and in particular to a method for simulating the growth morphology evolution of conductive filaments in CBRAM based on thermoelectric coupling. Background Technology

[0002] CBRAM (Containerless Circuit Resistive RAM) has emerged as a promising candidate for next-generation non-volatile memory due to its low power consumption, high speed, and excellent scalability. However, the performance of CBRAM devices is highly dependent on the growth and fracture behavior of the conductive filaments, a process influenced by a combination of factors including potential barriers, electric fields, and thermal effects. During the formation of the conductive filaments, the current flowing through the conductive channels generates localized Joule heating, leading to changes in ion migration rates and other parameters, thereby affecting the morphological evolution of the conductive filaments and the electrical performance of the device.

[0003] Existing methods for CBRAM conductive wire evolution generally employ electric field-driven ion migration models to simulate the CBRAM conductive wire growth mechanism. While these models can describe the formation and breakage processes of the conductive wires, they mostly focus only on electric field-driven ion migration and interfacial reactions, often neglecting or oversimplifying the crucial role of thermal effects and heat conduction in device operation. This makes it difficult to accurately simulate the dynamic evolution behavior of the conductive wires under actual operating conditions. This, to some extent, limits the accuracy of device performance optimization and reliability analysis. Summary of the Invention

[0004] This invention provides a method for simulating the growth morphology evolution of CBRAM conductive wires based on thermoelectric coupling, which addresses the shortcomings of existing technologies in accurately simulating the dynamic evolution behavior of conductive wires under actual working conditions. It enables accurate simulation of the impact of local heating on the growth and morphology evolution of conductive wires, thereby improving the accuracy of CBRAM device performance simulation.

[0005] This invention provides a method for simulating the growth morphology evolution of CBRAM conductive filaments based on thermoelectric coupling, comprising:

[0006] A two-dimensional matrix is ​​used to simulate the CBRAM device. Based on the resistor network model, a sparse matrix solver is used to process the nodal voltage equations of the grid points to obtain the potential distribution inside the CBRAM device.

[0007] The temperature field during the growth of the conductive filament inside the CBRAM device is calculated using the finite difference method. Based on the potential distribution and temperature field, the rate of various reactions of each particle in the CBRAM device is calculated. The reaction of each particle is determined by the rate based on the KMC method.

[0008] Each particle updates the two-dimensional matrix, local resistance network, potential, and temperature field of the CBRAM device after performing a reaction, and triggers a resistive switching event when a connection is detected between the top and bottom electrodes of the CBRAM device via a conductive filament.

[0009] According to the present invention, a simulation method for the growth morphology evolution of CBRAM conductive filaments based on thermoelectric coupling is provided, which uses the finite difference method to calculate the temperature field during the growth process of the conductive filaments inside the CBRAM device, including:

[0010] The thermal diffusivity of the CBRAM device is calculated based on the material thermal properties of the CBRAM device, and the stability factor of the CBRAM device is calculated based on the thermal diffusivity.

[0011] The local Joule heat source is obtained based on the cell current and the potential difference between adjacent grids of the CBRAM device.

[0012] The temperature field during the growth of the conductive wire is calculated based on the stability factor of the CBRAM device, the local Joule heat source, and the material thermophysical parameters.

[0013] According to the present invention, a method for simulating the growth morphology evolution of CBRAM conductive filaments based on thermoelectric coupling is used to calculate the thermal diffusivity of the CBRAM device based on the material thermal properties of the CBRAM device using the following formula:

[0014]

[0015] in, Where is the thermal diffusivity, Thermal conductivity, For material density, Specific heat capacity.

[0016] According to the present invention, a method for simulating the growth morphology evolution of CBRAM conductive filaments based on thermoelectric coupling is provided, and the stability factor of the CBRAM device is calculated based on the thermal diffusivity using the following formula:

[0017]

[0018] Where F is the stability factor, Where is the thermal diffusivity, For temperature evolution step size, The lattice spacing is F, which is less than the first preset threshold.

[0019] According to the present invention, a method for simulating the growth morphology evolution of CBRAM conductive filaments based on thermoelectric coupling is used to obtain the local Joule heat source based on the cell current and the potential difference between adjacent grids of the CBRAM device using the following formula:

[0020]

[0021]

[0022]

[0023] in, This indicates that once the conductive path is completed, the total current is evenly distributed among the visible metal atoms in each row according to their quantity. The obtained coordinates Unit current at that point The number of visible metal atoms in the j-th row. Indicates the potential difference between adjacent grid cells. Representing coordinates Voltage at that point Representing coordinates Voltage at that point Representing coordinates The local Joule heat source at that location, Let Q(i,j) be the temperature evolution step size. In the non-metallic region, the current is zero, and Q(i,j)=0.

[0024] According to the present invention, a method for simulating the growth morphology evolution of CBRAM conductive wires based on thermoelectric coupling is provided, with a temperature evolution step size. The calculation formula is:

[0025] = +( ) / total rate

[0026] in, It is a random number within the interval [0,1].

[0027] According to the present invention, a method for simulating the growth morphology evolution of CBRAM conductive wires based on thermoelectric coupling is provided. For non-bottom nodes inside the CBRAM device, the temperature field during the growth of the conductive wire is calculated using the following formula based on the stability factor of the CBRAM device, the local Joule heat source, and the material thermophysical parameters:

[0028]

[0029] For the bottom-level nodes inside the CBRAM device, the temperature field during the growth of the conductive filament is calculated using the following formula based on the stability factor of the CBRAM device, the local Joule heat source, and the material thermophysical parameters:

[0030]

[0031] in, and Representing the coordinates at time t+1 and t respectively. The temperature at that location Indicates the stabilizing factor. The unit conversion scaling factor.

[0032] According to the CBRAM conductive filament growth morphology evolution simulation method based on thermoelectric coupling provided by the present invention, after calculating the temperature field during the conductive filament growth process inside the CBRAM device using the finite difference algorithm, the method further includes:

[0033] If the temperature of any grid in the CBRAM device is greater than the second preset threshold and the grid is not located in the bottom electrode layer, then high-temperature induced material instability or atomic desorption occurs at that grid. The corresponding metal atoms and adsorbed ions are marked to zero to reflect structural damage or local melting behavior.

[0034] According to the present invention, a method for simulating the growth morphology evolution of CBRAM conductive filaments based on thermoelectric coupling is provided, which calculates the rate of various reactions occurring in each particle of the CBRAM device based on the potential distribution and temperature field using the following formula:

[0035]

[0036] Where Γ is the rate at which each particle undergoes various reactions, Ea is the activation energy required for each reaction to occur, and f is the vibrational frequency of the electron in the solid. Δφ is the correction factor, Δφ is the potential energy change before and after each reaction, k is the Boltzmann constant, q is the charge of a single electron, and T is the temperature.

[0037] The present invention provides a CBRAM conductive filament growth morphology evolution simulation method based on thermoelectric coupling. By introducing a local heat source calculation model and an explicit finite difference heat conduction solution module, it achieves a strong coupling simulation of electric field, thermal effect and conductive filament morphology evolution. This allows the temperature distribution to dynamically reflect the spatial accumulation and diffusion of Joule heat, and adjusts ion migration, interface reaction and material state update through temperature rise feedback. This enables accurate simulation of the influence of local heating on conductive filament growth and morphology evolution, improves the accuracy of CBRAM device performance simulation, and effectively makes up for the limitations of existing technologies in terms of insufficient consideration of thermal effects and lack of electro-thermal coupling. Attached Figure Description

[0038] To more clearly illustrate the technical solutions in this invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0039] Figure 1 This is a flowchart illustrating the simulation method for the growth morphology evolution of CBRAM conductive wires based on thermoelectric coupling provided by the present invention.

[0040] Figure 2 This is a schematic diagram of the simulation model of the CBRAM device in the CBRAM conductive filament growth morphology evolution simulation method based on thermoelectric coupling provided by the present invention.

[0041] Figure 3 This is a schematic diagram of the overall model simulation process of CBRAM device in the CBRAM conductive filament growth morphology evolution simulation method based on thermoelectric coupling provided by the present invention.

[0042] Figure 4 This is a schematic diagram of the CBRAM device resistor network model in the CBRAM conductive filament growth morphology evolution simulation method based on thermoelectric coupling provided by the present invention.

[0043] Figure 5 This invention provides a simulation method for the evolution of the growth morphology of conductive filaments in CBRAM based on thermoelectric coupling, which includes growth diagrams of conductive filaments at different temperatures and schematic diagrams of forming voltages. Detailed Implementation

[0044] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0045] The following is combined Figure 1 This invention describes a method for simulating the growth morphology evolution of CBRAM conductive filaments based on thermoelectric coupling, comprising:

[0046] Step 101: The CBRAM device is simulated using a two-dimensional matrix. Based on the resistor network model, a sparse matrix solver is used to process the nodal voltage equations of the grid points to obtain the potential distribution inside the CBRAM device.

[0047] Step 102: The temperature field during the growth of the conductive filament inside the CBRAM device is calculated using the finite difference method. Based on the potential distribution and temperature field, the rate of each particle undergoing various reactions in the CBRAM device is calculated. The reaction of each particle is determined by the rate based on the KMC method (kinetic Monte Carlo method).

[0048] Step 103: After each particle performs the reaction, it updates the two-dimensional matrix, local resistance network, potential and temperature field of the CBRAM device, and triggers a resistive switching event when it detects that the top electrode and bottom electrode of the CBRAM device are connected through a conductive filament.

[0049] This embodiment designs a CBRAM device simulation model and temperature module based on the KMC method. The simulation platform uses Matlab software. A schematic diagram of the simulation model is shown below. Figure 2 The overall model simulation process is as follows: Figure 3 As shown. The simulation model consists of a main program module, a current calculation module, a particle reaction rate calculation module, a reaction event selection module, a filament connection determination module, a heat conduction module, and a potential calculation module. The specific setup steps are as follows:

[0050] 1. First, a two-dimensional matrix with a thickness of m and a width of n is used to simulate the CBRAM device. The bottom represents an inert electrode, the top is considered to be an active electrode that can be infinitely replenished with metal atoms, and the middle is the insulating dielectric region.

[0051] 2. To simulate the distribution of resistance, voltage, and potential within the device, a discretized resistor network model was adopted. Figure 4 The system calculates resistance and writes a voltage and potential calculation module based on Ohm's law and Kirchhoff's law to update the potential distribution within the device.

[0052] 3. A particle reaction rate calculation module is developed based on the Boltzmann distribution formula. The rate of each particle to undergo various reactions is updated according to the potential distribution. The probability of each reaction to occur is determined by the rate magnitude based on the KMC method. Then, the occurrence of a certain reaction is determined by random sampling.

[0053] 4. Using the finite difference algorithm, a thermo-electric coupling solution module was developed to realize the real-time calculation and dynamic visualization of the two-dimensional temperature field during the growth of the conductive wire. This module achieves accurate thermal management simulation through the following techniques: (1) Establishing a solution framework for the unsteady heat conduction equation and coupling the Joule heat source term Q=J 2 ρ accurately reflects the influence of current density distribution on temperature field; (2) Develop dynamic boundary condition processing algorithm to optimize the calculation format for electrode boundary (three-dimensional difference) and resistive switching layer region (four-dimensional difference); (3) Integrate temperature-sensitive reaction rate correction mechanism to automatically adjust ion migration rate when local temperature exceeds threshold.

[0054] 5. A dynamic parameter update mechanism is adopted, and the following parameters are updated after each KMC step: matrix (recording the distribution of atoms and ions), local resistance network (adjusting the grid resistance value according to the filament morphology), potential and temperature distribution.

[0055] 6. Based on the principles of digital image processing and graph theory, a binary image connected component analysis algorithm is used to develop a module to determine whether the conductive wire is connected. By scanning the atomic distribution matrix formed by the conductive wire, the module detects whether there is a connected atomic path between the top and bottom electrodes, thereby determining whether the device has undergone resistance switching.

[0056] First, in the initialization phase, Cu / HfO is constructed. x The Pt sandwich structure features a two-dimensional discretized grid (1nm precision) with an 18nm thick resistive switching layer. The top Cu electrode is defined as an infinite atomic Cu source, the bottom Pt electrode is grounded, and a material parameter library (including various particle reaction barriers Ea, material dielectric constants, etc.) is loaded.

[0057] Then the main simulation loop begins, with each time step (adaptive Δt, typical value 10)... -12 ~10 -9 s) Execute sequentially: 1) Potential field calculation module, based on dynamic resistance network (Cu filament region ρ=1680Ω·nm, HfO x Zone ρ=10 8 Ω·nm) is processed using a sparse matrix solver to handle 10 4 1) Nodal voltage equations for grid points of magnitude; 2) Reaction kinetics module, calculating the rates of various particle reactions, tip reduction, etc.; 3) Conductive filament morphology update module, dynamically adjusting grid parameters (ΔR up to 10) after executing selected events. 6 (Times); Simultaneously, a multiphysics solver is coupled, and the thermal field evolution is calculated using the explicit finite difference method. When the local temperature exceeds 800K, a resistivity temperature correction (α=0.0017 / K) is triggered. Every 100 steps, the connection status of the conductive wire is detected through 8-neighborhood connected domain analysis (bwlabel algorithm). When any top-level mesh is connected to the bottom layer (resistance drop > orders of magnitude), a resistivity change event is triggered.

[0058] like Figure 2 As shown, the bridged random access memory (CBRAM) is mainly composed of three layers: bottom electrode, resistive switching layer, and top electrode. Figure 2 Behaviors 1-10 in the KMC simulation model represent the reactions that occur to particles: (1) Metal atoms on the anode surface are oxidized and become adsorbed ions; (2) Adsorbed ions become desorbed ions under the action of an electric field; (3) Ions near the electrode are adsorbed on the electrode surface; (4) Desorbed ions migrate in the functional layer; (5) Metal ions migrate in the resistive switching layer; (6) Adsorbed ions are reduced to atoms at the cathode platform; (7) Adsorbed ions are reduced to atoms at the cathode steps; (8) Adsorbed ions are reduced to atoms at the cathode pores; (9) Metal ions are adsorbed at the cathode; (10) Metal ions are desorbed at the cathode.

[0059] Traditional models typically treat temperature as a fixed parameter or use a global average temperature, failing to reflect the local Joule heating generated by current in the conductive channel and neglecting to accurately solve for the heat conduction process between different material layers. Therefore, they struggle to characterize potential local hotspots near the conductive filament and their actual impact on migration rates, interfacial reaction rates, and material states. Furthermore, existing methods lack a multi-physics coupling relationship between potential distribution, Joule heating generation, temperature field changes, and the evolution of the conductive filament morphology. This results in inaccurate reproduction of key filament morphologies and a lack of temperature threshold-triggered material damage or phase transition mechanisms, leading to insufficient predictive power for device reliability and thermally induced failure. Therefore, establishing a method that comprehensively considers thermal effects and heat conduction mechanisms to analyze the evolution of conductive filament growth morphology is crucial for revealing the internal physical mechanisms of CBRAM devices and guiding material design and process optimization.

[0060] This embodiment achieves a strong coupling simulation of electric field, thermal effect, and conductive filament morphology evolution by introducing a local heat source calculation model and an explicit finite difference heat conduction solution module. This enables the temperature distribution to dynamically reflect the spatial accumulation and diffusion of Joule heat, and regulates ion migration, interface reaction, and material state renewal through temperature rise feedback. As a result, it can accurately simulate the influence of local heating on the growth and morphology evolution of conductive filaments, improve the accuracy of CBRAM device performance simulation, and effectively make up for the limitations of existing technologies in terms of insufficient consideration of thermal effects and lack of electro-thermal coupling.

[0061] Based on the above embodiments, this embodiment uses a finite difference algorithm to calculate the temperature field during the growth of the conductive filament inside the CBRAM device, including:

[0062] The thermal diffusivity of the CBRAM device is calculated based on the material thermal properties of the CBRAM device, and the stability factor of the CBRAM device is calculated based on the thermal diffusivity.

[0063] The local Joule heat source is obtained based on the cell current and the potential difference between adjacent grids of the CBRAM device.

[0064] The temperature field during the growth of the conductive wire is calculated based on the stability factor of the CBRAM device, the local Joule heat source, and the material thermophysical parameters.

[0065] This embodiment proposes a thermal effect solution module. By establishing a two-dimensional discrete temperature field, this module comprehensively considers local Joule heating generation, thermal diffusion, and temperature threshold-triggered structural evolution, thereby achieving accurate simulation of transient temperature and morphological changes of conductive paths under the action of an applied electric field.

[0066] Based on the above embodiments, this embodiment calculates the thermal diffusivity of the CBRAM device using the following formula based on the material thermal properties of the CBRAM device:

[0067]

[0068] in, Where is the thermal diffusivity, Thermal conductivity, For material density, Specific heat capacity.

[0069] Based on the above embodiments, this embodiment calculates the stability factor of the explicit differential scheme of the CBRAM device using the following formula based on the thermal diffusivity:

[0070]

[0071] Where F is the stability factor, Where is the thermal diffusivity, For temperature evolution step size, For the lattice spacing, F is less than a first preset threshold, such as... It meets the numerical stability condition.

[0072] Based on the above embodiments, this embodiment obtains the local Joule heat source using the following formula based on the cell current and the potential difference between adjacent grids of the CBRAM device:

[0073]

[0074]

[0075]

[0076] in, This indicates that once the conductive path is completed, the total current is evenly distributed among the visible metal atoms in each row according to their quantity. The obtained coordinates Unit current at that point The number of visible metal atoms in the j-th row. Indicates the potential difference between adjacent grid cells. Representing coordinates Voltage at that point Representing coordinates Voltage at that point Representing coordinates The local Joule heat source at that location, Let Q(i,j) be the temperature evolution step size. In the non-metallic region, the current is zero, and Q(i,j)=0.

[0077] Based on the above embodiments, the temperature evolution step size in this embodiment The calculation formula is:

[0078] = +( ) / total rate

[0079] in, It is a random number within the interval [0,1].

[0080] Based on the above embodiments, this embodiment uses a two-dimensional explicit finite difference scheme to solve the heat conduction equation for the temperature field. For the non-bottom nodes inside the CBRAM device, the temperature field during the conductive filament growth process is calculated using the following formula based on the stability factor of the CBRAM device, the local Joule heat source, and the material thermophysical parameters:

[0081]

[0082] For the bottom-level nodes inside the CBRAM device, due to the lack of a lower neighborhood, the temperature field during the conductive wire growth process is calculated using the following formula based on the stability factor of the CBRAM device, the local Joule heat source, and the material thermophysical parameters:

[0083]

[0084] in, and Representing the coordinates at time t+1 and t respectively. The temperature at that location Indicates the stabilizing factor. The unit conversion scaling factor.

[0085] Based on the above embodiments, this embodiment, after calculating the temperature field during the growth of the conductive filament inside the CBRAM device using the finite difference algorithm, further includes:

[0086] If the temperature of any grid in the CBRAM device is greater than the second preset threshold and the grid is not located in the bottom electrode layer, then high-temperature induced material instability or atomic desorption occurs at that grid. The corresponding metal atoms and adsorbed ions are marked to zero to reflect structural damage or local melting behavior.

[0087] After temperature evolution, when the temperature of a certain grid satisfies: (For example If the temperature is 1000K and the location is not part of the bottom electrode layer, it is considered that high-temperature induced material instability or atomic desorption has occurred at that location. Therefore, the corresponding metal atoms and adsorbed ions are marked as zero to reflect structural damage or localized melting behavior. Isothermal boundary conditions are applied to the left and right boundaries of the system. To simulate lateral heat dissipation.

[0088] Through the above calculation process, this module realizes the generation of local Joule heat, the diffusion of heat in the conductive wire structure, and the material state changes caused by temperature rise, forming a closed-loop coupled model that can reproduce the real behavior of conductive wires under working conditions, such as growth, fracture, and reconstruction caused by thermal effects.

[0089] Based on the above embodiments, this embodiment calculates the rate of various reactions occurring in each particle of the CBRAM device using the following formula, based on the potential distribution and temperature field:

[0090]

[0091] Where Γ is the rate at which each particle undergoes various reactions, Ea is the activation energy required for each reaction to occur, and f is the vibrational frequency of the electron in the solid. Δφ is the correction factor, Δφ is the potential energy change before and after each reaction, k is the Boltzmann constant, q is the charge of a single electron, and T is the temperature field.

[0092] Based on the definition of dynamic Monte Carlo (KMC), the rates during the device model transition process all conform to the Boltzmann distribution: Based on the principles of non-equilibrium statistical mechanics and electrochemical kinetics, a refined multiphysics coupling algorithm was used to achieve atomic-level precision reaction kinetic modeling. This module first establishes a 1nm resolution spatial discrete grid and then uses the central difference method to calculate the potential gradient at each grid point. This allows for the acquisition of local electric field distributions; for 10 key reactions (including anodic oxidation, ion migration, and cathode step-site reduction), a modified Boltzmann formula is used. Rate calculation models were constructed, with the modifier α=0.5 calibrated using first-principles calculations, and the temperature field T(x,y) coupled with the Joule heating effect in real time. Table 1 shows the main parameter values ​​in the simulation, and the activation energy parameters required for the reaction process were obtained from first-principles calculations and experimental fitting.

[0093] Table 1. Reaction types, rate formulas, various potential barriers, and their values ​​in the KMC model.

[0094]

[0095] The thermal effects module implements a multiphysics simulator coupled with electrothermal effects. Its core principle is to simulate the dynamic influence of Joule heating on the temperature field and phase state of a material system by solving the heat conduction equation. The specific implementation process is as follows: First, the function extracts key physical quantities from the input parameters, including grid size, material thermal properties (heat capacity, thermal conductivity, ion density), and temperature limit threshold. Based on these parameters, it calculates the thermal diffusivity coefficient alpha = thermal conductivity / (density × heat capacity) and the stability condition coefficient F. This F coefficient is a key parameter for the stability of the explicit finite difference method and needs to satisfy the numerical stability condition F ≤ 0.25. In the heat source calculation stage, the function calculates the Joule heat Q at each grid point based on the current distribution and potential difference: for the location of copper atoms (visibility_atom(ii,jj)==1), the heat generation is calculated by multiplying the potential difference by the distributed current current_elem and then by the time step; for oxide regions, heat generation is not considered. Here, current_elem is obtained by evenly distributing the total current current according to the number of copper atoms, reflecting the distribution characteristics of the current in the conductive channel. The heat conduction equation is solved using an explicit finite difference method: for internal nodes, the temperature update includes the heat conduction contribution from four adjacent nodes and a local heat source term; for the bottom boundary node (ii==r), it simplifies to a three-neighbor conduction mode. This discretization method transforms the continuous heat conduction partial differential equation into an algebraic iterative form, simulating the temperature field evolution through time stepping. After the temperature field is updated, the function performs phase state detection: when the temperature at a certain location exceeds a threshold (1000K), the corresponding copper atoms and adsorbed ions are set to zero, simulating material decomposition or phase transition caused by high temperature. The boundary temperature is kept constant at 300K to simulate heat dissipation conditions. The entire algorithm achieves multi-field coupling of electro-thermal-structure: the current affects the temperature field through Joule heating, and the temperature field reacts to the electrical conductivity by changing the material state. This bidirectional coupling mechanism can simulate the thermally induced resistance change phenomenon in devices, providing a numerical experimental platform for studying electrothermal synergistic effects. Finally, the function outputs the updated temperature field T, time step delta_t, and ion / atom distribution matrix, providing updated material parameters for subsequent electrical calculations.

[0096] The core of the main program module is to simulate the resistive switching behavior of CBRAM by coupling atomic migration and thermal effects through a multi-scale modeling method. The program first initializes physical constants and material parameters, including the fundamental charge q, Boltzmann constant k, and reaction barrier, and constructs a sandwich structure device model consisting of a top electrode, a bottom electrode, and a resistive switching layer. For spatial discretization, the device is divided into an m×n grid, with each grid representing a 1nm physical dimension. The positions of copper atoms are marked using the `visibility_atom` matrix, and the `visibility_ion` matrix tracks the ion distribution. The core loop of the program consists of two stages: a SET process (positive voltage scan) and a RESET process (negative voltage scan). In the SET stage, the potential distribution is obtained by solving the Poisson equation, and various reaction rates (atomic oxidation / reduction, ion adsorption / desorption, ion migration, etc.) are calculated using the Arrhenius equation. A kinetic Monte Carlo method is then used to select the occurring reactions and update the atomic / ion distribution. When the conductive filament (CF) forms and causes the device to conduct, the system switches to the RESET stage. At this stage, metal atom oxidation and Joule heating become dominant. The temperature field evolution is simulated using the heat conduction equation. When the local temperature exceeds 1000K, copper atoms decompose or oxidize into ions, causing the filament to break and achieving a resistance change. Current calculation employs an improved nodal analysis method, considering the parallel resistance effect of the conductive filament. The nodal potential is solved using the potential calculation module, and the current is calculated based on Ohm's law. The thermal management module couples the Joule heating generation and heat conduction processes. The heat source term Q is calculated from the local current density and potential difference. The temperature field update uses an explicit finite difference method, considering the effect of temperature on the copper resistivity. Figure 5 As shown, the program records the IV characteristic curves and the microstructure of the conductive filaments in real time during runtime, saving the visualization results of atomic configuration, potential distribution, and temperature field every 50 steps (generated as a pseudo-color image using the pcolor function). The IV characteristic curves reveal that as temperature increases, both the forming time and voltage of the device decrease, while the abrupt change in current increases. The conductive filament growth diagram shows that the height and number of weak filaments increase with increasing temperature. Furthermore, the system innovatively introduces an electronic dynamic model, tracking electron distribution through an Electron matrix to simulate the influence of electron-ion interactions on the reduction reaction. The entire simulation is controlled by setting the number of loops, typically running 10 complete SET-RESET cycles, ultimately outputting data including time, voltage, current, and resistance. In the core loop of the main program, each time a KMC reaction event selection and structure update operation (including ion migration, reduction deposition, etc.) is executed, the program synchronously updates the atomic distribution matrix (visibility_atom) and the local potential distribution (potential).

[0097] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for simulating the growth morphology evolution of CBRAM conductive filaments based on thermoelectric coupling modeling, characterized in that, include: A two-dimensional matrix is ​​used to simulate the CBRAM device. Based on the resistor network model, a sparse matrix solver is used to process the nodal voltage equations of the grid points to obtain the potential distribution inside the CBRAM device. The temperature field during the growth of the conductive filament inside the CBRAM device is calculated using the finite difference method. Based on the potential distribution and temperature field, the rate of various reactions of each particle in the CBRAM device is calculated. The reaction of each particle is determined by the rate based on the KMC method. Each particle updates the two-dimensional matrix, local resistance network, potential, and temperature field of the CBRAM device after performing a reaction, and triggers a resistive switching event when a connection is detected between the top and bottom electrodes of the CBRAM device via a conductive filament.

2. The method for simulating the growth morphology evolution of CBRAM conductive wires based on thermoelectric coupling according to claim 1, characterized in that, The temperature field during the growth of the conductive filament inside the CBRAM device was calculated using the finite difference method, including: The thermal diffusivity of the CBRAM device is calculated based on the material thermal properties of the CBRAM device, and the stability factor of the CBRAM device is calculated based on the thermal diffusivity. The local Joule heat source is obtained based on the cell current and the potential difference between adjacent grids of the CBRAM device. The temperature field during the growth of the conductive wire is calculated based on the stability factor of the CBRAM device, the local Joule heat source, and the material thermophysical parameters.

3. The method for simulating the growth morphology evolution of CBRAM conductive wires based on thermoelectric coupling according to claim 2, characterized in that, The thermal diffusivity of the CBRAM device is calculated using the following formula based on the material's thermal properties: ； in, Where is the thermal diffusivity, Thermal conductivity, For material density, Specific heat capacity.

4. The method for simulating the growth morphology evolution of CBRAM conductive wires based on thermoelectric coupling according to claim 2, characterized in that, The stability factor of the CBRAM device is calculated based on the thermal diffusivity using the following formula: ； Where F is the stability factor, Where is the thermal diffusivity, For temperature evolution step size, The lattice spacing is F, which is less than the first preset threshold.

5. The method for simulating the growth morphology evolution of CBRAM conductive wires based on thermoelectric coupling according to claim 2, characterized in that, The local Joule heat source is obtained using the following formula based on the cell current and the potential difference between adjacent grids of the CBRAM device: ； ； ； in, This indicates that once the conductive path is completed, the total current is evenly distributed among the visible metal atoms in each row according to their quantity. The obtained coordinates Unit current at that point The number of visible metal atoms in the j-th row. Indicates the potential difference between adjacent grid cells. Representing coordinates Voltage at that point Representing coordinates Voltage at that point Representing coordinates The local Joule heat source at that location, Let Q(i,j) be the temperature evolution step size. In the non-metallic region, the current is zero, and Q(i,j)=0.

6. The method for simulating the growth morphology evolution of CBRAM conductive wires based on thermoelectric coupling according to claim 4 or 5, characterized in that, Temperature evolution step The calculation formula is: = +( ) / total rate; in, It is a random number within the interval [0,1].

7. The method for simulating the growth morphology evolution of CBRAM conductive wires based on thermoelectric coupling according to claim 2, characterized in that, For the non-bottom nodes inside the CBRAM device, the temperature field during the conductive filament growth process is calculated using the following formula based on the stability factor of the CBRAM device, the local Joule heat source, and the material thermophysical parameters: ； For the bottom-level nodes inside the CBRAM device, the temperature field during the growth of the conductive filament is calculated using the following formula based on the stability factor of the CBRAM device, the local Joule heat source, and the material thermophysical parameters: ； in, and Representing the coordinates at time t+1 and t respectively. The temperature at that location Indicates the stabilizing factor. The unit conversion scaling factor.

8. The method for simulating the growth morphology evolution of CBRAM conductive wires based on thermoelectric coupling according to claim 1, characterized in that, After calculating the temperature field during the growth of the conductive filament inside the CBRAM device using the finite difference algorithm, the method further includes: If the temperature of any grid in the CBRAM device is greater than the second preset threshold and the grid is not located in the bottom electrode layer, then high-temperature induced material instability or atomic desorption occurs at that grid. The corresponding metal atoms and adsorbed ions are marked to zero to reflect structural damage or local melting behavior.

9. The method for simulating the growth morphology evolution of CBRAM conductive wires based on thermoelectric coupling according to claim 1, characterized in that, The rates of various reactions occurring in each particle of the CBRAM device are calculated using the following formulas based on the potential distribution and temperature field: ； Where Γ is the rate at which each particle undergoes various reactions, and E a The activation energy required for various reactions to occur, where f is the vibrational frequency of electrons in the solid. Δφ is the correction factor, Δφ is the potential energy change before and after each reaction, k is the Boltzmann constant, q is the charge of a single electron, and T is the temperature.

Images