Modeling Defect-Level Switching for Nonlinear and Hysteretic Electronic Devices

Previously, we demonstrated hysteretic and persistent changes of resistivity in two-terminal electronic devices based on charge trapping and detrapping at immobile metastable defects [H. Yin, A. Kumar, J.M. LeBeau, and R. Jaramillo, Phys. Rev. Applied 15 , 014014 (2021)]; we termed these devices as defect-level switching (DLS) devices. DLS devices feature all-electronic resistive switching and thus are volatile because of the “voltage - time” dilemma . However, the dynamics of volatile resistive switches may be valuable for emerging applications such as selectors in crosspoint memory, and neuromorphic computing concepts. To design memory and computing circuits using these volatile resistive switches, accurate modeling is essential. In this work we develop an accurate and analytical model to describe the switching physics in DLS devices, based on the established theories of point defect metastability in Cu(In,Ga)Se 2 (CIGS) and II-VI semiconductors. The analytical nature of our model allows for time-efficient simulations of dynamic behavior of DLS devices. We model the time durations of SET and RESET programming pulses, which can be exponentially shortened with respect to the pulse amplitude. We also demonstrate the concept of inverse design: given desired resistance states, the width and amplitude of the programming signal can be chosen accordingly.


Introduction
The trapping and detrapping of electronic charge at immobile point defects can change the conductivity of metal/insulator/metal (MIM) stacks.][4][5] This class of devices are known as all-electronic resistive switches.We recently proposed and demonstrated a mechanism for all-electronic resistive switching based on capture and release of electronic charge at metastable point defects at MoO3/CdS heterojunctions. 6These metastable defects can exist in two distinct lattice configurations: a shallow-and a deep-double-donor state, and they will switch between these two configurations upon stimulation with light or voltage bias; we termed this effect as defect level switching (DLS).In the case of MoO3/CdS heterojunctions, the DLS-active defects are the sulfur vacancies V S in CdS. 72][13][14] In DLS devices, the state retention time is determined by point defect metastability; compared this to MIM stacks, in which the designed potential wells determine the retention time.
Because of the "voltage-time" dilemma analyzed by Schroeder et al., all-electronic resistive switching is not suitable for nonvolatile memory: the high current densities required for short READ/WRITE pulse (~ns) will result in the retention time many orders of magnitude less than the benchmark value (~years). 15However, for emerging applications such as selectors in a crosspoint array, 16 and neuromorphic computing concepts that operate in the temporal domain, such as reservoir computing, 17 the dynamical behavior of volatile resistive random access memory (RRAM) is valuable.For instance, Wang et al. utilized the dynamics of silver nanoparticles to functionally resemble the synaptic Ca 2+ behavior. 18Using the volatile switching effect of silver filaments, Zhang et al. achieved the emulation of multiple essential features of artificial neuron in a unified manner without auxiliary circuits. 191][22][23][24][25] All-electronic volatile RRAM may prove particularly useful if the switching dynamics can be predictively modelled based on established semiconductor physics concepts, and if devices not relying on mass transport prove more reliable in operation than devices relying on ion motion, nanoscale redox chemistry, and other complex coupled electronic-ionic phenomena.In this context, DLS devices occupy a unique niche: they feature dynamics that can be quantitatively modeled and tuned through materials selection and device design, and they operate on the basis of charge state transitions at immobile defects.
Here we develop an analytical and highly efficient model to describe the dynamics of DLS devices, starting from the statistics of point defect charge state transition.In §1, we describe the statistics of DLS-active point defects, following established theories of metastable point defects in CIGS solar cells, [11][12][13][14] and II-VI semiconductors. 7,26,27In §2, we model a metal-semiconductor heterojunction with DLS-active defects, as we have demonstrated experimentally. 6The resistive switching is caused by DLS defects that remain in metastable charge states after voltage stimulus, resulting in an out-of-equilibrium band diagram.Using Kimerling's model of deep traps, we derive an analytical expression of the interface band diagram by approximating the space charge profile using a step function. 28The effects of bias treatment on the junction conductivity are consistent with our experimental results, and the heterojunction band diagram and metastable defect distributions computed from our model are in great agreement with those computed from numerical iteration (as proposed by Decock et al.). 14Also, we estimate that our model is at least one million times more efficient compared to the numerical approach.
The analytical nature of our model allows for simulations of the device dynamics in a timeefficient manner.In §3, we illustrate the nonlinear and hysteretic resistive switching of a DLS heterojunction device upon a sinusoidal driving voltage.We model and discuss the hysteretic behavior at different frequencies of the driving voltage.We model the time durations of SET and RESET programming pulses, which can be exponentially shortened with respect to the pulse amplitude.These simulations demonstrate the concept of inverse design: given desired resistance states, the time width and voltage amplitude of the programming signal can be chosen accordingly.The concept of inverse design may extend to materials selection led by circuit design and system optimization, provided that a larger and accurate database of DLS-active defects becomes available.DLS-active defects are found in a number of compound semiconductors, including DX centers in n-AlGaAs and III-nitrides, [8][9][10] the oxygen vacancy V O in n-ZnO, 26,27 the sulfur vacancy V S in n-CdS. 7Though in this section we use specific features of anion vacancies in wide-bandgap II-VI semiconductors to discuss configuration transition kinetics, the model can be generalized to any other semiconductor system containing DLS-active defects.

Configuration transition kinetics of DLS-active defects
In Fig. 1 we present the configuration coordinate diagram for DLS-active defects, using for illustrative example the case of anion vacancies that are double donors.We represent the deep, neutral state of the anion vacancy as V An × , which is occupied by two electrons, corresponding to V O × in ZnO and V S × in CdS; upon ionization, V An × can transform into the doubly-ionized, shallow donor state (V An ⦁⦁ ) * ; the asterisk indicates a metastable lattice configuration.In ZnO, due to the negative-U effects the singly-ionized state V O ⦁ is never thermodynamically stable. 26,29In §A1 we plot the equilibrium distributions of different charge states for the oxygen vacancy, and the concentration of singly-ionized state is at least ten orders of magnitude smaller than V O × or (V O ⦁⦁ ) * .Therefore, the transitions V An × /V An ⦁ and V An ⦁ /(V An ⦁⦁ ) * are experimentally inaccessible, and we only address the configuration conversion between the neutral and doubly-ionized metastable states V An × and (V An ⦁⦁ ) * .
In n-type semiconductors, the transition from deep to shallow donor configuration can happen via electron emission process, and the reaction path is given by: The transition in Eq. ( 1) requires simultaneous emission of two electrons as well as thermal activation over an energy barrier Δ 2EE , as indicated in the configuration coordinate diagram in Fig. 1.Accompanying the deep-to-shallow configuration transition, the cation-cation distances around the anion vacancy may increase by upwards of 30%. 30An energy of Δ 2EE is required to distort the lattice configuration, and shift the defect level to be resonant with the conduction band.The optical transition energy is higher than the thermal energy barrier, because crystal-momentumconserving optical transitions can occur only at fixed configuration coordinate. 8In addition to the reaction path described in Eq. ( 1), when valence band holes are injected through bias treatment in a junction, the deep-to-shallow transition is possible via simultaneous capture of one hole and emission of one electron, or simultaneous capture of two holes, with thermal energy barriers Δ 1HC , and Δ 2HC , respectively, and the corresponding reaction pathways are given by: As shown in Fig. 1, the reverse transition from the metastable state to the neutral state happens through simultaneous capture of two electrons from the conduction band as well as thermal activation over an energy barrier Δ 2EC , and the reaction path is given by: In principle, the transition from the doubly-ionized state to the neutral state can happen through hole emission processes, i.e., simultaneous capture of one electron and emission of one hole, or simultaneous emission of two holes.To emit holes to the valence band, the defect needs to be in acceptor states.Because the defects discussed in this work (anion vacancies in II-VI semiconductors) are either in shallow donor states or deep donor states, hole emission processes are unlikely to happen.However, in some other semiconductor systems, the contribution from hole emission needs to accounted for.For instance, the defect complexes in CIGS are amphoteric and can transform between donor and acceptor states. 13In this case, the donor-to-acceptor configuration transition pathways should include those through hole emission processes.
The energy barriers related to these configuration transitions can be calculated by ab initio methods or inferred from experimental data.Using the local density approximation (LDA) of density function theory (DFT), Lany and Zunger calculated the configuration coordinate diagram for the oxygen vacancy in ZnO, which gives Δ 2EE = 3.7 eV, Δ 2EC = 0.2 eV and Δ 1HC = 1.3 eV. 27For V S in CdS, using the experimental data on photoconductivity, we infer that the thermal activation energy for electron recombination is Δ 2EC = 0.6 eV; using the experimental data on resistivity of MoO3/CdS heterojunction devices, we infer Δ 2EE = 0.9 eV . 6,7In the presence of free holes, the deep donor DX state is immediately ionized without thermal activation, i.e., Δ 2HC = 0 for both defects. 31An important future work for us is to use DFT calculations to accurately compute the configuration coordinate diagram for V S as well as other DLS-active defects, and compare the calculations results to the experimentally determined activation energy barriers.Now we discuss the kinetics of transitions described in Eqs. ( 1)- (4).According to the Shockley-Read-Hall (SRH) mechanism, we denote one-electron emission/capture rates as  ee/ec − . 13Hence, the transition rate  2EE −1 is given by the following expression: Here, the assumption is that the one-electron capture cross sections are the same in the first and second ionization processes.We acknowledge that it's a crude assumption, and it may lower the accuracy of simulation results.According to the analysis by Chicot et al., the cross section is larger in the second ionization than the first one, based on an argument of the Huang and Rhys parameter. 32However, more rigorous analysis and calculations are required to accurately determine the capture cross sections in the two ionization processes.As we will show in §3, Eq. ( 25) determines the device dynamics and the capture cross sections only enter the prefactor of Eq. ( 25), without affecting the hysteretic behavior and voltage dependence and of the device dynamics.Therefore, in this work we are not bothered unduly by this assumption.
Following our discussion on  2EE −1 ,  2EC −1 ,  1HC −1 and  2HC −1 are given by: where  V is the effective density of states of the valence band.
We should note that in this work the charge transition rates are basically calculated by the SRH model, and we account for the lattice configuration changes by introducing Boltzmann factors.The SRH theory is a simplified model and it was shown to be invalid to describe the charge exchange rates between two materials. 334][35][36] Integration of these advanced defect models into our model will effectively improve the accuracy of simulation results, that we leave in our future work.
We denote the shallow-to-deep and deep-to-shallow configuration transition rates as  S→D −1 and  D→S −1 , respectively.For the shallow-to-deep transition, Eq. ( 4) is the only pathway considered in this work.For the deep-to-shallow transition, the two-hole capture process in Eq. ( 3) is the dominant mechanism because it has the lowest thermal activation energy.The other processes have higher energy barriers and can gain importance at high temperatures.In §A1 we analyze  2HC −1 ,  1HC −1 and  2EE −1 as a function of Fermi energy and demonstrate that  2HC −1 is the most important contribution upon junction reverse biasing and hole injection.Therefore, −1 , which also preserves the symmetry of rate equations: We denote the Fermi level position where the deep donor state V An × and the metastable shallow donor state (V An ⦁⦁ ) * are equally likely as the transition energy  trans , i.e.,  D→S −1 =  S→D −1 when  F =  trans .Using Eq. ( 9) and Eq. ( 10), we have: where  g is the semiconductor band gap.For  F >  trans ,  D→S −1 <  S→D −1 and the deep donor state V An × is thermodynamically favored; for  F <  trans ,  D→S −1 >  S→D −1 and the metastable shallow donor state (V An ⦁⦁ ) * is thermodynamically favored.
If we denote the concentration of DLS-active defects as  DLS , and the concentration of deep and shallow donor states as  deep and  shallow , respectively, the time variation of  shallow () and  deep () are given by: At steady states, the time variation of  deep () and  shallow () are both equal to zero, yielding steady state concentrations  deep SS and  shallow SS : DLS (15)   The equations in this section connect the dynamics of DLS devices to defect and material parameters, via fundamental semiconductor statistics.In Table 1 32 The order of magnitude of 10 −13 cm −2 is in agreement with a positively charged defect center and we assign  = 10 −13 cm −2 for both V S in CdS and V O in ZnO.We take the phonon frequency  ph = 1 THz for both defects. DLS and  d are the concentrations of DLS-active defects and non-DLS-active shallow donors.We keep  DLS the same as the metastable defect concentration that we use in device simulations in our previous work. 6The layer thickness is used in the simulations below.

Resistive switching physics in DLS-heterojunction devices
Configuration conversion between deep neutral and ionized states of DLS-active defects leads to persistent photoconductivity, and is useful in high-responsivity photoconductive sensors.However, it can also be achieved without illumination by applied bias in a metal-semiconductor Schottky junction, where the metal serves as a hole injection layer to induce the configuration transition defined in Eq. (3).Though operating at an interface, DLS devices are not based on the valence-change resistive switching mechanism in which crystal defects drift under an electric field, and there is no mass transport in DLS devices. 37In case of the MoO3/CdS heterojunction that we have experimentally demonstrated, MoO3 has a higher work function than CdS and thus can inject holes into CdS, leading to a space charge layer (SCL) within CdS. 6 Within the SCL, most metastable defects are in the shallow donor state (V S ⦁⦁ ) * .In Fig. 2a we present the energy levels of materials used to make a DLS-heterojunction device: an n-type semiconductor with DLS-active defects, and a metal or highly-doped semiconductor. m and  s refer to the metal work function and the semiconductor electron affinity, respectively.In the semiconductor,  deep and  shallow refer to the energy levels of deep and shallow donor states, respectively.Because of the negative-U effects, the first ionization energy  deep (corresponding to the V An × /V An ⦁ transition) is deeper than  trans (corresponding to the V An × / ( V An ⦁⦁ ) * transition). 38The Fermi level  F is above the transition energy  trans .We denote by  n the energy difference between  C and the Fermi level  F .In Fig. 2b we present the schematic band diagram of an unbiased DLS-heterojunction.The formation of Schottky barrier  S is due to the metal work function  m being larger than the semiconductor electron affinity  s .In the semiconductor side, the conduction band edge energy  C () and the metastable state transition energy  trans () vary with distance  from the heterojunction.We denote by  trans the location where  trans () intersects the Fermi energy: In Fig. 2b    trans defines where V An × and (V An ⦁⦁ ) * are equally likely at equilibrium.Under time-varying external stimulus, the distribution of metastable states may be driven out-of-equilibrium, and  trans may no longer describe the position where  deep () and  shallow () are equal to  DLS /2.We define  DLS as the instantaneous position where the deep and shallow donor state concentrations are the same: At  <  DLS ,  shallow () >  deep () ; while at  >  DLS ,  shallow () <  deep () .At thermodynamic equilibrium,  DLS =  trans .
Under positive bias applied to the metal side of the heterojunction in Fig. 2b, the concentration of electrons in the space charge region of the semiconductor will increase.Consequently, around  DLS where (V An ⦁⦁ ) * and V An × coexist, electron recombination will convert (V An ⦁⦁ ) * to V An × as described in Eq. ( 4).As a result, fewer defects are in the ionized state (V An ⦁⦁ ) * and more defects are in the deep neutral state V An × , and  DLS will move closer to  = 0 according to Eq. ( 17).This process requires thermal activation, and it's thus slow to follow the applied bias, leading to hysteretic redistribution of  deep (, ) and  shallow (, ).Hence,  deep (, ),  shallow (, ) and  DLS () are dynamical, history-dependent variables.
To track the evolution of  deep (, ) and  shallow (, ), the Poisson equation must be solved at each time to compute  C (, ) and  V (, ), which are used to calculate the configuration transition rates defined in Eq. ( 12) and Eq. ( 13), and then update  deep (, ) and  shallow (, ) at next time step.Solving the Poisson equation at each time is computationally demanding because it requires numerical iteration.For faster computational simulation of the device dynamics, we simplify the profile of space charge density to a step function, as was proposed by Kimerling for p + n junctions that contain deep traps. 28Kimerling assumed that deep traps contribute to the space charge only within a region of width  from the interface, which bears the same meaning as the  DLS in Eq. (17).Hence, the space charge density is given by the following expression: For the cases simulated here, DLS-active defects are chosen to be the dominant type of defects in the materials discussed, i.e.,  DLS is much larger than the shallow (non-DLS) donor concentration  d (Table 1), and the contribution from regular shallow donors to the space charge profile is negligible.For cases in which both metastable and shallow defects make meaningful contributions to the space charge profile, as in CIGS thin films, Decock et al. proposed to use a staircase space charge density profile (this multi-step model can be seen as an extension of the one-step model). 14n §A2 we explicitly compare the results of this staircase approximate model to direct numerical simulation and find good agreement.Therefore, the assumption of step-function changes in space charge density profile is generally applicable to junctions with deep and shallow levels.
We use Gauss's law to calculate the electric potential profile variation through the junction.The boundary conditions of the conduction band energy  C (, ) (referenced to the Fermi energy) on the two sides of the semiconductor are: s and  s are the dielectric constant and width of the semiconductor, respectively. bias is the applied bias at the given time . C (, ) is then given by: where ( bias ) is given by: Taking the parameters of the MoO3/CdS heterojunction (that we list in Table 1, and for MoO3 we take its work function to be 6.7 eV), we illustrate the effect of bias on the interface band diagram and metastable defect states in Fig. 3.As shown in Fig. 3a and Fig. 3b are the transition energy level  trans () and the concentration of shallow donor states  shallow () in the CdS layer. trans () is computed by combining Eq. ( 11) and Eq. ( 20) while  shallow () is computed from Eq. ( 15).In Fig. 3: the dashed blue line indicates  trans () or  shallow () at zero bias, while the orange and green solid lines indicate those after forward bias treatment (+0.8 V) and reverse bias treatment (-0.8 V), respectively; in Fig. 3b, the vertical dashed lines indicate positions of  DLS .After forward bias treatment,  DLS moves from 16.1 nm to 9.5 nm.According to Eq. ( 17), the decrease of  DLS means fewer defects are in the ionized shallow state, which reduces the effective donor doping concentration and the electric field at the junction; as a result, the interface barrier is thicker than before, switching the junction to high resistance state (HRS).After reverse bias treatment,  DLS moves from 16.1 nm to 20.7 nm, meaning that more defects are in the ionized state, increasing the effective donor doping concentration and the electric field at the junction; therefore, the interface barrier is thinner than before, switching the junction to low resistance state (LRS).These simulation results are consistent with the bipolar resistance switching that we demonstrated experimentally: DLS devices switch to HRS under forward bias, and switch to LRS under reverse bias. 6We note that the changes of defect charge states and depletion region upon bias treatment shown in Fig. 3 are in analogy to those in deep level transient spectroscopy (DLTS) measurements.In DLTS, a forward bias treatment can neutralize defects and widen the depletion region, causing a temporary decrease of the junction capacitance; in this work, the thickening of interface barrier induced by forward bias treatment causes a decrease of the junction conductivity.11) and Eq. ( 20): without bias treatment (blue dashed line), after forward bias (orange solid line) and reverse bias treatments (green solid line).The square symbols represent the results that we compute using the numerical iteration method developed by Decock et al. 14 Energy is referenced to the Fermi level (grey dot-dash line).(b) Distribution of shallow donor states  shallow () that is computed from Eq. ( 15): without bias treatment (blue dashed line), after forward bias (orange solid line) and reverse bias treatments (green solid line).The vertical dashed lines indicate positions of  DLS .The heterojunction is in HRS and LRS after forward and reverse bias treatments, respectively, and the corresponding positions of  DLS are labeled as  DLS,HRS and  DLS,LRS , respectively.
To examine the accuracy and efficiency of our model, we use the numerical iteration method developed by Decock et al. to compute  trans () and  shallow (). 14The results are shown as square symbols in Fig. 3 and they are in good agreement with those from our model which are represented by dashed and solid lines.Also, we estimate that our analytical model is at least one million times more efficient than the numerical iteration approach: to compute one band diagram in Fig. 3a, it takes less than 60 micro-seconds with Eq. ( 20), while it takes about 60 seconds to achieve the self-consistency with the numerical iteration method.All our simulations are executed with Python code and NumPy package.
Using the interface band diagram that we compute from Eq. ( 20), we estimate that after forward bias treatment the junction resistivity is 783 times under zero bias (conductivity modeling described in §A3).At the doping concentration  DLS = 2 × 10 18 cm −3 , as in Table 1, tunneling is the dominant transport mechanism, and it is exponentially dependent on both of the height and width of the interface barrier.In experiments, we found that the MoO3/CdS heterojunction resistance could increase by up to 200 times after a forward bias treatment of +0.8 V at room temperature. 6When the concentration of DLS-active defects is reduced  DLS to 1.5 × 10 18 cm −3 , we calculate the resistivity ratio of the forward-biased to the zero-biased junction to be 187, which is close to the experimental value.
As shown in Fig. 3b, in the HRS more metastable defects are in the neutral deep state (more electrons are trapped at the interface) while in the LRS more are in the ionized shallow donor configuration (fewer electrons are trapped at the interface).The physical insight is the same as presented by Simmons and Verderber: trapped charge distorts the interface band structure, and as a consequence changes the electric field and tunneling currents passing through the interface. 1.( 20) and Fig. 3 indicates that  DLS solely determines the interface band diagram as well as the resistivity.Therefore,  DLS is a state variable of the DLS device and the device dynamics are represented by the movement of  DLS ().Its velocity,  DLS (), can be determined by the space and time variation of  shallow (, ) at  DLS (): shallow ( DLS , )   shallow ( DLS , )  (22)   Combining Eqs. ( 9)-( 12) and ( 17), and estimating  C =  V , we have the numerator of Eq. ( 22) given by the following expression: Given the boundary conditions of  C (, ) in Eq. ( 19), we have  F =  bias , which indicates that the numerator of Eq. ( 22) can be increased exponentially under applied bias.
We rewrite the denominator of Eq. ( 22) as: shallow ( DLS , )  = −  DLS  C (24)    C represents the characteristic length of the truncation of the shallow donor distribution  shallow () around  DLS .For the  shallow () shown in Fig. 3b, we calculate  C to be 3.1 nm for the junction under zero bias, and 3.6 nm and 2.8 nm for that after forward and reverse bias treatment, respectively.Since  C only changes slightly after bias treatment, we take it to be constant and thus Eq. ( 24) is constant.Combining Eqs. ( 22)-( 24) yields: where  0 is a constant velocity, that is given by: In deriving Eq. ( 25), we have simplified the complex process of defect ionization and electron recapture in DLS devices into the motion of  DLS , while preserving the essential switching physics.The expression of the velocity of  DLS has a similar form to that describing the growth velocity of conductive filaments in metal oxide resistive switching devices: both have a hyperbolic-sinusoidal dependence on the applied bias, which means that in both cases the resistive switching is bipolar and the switching speed can be increased exponentially with the drive voltage amplitude. 39owever, the switching physics is different: in the case of DLS-heterojunction devices, the exponential dependence stems from electron or hole injection whose densities are exponentially dependent upon applied bias, while in the case of metal oxide resistive switching devices, the exponential dependence originates from electric field driven ion hopping.Also, in Eq. (25) all of the parameters are constructed from point defect physics, and there is no mass transport.

Dynamic behavior of DLS-heterojunction devices upon voltage stimulus
As an example, in Fig. 4a, we show the response of  DLS (red solid line) upon a sinusoidal driving voltage (blue solid line) at a frequency of 1 MHz and an amplitude of 1.8 V. We also show the steady-state (i.e., equilibrium) position of  DLS upon the driving voltage with the gray dashed line, and we see that the device operates in out-of-equilibrium states of  DLS .As illustrated in Fig. 4b,  DLS is hysteretic with respect to the driving voltage, and it changes in the interval of 13.8 -17.4 nm.When the device is at HRS,  DLS =  DLS,HRS = 13.8 nm, and when the device is at LRS,  DLS =  DLS,LRS = 17.4 nm.Under this sinusoidal driving voltage, we estimate the resistivity ratio of HRS to LRS to be 50.Since the hysteretic resistive changes are based on non-equilibrium charge states of metastable point defects, DLS devices are predicted to be volatile and not suitable for nonvolatile memories. 15t lower driving frequency, we expect reduced resistive switching because the metastable defects have more time to approach their thermodynamic equilibrium states.For sufficiently low frequency (smaller than the order of magnitude of the inverse of the device retention time), the hysteretic loop shown in Fig. 4 will collapses into a line, i.e., the hysteresis vanishes and  DLS closely follows the applied bias.The room-temperature retention time of the HRS at MoO3/CdS heterojunctions has been found to be about 90 s. 6 At the HRS,  DLS,HRS <  DLS,0 (as previously demonstrated in Fig. 3b), and the DLS-active defects will gradually recover their thermodynamic equilibrium states, by emptying the occupied electrons through electron emission.Using Eq. ( 5) and the experimental fact that  2EE = 90 s, we estimate Δ 2EE = 0.9 eV. 2EE −1 impacts the device dynamics at low driving frequencies.Adding Eq. ( 5) into Eq.( 10) and Eq. ( 12), we can study the hysteretic behavior of  DLS at different frequencies and compute the corresponding resistivity ratio of HRS to LRS.In Fig. 5 we show the resistivity ratio  HRS / LRS upon sinusoidal driving voltage of different frequencies; insets show hysteretic loops of  DLS .At frequency  = 10 −4 Hz, which is much lower than the inverse of the retention time, i.e., 1/90 s −1 , the hysteretic loop of  DLS collapses into a line as expected, and  HRS / LRS = 1.With increasing frequency,  HRS / LRS increases to a maximum and then decreases.The decrease is because, at very high frequency, defects have insufficient time to switch between configurations.As indicated by Eq. ( 25), the equation of motion of  DLS , the resistivity ratio is limited by the single electron/hole emission rates, the configuration transition barriers, the driving voltage amplitude, as well as other material properties.Using higher voltage amplitude will increase the injected carrier concentration and thus increase the resistivity ratio  HRS / LRS .Using Eq. ( 25), we can study the switching time between HRS and LRS of DLS-devices under voltage pulses, which are typically used to program resistive states in RRAM devices, and are relevant to neuromorphic computing strategies.To understand the switching process intuitively, we derive analytical expressions of the switching time.Eq. ( 25) is analytically solvable if we expand  trans ( DLS , ) linearly with respect to  DLS , yielding: trans ( DLS , ) = (  trans  DLS ) 0 ( DLS −  DLS,0 ) +  DLS  s  bias (27)   where  DLS,0 is the equilibrium position of  DLS at zero bias (as previously demonstrated in Fig. 3b).Eq. ( 27) approximates  trans ( DLS , ) fairly well in the range of  DLS and q bias that we are interested in ( §A3).
Under reverse voltage pulse  bias < 0 , the device changes from HRS to LRS (the SET operation):  DLS moves from  DLS,HRS to  DLS,LRS .Combining Eq. ( 25) and Eq. ( 27), we solve for the switching time  SET : Under forward voltage pulse  bias > 0, the device changes from LRS to HRS (the RESET operation):  DLS moves from  DLS,LRS to  DLS,HRS .Combining Eq. ( 25) and Eq. ( 27), we solve for the switching time  RESET : We see in Eq. ( 28) and Eq. ( 29) that the SET/RESET switching times can be exponentially reduced by increasing the pulse amplitude | bias |.So, using voltage pulses of larger amplitude is preferred to achieve faster switching.We note that  SET is only dependent on  DLS,LRS and independent of  DLS,HRS ; vice versa,  SET is solely dependent on  DLS,HRS , regardless of the position of  DLS,LRS .
Eq. and Eq. ( 29) enable predictive device design and selection of operating parameters.Provided that we are required to design a device with a resistance ratio of 50.From Fig. 4 we learn that using a combination of  DLS,HRS = 13.8 nm and  DLS,LRS = 17.4 nm can achieve a resistance ratio of 50.Then, to switch the  DLS between 13.8 nm and 17.4 nm, we need to compute the time widths ( SET and  RESET ) and amplitude (| bias |) of programming pulses.
In Fig. 6a we compute and plot  SET and  RESET as a function of | bias |, using Eq. ( 28) and Eq. ( 29); the pulse time widths and amplitude should be chosen accordingly.For instance, we use a pulse amplitude of | bias | = 1.8 V; then the pulse time widths need to be  SET = 50.0ns and  RESET = 47.2 ns, respectively.In Fig. 6b we demonstrate the dynamic evolution of  DLS under a train of such SET and RESET voltage pulses, and the  DLS starts from the equilibrium position  DLS,0 , and switches between  DLS,LRS and  DLS,HRS as designed.Here, the evolution of  DLS is computed with Eq. (25).We compute and plot the dynamic evolution of  DLS (red dashed line) using Eq.(25). DLS starts from  DLS,0 = 16.1 nm, which is the equilibrium position of  DLS at zero voltage, and it then switches between  DLS,HRS and  DLS,LRS as designed.

Discussion and conclusion
In this work, we use the statistics of metastable point defects and an approximation of stepchanges in charge density to develop a physics-based model of resistive switching dynamics in DLS-based heterojunction devices, as we previously demonstrated experimentally. 6The analytical nature of our model allows for highly efficient device simulations, compared to iterative methods, and we validate our model approximations by reference to the numerical iteration.We explore how devices respond to voltage pulses and variable-frequency drive.
DLS devices are based on charge trapping/detrapping at immobile point defects, and thus the relevant resistive switching mechanism is purely electronic and inherently volatile.Therefore, in DLS devices the programmed resistive state will relax into a thermodynamically stable state after removing the programming signal, but the relaxation may offer desirable dynamics for the emulation of biological neurons and synapses. 24A notable feature of DLS devices is that they don't suffer from stochastic ion migration events which occur in devices based on mass transport, such as conductive filament (CF) based RRAM. 40As a result, DLS devices are more amenable to physics-based design and are likely less affected by cycle-to-cycle variation and offer long-term stability.
Predictive modeling of these volatile resistive switching devices will help in designing computing circuits that leverage their dynamics.For future work, a computational database of DLS-active defects is required.The database should include defect parameters which are necessary to guide the design of DLS devices, and these parameters can be computed by DFT calculations and calibrated by experiments.Such a computational database of defects, combined with the physics-based model presented here, suggests a future in which circuits employing DLS devices could be optimized through an inverse design process: first by identifying the circuit dynamics that best suit the proposed algorithm, then resolving dynamics to the level of individual devices, and then achieving the desired device dynamics by materials selection.Volatile DLS devices with designer dynamics could be useful for circuits of coupled, hysteretic, and nonlinear oscillators, such as concepts in physical reservoir computing. 41Our model is also useful in other junction devices containing metastable defects, such as solar cells, in which metastable defects have direct impact on device performance. 42

A1. Equilibrium distributions of different charge states and configuration transition rates for the oxygen vacancy in ZnO
The singly-ionized state V An ⦁ is thermodynamically unstable, i.e., the energy of the singlyionized defect V An ⦁ is always higher than the neutral state V An × and the metastable state (V An ⦁⦁ ) * , as shown in Fig. A1a. 26,27According to the formation energies, we compute the equilibrium concentrations of different charge states of the oxygen vacancy as a function of Fermi energy, as shown in Fig. A1b.The dominant charge state is either V An × or (V An ⦁⦁ ) * , and the concentration of the singly-ionized state V An ⦁ is at least 10 orders of magnitude smaller than the dominant charge state.

A2. Calculation of apparent defect density profiles using staircase space charge density profile
We reexamine the calculation of apparent defect density profiles using staircase charge density profile, as proposed by Kimerling and Decock et al. 14,28 Decock et al. proposed to use this approximation to explain the measurement results of apparent doping density profiles in CIGS solar cells (Eqs.( 13), ( 16) and (17) in ref. 14 ), but they did not directly compare modeling results of this approximate model to their experimental data and to iterative numerical simulations.
The DLS-active defects in CIGS solar cells are V Se − V Cu complexes.We optimize the following four parameters: (1) the built-in barrier of the CdS/CIGS heterojunction, (2) the shallow doping density  A , (3) the concentration of DLS-active defects  DLS (the transition energy level is  trans −  V = 0.19 eV), and (4) the density of the additional acceptor defect  t (the energy level is  t −  V = 0.33 eV).The optimized parameters are listed in Table A1, together with all the other parameters that are used for the calculation of apparent defect density profiles.

Materials
ZnO/CdS/CIGS Temperature (K) 200 Thickness of the CIGS layer (um) 1.0 Hole effective mass of the CIGS 0.087 e Built-in barrier of the CdS/CIGS heterojunction (V) 0.5  A (cm −3 ) 0.4 × 10 15 DLS-active defects V Se − V Cu complexes  DLS (cm −3 ) 0.8 × 10 15  trans −  V (eV) 0.19  t (cm/s) 1 × 10 16  t −  V (eV) 0.33 Table A1: Overview of the material parameters of CIGS solar cells that we use for the calculation of apparent defect density profiles.
We present in Fig. A3 the calculated apparent doping density profiles for the CdS/CIGS heterojunction without bias treatment (blue solid line) and after a reverse bias treatment of -2 V (red solid line).The square and triangle symbols represent the numerical simulation and experimental data from Decock et al., respectively; the blue and red colors represent those with and without reverse bias treatment, respectively.The calculated apparent doping density profiles are in satisfactory agreement with the results from Decock et al.
For the heterojunction after reverse bias treatment, there is a sharp increase of the apparent doping density (red solid line) when the depletion width reduces below 0.63 um, creating a local maximum of the apparent defect density.The reason is that, in the heterojunction, the length of the distribution of metastable states, that we define as  DLS in the main text, is 0.63 um.When the depletion width is larger than 0.63 um,  DLS is fixed at 0.63 um and DLS-active defects don't respond to the capacitance-voltage testing.When the depletion width shrinks below 0.63 um,

Figure 1 :
Figure 1: Configuration coordinate diagram for DLS-active point defects that exhibit metastable behavior -here, anion vacancies that are double donors. 0 and  * represent stable configuration coordinates for the neutral deep donor state V An × and doubly-ionized shallow donor state ( V An •• ) * , respectively; the asterisk indicates a metastable lattice configuration.System enthalpy vs. configuration coordinate V An × and (V An •• ) * are represented by orange and green solid lines, respectively.Thermal activation energies Δ 2EE and Δ 2EC necessary for transitions in Eq. (1) and Eq.(4) are indicated by orange and green arrows, respectively.

Figure 2 :
Figure 2: Schematic of a DLS-heterojunction.(a) Illustration of energy levels of the materials used: an n-type semiconductor with DLS-active defects, and a metal or highly-doped semiconductor.Terms are described in the text.(b) Illustration of the energy band diagram after junction formation.The black arrow indicates the height of the Schottky barrier  S .The gray dash-dotted line indicates the location of the interface where  = 0.The black star represents  =  trans , i.e., the intersection of the transition energy  trans () and Fermi level  F .For  <  trans ,  trans () >  F and the shallow donor state is thermodynamically favored; while for  >  trans ,  trans () <  F and the deep donor state is thermodynamically favored.

Figure 3 :
Figure 3: Relationship between interface band diagram and metastable defect states.(a) Transition energy level  trans () that is computed by combining Eq. (11) and Eq.(20): without bias treatment (blue dashed line), after forward bias (orange solid line) and reverse bias treatments (green solid line).The square symbols represent the results that we compute using the numerical iteration method developed by Decock et al.14 Energy is referenced to the Fermi level (grey dot-dash line).(b) Distribution of shallow donor states  shallow () that is computed from Eq. (15): without bias treatment (blue dashed line), after forward bias (orange solid line) and reverse bias treatments (green solid line).The vertical dashed lines indicate positions of  DLS .The heterojunction is in HRS and LRS after forward and reverse bias treatments, respectively, and the corresponding positions of  DLS are labeled as  DLS,HRS and  DLS,LRS , respectively.

Figure 4 :
Figure 4: Resistive switching and hysteresis behavior of a DLS-heterojunction device under a sinusoidal driving voltage whose frequency is 1 MHz.In (a) we plot the voltage stimulus (blue solid line) and the corresponding change of  DLS (red solid line) versus time.We also plot the equilibrium position of  DLS (gray dashed line) versus time.The device operates in out-of-equilibrium states of  DLS .In (b) we plot  DLS with respect to the voltage stimulus.We label the processes of SET (from HRS to LRS) and RESET (from LRS to HRS) with black arrows.

Figure 5 :
Figure 5: The resistivity ratio of HRS to LRS upon sinusoidal driving voltage of frequencies from 0.1 mHz to 1 MHz.In the three insets we plot the hysteretic loop of  DLS (i.e., the change of  DLS with respect to the voltage stimulus) at three different frequencies: 0.1 mHz, 10 mHz and 1 MHz.

Figure 6 :
Figure 6: SET/RESET of a DLS-heterojunction device with a resistance ratio of 50, triggered by voltage pulses, and  DLS of HRS and LRS are set to be 13.8 nm and 17.4 nm, respectively, i.e.,  DLS,HRS = 13.8 nm and  DLS,LRS = 17.4 nm.(a) Time widths of voltage pulse required to trigger the SET (green solid line) and RESET (orange solid line) versus amplitudes of voltage pulse.(b) A train of SET and RESET voltage pulses (blue solid line) that switch the device between HRS and LRS.The voltage pulses have

Figure A1 :
Figure A1: (a) Formation energies and (b) equilibrium concentrations of charge states V O × , V O ⦁ and (V O ⦁⦁ ) * for the oxygen vacancy in ZnO as a function of Fermi energy  F at temperature T=300 K.  V is the valence band energy level.Formation energies were calculated by Lany and Zunger. 27Here, the total concentration of oxygen vacancies is 2 × 10 18 cm −3 .We use black solid line, black dashed line and grey dash-dotted line to represent the formation energy or equilibrium concentration of charge state V O × , V O ⦁ and (V O ⦁⦁ ) * , respectively.
Fig. A2a and Fig. A2b, respectively. 2HC −1 ,  1HC −1 and  2EE −1 can contribute to the deep-to-shallow transition rate  D→S −1 .As shown in Fig. A2, at a high concentration of injected holes,  2HC −1 is the dominant contribution.Also, the configuration transition rates need to be comparable to the frequency of driving voltage (i.e., 1 MHz = 10 6 s −1 , which is indicated by a black dashed line in Fig. A2) to have an impact on the device dynamics.At T=300 K,  1HC −1 is way lower than 10 6 s −1 .

Figure A2 :
Figure A2: Configuration transition rates  2HC −1 (green solid line),  1HC −1 (grey solid line),  2EE −1 (blue solid line) and  2EC −1 (orange solid line) as a function of Fermi energy  F at temperature (a) T=300 K and (b) T=1000 K, respectively. V is the valence band 1, which are given by  ee −1 =  th  C  and  ec −1 =  th , respectively, where  th is the carrier thermal velocity,  C the effective density of states of the conduction band,  the electron density in the conduction band and  the carrier capture cross section.Since  is given by  =  C exp (−( C −  F )/), where  C is the energy level of the conduction band and  F the Fermi level, we can rewrite the electron capture rate as  ec −1 =  th  C exp (−( C −  F )/).

Table 1 :
, we list the properties of DLSactive defects in II-VI semiconductors, i.e., V O × in ZnO and V S × in CdS, and other material parameters, which we will use in simulations in §2 and §3.For future work, a computational database of DLS-active defects is required.The database should include defect parameters which are necessary to guide the design of DLS devices, such as thermal activation energies of configuration transitions.These properties could differ greatly from one defect to another, and should computed carefully and accurately by DFT calculations and calibrated by experiments.Such a computational database, combined with the model that we develop in this work, suggests promising opportunities for technology computer-aided design (TCAD) of resistive switching devices based on DLS phenomena.Properties of two wide-band gap II-VI semiconductors that host DLS-active defects: V O × in ZnO and V S × in CdS. C −  trans is calculated according to Eq. (11).For the carrier capture cross section , Chicot et al. identified a deep donor of a large electron capture cross section of 1.6 × 10 −13 cm −2 to be related to V O in ZnO.