2 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 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 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. 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 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. R$_{T}$ = c N ∑ i =1 ρ$_{i}$ = ∑ α ∑ i ∈ α A$_{α}$δ$_{i}$, (2) p$_{ab}$ = δ$_{a}$ (1 − δ$_{b}$ ) exp( − V$_{0}$ + Δ V$_{a}$ ) . (1) FIG. 1: a) Schematic model with a single conductive channel within the dielectric. The three regions L , R and B corre- spond to the two high resistance interfaces and the more con- ductive central bulk, respectively. The small boxes indicate the domains. b) Detailed scheme of the conductive path. The grayscale qualitatively depicts the variation in the oxygen va- cancy concentration through the channel (darker corresponds to higher concentration). The top figure shows the initial state with uniformly distributed oxygen vacancies, and the bottom one shows the inhomogeneous distribution after the first few “forming” voltage cycles (see text).