Finite-element simulation of interfacial resistive switching by Schottky barrier height modulation

This study demonstrates a numerical model for interfacial switching memristors based on the Schottky barrier height modulation mechanism. A resistive Schottky contact is formed for an n-type semiconductor and a high work-function metal (e.g., strontium titanate and platinum). The contact resistance is determined by the Schottky barrier height, which is influenced by the concentration of oxygen vacancies serving as space charges. Accordingly, the spatial distribution of vacancies and cell conductance can be controlled by applying a bias voltage. This interfacial switching is advantageous over filamentary switching, owing to the conductance change being more gradual in interfacial switching. In this study, a two-step numerical analysis was performed to model the conductance change in an interfacial switching memristor having a metal–oxide–metal structure of Pt/SrTiO3/Nb-SrTiO3, where Pt and SrTiO3 form a Schottky contact. In the first step, the change in the spatial distribution of vacancies by an applied switching voltage was obtained by solving the drift and diffusion equations for vacancies. In the second step, after setting the Schottky barrier height according to the vacancy concentration near the contact, the cell conductance was obtained by calculating the current value by applying a small read voltage. Consequently, our simulation successfully reproduced the experimental results for the SrTiO3-based memristor. Through this study, our device simulation for interfacial switching was successfully established, and it can be utilized in the computational design of various device architectures.


Introduction
In recent years, resistance-based nonvolatile memory devices have emerged rapidly and have exhibited enormous growth in the field of resistive random-access memory (ReRAM or RRAM).The growing demands in the semiconductor industry necessitate computations beyond the standards of CMOS-based platforms and new technologies for memory.Currently, RRAM is crucial to the fields of memory devices and storage technologies.A memristor is a two-terminal electric element that can change its conductance or resistance based on the state of the active material [1,2].Various materials have contributed their best performances to the development of memristive devices, including oxides [3][4][5], polymers [6], ferrites [7], and biomaterials [8], owing to their functional electronic, magnetic, and physical properties.Specifically, in terms of switching uniformity, cycle-to-cycle or device-to-device reliability, and low-power consumption, metal-oxide-based resistive memory devices perform comparatively well.
For several years, oxide-based memristors have been considered a promising candidate for neuromorphic computing [3,4], with the most probable reality being their continued impact with attractive advantages and rising opportunities in the field of memory devices.Over the past decade, studies on Schottky barrier height (SBH) modulation have exhibited a sharp escalation for various material systems, such as InP [9], BiFeO 3 [10], HfO 2 [11], and TaO x [12], including SrTiO 3 (STO) used in our model system as an example [13][14][15][16].In the case of STO that is n-type, a Schottky contact is formed at the contact with Pt, a metal with high work function.The SBH depends upon the concentration Sagar Khot and Dongmyung Jung have contributed equally to this study.
of oxygen vacancies ( V O ⋅⋅ ), that is positive space charges, in the STO.The local concentration of vacancies near the Schottky contact that determines the SBH can be increased or decreased via the application of bias voltages that disturb the spatial distribution of vacancies.Finally, the SBH change results in the conductance change in the memristor.
The resistive switching mechanism can be divided into two different types: filamentary and interfacial.The current flow can be explained differently based on the type of switching.In the case of filamentary switching, the resistance state is determined by the formation and rupture of a conductive filament (CF).The current in the low-resistance state (LRS) flows through the localized conduction path that is the CF in the insulating matrix, whereas that in the high-resistance state (HRS) where the CF is ruptured flows through the films homogenously [17].Inevitably, the current change between LRS and HRS is abrupt.Conversely, the current in interfacial switching devices is determined via the SBH at the interface between the electrode and oxide.Compared to the filament change, the SBH change is more gradual.Therefore, interfacial switching can be used for applications in analog or neuromorphic computing because it exhibits better performance in terms of controllability and variability compared to filamentary switching [18].In the case of storage applications, interface switching is plagued by certain limitations owing to its poor retention characteristics.However, filamentary switching can be considered a better candidate for storage and memory applications owing to its superior retention characteristics.
This study used a commercial finite-element solver, COMSOL Multiphysics, where the migration of charged species by an applied electric field can be calculated.A twostep simulation was performed for the interfacial switching by the SBH modulation described earlier.The first step is to obtain the spatial distribution of vacancies disturbed by an applied bias voltage, which corresponds to the "switching simulation."The second step is to measure the conductance value resulting from the modulated SBH, a function of the vacancy concentration near the contact.The conductance value is obtained by calculating the current value for a small positive voltage, which corresponds to the "read" simulation.We compared our modeling results with those through experiments reported in ref. [13] for better clarification of the suitability of our model.We believe that the model demonstrated in this study can be utilized in investigating the engineering directions of device architectures and materials for synaptic devices.

Model system and workflow
The model memristor system considered here as an example was Pt/STO/Nb-STO in Ref. [13], where Nb-STO denotes Nb-doped STO, as illustrated in Fig. 1a.The Pt/STO contact is Schottky whereas the STO/Nb-STO contact is Ohmic.All three layers have a thickness of 5 nm.Most importantly, a 0.5-nm-thick interfacial layer is placed between Pt and STO to depict the SBH modulation.The interfacial layer is designated as an STO material, and the average vacancy concentration value within this layer is used to determine the SBH and contact resistivity (Ωm 2 ) or conductivity (S/ m 2 ).The thickness of 0.5 nm is selected, considering both accuracy and numerical stability.The contact area was set to 25 by 25 nm.In the actual 2D simulation, both device width and area factor were 25 nm.In addition, a scaling factor was used to fit the experimental results where much larger electrodes were used, which will be explained in Sect.2.3.This system can be fully described one-dimensionally (1D).Nevertheless, the two-dimensional (2D) model was developed by considering its applicability to various architectures in the future.Moreover, the developed 2D model can be easily converted to a 3D model.
The SBH, Φ b , is determined from the average vacancy concentration, C , in the interfacial layer using Eq. ( 1).
(1) 1) is an empirical equation that we designated and is plotted in Fig. 1b.Φ bM and Φ bm are the maximum and minimum values of SBH, respectively, and C M and C m are the maximum and minimum vacancy concentrations (mol/ m 3 ), respectively.Φ bM , Φ bm , C M , and C m are model param- eters listed in Table 1.The reason why the adoption of the logarithm term is simply because it resembles the relation between doping concentration and Fermi level in the device physics textbook.This equation needs to be verified more.However, it yields the best match of switching curves to experimental ones at this moment.
The workflow, splitting it into four different parts involving definitions, two-step simulations, and outputs, is summarized in Fig. 2. In the first-step simulation, the drift and diffusion equations for oxygen vacancies were solved for the 5.5 nm STO including the interfacial layer, which will be explained in Sect.2.2.In the second-step simulation, the SBH value was determined from the average oxygen vacancy concentration in the interfacial layer, C (mol/m 3 ), using Eq. ( 1).Subsequently, it was used to estimate the conductance of the device, which is explained in Sect.2.3.
A Schottky contact is associated with the Pt/STO interface, whereas an ohmic contact exists at the STO/Nb-STO interface in the Pt/STO/Nb-STO memristor.The resistance of the device comprises the insulator and the Schottky contact, which are in series, the bulk and interface parts, respectively ( R total = R bulk + R interface ).The bulk resistance is a function of vacancy concentration in the STO.The interfacial resistance is  a function of SBH that is also a function of vacancy concentration near the interface.We assume that only R interface changes while R bulk is constant.The migration distance of vacancies is the order of picometer, which will be shown in Sect.2.3, and the drift-diffusion process of vacancies does not change R bulk .
The interfacial conductance, G interface = 1∕R interface is primarily owing to the thermionic emission current across the Schottky barrier.The range of the SBH value ( Φ b ) is set to 0.6 ~ 0.7 eV, as listed in Table 1, referring to an experiment reported by Gwon et al. [15].This explains the variation in the SBH based on the doping level and other papers also show the doping dependence [19,20].Robertson et al. reported a calculated barrier height for SrTiO 3 on Pt of 0.89 eV [21].
As illustrated in the workflow in Fig. 2, defining the device architecture, that is the geometry and constituent materials, and assigning material and interface properties is the starting point.Numerical calculation for a single operation comprises two steps "switching" and "reading."The first step is to obtain the spatial distribution of the oxygen vacancy concentration after applying a single-voltage pulse by performing a transient analysis of the drift and diffusion equations.The second step is to obtain the conductance caused by the vacancy-concentration-dependent SBH by calculating the current for a small applied voltage.By repeating a pair of the first-and secondstep calculations, the conductance change can be obtained by the successive application of switching pulses.In other words, the potentiation and depression curves for synapse devices can be simulated.The switching pulse can be both positive and negative while the reading pulse is always positive.Two types of switching pulses were applied, namely constant and incremental-step pulses (CP and ISP).Furthermore, to check the uniformity and stability of our model, we repeated the same cycle four times for the CP; in the case of the ISP, we repeated the same cycle five times.

Switching simulation
Transient analysis was performed for the drift and diffusion of vacancies according to the applied voltage pulse, V a (t) .Initially, the vacancy concentration (mol/m 3 ) is uniform as C ⃗ r, 0 = C i , where C i denotes the initial vacancy concentra- tion.At a time step, the electric potential distribution, V ⃗ r , where ⃗ r denotes position vector, in the device was obtained by solving the Poisson equation for given boundary conditions for the applied voltages of V a (t) and 0 at the top and bottom electrodes.
where denotes electrical conductivity.In our device, the potential gradient is constant because all -values of all materials are constant.The electromigration or drift of the charged oxygen vacancies, V O ⋅⋅ , occurs under an electric field and consequently concentration gradient is developed.Diffusion of the vacancy then occurs owing to the concentration gradient.The flux of vacancies, J , can be described using the classic Nernst-Planck transport equation in our model, as follows: where D , C , z , and denote diffusivity, concentration, charge number, and mobility of vacancies; F is the Fara- day constant (96,485.3sA/mol), and z denotes the effective valence for oxygen vacancy as + 2 based on a previous study [22].In our model, no generation and annihilation of vacancies are considered, meaning that the total amount of vacancies is conserved.Therefore, the continuity equation for vacancy also holds.These Eqs. ( 2)-(4) were solved using AC/DC and Chemical Reaction Engineering modules provided in COMSOL Multiphysics.

Read simulation
In this step, the conductance of the memristor is calculated for the spatial distribution of the vacancy concentration obtained from the switching simulation.The SBH is defined as a function of the average vacancy concentration in the interfacial layer, Φ b = Φ b C , using Eq. ( 1), as mentioned in Sect.2.1.The thermionic emission (TE) current density, J TE , for an applied bias voltage, V , is then expressed as: where A * denotes the Richardson constant (1.56 × 10 6 A/m 2 / K 2 ), T is the absolute temperature (300 K), q is the elementary charge (1.6 × 10 -19 C), and k is the Boltzmann constant (8.617 × 10 -5 eV/K).The specific contact conductivity, S C (S/m 2 ), is then expressed as: Finally, the conductance of the simulated device in Fig. 1a, G sim , is the product of the contact conductivity for (3) the read voltage of 0.1 V, and the effective electrode area, A electrode .
In our simulation, the S C value is calculated using the C value obtained in the first step and is assigned to the Pt/ STO interface.The G sim value is automatically calculated reflecting the device geometry.Finally, a scaling factor A f is multiplied to compensate for the difference in the electrode area of the experimental and simulated devices.
The A f value is calibrated to 5000 to match the experi- mental data.In Ref. [13], the electrode area is not presented, while our simulated device has 625 nm 2 (25 nm by 25 nm).The A f value corresponds to an electrode area of approxi- mately 3.125 μm 2 for the experimental device.
All model parameters used in our simulation are listed in Table 1.The diffusion coefficient of oxygen vacancies in the STO was chosen to be 6 × 10 -17 m 2 /s, which yielded conductance change curves close to the experimental ones.This value is somewhat higher than reported values of 4 × 10 -18 m 2 /s, 2 × 10 -21 m 2 /s, and 3.2 × 10 -22 m 2 /s at 300 K [23][24][25].The electrical conductivity of STO, , is assumed to be uniform.When a voltage and pulsing time of 1 V and 10 μs are applied to a 5-nm-thick device at 300 K, the migration distance, d , of a vacancy by a switching pulse is to be only a few picometers, using the equation, d ∼ Et = zqD kT Et .The numerical calculation is , where e denotes the elementary charge.The mobility of vacancies, , was calculated using the Einstein relation, = zqD∕kT .This short migration distance implies that the vacancy concentration profile changes only near the interface, which will be shown in Sect.3.3, and the bulk resistance, R bulk , does not change by the switching operations.Therefore, the assumption of the uniform constant conductivity of STO is valid.The conductance of the device is controlled only by the interface portion, R interface , which changes with the vacancy concentration.

Constant voltage pulses
Conductance measurement tests were performed to verify the synapse operations for the reported device.The device performance in the form of potentiation and depression after repetitive constant potentiating and depressing (positive and negative bias) pulses is shown in Fig. 3.The pulse-potential profile is depicted in Fig. 3a.The pulse function used here was identical to that used by experimentalists for their operation.Four cycles were repeated to analyze the behavior of the conductance for consecutive application of 30 pulses for the positive and another 30 pulses for the negative half to perform one cycle.All the pulse widths and intervals were maintained at 10 µs, with rising and falling time of 0.1 µs Fig. 3 Conductance variation on the application of potentiation and depression pulses for constant voltage amplitude.a Operated pulse with the constant potential for modulation of conductance.The amplitude of the pulses for potentiation and depression are V p = +1.62Vand V d = −1.80V, respectively.b Combined simula-tion (upper) and experimental graph (lower).The decreasing peak conductance values are observed in the simulation.This is the initial transient effect, and the peak conductance converges as the cycling is repeated and pulse amplitude of V p = +1.62V and V d = −1.80Vfor potentiating and depressing pulses, respectively.
When a positive potential was applied to the top (Pt) electrode, oxygen vacancies migrated from the top (Pt) electrode to the bottom (Nb-STO) electrode, which lowered the Schottky barrier height ( Φ b ) and resulted in a low resist- ance state or an increase in the conductance.Reciprocally, compensation of these positive charges at the interface layer was achieved by applying a negative potential to the top (Pt) electrode resulting in vacancies migrating back to the top electrode side.Consequently, the Schottky barrier height ( Φ b ) increased, the resistance evolved to the HRS, and the conductance decreased reasonably.Thus, a quasi-continuous change in conductance was observed.Conductance modulation can be achieved for positive and negative pulses, indicating a stable and uniform behavior.The conductance change as a function of the number of pulses at a fixed pulse height for four consecutive cycles confirmed the consistent switching of potentiation and depression.When considering the depressing pulses after potentiation revealed an abrupt decay in conductance, and a further gradual reduction in the conductance resulted in a closing value by subsequently repeating the depression pulses.It is evident from the graph that depressing pulses accompany a faster and more abrupt decay in conductance as compared to the increase through potentiating pulses.To signify the comparison and confirm the trend, behavior, and candidness of our results, a clear presentation of both the simulation and experimental results in the collusion is presented in Fig. 3b.The conductance began to increase from 1 × 10 -9 S with the highest conductance values attained for the simulation and experiment being similar and closer to 2.7 × 10 -9 S. The simulation results were consistent with the experimental data, thus proving that the model developed for the interfacial switching mechanism functions smoothly.

Incremental step pulse (ISP)
Similar to the potentiation and depression characteristics explained in Fig. 3 with constant pulse amplitude, we also checked with linearly increasing potentiating and depressing (positive and negative) bias shown in Fig. 4a pulses starting from + 1.5 up to − 2.45 V for positive and + 1.0 to − 1.95 V for negative pulse operation at intervals of ± 0.05 V.In total, 20 consecutive potentiating and depressing pulses were performed for the operation of one cycle, and the observed behavior along with a comparison with experimental data is shown in Fig. 4b.Five consecutive cycles were performed, resulting in a uniform and stable switching.A well-matching behavior, closely matching conductance values, and a comparable trend confirmed the correctness of our model.However, a considerable difference was obtained for the summit value of the conductance in case of the simulation compared to the experimental value.In contrast to the constant pulse operation depicted in Fig. 3, an increasing pulse potential results in a slightly higher value, which is approximately 3.5 × 10 -9 S. Initial conductance depicted a value of approximately 1 × 10 -9 S, which is similar to the one observed in the constant pulse.
Fig. 4 Conductance variation upon consecutive application of potentiation and depressing pulses with successively increasing amplitude.a Operated nonlinear pulse with increasing potential for modulation of conductance.The amplitude increases from + 1.5 to + 2.45 V for positive and from − 1.0 to − 1.95 V at a step of ± 0.05 V. Combined simulation (upper) and experimental graph (lower).The increasing peak conductance values are observed in the simulation.This is the initial transient effect, and the peak conductance converges as the cycling is repeated When short pulses were applied, analog-mannered conductance modulation was observed.The determination of the overall memristor conductance is viable owing to the history of the applied voltages; specifically, the time interval between the applied voltage pulses.Chang et al. reported a similar behavior in their study [1].For instance, the application of a train of positive (negative) pulses strengthens (weakens) the conductance through the determination of whether the potentiation (depression) characteristics of the device is feasible.Further repetition of the depressing pulses propels the conductance to convey the initial value.The conductance of our device exhibited a gradual change after every subsequent pulse.The experimental work used for the comparison with the simulation results was based on an article by Yang et al., which explained the synaptic properties and interfacial switching mechanism of STO [13].
An obvious consideration can be made that depends on the number of stimulating pulses, and device performance can be adjusted from cycle to cycle.The observed stable and uniform potentiation and depression beyond the polarity of the applied bias are immensely profitable for memory enhancement in electronic synapse applications.

Concentration profile
The evolution of the vacancy concentration profiles according to the application of repetitive pulses is shown in Fig. 5.A pink dashed line indicates the boundary between the bulk and interfacial regions of STO.All concentration changes only occur at both top and bottom contacts because the migration distance is short as picometer-order, as explained in Sect.2.3.The initial concentration profile is set uniformly.As evident in Fig. 5a and c, the positive bias pushed the oxygen vacancies down, which ultimately lowered the interfacial concentration, C .In contrast, the negative bias drove back the oxygen vacancies to the original position, thereby resulting in an increase in the interfacial concentration, as shown in Fig. 5b and d.
STO thin films are n-type semiconductors with oxygen vacancies as the donors [13].The modification of the SBH at the Pt/STO interface conveys resistive switching because of oxygen vacancy electromigration [26].When a positive potential was provided, positively charged oxygen vacancies moved toward the bottom electrode.This lowered the concentration at the junction between the semiconductor and metal, where our Schottky barrier function by Eqs.
(1) and ( 6) was applied to regulate the higher conductivity.The concentration in the interface layer, C , accounts for the change in conductance or resistance throughout the device.In contrast, the negative potential at the top electrode brought back oxygen vacancies to the top electrode, reshaping the concentration of the interface layer, which ultimately resulted in lower conductivity.
Figure 5a and b shows the concentration profiles during the potentiating and depressing operations by CP, which were obtained initially and after the 15th and 30th pulses, that is, exactly at the mid and end of the cycling operation.Similarly, Fig. 5a and b show the concentration profiles by ISP which were obtained initially and after the 10th and 20th pulses, that is, exactly at the mid and end of the cycling operation.These concentration plots were only used to observe the shift of concentration in the device from top to bottom and bottom to top upon positive and negative pulses.

Validation for concentration profile
The evolution of the concentration in the interfacial layer is crucial for investigating and presenting explicit validation.A clear interpretation of the conductance changes can be perceived by focusing on the concentration profile associated with the interfacial layer owing to the increase and decrease in the SBH.As explained in Sect.2.3, we designated the SBH as a function of the concentration.The device structure represented in Fig. 1a, which highlights an interfacial layer, was considered for measuring the concentration, and the evolution of the concentration profile in this layer for both CP-and ISP-operated systems is shown in Fig. 6a and  b, respectively, along with their operating pulses depicted in Fig. 6c and d, respectively.
Four cycles were conducted for the CP and five cycles for the ISP for pulsing such that the concentration profiles observed for similar cycles are as shown in Fig. 6a and b, with continuous change over every next pulse.A positive bias lowered the concentration, which reasonably increased the conductance, leading to potentiation.However, gradually decreasing conductance for a higher concentration for negatively biased cycles disclosed depression behavior.In the case of the CP application, repeating over the cycles increased the minimum and maximum concentration points, which ultimately reduced the minimum and maximum values of conductance.In contrast, as evident in Fig. 6b, for the ISP application, these minimum and maximum values were found to decrease and eventually increase the conductance values.The observed concentration profile validated the concentration shown in Fig. 5.

Conclusions
This study provided a brief view of computational modeling for the modulation of SBH by using oxygen vacancies as a source for investigating the memory and synaptic characteristics of an oxide-based memristive device.Furthermore, our model system was compared with strontium titanate material as an example.An attractive approach for simulating interfacial memristive devices was successfully implemented using this model system.In particular, a two-step simulation method was executed, and the basic parameters and properties were estimated for the advancement of the model.Interfacial switching, diffusion profile, and modulation of Schottky barrier height are the key attributes of this study.Moreover, the synaptic properties were demonstrated reliably, and the theoretical model proposed for neuromorphic simulation was consistent with the experimental results for a similar switching mechanism.In addition to the pulse profiles for synapse operations, two distinct pulse profiles were employed (viz., constant and incremental step pulse).These pulses for the synaptic properties can be attributed to behavioral changes in the conductance graph.The concentration profile validated the change in the resistance/conductance values.

Fig. 1 a
Fig. 1 a Schematic model structure of the developed 25 × 25 nm 2 size, b Logarithmic function, given as Eq.(1), used in the electrical contact to modulate the Schottky barrier height

Fig. 2
Fig. 2 Workflow of the interfacial switching model for Schottky barrier height modulation splitting up into four different sections

Fig. 5
Fig. 5 Evolution of vacancy concentration profile along the z-axis for CPs with a positive (potentiating) and b negative (depressing) bias, respectively, and for ISPs with c positive (potentiating) and d

Fig. 6 a
Fig. 6 a & b Concentration evolution in the interfacial layer, C , for synapse characteristics.Inset discloses the concentration profile for single potentiation and depression cycle.c & d Operated pulse for the synapse characteristics

Table 1
Model parameters used in the development of the model STO 7.8 × 10 -4 S/m Electrical conductivity of SrTiO 3 z + 2 dimensionless Charge number of oxygen vacancy A * 1.56 × 10 6 A/m 2 /K 2 f 5000 dimensionless Electrode area scaling factor