FIG. 4: Top (semi-log scale): Snapshots of the density concen- tration profiles for the symmetric configuration, normalized to the initial uniform density value δ$_{i}$/δ$_{o}$ , during the hysteresis cycle (first half). Bottom: The corresponding profiles of local resistance ρ$_{i}$ for the same snapshots. Inset: Position of the snapshots in the (first half) hysteresis loop. Once in the bulk, the contribution to the total resistance of these migrating vacancies is reduced ( A$_{B}$ << A$_{L}$ ), thus the system reaches the R LO T state (leg of the table). Notice that in this state, the largest fields occur at the boundaries of the bulk and the interfacial regions. Thus, as V is further increased, the left interfacial region de- pletes further and therefore the voltage drop gets lower there. In contrast, on the right side the electric fields are further enhanced and there is now a migration of vacan- cies from the accumulation peak of the bulk towards the right interfacial region. This leads to an increase of the total resistance and, at the maximal voltage V/V$_{max}$ =1, we find that the vacancies have entered the interfacial region and already piled-up at the right end. Thus, the system is back to R HI T . The voltage protocol continues with the decrease of V back to zero but keeping the same (positive) polarity. Therefore no significant change in the δ$_{i}$ and ρ$_{i}$ profiles occurs, and half of the table with legs is already formed. When the negative polarity part of the cycle begins a similar analysis follows, since the distribu- tion of vacancy concentrations is a mirror image of the initial one. This forms the other half of the table. To complete our study we consider the effect of intro- ducing an asymmetry in the model parameters A$_{R}$ and A$_{L}$ . In Fig. 5 we display the results for the hystere- sis loops upon increasing the asymmetry. The results of the simulations show the gradual evolution of the “table with legs” towards a conventional rectangular hysteresis cycle $^{24}$. Asymmetry in experimental samples can either be due to the fabrication process, as for instance, in the obvious case of using different type of metal for the elec- trodes. Though our samples were fabricated in a sym- metric configuration, we induced a significant asymme- try by vigorous pulsing with a given polarity (ie, poling) before measuring the resistance hysteresis characteristic. Significantly, as shown in Fig. 5, we find that the ex- perimental data obtained in an asymmetric manganite sample and the cuprate, induced by intensive same po- larity pulsing, are in good qualitative agreement with our simulations. These results provides additional validation to our model. side of the interfacial/bulk boundary. Significantly this accumulation, as seen in Fig. 3, is quite large and narrow leading to a substantial increase in the local resistance. This feature will be a key to understand the origin of the legs of the hysteresis loop, which we shall consider next. In Fig. 4 we show snapshots, during the first half of the voltage cycle, of the vacancy concentration and lo- cal resistance profiles in the interfacial regions and their boundaries with the bulk. These are the active regions of the system where the electric-field-enhanced migration takes place. Notice that since the current I ( t ) is uniform along the conductive path, the local electric fields are di- rectly proportional to the local resistance. Let’s start at V =0, from the state at the beginning of a cycle. Interest- ingly, the ρ$_{i}$ profile indicates significant electric fields in both interfacial regions which extends into the neighbor- ing bulk. Thus the depletion of vacancies in the inter- faces and the accumulation at the boundaries are such that compensate for the difference between the respec- tive A coefficients, to yield similar electric fields across the boundaries of the regions. Yet, at the very left end of the system there is a small pile-up of vacancies which, as soon as the voltage is ramped up, will translate into the largest local fields and initiate the ionic migration. At V/V$_{max}$ =0.25, we observe a snapshot of the migration of the vacancies across the left interfacial region towards the bulk. From the resistance expression Eq. (2), so long the vacancies remain within the interfacial region, the system remains in R HI T state. At V/V$_{max}$ =0.62, the va- cancies have moved out the left interfacial region and entered the bulk, where their migration suddenly stops. To conclude, our results put on solid theoretical grounds the key role played by oxygen vacancies in the mechanism of resistive switching in TMO. They also pro- vide valuable insights, predicting a non-trivial spatial profile of the oxygen vacancy distribution which may be of help for device design. An exciting idea for future work is to explore the possibility of using atomistic, first principles, calculations to study properties of electrode - transition metal oxide interfaces to estimate the param- eters of the model and provide guidance in the material choice for actual memory devices. IV. CONCLUSIONS 4