## Abstract

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 MO_{2−x} systems such as HfO_{2−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 next-generation 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 to ~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} ~ 10^{−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 oxide-based 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_{\ddot O}} = 1.2 \times 10^{27}\,{\mathrm{m}}^{ - 3}\) and vanish in the remaining region in the initial state. Then a triangular voltage sweep with rate d*V*/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_{\ddot 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 oxide-based RRAM. The local strain \(\varepsilon _{ij}^0\) induced by the change of *V*_{ö} density can be calculated as \(\varepsilon _{ij}^0 = V_{ij}{\Delta}N_{V_{\ddot O}}\delta _{ij}\), where *V*_{ij} is the Vegard strain coefficient or chemical expansion coefficient, *δ*_{ij} is the Kronecker operator, \({\Delta}N_{V_{\ddot O}} = \left( {N_{V_{\ddot O}} - N_{V_{\ddot Oi}}} \right)\) is the variation of oxygen vacancy density. The constant value \(N_{V_{\ddot 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\varepsilon _{ij}^0/dz\) and \(d\varepsilon _{ij}^0/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_{\ddot 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_{\ddot 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} ~15 at *K*_{1} = 7.5 shown in Fig. 3e. To further understand this behavior, we plot the \(N_{V_{\ddot 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.

### 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 \(\left| {V_{ij}} \right|\)), 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} = {\mathrm{18}}.{\mathrm{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. 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 oxide-based 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) and (2).

From Eq. (1), by choosing *V*_{ij} = 0–0.24 which covers most metal oxide, log(*t*_{switch}) monotonously increases with increasing \(\frac{{k_{th}}}{\sigma }\) (\(\frac{{{\mathrm{W}}\,{{m^{-1}}K^{-1}}}}{{{\mathrm{10}}^3{\mathrm{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 \(\frac{{k_{th}}}{\sigma }\) (green lines in Supplementary Fig. S6). Equation (2) reveals that *I*_{on}/*I*_{off} generally decreases with increasing \(\frac{{k_{th}}}{\sigma }\). At very small \(\frac{{k_{th}}}{\sigma }\)(<2), *I*_{on}/*I*_{off} increases with increasing \(\frac{{k_{th}}}{\sigma }\) (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 \(\frac{{k_{th}}}{\sigma }\). 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 \(\frac{{k_{th}}}{\sigma }\) 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., \(\frac{{k_{{\mathrm{th}}}}}{\sigma } \sim LT\), where *L* is the Lorenz number. Theoretically, \(\frac{{k_{{\mathrm{th}}}}}{\sigma }\) 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. High-throughput 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 \(\frac{{k_{th}}}{\sigma }\) 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.

## Methods

### 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 phase-field model, we choose *V*_{ö} density (\(N_{V_{\ddot 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_{\ddot 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 \(\varepsilon _{ij}^0\) is the local eigenstrain induced by the variation of oxygen vacancy density. In this model, we assume that *ε*_{ij} = 0. The local eigenstrain \(\varepsilon _{ij}^0\) induced by the variation of oxygen vacancy density is determined based on the converse Vegard’s law, which is written as:

where \(V_{ij} = \frac{1}{a} \cdot \frac{{da}}{{dN_{V_{\ddot o}}}}\) is the Vegard strain coefficient, which describes the lattice parameter (*a*) change with the oxygen vacancy density (\(N_{V_{\ddot O}}\)), *δ*_{ij} is the Kronecker operator; \({\Delta}N_{V_{\ddot O}} = (N_{V_{\ddot O}} - N_{V_{\ddot Oi}})\) is the variation of the oxygen vacancy density. The constant value of \(N_{V_{\ddot 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_{\ddot Oi}}\), i.e.,

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 \(D = D_0e^{ - \frac{{E_A}}{{kT}}}\). Here *D*_{0} is the pre-exponential 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_{\ddot O}}\), the electrical potential *φ* and the temperature *T* using finite element method based on the platform of COMSOL Multi-physics 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_{\ddot 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}, \(\sigma = \sigma _0e^{ - \frac{{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_{\ddot O}}\)with a slope of *K*_{1} = 2.38. The *E*_{AC} is chosen to be 0.05 eV for \(N_{V_{\ddot 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_{\ddot 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 over-simplified, 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_{\ddot 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.

## References

- 1.
Yang, J. J., Strukov, D. B. & Stewart, D. R. Memristive devices for computing.

*Nat. Nanotechnol.***8**, 13–24 (2013). - 2.
Waser, R. & Aono, M. Nanoionics-based resistive switching memories.

*Nat. Mater.***6**, 833–840 (2007). - 3.
Ielmini, D. et al. Resistive switching memories based on metal oxides: mechanisms, reliability and scaling.

*Semiconduct. Sci. Technol***31**, 063002 (2016). - 4.
Waser, R., Dittmann, R., Staikov, G. & Szot, K. Redox-based resistive switching memories—nanoionic mechanisms, prospects, and challenges.

*Adv. Mater.***21**, 2632–2663 (2009). - 5.
Nardi, F., Larentis, S., Balatti, S., Gilmer, D. C. & Ielmini, D. Resistive switching by voltage-driven ion migration in bipolar RRAM-part I: experimental study.

*IEEE Trans. Electron Devices***59**, 2461–2467 (2012). - 6.
Kozicki, M. N., Park, M. & Mitkova, M. Nanoscale memory elements based on solid-State electrolytes.

*IEEE Trans. Nanotechnol.***4**, 331–338 (2005). - 7.
Guo, X., Schindler, C., Menzel, S. & Waser, R. Understanding the switching-off mechanism in Ag

^{+}migration based resistively switching model systems.*Appl. Phys. Lett.***91**, 133513 (2007). - 8.
Jeong, D. S., Schroeder, H., Breuer, U. & Waser, R. Characteristic electroforming behavior in Pt/TiO

_{2}/Pt resistive switching cells depending on atmosphere.*J. Appl. Phys.***104**, 123716 (2008). - 9.
Ielmini, D. et al. Scaling analysis of submicrometer nickel-oxide-based resistive switching memory devices.

*J. Appl. Phys.***109**, 034506 (2011). - 10.
Ielmini, D., Nardi, F. & Cagli, C. Physical models of size-dependent nanofilament formation and rupture in NiO resistive switching memories.

*Nanotechnology***22**, 254022 (2011). - 11.
Bersuker, G. et al. Metal oxide resistive memory switching mechanism based on conductive filament properties.

*J. Appl. Phys.***110**, 124518 (2011). - 12.
Zhuang, H. L., Ganesh, P., Cooper, V. R., Xu, H. & Kent, P. R. C. Understanding the interactions between oxygen vacancies at SrTiO

_{3}(001) surfaces.*Phys. Rev. B***90**, 064106 (2014). - 13.
Yang, Y. C. et al. Observation of conducting filament growth in nanoscale resistive memories.

*Nat. Commun.***3**, 732 (2012). - 14.
Yang, Y. et al. Electrochemical dynamics of nanoscale metallic inclusions in dielectrics.

*Nat. Commun.***5**, 4232 (2014). - 15.
Kwon, D. H. et al. Atomic structure of conducting nanofilaments in TiO

_{2}resistive switching memory.*Nat. Nanotechnol.***5**, 148–153 (2010). - 16.
Chen, J. Y. et al. Dynamic evolution of conducting nanofilament in resistive switching memories.

*Nano Lett.***13**, 3671–3677 (2013). - 17.
Chen, J. Y., Huang, C. W., Chiu, C. H., Huang, Y. T. & Wu, W. W. Switching kinetic of VCM-based memristor: evolution and positioning of nanofilament.

*Adv. Mater.***27**, 5028–5033 (2015). - 18.
Celano, U. et al. Three-dimensional observation of the conductive filament in nanoscaled resistive memory devices.

*Nano Lett.***14**, 2401–2406 (2014). - 19.
Celano, U. et al. Imaging the three-dimensional conductive channel in filamentary-based oxide resistive switching memory.

*Nano Lett.***15**, 7970–7975 (2015). - 20.
Wang, J. Y., Li, L. Z., Huyan, H. X., Pan, X. Q. & Nonnenmann, S. S. Highly uniform resistive switching in HfO

_{2}films embedded with ordered metal nanoisland arrays.*Adv. Funct. Mater.***29**, 1808430 (2019). - 21.
Sankaran, K. et al. Modeling of copper diffusion in amorphous aluminum oxide in CBRAM stack.

*ECS Transactions***45**, 317–330 (2012). - 22.
Kamiya, K. et al. ON-OFF switching mechanism of resistive-random-access-memories based on the formation and disruption of oxygen vacancy conducting channels.

*Appl. Phys. Lett.***100**, 073502 (2012). - 23.
Clima, S. et al. First-principles simulation of oxygen diffusion in HfOx: role in the resistive switching mechanism.

*Appl. Phys. Lett.***100**, 133102 (2012). - 24.
Jiang, H. & Stewart, D. A. Enhanced oxygen vacancy diffusion in Ta

_{2}O_{5}resistive memory devices due to infinitely adaptive crystal structure.*J. Appl. Phys.***119**, 134502 (2016). - 25.
Bondi, R. J., Fox, B. P. & Marinella, M. J. Role of atomistic structure in the stochastic nature of conductivity in substoichiometric tantalum pentoxide.

*J. Appl. Phys.***119**, 124101 (2016). - 26.
Onofrio, N., Guzman, D. & Strachan, A. Atomic origin of ultrafast resistance switching in nanoscale electrometallization cells.

*Nat. Mater.***14**, 440–446 (2015). - 27.
Onofrio, N., Guzman, D. & Strachan, A. Atomistic simulations of electrochemical metallization cells: mechanisms of ultra-fast resistance switching in nanoscale devices.

*Nanoscale***8**, 14037–14047 (2016). - 28.
Larentis, S., Nardi, F., Balatti, S., Gilmer, D. C. & Ielmini, D. Resistive switching by voltage-driven ion migration in bipolar RRAM-part II: modeling.

*IEEE Trans. Electron Devices***59**, 2468–2475 (2012). - 29.
Kim, S. et al. Physical electro-thermal model of resistive switching in bi-layered resistance-change memory.

*Sci. Rep.***3**, 1680 (2013). - 30.
Kim, S., Choi, S. & Lu, W. Comprehensive physical model of dynamic resistive switching in an oxide memristor.

*Acs Nano***8**, 2369–2376 (2014). - 31.
Menzel, S. et al. Origin of the ultra-nonlinear switching kinetics in oxide-based resistive switches.

*Adv. Funct. Mater.***21**, 4487–4492 (2011). - 32.
Marchewka, A. et al. Nanoionic resistive switching memories: on the physical nature of the dynamic reset process. Adv. Electron.

*Mater.***2**, 1500233 (2016). - 33.
Ambrogio, S., Balatti, S., Choi, S. & Ielmini, D. Impact of the mechanical stress on switching characteristics of electrochemical resistive memory.

*Adv. Mater.***26**, 3885–3892 (2014). - 34.
Nicholas, J. D., Qi, Y., Bishop, S. R. & Mukherjee, P. P. Introduction to mechano-electro-chemical coupling in energy related materials and devices.

*J. Electrochem. Soc.***161**, Y11–Y12 (2014). - 35.
Billah, M. M., Hasan, M. M. & Jang, J. Effect of tensile and compressive bending stress on electrical performance of flexible a-IGZO TFTs.

*IEEE Electron Device Lett.***38**, 890–893 (2017). - 36.
Korobko, R. et al. Giant electrostriction in Gd-doped ceria.

*Adv. Mater.***24**, 5857–5861 (2012). - 37.
Schmitt, R., Spring, J., Korobko, R. & Rupp, J. L. M. Design of oxygen vacancy configuration for memristive systems.

*ACS Nano***11**, 8881–8891 (2017). - 38.
Esposito, V. & Traversa, E. Design of electroceramics for solid oxides fuel cell applications: playing with ceria.

*J. Am. Ceram. Soc.***91**, 1037–1051 (2008). - 39.
Ghicov, A. & Schmuki, P. Self-ordering electrochemistry: a review on growth and functionality of TiO

_{2}nanotubes and other self-aligned MO(x) structures.*Chem. Commun.***28**, 2791–2808 (2009). - 40.
Lee, M. J. et al. A fast, high-endurance and scalable non-volatile memory device made from asymmetric Ta

_{2}O_{5-x}/TaO_{2-x}bilayer structures.*Nat. Mater.***10**, 625–630 (2011). - 41.
Balachandran, J., Lin, L., Anchell, J. S., Bridges, C. A. & Ganesh, P. Defect genome of cubic perovskites for fuel cell applications.

*J. Phys. Chem. C.***121**, 26637–26647 (2017). - 42.
Hu, G., Fung, V., Sang, X., Unocic, R. R. & Ganesh, P. Predicting synthesizable multi-functional edge reconstructions in two-dimensional transition metal dichalcogenides.

*npj Comput. Mater.***6**, 44 (2020). - 43.
Shen, Z. H. et al. Phase-field modeling and machine learning of electric-thermal-mechanical breakdown of polymer-based dielectrics.

*Nat. Commun.***10**, 1843 (2019). - 44.
Zhang, Y. et al. Unsupervised discovery of solid-state lithium ion conductors.

*Nat. Commun.***10**, 5260 (2019). - 45.
Chimata, R., Shin, H., Benali, A. & Heinonen, O. Defect energetics of cubic hafnia from quantum Monte Carlo simulations.

*Phys. Rev. Mater.***3**, 075005 (2019). - 46.
Ouyang, R., Curtarolo, S., Ahmetcik, E., Scheffler, M. & Ghiringhelli, L. M. SISSO: a compressed-sensing method for identifying the best low-dimensional descriptor in an immensity of offered candidates.

*Phys. Rev. Mater.***2**, 083802 (2018). - 47.
JC & J. L., F. Sure independence screening for ultrahigh dimensional feature space.

*J. R. Stat. Soc. B***70**, 849–911 (2008). - 48.
Fan, J. Q. et al. Ultrahigh dimensional feature selection: beyond the linear model.

*J. Mach. Learn Res.***10**, 2013–2038 (2009). - 49.
Lee, S. et al. Anomalously low electronic thermal conductivity in metallic vanadium dioxide.

*Science***355**, 371–374 (2017). - 50.
Ganesh, P. et al. Doping a bad metal: origin of suppression of metal-insulator transition in non-stoichiometric VO

_{2}.*Phys. Rev. B***101**, 155129 (2020). - 51.
Cao, Y., Morozovska, A. & Kalinin, S. V. Pressure-induced switching in ferroelectrics: Phase-field modeling, electrochemistry, flexoelectric effect, and bulk vacancy dynamics.

*Phys. Rev. B***96**, 184109 (2017). - 52.
Ielmini, D., Nardi, F. & Cagli, C. Universal reset characteristics of unipolar and bipolar metal-oxide RRAM.

*IEEE Trans. Electron Devices***58**, 3246–3253 (2011). - 53.
Ielmini, D. E. A. Physical modeling of voltage-driven reistive swithcing in oxide RRAM. 2012 IIRW final report, 9–15.

- 54.
Zhang, C., Liu, G., Geng, X., Wu, K. & Debliquy, M. Metal oxide semiconductors with highly concentrated oxygen vacancies for gas sensing materials: a review.

*Sens. Actuators A: Phys.***309**, 112026 (2020). - 55.
Zheng, X. D. The influence of ion implantation-induced oxygen vacancy on electrical conductivity of WO

_{3}thin films.*Vacuum***165**, 46–50 (2019). - 56.
Blochl, P. E. Projector augmented-wave method.

*Phys. Rev. B Condens Matter***50**, 17953–17979 (1994). - 57.
Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple.

*Phys. Rev. Lett.***77**, 3865–3868 (1996). - 58.
Dudarev, S. L., Botton, G. A., Savrasov, S. Y., Humphreys, C. J. & Sutton, A. P. Electron-energy-loss spectra and the structural stability of nickel oxide: an LSDA+U study.

*Phys. Rev. B***57**, 1505–1509 (1998).

## Acknowledgements

Y.C. acknowledges supports from the Faculty Science and Technology Acquisition and Retention (STARs) Program in the University of Texas System, and the startup funding from the University of Texas at Arlington. K.Z. and Y.C. acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper (http://www.tacc.utexas.edu). P.G. was supported by the Center for Nanophase Materials Sciences, which is a DOE Office of Science User Facility. J.-J.W. acknowledges the partial support from the Army Research Office under grant number W911NF-17-1-0462. L.Q. Chen acknowledges partial support from the Computational Materials Sciences Program funded by the US Department of Energy, Office of Science, Basic Energy Sciences, under Award Number DE-SC0020145. J.-J.W. and L.-Q.C. also acknowledges the partial support from the Donald W. Hamer Foundation through the Hamer Professorship at Penn State. Y.H.H. acknowledges support from National Natural Science Foundation of China under grant number 51802280. This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

## Author information

### Affiliations

### Contributions

Y.C, J.W., and P.G. conceived and designed the study. P.G. implemented ab initio approaches to benchmark certain phase-field parameters and applied machine-learning methods to high-throughput phase-field datasets. K.Z. developed the phase-field models and performed all phase-field simulations under the guidance of Y.C., J.W. and Y.C. supervises the progress of the entire work. All authors contributed to interpreting the results and writing of the paper.

### Corresponding authors

## Ethics declarations

### Competing interests

The authors declare no competing interests.

## Additional information

**Publisher’s note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Supplementary information

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

## About this article

### Cite this article

Zhang, K., Wang, J., Huang, Y. *et al.* High-throughput phase-field simulations and machine learning of resistive switching in resistive random-access memory.
*npj Comput Mater* **6, **198 (2020). https://doi.org/10.1038/s41524-020-00455-8

Received:

Accepted:

Published: