High-throughput phase-field simulations and machine learning of resistive switching in resistive random-access memory

Metal oxide-based Resistive Random-Access Memory (RRAM) exhibits multiple resistance states, arising from the activation/deactivation of a conductive filament (CF) inside a switching layer. Understanding CF formation kinetics is critical to achieving optimal functionality of RRAM. Here a phase-field model is developed, based on materials properties determined by ab initio calculations, to investigate the role of electrical bias, heat transport and defect-induced Vegard strain in the resistive switching behavior, using MO2−x systems such as HfO2−x as a prototypical model system. It successfully captures the CF formation and resultant bipolar resistive switching characteristics. High-throughput simulations are performed for RRAMs with different material parameters to establish a dataset, based on which a compressed-sensing machine learning is conducted to derive interpretable analytical models for device performance (current on/off ratio and switching time) metrics in terms of key material parameters (electrical and thermal conductivities, Vegard strain coefficients). These analytical models reveal that optimal performance (i.e., high current on/off ratio and low switching time) can be achieved in materials with a low Lorenz number, a fundamental material constant. This work provides a fundamental understanding to the resistive switching in RRAM and demonstrates a computational data-driven methodology of materials selection for improved RRAM performance, which can also be applied to other electro-thermo-mechanical systems.


INTRODUCTION
Oxide-based Resistive-Random Access Memory (RRAM) has attracted broad attention as a potential candidate for nextgeneration nonvolatile memories, due to its fast switching speed, small programing current 1 , controllable resistance states 2 and ease of fabrication 3 . The key component of a RRAM is a memristor, a two-terminal "memory resistor" that can maintain (memorize) its electrical resistance based on the history of an applied voltage and/or current. Thus, its resistance is tunable on demand and can be reversely changed between a high resistance/reset state (HRS) and a low resistance/set state (LRS) 4,5 . The typical structure of RRAM consists of an electrically insulating layer (oxide) covered by top and bottom electrodes. Under the applied voltages, electrically conductive filaments (CFs) can either be formed to connect both electrodes resulting in LRS in the RRAM, or ruptured inside the oxide layer to switch the device back to HRS. Therefore, the resistance switching mechanisms are generally attributed to the filamentary modification of the conduction properties. Depending on the type of atoms (ions) that form the CFs, RRAM can be classified into (a) electrochemical metallization devices (ECM) where the CF is made of extrinsic metal ions from the electrochemical active electrode such as Ag 6,7 , and (b) valance change memory devices (VCM) where the active species are defects such as oxygen vacancies that are present internally in the CFs [8][9][10][11][12] .
One of the critical challenges that inhibit the applications of RRAM is the randomness of the formation and dissolution of the CFs, which results in large variance from one device to another and from one switching cycle to another in the same device. Although the physical reconfigurations of the CFs have been captured by direct imaging methods such as TEM [13][14][15][16][17] and conductive AFM (c-AFM) [18][19][20] , the coupling of the coexisting and concurrent chemical, electrical, mechanical and thermal processes during growth of CFs make a resistive switching process extremely complicated and unpredictable, and its underlying mechanism far from being fully understood.
Several theoretical/physical models have been used to study the CF formation dynamics during resistance switching. For example, atomic-scale simulations based on density functional theories and molecular dynamics have been employed to calculate the oxygen vacancy formation energy and the migration barriers in various oxide switching layers such as Al 2 O 3 21 , TiO 2 22 , HfO 2 23 , and Ta 2 O 5 24,25 , and demonstrate the formation/rupture of metallic CFs in Cu/SiO 2 cells 26,27 . This atomistic modeling is limited by its spatial and temporal scales that are not typically accessible in experiments. Mesoscale models, based on the solutions to partial differential equations for mass transport, electrical conduction and heat generation/dissipation, have been applied to quantitatively describe the resistive switching behaviors in both VCM and ECM devices [28][29][30][31][32] . However, most existing models do not take into consideration the possible impact of mechanical stress during CF evolution 33 . It is found that the internal point defect concentrations can alter the local mechanical strains and electrical potential [34][35][36][37] . Recently, chemo-mechanical coupling in nonstoichiometric metal oxides has been studied and is shown to have strong influence on the properties of many electrical and energy devices such as solid-oxide fuel cells and catalysts 38,39 . Therefore, it is necessary to carefully assess the mechanical effect induced by the defect distribution and transport on the resistance switching behavior of RRAM.
The large current on/off ratio and fast switching speed are the two key characteristics for RRAM. For metal oxide based VCM devices, both characteristics are highly dependent on the kinetics of oxygen vacancy migration that accounts for the formation/ dissolution/reconnection of the CFs. While reliable switching up tõ 10 9 cycles have been demonstrated in TaO x 40 and HfO x 5 based systems recently, other oxides with similar electrical performances have not been explored as potential candidates for RRAM. Exploration of new materials and structures through conventional experimental trial-and-error approaches is highly inefficient. Here, we propose to employ a combination of high throughput computation and machine learning to guide the materials discovery and design for RRAM. Such an approach has recently been demonstrated for defective oxides and nanomaterials 41,42 . Shen et al. also employed such a methodology to screen and identify potential oxide nanofillers in polymer-based dielectrics for improved dielectric breakdown strength 43 . Zhang et al. demonstrated an unsupervised approach for discovering new candidates of solid state Li-ion conductors with conductivities of 10 −4~1 0 −1 S cm −1 predicted in ab initio molecular dynamics simulations 44 .
In this work, we first developed a phase-field model to incorporate the impact of the mechanical stress on the resistive switching behavior and its interaction with mass diffusion, thermal transport, and electrical conduction dynamics in metal oxidebased RRAM. Thin-film metal-oxides, such as HfO 2 , are chosen as prototypes in our simulations. By parameterizing the three key material constants of metal oxides, i.e., the Vegard strain coefficient (V ij ), the electrical conductivity (σ), and the thermal conductivity (k th ) to represent different insulating oxide layers in RRAM, we perform high-throughput phase-field simulations to independently and systematically explore each of these intrinsic material constants on the performance of metal oxide-based RRAM and generate a material-property database. Based on this we apply compressed-sensing based machine learning to derive analytical prediction models for the device performance including the current on/off ratio (I on /I off ) and resistance switching time (t switch ) of RRAM as a function of the aforementioned key material constants. These analytical models provide fundamental physical insights for accelerating future discovery of RRAM materials.

RESULTS AND DISCUSSION
In this work, we employ HfO 2 metal oxide as an example and describe resistive switching process through defect migration driven by a local electric field and a temperature gradient. A two-dimensional (2D) phase-field simulations model is used to describe the resistance switching process by solving the Nernst-Planck equations for defect migration, current continuity equation for electronic conduction, and thermal transport equation for Joule heating (see "Methods" for details). Figure 1a illustrates the 2D simulation geometry of a 3D axis-symmetric HfO 2 memristor device, which consists of a highly conductive filamentary region (HfO 2−x ) with high oxygen deficiency and the remaining insulating region of stoichiometric HfO 2 . The simulation starts after the completion of the electroforming process where the CF is formed to connect the top (TE) and bottom electrodes (BE), which acts as n-type doping sites for electrical and thermal conduction, as shown in Fig. 1a. Here the oxygen vacancies are assumed to be uniformly distributed in the CF region with a concentration of N V€ O ¼ 1:2 10 27 m À3 and vanish in the remaining region in the initial state. Then a triangular voltage sweep with rate dV/dt = 0.1 V/s (inset of Fig. 1b) is applied on the TE while the BE is grounded.
A hysteresis-like current-voltage behavior is clearly seen, indicating the bipolar resistive switching characteristics. The reset voltage is around 0.4 V where the resistance starts to increase, and the device is switched from LRS to HRS. On the other hand, the set transition occurs at negative voltage (−0.57 V) and switches the HRS to LRS. Our simulation results agree well with experimental measurements and theoretical calculations of HfO 2 based memristor device from existing literature 5 , indicating that the current phase-field model is able to capture the resistive switching behavior in metal oxide-based RRAM.
To further understand the physical nature of the reset and set processes, we studied the evolutions of the oxygen vacancies N V€ O , temperature T and electric potential φ. The results are shown in Fig. 1c-k. During the reset process, V ö migrates to the BE and forms a local V ö deficiency gap along the CF accompanied by the local electric field enhancement and temperature increase due to the Joule heating effect. Therefore, the resistance increases with the increasing applied voltage. During the set process, the gap formed by reset process is filled via the V ö migration towards TE under negatively applied voltage, accompanied by the decrease of local temperature along the CF. The reconnection of the CF allows fast electrical conduction along the CF and reduces the overall resistance of the device. It is also seen that an additional V ö depletion region is formed near the BE during the set process, as depicted in Fig. 1i. This is because the BE is unable to provide sufficient V ö continuously. In this case, a relatively higher overall resistance after the set process is induced compared with the initial value.
Mechanical strain effect in RRAM In metal oxide-based RRAM, the oxygen vacancy concentration can reach up to~10 27 m −3 . This creates a lattice shrinkage strain up to 1%, which could influence the vacancy transport and the resistance switching dynamics. Here we study the effect of V ö induced local strain on resistive switching behavior in metal oxidebased RRAM. The local strain ε 0 ij induced by the change of V ö density can be calculated as ε 0 is the variation of oxygen vacancy density. The constant value N V€ Oi is the V ö density in the reference stress-free state. Here the initial electroforming state is chosen as the reference stress-free state. This mechanical strain causes an additional elastic driving force for the V ö migration based on Eq. (10). We assume V ij to be −0.064 based on our density functional theory (DFT) calculations, benchmarked to more accurate many-body theory based methods 45 , indicating that the oxygen vacancy will create a lattice shrinkage in HfO 2 (see "Methods" for details about DFT calculations). Figure 2a illustrates the local strain distribution after reset process under 1.1 V applied voltage. Due to the migration of oxygen vacancies toward the BE, a local tensile strain up to 0.5% is seen in the V ö -rich region inside the CF near the BE (blue region), and a local compressive strain arises from the CF gap region where oxygen vacancies deplete (red region). In addition, a tensile strain is also seen aside the CF due to the lateral diffusion of oxygen vacancies. This strain dipole creates strain gradients (dε 0 ij =dz and dε 0 ij =dx) along both vertical and lateral directions, which act as an additional driving force for the oxygen vacancy transport. We compare the elastic (μ elastic ), electrical (μ electric ), and chemical potentials (μ chem ) for the initial electroforming state, the intermediate state and the final reset state as shown in Supplementary Fig. S3. The μ elastic map ( Supplementary Fig.  S3 a, d, g) seems to reflect the oxygen vacancy transport, i.e., it is homogeneous inside the CF in the initial state and becomes higher/ lower near the BE/TE during the reset process. This is opposite to the electrical potential distribution which reduces from TE to BE (Supplementary Fig. S3 c, f, i), indicating that the elastic effect counterbalances the electrical effect. Meanwhile, the μ elastic map ( Supplementary Fig. S3 d, g) also induces an additional lateral elastic potential gradient which promotes the lateral diffusion of oxygen vacancies. We analyze the 1D profiles of each of the energy potentials along z direction across the center of the CF, as shown in Fig. 2b. At z = 17 nm (gap region), the chemical, electrical and elastic potentials all exhibit significant positive gradients, among which the electrical potential shows the largest drop (~10 5 J mol −1 ), acting as the major driving force for the ion migration from TE to BE which makes a gap nearby the TE. On the other hand, at z = 5-15 nm where V ö locally segregates/depletes, the electric potential gradually increases, while the chemical and elastic potentials decrease. The negative elastic potential gradient partially offsets the positive electrical potential gradient in this region. Based on Fig. S3 and Fig. 2b, the elastic potential inhibits the vertical V ö migration in the CF, and promotes the lateral migration of V ö into the insulating region. This results in a reduced area of local V ö segregation inside the CF near the BE (Fig.  2d), as compared to the case when μ elastic is not taken into consideration (Fig. 2c). This is further illustrated by the comparison of 1D profiles of N V€ O with and without considering μ elastic , along both vertical and horizontal directions (Fig. 2e).
The effect of mechanical strain on the temporal evolution of the overall resistance of the entire simulation system during the reset process is shown in Fig. 2f. Here we use a constant voltage pulse (rectangular pulse) of 1.1 V with pulse width from 10 ps to 10 us. For both cases, the overall resistance in the initial set/current on state is R on~0 .335 kΩ (or I on~3 mA) and increases with time. However, when the mechanical effect is considered, the final resistance in reset/current off state is calculated to be R off = 0.969 kΩ, nearly 12% lower than that without considering the mechanical effect (1.1 kΩ). This causes a decrease of current on/off ratio (I on /I off ) from 3.28 to 2.89, and a slight increase in the resistive switching time (t switch ) from 18 to 23 ns (here t switch is defined when the cell resistance increase by more than 50% from its initial value). This is because the elastic potential inhibits the V ö transport and local segregation, resulting in a longer time for the CF gap formation and a reduced value of R off . Our results indicate that the mechanical effect could influence the ion transport dynamics and the resistive switching characteristics.
Effect of electrical, thermal, and elastic properties on the performance of memristor From previous studies, it is clear that the resistive switching behaviors, such as I on /I off and t switch , are heavily dependent on the oxygen vacancy transport driven by the combined chemical, electrical, thermal, and elastic effects. These effects scale with their relevant kinetic coefficients, such as the σ and k th of the insulating metal oxides and their corresponding metallic states (of high oxygen deficiency), and the Vegard strain coefficients V ij . While V ij is assumed to be a constant, σ and k th are spatially dependent on the local V ö density. This significantly increases the complexity of the vacancy transport dynamics and the resistive switching behavior. To untangle all these couplings and understand their independent roles in the resistive switching, we performed systematic phase-field simulations by tuning one of the forementioned three material constants while fixing the other two. Since σ and k th for stoichiometric metal oxides (such as HfO 2 , TaO 2 , etc.) are much lower than their corresponding metallic states (such as Hf, Ta, etc.), for simplicity we assume σ and k th remain constant for different oxides in the absence of V ö , and are linearly extrapolated to their metallic values with increasing N V€ O . In this way σ and k th for different metals oxides can simply be represented by the slopes (K 1 and K 2 ) of the positive linear extrapolations. (see "Methods" and "Supplementary Information Note 1" for more details). Figure 3 shows the temporal evolutions of the overall resistance and the dependence of I on /I off on V ij , K 1 , and K 2 , respectively. When V ij increases, the initial resistance (on state) remains constant for all V ij 's, while the final resistance (off state) decreases with increasing V ij , as shown in Fig. 3a, which results in a decreasing I on /I off with increasing V ij as seen in Fig. 3d. This agrees with our previous study that the additional mechanical driving force tend to reduce the current on/off ratio.
The effect of electrical conductivity σ (by its slope K 1 ) on the overall resistance is shown in Fig. 3b. The initial resistance decreases with increasing K 1 , while the final resistance experiences a decrease-increase-decrease behavior. This results in a peak value of I on /I off~1 5 at K 1 = 7.5 shown in Fig. 3e. To further understand this behavior, we plot the N V€ O , T and φ distributions under three typical electrical conductivities corresponding to K 1 = 0.5, 7.5, and 15, as shown in Supplementary Fig. S4. Under low electrical conductivity (K 1 = 0.5), oxygen vacancies cannot migrate from TE to BE and no gap is formed inside the CF (Supplementary Fig. S4a), no resistive switching occurs and I on /I off remains low. When σ further increases (K 1 = 7.5), a gap is formed and enlarged inside the CF by the V ö migration ( Supplementary  Fig. S4d). This causes a significant increase in the resistance of the final reset state (Fig. 3b, cyan curve), as well as I on /I off increase. When K 1 = 15, large heat generation occurs by V ö migration due to the Joule heating effect (Supplementary Fig. S4h), which gives rise to a significant oxygen vacancy transport in the lateral direction and reduces the local V ö segregation/depletion inside the CF, as shown in Supplementary Fig. S4g. Thus, the CF gap disappears and I on /I off reduces with increasing electrical conductivity.
Finally, the effect of thermal conductivity k th (by its slope K 2 ) on the switching behavior is shown in Fig. 3c. It is seen that the initial resistance of RRAM remains the same, while the final resistance in reset state decreases with increasing k th . This results in a monotonous decrease of I on /I off with increasing k th . Supplementary Fig. S5a-c demonstrates that lower k th favors the V ö migration along the electrical field direction and the formation of larger CF gap, giving rise to a higher I on /I off , as shown in Fig. 3f. Higher k th inhibits the vertical transport of the V ö . When k th reaches a very high value, the heat generated by Joule heating can easily dissipate through the CF (Supplementary Fig. S5g-i). The reduction of the thermal mobility of V ö increases the difficulty of the gap formation along the CF, resulting in the decrease of I on /I off . When K 2 > 0.8, I on /I off is around 1.0 indicating that resistive switching no longer occurs. a b c d e f Fig. 3 Effects of material parameters on the device performance. a-c Temporal evolutions of overall resistance with different Vegard strain coefficient V ij , electrical conductivity (by K 1 ), and thermal conductivity (by K 2 ), d-f dependence of current on/off ratio I on /I off on V ij , K 1 and K 2 .
K. Zhang et al.
High-throughput phase-field simulations and machine learning Based on the general trends obtained by isolating and studying each of the three characteristic material constants, we performed high-throughput phase-field simulations by parameterizing V ij , K 1 , and K 2 to calculate the corresponding I on /I off and t switch , and establish a "material constant-resistive switching property" database in metal oxides based RRAM. The ranges of V ij , σ and k th are chosen to be 0-0.24 (in absolute value V ij ), 0.06-18 (10 3 S cm −1 ), and 6.5-120.5 (W m −1 K −1 ), respectively, which include most of the metal oxides used in RRAM. The results are shown in Fig. 4a, b. It is clearly seen that most I on /I off ratios are in the range of 1-10, with a few exceptions reaching up to 20. On the other hand, from the database, it shows a general trend that metal oxides with larger electrical conductivity and lower thermal conductivity exhibit relatively larger I on /I off and smaller t switch , as indicated by the red region in Fig. 4a and blue region in Fig. 4b. By comparison, the Vegard strain coefficient has relatively smaller influence on the performance of metal oxide-based RRAM.
To clearly identify the trend from 3D database, we plotted the 2D mapped I on /I off and t switch as a function of two out of the three materials parameters with a fixed remaining parameter (Fig. 4c-h). Figure 4c shows the dependence of I on /I off on σ and V ij at fixed k th ¼ 18:5 W m −1 K −1 . These results indicate that the electrical and thermal conductivities play dominant roles than the Vegard strain coefficient during resistive switching in metal oxide materials, as shown in Fig. 4c-e. It is clearly seen that higher σ and lower k th give rise to higher I on /I off , as marked by the black circle shown in Fig. 4e. When σ is even higher than 12 × 10 3 S cm −1 while k th is lower than 20 W m −1 K −1 it causes further drop of I on /I off , as shown in the up-left corner. On the other hand, t switch demonstrates generally opposite trend to I on /I off , i.e., smaller t switch (or faster switching) is realized where I on /I off is large, as illustrated in  Fig. 4 High-throughput phase-field simulations. a The current on/off ratios (I on /I off ) and b the resistive switching time (t switch ) for metal oxides-based RRAM by parameterizing the three characteristic materials constant, i.e., the Vegard strain coefficient (V ij ), electrical conductivity (σ) and thermal conductivity (k th ). 2D mapped I on /I off c-e and t switch f-h as a function of two out of the three materials parameters (V ij , σ, k th ) at a fixed remaining parameter. Fig. 4f-h. In addition, several data points marked with colored solid symbols in Fig. 4e represent the metal oxides of VO 2 , Ta 2 O 5 , SnO 2 , NiO, HfO 2 , respectively, which are commonly used as switching layer in metal oxide-based RRAM 1,3 . It shows that the performance of metal oxide-based RRAM employed with these materials are not optimal. To further improve the performance of metal oxidebased RRAM, increasing the electrical conductivity and decreasing thermal conductivity, as well as avoiding large Vegard strain is necessary.
Based on the database from high-throughput simulations, a recently developed compressed-sensing based machine learning (ML) approach 46 is employed to further elucidate the correlation between (V ij , σ, k th ) and (I on /I off , t switch ). First, linear correlations among V ij , σ, k th , I on /I off , t switch were investigated using the Pearson correlation coefficient, as shown in Fig. 5a. Clearly V ij , σ, k th are positively correlated with themselves, and non-correlated between each other. On the other hand, σ positively correlates with I on /I off , but negatively correlates with t switch , indicating that I on /I off (t switch ) increases (decreases) with increasing σ. Meanwhile, k th exhibits negative (positive) correlations I on /I off (t switch ). These results agree with our previous calculations. Finally V ij shows slight negative correlation with I on /I off , and almost no correlation with t switch , indicating that V ij bearly influences t switch .
We choose Vegard strain coefficient V ij , electrical conductivity σ (10 3 S cm −1 ), and thermal conductivity k th (W m −1 K −1 ) of metal oxides as the fingerprints to build an interpretable machine learning model. Since some of the properties could be correlated (e.g., Wiedemann-Franz law that relates electrical conductivity (σ) and thermal conductivity (k th )), we use the Sure Independence Screening and Sparsity Operator (SISSO) approach 46 , to predict I on /I off as well as t switch . SISSO not only allows dimensional reduction of the size of the (correlated) feature-space using the more numerically stable sure independence screening (SIS) approach 47,48 , but also allows us to progressively expand the dimensionality of the descriptor (i.e., construct a nD-descriptor), even in the presence of correlated features, using sparsifying operators (SO) iteratively on the residuals, thereby allowing us to build optimally interpretable physics-based models. The training error root-mean-squared-error (RMSE) is used as the criterion to select the suitable nD-descriptors. After performing machine learning on the training datasets, 2D-descriptor based predictive expressions of I on /I off with RMSE = 0.453 and t switch with RMSE = 0.385 were found to be optimal (i.e., with relatively low and high interpretability), as shown in Eqs. (1) From Eq. (1), by choosing V ij = 0-0.24 which covers most metal oxide, log(t switch ) monotonously increases with increasing k th σ ( W m À1 K À1 10 3 S cm À1 ) from 0.5 to 25, a range where resistive switching occurs. On the other hand, increasing V ij (magnitude only) results in increasing t switch at given k th σ (green lines in Supplementary Fig. S6). Equation (2) reveals that I on /I off generally decreases with increasing kth σ . At very small kth σ (<2), I on /I off increases with increasing k th σ (red line in Supplementary  Fig. S6). It is noteworthy that I on /I off is independent of V ij based on Eq. (2). From the predicted machine learning expressions, it is clearly seen that the properties of memristor can be enhanced by decreasing k th σ . Meanwhile the mechanical strain effect plays an important role in the resistive switching kinetics, which is dependent on the dynamic evolution of V ö . However, the mechanical effect becomes less important in determining the I on /I off , i.e., the change of the electrical current, governed by the initial concentration and final distribution of V ö during reset process are less affected by V ij . This trend obtained from the machine learning predictive models agree generally with the earlier shown phase-field calculations in Fig. 4e and h. Figure 5b, c shows the comparisons of t switch and I on /I off calculated from the phase-field simulation and predicted from the machine learning model. It is seen that the training datasets are clustered near the orange straight line where the phase-field calculated and machine learning predicted properties are equal.
In addition, six testing datapoints of NiO, TiO 2 , SnO 2 , VO 2 , HfO 2 , Ta 2 O 5 are also plotted in Fig. 5b, c (with solid colored symbols). All these six data points are located around the solid orange lines, which further validates our machine learning as a predictive model. It should be noted that none of these six data points was used as training datasets. The materials parameters and the calculated/predicted properties of these metal oxides are listed in Supplementary Table 1. From the table, it is found that when kth σ ratio decreases, t switch decreases while I on /I off increases, from both phase-field calculations and machine learning predictions. According to the Wiedemann-Franz law, for nearly-free electron metals (and oxide metals), the ratio of the electronic contribution of the thermal conductivity (k th ) to the electrical conductivity (σ) of a pure metal is proportional to the temperature (T), i.e., kth σ $ LT, where L is the Lorenz number. Theoretically, k th σ has approximately the same value for different pure metals at the same temperature. But for binary transition-metal oxides, the strong electron-electron correlations have been shown to drastically lower the Lorenz number (L) by several orders of magnitude 49,50 , due to an anomalously low electronic thermal conductivity. Therefore, Lorenz number would be at lower values for metal oxides with stronger electron-electron correlations at given temperatures. In short, a key design principle for future materials-discovery to achieve optimally performing memristor materials would be to select oxides with a relative low Lorenz number at the same temperature, such as those with strong electron-electron correlations. In summary, a comprehensive computational model based on fundamental thermodynamics and kinetics, including the ion and thermal transport, electrical conduction, and chemical expansion strain effect is developed to investigate the resistance switching process in metal oxide-based RRAM. It can successfully capture the resistive switching process. The local oxygen vacancy distribution induces a lattice expansion strain and a strain gradient, which act as an additional driving force that inhibits the V ö migration in the conductive filament, and reduces the current on/off ratios. Highthroughput phase-field simulations are performed to construct a "materials fingerprints-targeted performance" database, and machine learning approaches are employed to establish interpretable analytical predictive expressions for I on /I off and t switch of memristor in different metal oxides. The machine learning model is further verified by additional phase-field calculations of real metal oxides. This high-throughput/machine-learning approach reveals that metal oxides with relatively small k th σ ratios, found in bad-metals with a low Lorenz number, yield high memristor performance, thereby establishing a key materials-design principle for designing new memristor materials by relating device performance metrics to fundamental material constants. This work thus establishes a strategy to select the metal oxides to optimize the performance and provide guidance to experimentalists in designing high performance metal oxide-based RRAM device.

Phase-field model of electro-thermo-mechanical coupling in resistive switching
The resistance switching process is determined by the growth and dissolution of oxygen vacancy (V ö ) rich conductive filaments. In the phasefield model, we choose V ö density (N V€ O ) as the order parameter. The total free energy of the system (F total ) includes the synergistic contributions from the chemical, electrical, and mechanical effects, which is written as: where f chem , f electric , and f elastic represent the chemical energy density, the electrical energy density, and the elastic energy density, respectively. In this work, the metal oxide with oxygen vacancies is assumed as a dilute solution, such that the chemical energy density depends on the configuration entropy of oxygen vacancy and is described as: where k B is the Boltzmann constant, N o , N V€ O are the total number of lattice sites of oxygen atoms and number of oxygen vacancies, respectively. The electric energy density is described as: where φ is the electrical potential, ε 0 and ε r are the vacuum permittivity and relative permittivity, respectively. The elastic energy density is written as 51 : where C ijkl is the elastic constant tensor, ε ij is the total strain, and ε 0 ij is the local eigenstrain induced by the variation of oxygen vacancy density. In this model, we assume that ε ij = 0. The local eigenstrain ε 0 ij induced by the variation of oxygen vacancy density is determined based on the converse Vegard's law, which is written as: where is the Vegard strain coefficient, which describes the lattice parameter (a) change with the oxygen vacancy density (N V€ O ), δ ij is the Kronecker operator; ΔN V€ O ¼ ðN V€ O À N V€ Oi Þ is the variation of the oxygen vacancy density. The constant value of N V€ Oi is the oxygen vacancy density in the reference stress-free state.
The chemical, electrical and elastic potentials of the system are obtained by taking the variational derivative of the corresponding free energy densities with respect to N V€ Oi , i.e., μ chem Thus, the total potential of the system is described as: The flux of oxygen vacancy is proportional to the gradient of the potential: where the diffusion coefficient D is given by Here D 0 is the preexponential factor of diffusivity and E A is the activation energy of migration. In this work, D 0 = 2 × 10 −3 cm 2 s −1 and E A = 1 eV are chosen 52 .
The oxygen vacancy transport process can be described by the Nernst-Planck equation. To investigate the defect-induced mechanical effect on the resistive switching process, we either include or exclude the elastic potential in the Nernst-Planck equation, as shown in Eqs. (13) and (14): Equations (13) and (14) are coupled with the current continuity equation for electrical conduction (Eq. (15)) and the thermal transport equation for Joule heating (Eq. (16)), where σ and k th are the electrical and thermal conductivity, respectively. Equations (13)-(16) are self-consistently solved to obtain the oxygen vacancies density N V€ O , the electrical potential φ and the temperature T using finite element method based on the platform of COMSOL Multiphysics 5.4. The total system size is 35 × 20 nm 2 and the mesh size is selected as 0.5 nm. The two electrodes are assumed to act as heat sinks with fixed temperature T = 300 K. We consider that both electrodes are blocking the oxygen transfer, i.e., there is no oxygen vacancy flux at the interfaces between the oxide layer and the electrodes. Also the oxide/ electrode interfaces are assumed to be ohmic contact, where the electrons transport barrier is very small and can be ignored 5,[28][29][30]53 . Therefore we simulate the resistance change by solving the current continuity equation (Eq. (15)), assuming that electrons are moving much faster than oxygen vacancies and are always able to compensate the excessive charges induced by the local segregation and depletion of the oxygen vacancies. Since electron transport is not directly simulated, we ignore the contact resistance at the resistive switching layer (RSL)/electrode interface, and only consider the intrinsic resistance and the redox states inside the RSL. This also allows us to focus on the main driving forces on the oxygen vacancy transport, and to study the effects of electrical and thermal conductivities of difference materials on the resistive switching properties of metal oxide-based RRAM.
To solve Eqs. (15) and (16), a model for electrical and thermal conductivities is needed. For most metal oxide-based RRAM, the CF is assumed to consist of a metallic phase with locally enhanced oxygen vacancy density, which acts as n-type donors 10,54,55 . Meanwhile electrons are of much higher mobility than oxygen vacancies and are assumed to be always compensating the excessive charges induced by the local oxygen vacancy accumulation on the CF almost instantaneously. Since the electrical and thermal conductivity are strongly dependent on the electronic defect concentration (here it is mainly electrons), given the instantaneous response of electrons to the oxygen vacancy, it is reasonable to assume that the oxygen vacancy density N V€ O controls the local electrical and thermal conductivity. Therefore we assume that the electrical and thermal conductivity are linearly dependent on the oxygen vacancy density, based on similar descriptions from previous literature [28][29][30] . The electrical conductivity (σ) is given by the Arrhenius equation 28 , σ ¼ σ 0 e À E AC kT , where σ 0 is a pre-exponential factor and E AC is the activation energy for conduction. As shown in Fig. S1a, σ 0 is assumed to linearly increase with N V€ O with a slope of K 1 = 2.38. The E AC is chosen to be 0.05 eV for N V€ O ¼ 0, which is consistent with the activation energy in pure HfO 2 28 . The activation energy is zero for oxygen vacancy density above 0.3 × 10 27 m −3 , which corresponds to the metallic CF nature in set state. As shown in Fig. S1c, a linear dependence of k th on N V€ O with a slope of K 2 = 1.875 is assumed based on the Wiedemann-Franz Law. Although the linear approximations of σ 0 , E AC , and k th on the oxygen vacancy density might appear to be oversimplified, the calculation results based on these assumptions show good consistency with the experimental data, the detailed description of parameters in the model can be found in Supplementary Note 1.
Based on the above assumptions, K 1 and K 2 can be used to describe the change in electrical and thermal conductivity as a function of oxygen vacancy density in metal oxide. For high-throughput calculations, σ and k th remain constant for different defect-free oxides and are linearly proportional to N V€ O . Therefore the values of σ and k th in metallic state of high oxygen vacancy concentration (~1.2 × 10 27 m −3 ) for different metal oxides can be represented by the slopes (K 1 and K 2 ) of the linear relations. This enables us to obtain a dataset of metal oxide materials with different σ and k th combinations by simply tuning the magnitudes of K 1 and K 2 .

DFT calculations of Vegard strain coefficient
Density functional theory calculations of monoclinic (P2 1 /c) HfO 2 phase are performed by using the Vienna Ab initio Simulation Package (VASP) with PAW atomic potentials 12,56 to calculate the Vegard strain coefficient. Based on the earlier electronic-structure calculations 45 , it was determined to use the generalize gradient approximation (GGA) with the Perdew-Burke-Ernzerhof 57 functional and an additional Hubbard "U" potential of 2.2 eV applied to the "f" orbitals of Hafnia for treating the on-site Coulomb interaction, following the Dudarev's rotationally invariant approach 58 .
To obtain the Vegard strain at high concentrations of~8.3% (one oxygen-vacancy in a 12-atom unit-cell), the DFT calculations were calculated with a 6 × 6 × 6 Gamma-centered k-mesh to integrate the Brillion Zone with the atomic forces relaxed by using a conjugate-gradient algorithm till the forces on all atoms converged to at least 0.01 eV/Å. On the other hand, to obtain the Vegard strain at low concentrations, a 3 × 3 × 3 supercell was used to perform oxygen defect cluster calculations. A single k-point (the Gamma point) was used to integrate the Brillion Zone and the atomic forces were converged down to 0.1 eV/Å with concomitant relaxation of lattice parameters. At the same time, the lattice parameters were simultaneously optimized to obtain the zero-pressure volume with and without oxygen vacancy. The results are summarized in Supplementary Note 2.

Machine learning approach
Three material properties of metal oxides including Vegard strain coefficient V ij , electrical conductivity σ, and thermal conductivity k th are selected as the materials fingerprints to identify the properties of metal oxides-based RRAM (I on /I off as well as the t switch ). For building an interpretable machine learning features space, we start with a set of features Φ 0 (V ij , σ, k th ) as the primary features where V ij is dimensionless, σ is in units of 10 3 S cm −1 and k th is in units of W m −1 K −1 . Then, additional features are constructed by using operators in the set H m ≡ (+, −, ×, /log, ∧ −1, ∧ 2, ∧ 3, sqrt, cbrt, | − |), where the superscript 'm' indicates that dimensional analysis is performed only for meaningful combinations (i.e., physically allowed combinations). Thus, the feature-space has been extended by using a combination of H m operators with three-rungs operation (i.e., Φ 0 → Φ 1 → Φ 2 → Φ 3 ), The number of elements in Φ 1 , Φ 2 , Φ 3 is~27,~619 and~448,330, respectively. Since some of the properties could be correlated, such as the electrical conductivity (σ) and thermal conductivity (k th )) related via the Lorenz factor based on the Wiedemann-Franz law, we use the recently developed SISSO approach 46 to identify a descriptor that predicts the memristor properties. SIS allows dimensional reduction of the size of the feature-space, while SO is used to pinpoint the optimal n-dimensional descriptor 47,48 . To obtain a more predictable machine learning model, we excluded very few discrete data points with I on /I off > 10 (Fig. 4a). A total of 1835 sets of phase-field simulations with different combinations of (V ij , σ, k th ) are performed to obtain the training dataset. The SISSO is then used to predict both I on /I off and t switch with two analytical equations (Eqs. (1) and (2)). To test the robustness of the ML predictive model, we performed additional phase-field simulations to obtain datasets that are not included for training the ML model. The log(t switch ) and I on /I off for these additional datasets, which are calculated from the phase-field simulations and predicted based on the machine learning predictive models are listed in Supplementary Information Table S1. More details of machine learning methods and results are shown in Supplementary Note 3.

DATA AVAILABILITY
The raw data of the phase-field simulation in this paper and its supplemental information files, as well as the data used for machine-learning are available from corresponding authors (ye.cao@uta.edu and ganeshp@ornl.gov) upon reasonable request.