The role of ion transport phenomena in memristive double barrier devices

In this work we report on the role of ion transport for the dynamic behavior of a double barrier quantum mechanical Al/Al2O3/NbxOy/Au memristive device based on numerical simulations in conjunction with experimental measurements. The device consists of an ultra-thin NbxOy solid state electrolyte between an Al2O3 tunnel barrier and a semiconductor metal interface at an Au electrode. It is shown that the device provides a number of interesting features such as an intrinsic current compliance, a relatively long retention time, and no need for an initialization step. Therefore, it is particularly attractive for applications in highly dense random access memories or neuromorphic mixed signal circuits. However, the underlying physical mechanisms of the resistive switching are still not completely understood yet. To investigate the interplay between the current transport mechanisms and the inner atomistic device structure a lumped element circuit model is consistently coupled with 3D kinetic Monte Carlo model for the ion transport. The simulation results indicate that the drift of charged point defects within the NbxOy is the key factor for the resistive switching behavior. It is shown in detail that the diffusion of oxygen modifies the local electronic interface states resulting in a change of the interface properties.

The research in the field of memristive devices dates back to the 1970s when Chua introduced his theory of memristors and memristive devices 1 . This idea has emerged a considerable amount of interest after 2008 when Strukov et al. linked their resistive switching device to Chua's theory 2 . Today, memristive (or synonymously resistive switching) devices have been identified as promising candidates for future non-volatile memory applications due to their distinct key features which are i) low power consumption, ii) passivity, and iii) scalability into the nanometer scale 3 . Beyond their applications as non-volatile memories, resistive switching devices turned out to be applicable as artificial synapses in neuromorphic circuits [4][5][6] .
Many of the available metal-oxide resistive switching devices rely on the stochastic phenomenon of creation and rupture of conductive metal filaments within a solid state electrolyte matrix. Inherent are initial electroforming procedures and distinct "on" and "off " states. The latter are in fact beneficial for memory applications. However, in the context of neuromorphic applications a wider dynamic range is required, rather than a purely "digital" behavior. To overcome these limitations various kinds of interface-based resistive switching devices have been developed. A profound review of the progress in the field with respect to materials, switching mechanisms, and performance has been provided by Pan et al. 7 . Predominantly the switching mechanisms rely either on the change of height of a tunnel barrier, and thus on the electron tunneling probability, or on the change of a Schottky contact [8][9][10][11] . Recently Hansen et al. presented a memristive device consisting of both, a tunnel barrier and a Schottky contact which embed a thin solid state electrolyte 12 . This device provides specific features such as a highly uniform current distribution, an intrinsic current compliance, an improved retention time, and no need for an initial electroforming procedure. These features make the device interesting for its application as artificial synapses in bio-inspired neural networks. However, the physical mechanisms which are responsible for the memristive behavior are still not completely understood. To explain the physics two hypotheses are proposed: The first assumes charging and de-charging of trap states within the solid state electrolyte and/or at the metal-semiconductor interface to be the main reason for the resistance change [13][14][15] . The second hypothesis is that the motion of charged point defects within a sufficiently high applied electric field controls the properties of the tunnel barrier as well as the Schottky barrier, and therefore the resistance of the device 16,17 .
The aim of this work is to provide a model which allows to investigate the role of the ion transport for the resistive switching of the double barrier memristive device. We utilize a 3D kinetic Monte Carlo code in order to describe the ion transport within the solid state electrolyte subject to both externally applied and Coulomb fields 18 . The kinetic Monte Carlo model is consistently incorporated into an efficient lumped element circuit model of the device. It is important to note that the simulation parameters are chosen in accordance with an experimental setup. We find that the ion transport is a key factor for the resistive switching of the device. Furthermore, we identify an adsorption mechanism of ions at the Au electrode to be a key factor for the retention characteristics of the device. The comparison of the simulation results with experimental findings shows excellent agreement.

Simulation Approach
The simulation scenario is shown in Fig. 1. It consists of a memristive device with the layer sequence Au/Nb x O y / Al 2 O 3 /Al. We assume the Nb x O y layer to act as an ionic/electronic mixed conductor. It represents a solid state electrolyte in which mobile oxygen ions (blue circles) drift under the influence of an externally applied electrical field. In order to guarantee for charge neutrality of the simulation domain, the negative oxygen ions are compensated by stationary positive ions (red circles) 19 . At the initial simulation stage the ions are randomly distributed within the Nb x O y assuming thermodynamic equilibrium. The Nb x O y /Au interface is assumed to be inert to prevent oxidation. It can be regarded as a Schottky contact. The Al 2 O 3 tunnel barriers are electrically high-quality barriers enabling elastic electron tunneling. They are nearly free of defects and highly stable. Both the tunnel barrier and the Au/Nb x O y interface represent chemical barriers which confine the ions within the Nb x O y layer. Details about the devices fabrication can be found in Hansen et al. 12 .
Before we describe the details of the model it might be helpful to firstly explain the overall idea of the simulation approach: The physical device is decomposed into distinct lumped elements which are connected in series, similar to an earlier approach 12 . The tunnel barrier is mimicked by a voltage controlled current source based on the famous Simmons formula 20 . The Nb x O y /Au interface is descibed by the Schottky contact model 21,22 . The solid state electrolyte is modeled by an ohmic resistance which depends on the inner atomic structure. As stated before, the transport of ions within the solid state electrolyte under the influence of the externally applied electric field and the Coulomb field due to the ions themselves is of particular importance. Therefore, the lumped element circuit model is consistently coupled with a 3D kinetic Monte Carlo (KMC) model for the ion transport in such a way that Kirchhoff 's voltage and current laws for the circuit model are satisfied at all instances of time. Within an inner loop the dynamical state of the solid state electrolyte is calculated subject to the electric field acting on the ions and subject to a set of appropriate boundary conditions. From this the ohmic resistance of and the voltage drop across the solid state electrolyte are calculated. Applying Kirchhoff 's current law which stated that the current through the device is constant, the individual voltage drops across the tunnel barrier and the Schottky contact as well as the current itself are calculated within an outer loop. Now we come to the details of the model: Using X-ray diffraction measurements the Nb x O y layer has found to be amorphous 12 . Due to the short range periodicity of amorphous material it is reasonable to model this layer by a 3D lattice with a lattice constant corresponding to the hopping distance of oxygen ions in Nb x O y 18,23 . We use lattice constants of 0.25 nm in z direction and and 0.33 nm in x and y direction. This is indicated in Fig. 1. It takes into account the dimensional lattice mismatch of Nb x O y . For an amorphous material it might be intuitive to chose equal lattice constants for all direction. However in order to have a better resolution in z direction, we apply a slightly reduced lattice constant, which is realistic and without any influence on the simulation results. The 3D simulation domain consists of a quadratic 9 nm × 9 nm basis area, a 2.5 nm thick solid state electrolyte and a 1.2 nm thick tunnel barrier (see Fig. 1).
For the KMC model which is described in detail in refs 18,24 and 25 three different hopping processes of the oxygen ions have been taken into account (see Fig. 1): i) diffusion within the solid state electrolyte, ii) adsorption Positively charged (red circles) and negatively charged (blue circles) defects are located within the solid state electrolyte. Positively charged defects are assumed to be immobile, negatively charged defects are allowed to move within the solid state electrolyte (i), adsorb at the Au electrode (ii), and desorb from the Au electrode (iii).
Scientific RepoRts | 6:35686 | DOI: 10.1038/srep35686 and iii) desorption at the respective interfaces of the Nb x O y layer. The related hopping rates k ij for the KMC model are given by an Arrhenius law, ν is the phonon frequency, T is the lattice temperature, k B is the Boltzmann constant, and E ij is the hopping (potential) barrier between the lattice sites i and j. For the hopping barrier we use represents a linear correction term of the hopping barrier due to an externally applied electric field and the local Coulomb field due to the ions themselves. The transport of ions within the solid state electrolyte evolves as follows: After inserting a certain number of stationary positive and mobile negative ions into the Nb x O y layer, the negative ions drift according to the local electric field which itself is a superposition of the field due to the externally applied voltage and the Coulomb field of the remaining ions in the system. When an oxygen ion reaches one of the interfaces a surface interaction takes place, depending on the local electric potential. The simulation parameters, which are mainly extracted from the experiment, are collected in Table 1.
The first step of the simulation loop is to calculate the electric potential within the Al 2 O 3 tunnel barrier. We assume that the Al 2 O 3 is a perfect insulator which contains no electric charge carriers. Therefore, the electric potential can be obtained from a simple Laplace equation subject to Dirichlet conditions, The permittivity of the Al 2 O 3 layer is assumed to be homogeneous. This leads of course to a spatially linear potential and a spatially constant electric field. With the voltage drop across the tunnel barrier the local electron tunneling current density can be calculated from the Simmons formula, Tunnel eff . The elementary charge, the free electron mass, and the Planck constant are given by q, m, and h respectively. d eff is the effective thickness of the tunnel barrier. Since the Simmons formula accounts only for elastic tunneling, the local ion concentration at the Al 2 O 3 is of particular importance for the resulting electron tunneling current. In order to take this into account the tunnel distance has been implemented as an effective distance d eff , which itself depends on the position of the negative ions within the Nb x O y layer, λ d is a fit parameter, d is the average relative change of the oxygen ion position from their inertial position, and d 0 is the initial effective thickness of the Al 2 O 3 layer (d 0 = 1.3 nm).
The electric potential within the Nb x O y layer cannot be obtained directly from Poisson's equation as the charge density within the layer is not known a priori. The assumption of an electric charge carrier-free insulator, as used for the tunnel barrier, is not valid in the Nb x O y layer due to an electric current which is present here. Instead, we calculate the local electric potential by solving the continuity equation ∇ ⋅ =  J 0 in conjunction with a constitutive law that couples locally the electric field E and the current density  J. The required relation between the electric field and the current density is given by a local Ohm's law, Here σ represents the electric conductivity of the Nb x O y layer, which is assumed to be homogeneous. It is important to note that the displacement current can be neglected because it scales with (L/Tc) 2 where c is the speed of light and L and T are the length and the time scales of the system. The displacement current is therefore small compared to the conduction current.
The electric field across the third device part, the Nb x O y /Au Schottky contact, depends mainly on the oxygen ion concentration close to the Nb x O y /Au interface. The ions influence the height of the Schottky barrier and the ideality factor. To calculate the electric field across the Schottky contact, the thermionic emission theory is employed. The Schottky contact is described by a set of analytical formulas 21,22 . Within this theory the Schottky diode current density is given by where k B and T are the Boltzmann constant and temperature, respectively. n is the ideality factor which describes the derivation from an ideal current. The reverse current J R at forward bias voltages is given by where Φ b is the Schottky barrier height and A * = 1.20173 × 10 6 A/(m 2 K 2 ) the effective Richardson constant. The reverse current in reset direction is dominated by the lowering of the Schottky barrier. If the apparent barrier height Φ b at the Schottky interface is significantly smaller than the conductive band gap of the insulator, the reverse current decreases gradually with the applied negative bias approximately as α r denotes a device dependent parameter which is used to describe the experimentally observed voltage dependence of the reverse current. Regarding (5)-(7) the effective current density is affected by both the ideality factor of the Schottky barrier n and the height of the Schottky barrier Φ b . The change of the ideality factor is calculated by taking the average distance of the oxygen ions to the Schottky barrier into account, n 0 λ n is a fit parameter and d is the average relative change of the oxygen ion position from their initial position. n 0 is the initial ideality factor which has been estimated to be 4.1 for the device 12 . The change of the Schottky barrier height is attributed to the spatial rearrangement of the mobile oxygen ions during external voltage applications, which leads to an image force adjustment. Taking this into account, the effective Schottky barrier height is calculated by with Φ b0 being the initial Schottky barrier height and Φ M being the average surface potential at the Au electrode. This involves in particular the Coulomb field of the ions within the solid state electrolyte which is superimposed to the local electric field produced by the external voltage. Therefore, for a complete description of the electric field the Coulomb field of the oxygen ions has to be taken into account. The Coulomb potential is calculated using Poisson's equation where the permittivity ε  r ( ) is either ε r ,Al O 2 3 or ε r ,Nb O x y , the relative permittivity of the Al 2 O 3 or the Nb x O y , respectively.

Results and Discussion
In Fig. 2 the calculated current-voltage characteristics (also referred to as I-V curve) is compared with a measured I-V curve. For both I-V curves the externally applied voltage is first ramped from 0 to 3 V using a sweep rate of 0.14 V/s in order to set the device resistance from the initial high resistance state (HRS) to a low resistance state (LRS). When the maximum voltage is reached, the sweep rate is switched to − 0.14 V/s until the applied voltage is zero again. Then a voltage sweep rate of − 0.1 V/s is applied until − 2 V is reached. The device is reseted back to the high resistance state. Now the voltage is increased with a sweep rate of 0.1 V/s back to 0 V. As it can be seen from Fig. 2 the calculated I-V curve is in a very good agreement with the experimentally obtained I-V curve. In line with the experimental data the most apparent feature of the memristive hysteresis is the asymmetry between positive and negative bias. This asymmetry can be attributed to the Nb x O y /Au Schottky contact. Furthermore, we are able to capture the nonlinear behavior of the I-V curve in which a high resistance at small voltages and a current saturation at higher voltages is obtained. The model shows moreover the gradual change of the device resistance by several orders of magnitudes. The fine structure for small measured currents at negative bias, which is due to the limitation of the current resolution of the experimental set up rather than physically relevant mechanism, is not captured by the model. To study the dynamics of oxygen ions within the solid state electrolyte under the influence of an external voltage and their impact on the device resistance the spatial ion distribution at four significant positions of the I-V curve are shown in a 3D view graph depicted in Fig. 3. At the initial state (Fig. 3a)) the negative and positive ions (blue and red markers, resp.) are concentrated in the center of the simulation box (marked by a vertical line representing the average position of the negative charged ions). By applying an external positive voltage at the left Au electrode the negative charged oxygen ions drift towards the Au interface (Fig. 3b)). When the mobile ions reach the Au interface they can adsorb with a certain probability, which prevents back diffusion. Their average position is then located at the Au interface (see vertical line of Fig. 3b)). In order to transport oxygen ions back from the Au interface a negative voltage is applied. It is shown in Fig. 3c) that for an external voltage of − 2 V the average negative charge position is close to the initial position, while for a complete resetting the external voltage has to be ramped back from − 2 V to 0 V (Fig. 3d)).
The numerical simulation indicates that the oxygen ion movement is strongly voltage dependent and therefore responsible for the gradual change of the device resistance. To study the voltage dependence of the resistance variation in some more detail, the change of the device resistance as a function of the applied voltage is depicted in Fig. 4a). For the data presented in Fig. 4a) the sweep rate is set to 0.14 V/s, while the device resistance in the LRS and HRS is determined at a voltage of 0.5 V. If the external voltage is ramped to voltages below 2 V almost no change of the device resistance can be observed. In contrast, if the external voltage is ramped up to 3.5 V the change of the resistance is of several orders of magnitude. Regarding the positions of the negatively charged mobile oxygen ions, as depicted in Fig. 4c) as 3D view graph, it can be clearly observed that they only drift in the case of a sufficiently high voltage. In particular, while for a voltage ramp up to 2 V the ion positions are nearly unaffected compared to the initial distribution. For a voltage ramp up to 3.5 V nearly all negative charged ions are located at the Au interface.
To get an idea of the physical mechanisms behind the strong voltage dependence, the voltage drops across the different device parts are given in Fig. 4b). These are the voltages across the tunnel barrier, the Nb x O y solid state electrolyte, and the Schottky contact. At low applied voltages the partial voltages across the Nb x O y layer and the tunnel layer (orange and green curves of Fig. 4b)) are almost zero since the current is blocked by the Schottky contact (blue curve of Fig. 4b)). In other words, the Schottky contact defines a threshold voltage for the device which has to be exceeded to change the resistance of the device. Here, the resulting electric field within the Nb x O y is too small to induce ionic diffusion. However, since the externally applied voltage exceeds 1.5 V, the voltage across the Nb x O y layer raises and excites voltage polarity directed ion transport. Moreover, taking Fig. 4b) into account, the reset process is mainly dominated by the voltage drop across the Schottky barrier. Due to the reverse direction of the diode like Schottky contact nearly the complete applied voltage is blocked by the Schottky contact during the reset process. The increasing voltage drop across the Schottky contact is able to initiate desorption of oxygen ions from the Au surface (E ij = 0.25 eV 26 ). Then the desorbed ions drift subject to the Coulomb potential of the remaining positive charges and subject to the concentration gradient back into the solid state electrolyte. This mechanism leads to a fast reset of the total device resistance.
The reason for the drastic resistance change of the device stems hypothetically from the variation of the two energy barriers which confine the solid state electrolyte 12 . These variations are induced by the change of the spatial oxygen concentration, as discussed above. This leads to a kind of interference of the both interfaces which embed the Nb x O y layer. In line with (4), (8), and (9) we assume for the model that the spatial distribution of negative charged oxygen ions in front of the barriers influences both the effective tunnel barrier length and the parameters of the Schottky barrier, such as the ideality factor and the barrier height. We find that the effective tunnel barrier width decreases from 1.3 nm to 1.2 nm if the voltage applied to the Au electrode raises from 0 V to 3.0 V. At the same time the ideality factor decreases from 4.0 to 3.4. The Schottky barrier height is lowered from 0.9 eV to 0.83 eV (see Fig. 5b)).
As given in (4) and (8) both the effective tunnel distance and the ideality factor are assumed to depend linearly on the average ion position, while the lowering of the Schottky barrier height depends on the Coulomb potential in front of the Au interface (see (9)). The variation of the Coulomb potential as a function of the applied external voltage is shown in Fig. 5a). For this study the local Coulomb potential is averaged along the x and y directions. We find in particular that the Coulomb potential at the Nb x O y /Au interface decreases from 0 mV to − 70 mV when ramping the external voltage from 0 V to 3.0 V. This means that the lowering of the Schottky barrier height (see Fig. 5b)) originates from the displacement of negatively charged ions towards the Au interface, which decreases the Coulomb potential at the interface.
In order to identify the rate determining mechanisms of the ion dynamics, the electric field within the Nb x O y layer is depicted in Fig. 6 at four distinct instances of time. To visualize the actual mechanism the electric field is projected onto the xz-plane. For t = 10 s the applied voltage is 0 V and the origin of the electric field are the positive and negative charges. This attractive Coulomb field has to be overcome to deflect charges from the stable position. Small setting voltages are blocked by the Schottky barrier (see Fig. 4a)). However, for setting voltages larger than approx. 1.5 V the voltage drop across the solid state electrolyte starts to increase. For an applied voltage of 2 V a small preferred orientation of the electric field in positive z direction can be observed (see Fig. 6b)) although the random structure of the Coulomb field can still be seen. This preferred orientation and the strength of the electric field increases for higher applied voltages up to 3 V (see Fig. 6c)). A positive voltage at the Au electrode on the device leads to a homogenization of the electric field. Although between 2 V and 3 V at certain positions the Coulomb field is still dominant and prevents ionic motion. The electric field due to the applied voltage is dominant so that an ionic drift occurs.
These observations indicate that the maximum applied voltage has a strong influence on the shape of the I-V curve in set direction. In particular, since for voltages smaller than approx. 1.5 V the directed ionic motion is suppressed. Furthermore, the hysteresis of the I-V curve is expected to disappear for maximum applied voltages smaller than approx. 1.5 V. To verify this hypothesis the positive branch of simulated and measured I-V curves for three different maximum applied voltages (1.8 V, 2.3 V and 3 V) are shown in Fig. 7. The simulated and measured I-V curves are obtained under comparable conditions. Both I-V curves show a nearly vanishing hysteresis for a maximum applied voltage of 1.8 V. However, for higher maximum applied voltages the shape of the hysteresis becomes significantly broader. These findings can be explained by ionic motion rather than by charging and de-charging of trap states. Therefore our results support ionic motion as a key feature for resistance change within the memristive double barrier device. In order to investigate the retention characteristics of the device, it is set to the LRS which is the initial state for this study. The change of the device resistance is then mapped over time. To track the evolution the resistance of the device, voltage pulses of 0.5 V are applied every 60 s in order to readout the resistance. Figure 8 shows the simulation results as well as experimental results for the same procedure. The results show a long resistance retention up to days without an externally applied voltage. Since almost the complete voltage during the reset process is blocked by the Schottky contact, the diffusion of the negative oxygen ions is similarly with or without an applied reset voltage. However, the desorption mechanism strongly depends on the Schottky contact and therefore on the applied reset voltage. The significantly different time scales of the change of the device resistance with or without an applied reset voltage indicate that the desorption mechanism is in fact the time limiting process of the resistance retention rather than the back diffusion of free oxygen ions into the bulk of the device. A fast resistance change during the first 60 s (see Fig. 8) which can be found in the experimental results but not in the simulation results is explained by the fast de-charging of trap states at the Nb x O y /Au interface. This mechanism is actually not included in the simulation 12 . Most probably both mechanisms, charge trapping and ion drift are present in the devices but leading to retention times on different time scales. It is important to note that the local electric field  and the resulting Coulomb potentials have to be analyzed in order to determine the rate determining mechanisms of the ion dynamics within the numerical investigations.

Conclusions
We report on numerical simulations of a quantum mechanical double barrier memristive device. The device consists of an Al/Al 2 O 3 /Nb x O y /Au structure where the Nb x O y solid state electrolyte is enclosed by a tunnel barrier and a Schottky contact. We present results from a kinetic Monte Carlo simulation of the (atomistic) ion transport, consistently coupled to a lumped circuit element model for the current through the device. The current voltage characteristics which is an important indicator for the dynamics of the device is numerically calculated and compared to experimentally obtained results. We find very good agreement between simulation and experiment. We claim that the model which is proposed captures most relevant physical processes for resistive switching on both, atomistic length scales and experimental time scales. The simulation results identify the transport of charged point defects within the solid state electrolyte as a key mechanism for the resistive switching of the double barrier device, although the charging and de-charging of trap states within the solid state electrolyte and/or at the metal semiconductor interface cannot be completely excluded from the physical picture. The results indicate that the local electric field is the rate determining quantity, particularly during the set process. Finally, we find that the threshold voltage for resistive switching in set direction and an adsorption mechanism is the reason for the long retention time. It is found that the Coulomb field and the concentration gradient are the main reasons for back diffusion of charged point defects. The calculated retention time is in excellent agreement with experimental findings. We believe that with the model we provide a powerful simulation tool that can be used to gain a deeper physical understanding of the resistive switching mechanisms of the double barrier memristive device.

Methods
Sample preparation. Memristive tunneling junctions are fabricated on 4-inch Si wafers with 400 nm of SiO 2 (thermally oxidized) using a standard optical lithography process. The junctions are arranged in 1 mm × 1 mm cells across the wafer, containing 6 different contact sizes ranging from 70 μm 2 to 2300 μm 2 . The devices are fabricated using the following procedure: First, the multilayer (including top-and bottom-electrode) is deposited without breaking the vacuum using DC magnetron sputtering. The Al 2 O 3 tunnel barrier is fabricated by depositing Al which was afterwards partially oxidized in-situ. The Nb x O y layer is deposited by reactive sputtering in an O 2 /Ar-atmosphere. Following the subsequent lift-off, the junction area is defined by etching the top electrode using wet etching (potassium iodide, for Au) and dry etching (reactive ion etching with SF 6 , for Nb). The etched parts are then covered with thermally evaporated SiO to insulate the bottom electrode from the subsequently deposited Nb-wiring to contact the top electrode. I-V measurements. All measurements are performed using an Agilent E5260 source measurement unit.
Current-voltage measurements are obtained by sweeping the applied voltage and measuring the device current simultaneously.