content_image/429699.jpg content_image/429697.jpg content_image/429696.jpg
the regions of the dielectric that are in physical proximity to the metallic electrodes. There is growing evidence that these are the regions where the RS takes place $^{13,17,20–23}$. 

Our model is schematically shown in Fig. 1 and consists of a single conductive channel within a more insulating dielectric, which is represented by a one dimensional re- sistive network on N links. The first and last N$_{I}$ links correspond to high resistance interfacial regions next to the external electrode, and the central N − 2 N$_{I}$ links de- scribe the bulk section. Each link is characterized by a certain concentration of oxygen vacancies, which deter- mines the resistivity of the link. They may be physically associated to small domains of nanoscopic dimensions which may actually correspond to grains of the polycrys- talline oxide. We take ρ$_{i}$ = A$_{α}$δ$_{i}$ , with α = B if i is in the bulk ( N$_{I}$ < i ≤ N − N$_{I}$ ), α = L if i is in the left interface ( i ≤ N$_{I}$ ) and α = R if i is in the right inter- face ( N − N$_{I}$ < i ≤ N ). In our study we set N = 100 and N$_{I}$ = 10. The following equation specifies how the vacancies diffuse through the network domains under an external voltage, 

It gives the probability for transfer of vacancies from do- main a to a nearest neighbor domain b . The probability 

<img src='content_image/429699.jpg'>

<img src='content_image/429697.jpg'>

is proportional to the concentration of vacancies present in domain a and to the concentration of “available va- cancy sites” at the target domain. The Arrhenius factor exp( − V$_{0}$ ), is controlled by a dimensionless constant, V$_{0}$ , related to the activation energy for vacancies diffusion. The important factor exp(Δ V$_{a}$ ) models the enhancement (or suppression) of the diffusive process due to the local electric field at domain a . 

From Eq.(1), a constant V$_{0}$ leads to an initially con- stant distribution δ$_{i}$ = δ$_{o}$ for all i . The value of the initial constant concentration, δ$_{o}$ must be much smaller than one, since it physically represents the concentration of defects (oxygen vacancies) within a domain. We adopt δ$_{o}$ = 10 − $^{4}$. 

Similarly to actual resistive switching experiments, we simulate the applied voltage protocol V ( t ) by lin- ear ramps that follow the sequence 0 → + V$_{max}$ → 0 → − V$_{max}$ → 0 (our convention is that the right electrode is grounded). The duration is of s time steps. The se- quence may be repeated a number of cycles n , for a to- tal duration τ$_{max}$ = ns . In our simulations we choose V$_{0}$/V$_{max}$ = 0 . 016, that provides a non-negligible but slow diffusive contribution to the evolution of δ$_{i}$ with respect to the total time duration of the simulations (ie, the to- tal number of time steps). We set V$_{max}$ = 1000 that provides for a sufficiently large electric stress. Our quali- tative results are rather robust with respect to the choice of model parameters, a detailed systematic study of their dependence is left for future work. 

The coefficients A$_{B}$ , A$_{R}$ and A$_{L}$ still remain to be spec- ified. With no loss of generality, we fix the value of the bulk coefficient to unity A$_{B}$ = 1 and leave the interfacial A$_{R}$ and A$_{L}$ free. In this initial study we shall concen- trate in the symmetric case, A$_{R}$ = A$_{L}$ , that corresponds to most common experimental devices. 

The numerical simulations are performed through the following steps: (i) at each simulation time step t (1 ≤ t ≤ τ$_{max}$ ) a given external voltage V ( t ) is applied to the resistive network (the electrodes are assumed perfect con- ductors). The current through the system is computed as I ( t ) = V ( t ) /R$_{T}$ , with R$_{T}$ the total (ie, two terminal) resistance 

where α = R, B, L denotes the three network regions, and c denotes an unessential geometrical constant related to the dimensions of the domains (which we set to unity). (ii) we compute the local voltage profile V$_{i}$ ( t ) = I ( t ) ρ$_{i}$ and the voltages drops Δ V$_{i}$ ( t ) = V$_{i}$$_{+1}$ ( t ) − V$_{i}$ ( t ). (iii) we use Eq. (1) to compute all the oxygen vacancy transfers between nearest neighboring domains, and update the values δ$_{i}$ ( t ) to a new set of concentrations δ$_{i}$ ( t + 1). (iv) we use these new values to recompute the current at τ$_{t}$$_{+1}$ under the applied voltage V ( t + 1) as indicated in the first step. 

<img src='content_image/429696.jpg'>

