The Role of Vacancy Dynamics in Two‐Dimensional Memristive Devices

Two‐dimensional layered transition metal dichalcogenides (TMDCs) are promising memristive materials for neuromorphic computing systems. Despite extensive experimental work, the underlying switching mechanisms are still not understood, impeding progress in material and device functionality. This study reveals the dominant role of defect dynamics in the switching process of 2D TMDC materials. The switching process is governed by the formation and annihilation dynamics of a local vacancy depletion zone. It explains the distinct features of the device characteristics observed experimentally, including fundamentally different device behavior previously thought to originate from multiple mechanisms. Key influence factors are identified and discussed with a fully coupled and dynamic charge transport model for electrons, holes, and ionic point defects, including image‐charge‐induced Schottky barrier lowering (SBL). Thermal effects and local Joule heating are considered by coupling the transient heat transfer equation to the electronic properties. The model is validated with hysteresis and pulse measurements for various lateral 2D MoS2‐based devices, strongly corroborating the relevance of vacancy dynamics in TMDC devices and offering a new perspective on the switching mechanisms. The insights gained from this study can be used to extend the functional behavior of 2D TMDC memristive devices in future neuromorphic computing applications.


Introduction
Two-dimensional (2D) and quasi-2D layered materials have attracted considerable attention for their application as memristive devices (memristors). [1,2]he current-voltage (I-V) characteristics of memristors are hysteretic, and therefore, their resistance can be modulated by applying voltage pulses, permitting switching between different resistance states. [3][5] They have been intensively investigated for neuromorphic computing [6,7] and were used in nonvolatile memory architectures, frequently referred to as resistive randomaccess memory (rRAM) or memristive random-access memory (mRAM). [8][14][15] They show excellent scaling potential, can easily be implemented in crossbar structures, and offer memristive switching at low voltages.For example, vertical electroforming-free MoS 2 memristors with switching voltage swings ΔU (ΔU = U set − U reset , where U set is the set voltage and U reset is the reset voltage) in the range 160 -550 mV have been reported, [12][13][14] while the lowest reported voltage swings for traditional vertical MIM (metal-insulator-metal) memristor voltages are around/above 1 V. [20,21] Lateral TMDC memristors (Figure 1), on the other hand, offer a higher degree of functionality since lateral structures can be equipped with additional terminals such as top or bottom gates. [22]18][19] While early lateral TMDC memristors were relatively large (channel length of several μm) and required high voltages for memristive switching (>15 V), [16,17,23] recently, devices with shorter channel lengths and much lower switching voltages have The three illustrated mechanisms have been suggested as the origin of hysteresis in lateral memristive devices based on 2D TMDC materials: a) the charging and discharging dynamics of trap states, [12][13][14][15][16][17][18] b) the migration of charged point defects in the electric field, [23,38] and c) the voltage-induced phase change from a stable, low-conductive semiconducting phase to a metastable, high-conductive metallic phase. [42,44]en demonstrated.For example, in [24] a gated lateral MoS 2 memristor with a channel length of ≈330 nm and a switching voltage as low as 0.42 V has been reported.[27] From both the physics and processing points of view, further scaling of the channel length of such devices can be expected with a scaling limit in the same range as that for the channel length of TMDC MOSFETs (metaloxide-semiconductor field-effect transistors).Unfortunately, the physics of lateral TMDC memristors is not yet understood, and even the mechanisms responsible for resistive switching in these devices are still debated.Indications for three different switching mechanisms were found experimentally.
30][31][32][33][34] These states could be introduced by defects far away from the electrodes (in the bulk) or at the semiconductor-electrode interfaces. [35,36]They can capture and release electrons or holes with a probability depending on the applied voltage, thereby altering the bulk and interface electrostatic potential, free charge carrier density, and contact barriers. [37]Experimental investigations based on current-voltage measurements and analytical approximations resulted in estimated trap area densities between ≈ 10 12 cm −2 [28] and ≈ 10 15 cm −2 [31] for MoS 2 .The second mechanism (Figure 1b) is based on charged point defects sufficiently mobile to drift in the applied electric field.Experimental indications were found for the accumulation and drift of sulfur vacancies along grain boundaries in the conducting channel of lateral monolayer MoS 2 devices. [38]This observation is consistent with other experimental work [35,39] and was correlated with hysteresis in the current-voltage characteristics. [38]Elsewhere, [23] measurements were performed on exfoliated MoS 2 multilayers, which were plasma-treated to increase the concentration of sulfur vacancies.The homogeneous drift of mobile sulfur vacancies in the conduction channel was confirmed by local Augerelectron spectroscopy and experimentally correlated with hysteresis in the I-V characteristics. [23]Sulfur vacancies could behave as n-type dopants in MoS 2 , [40,41] thereby changing the conductivity in the bulk material and modifying the potential energy barrier at the semiconductor-metal interfaces.The third mechanism (Figure 1c) is based on an electric-field-induced transition from a low-conductive semiconducting phase to a highly conductive metallic phase of the TMDC material.Such an effect has been demonstrated in lateral MoTe 2 devices, where a large electric field caused a transition from the semiconducting 2H phase to a metastable high-conductive 1T' metallic phase [42] and in 2H-Te 2 with a transition from the 2H phase to an unstable distorted 2H d phase. [43]In MoS 2, this effect was achieved indirectly by a preintercalation of Li + ions, which are highly mobile in between the monolayer sheets and induced the 2H to 1T' transition upon drift and agglomeration in the applied electric field. [44]Despite this experimental progress, the switching mechanisms in 2D TMDC devices lack a profound physical understanding.
In particular, the first two mechanisms seem promising for application [45] but are superficially similar because the change in conductance is based on a change in the charge density in both cases.Further, the two mechanisms could be coupled and occur in parallel since point defects can introduce additional electronic states. [46,47]These similarities make distinguishing the two mechanisms challenging or even impossible with experimental correlation studies.While trap-state-based hysteresis is assumed in most work on lateral devices to explain resistive switching, [12][13][14][15][16][17][18] it is unclear whether the drift of charged point defects could contribute as well or even dominate the hysteresis in some cases.Quantitative physical models could be powerful in separating such effects, which are experimentally difficult to access or distinguish.
However, strong simplifications are usually imposed on the physics, particularly on the charge carrier dynamics.The few reported models of lateral devices are based on analytical approximations of the I-V curve with many fitting parameters, [16] compact-model formulations. [48,49]or incompletely coupled approaches in adiabatic approximation, omitting dynamic effects. [50]While such models can be helpful in circuit simulations and device design, they do not capture the complex dynamics of the switching process crucial for device functionality.
Further, image-charge-induced Schottky barrier lowering (SBL) is frequently assumed to be highly relevant for interfacetype switching. [23,48,49]This effect describes the change of the Schottky barriers caused by an applied voltage and a dynamic redistribution of charge carriers in the semiconductor. [51,52]ence, SBL could lead to hysteresis by altering the interfacelimited conduction.Despite its fundamental character, up to now, there is no physically consistent mathematical formulation of this effect for charge-transport models, and the role of SBL and defect dynamics in the hysteresis remains unclear.
In this work, we analyze the role of mobile point defects in the switching mechanism of 2D memristive devices and demonstrate the dominant impact of their dynamics on the hysteresis and the device characteristics.A fully dynamic charge-transport model is introduced and validated with experimental data of lateral memristive devices based on exfoliated MoS 2 .The model self-consistently solves the dynamic drift-diffusion equations of electrons, holes, and ionic point defects together with the electrostatic potential.The presented formulation allows us to consider the nonlinear dynamics of mobile point defects, including their different behavior compared to electrons and holes.The effects of local Joule heating on the electronic properties and carrier dynamics are considered via the transient heat equation.Additionally, we develop a mathematical formulation of imagecharge-induced Schottky barrier lowering (SBL) to analyze its effect on the switching mechanism depending on the carrier dynamics.The implementation permits fully dynamic calculations over many voltage cycles and large time spans, including th first large-scale pulse simulations.
The paper is structured as follows.First, we present the system of model equations in general (Section 2), without SBL (Subsections 2.1 and 2.2), and with SBL (Subsection 2.3).Then, the mathematical model is applied to the geometry and material parameters of the simulated devices (Section 3).In Section 4, we validate the model with experimental data and present a detailed analysis of the switching mechanism.The results are discussed in Section 5.

Charge Transport Model
The charge transport model merges the semi-classical driftdiffusion equations for electrons, holes, and charged mobile point defects with Poisson's equation for the electrostatic po-tential in a self-consistent and fully coupled manner.The driftdiffusion equations are introduced in a quasi-Fermi potential formulation with different carrier statistics functions derived thermodynamically for the respective carrier species.This allows us to effectively consider the different behavior of mobile point defects compared to electrons and holes, such as their different nonlinear diffusion and limited maximum concentration prescribed by the density of available lattice sites in the crystal.To address the impact of applied voltage and charge carriers on the Schottky barriers, we develop a self-consistent formulation of imagecharge-induced Schottky barrier lowering (SBL).It requires modifying the equation system and solving for an additional potential.The equation system is resolved using a physics-preserving finite-volume discretization.
Additionally, we analyze the relevance of local self-heating effects to explore the dependency of the mobilities on the temperature and the electric field via an ion hopping model.This is done by coupling the electronic properties, i.e., the effective density of states, the mobilities, and the Fermi-Dirac integral, via the temperature to the transient heat equation.

Bulk Equations
Without the effects of SBL and Joule heating, the chargetransport model comprises a system of four fully coupled nonlinear partial differential equations: three drift-diffusion equations for the quasi-Fermi potentials   (x,t),  ∈ {n, p, x} of electrons n, holes p, and mobile defects x, as well as Poisson's equation for the electrostatic potential (x, t).The corresponding densities of the charge carriers are denoted by n  (x,t),  ∈ {n, p, x}.For the location vector x and the time t it is (x, t) ∈  × [0, t e ] for the spatial domain of the semiconductor  ⊂ ℝ d (d = 1, 2, 3) and some given end time t e > 0.
Charge carrier continuity equations.The drift and diffusion of electron, hole, and ionic defect densities n n , n p , and n x is described via the continuity equations where the charge q  of the respective carrier species  is given by with the positive elementary charge q.The charge numbers are z n = −1 and z p = +1 for electrons and holes, respectively, and z x ∈ ℤ for mobile defects.The source term q  r  on the right-hand side of Equation (1) includes the production/reduction rates r  , which can be used to consider the recombination and generation processes of charge carriers.The current densities j  , resulting from the drift and diffusion of  are proportional to the negative gradient − ∇  of their respective quasi-Fermi potentials   , and are given by with the charge carrier mobilities   .The charge carrier densities n  can be linked to the respective quasi-Fermi potentials   and the electrostatic potential  via the state equations [53] n with the temperature T, the Boltzmann constant k B , the effective densities of state of the conduction band and valence band N c = N n and N v = N p , and the statistics functions   .We use different statistics functions depending on the type of charge carriers to reflect their specific physical properties.For the statistics functions  n and  p of electrons and holes, we use the Fermi-Dirac integral of order 1/2, suitable for semiconductors with parabolic bands.This choice leads to nonlinear diffusion in the current densities (Equation ( 3)) as opposed to the exponential Boltzmann approximation. [53]In contrast, mobile point defects exhibit different behavior.Their concentration n x is limited by the density of available lattice sites in the crystal, resulting in nonlinear ionic diffusion. [54]The Boltzmann approximation often used for ionic carriers (e.g., in [50,55] ) does not account for this non-linearity, allowing non-physical behavior such as an excess of vacancies over available vacancy sites.To reflect experimentally observed ionic depletion accurately, we use the Fermi-Dirac integral of order −1, denoted as  x () = 1∕(exp(−) + 1), following a grand canonical formalism of an ideal lattice gas based on thermodynamic principles. [56]Here, the saturation limit for the defect density is represented by N x , such that n x < N x , and E x,0 represents the intrinsic energy level of the defects. [56]urther, E c,0 = E n,0 and E v,0 = E p,0 denote the intrinsic band edge energies of the conduction and valence bands for electrons and holes, respectively.These energies relate to the electron affinity  e and the bandgap E g via E c,0 ≔ −  e and E v,0 ≔ −  e − E g .Consistent with Eq. ( 4), the electric potential-dependent conduction band edge E c , valence band edge E v , and vacancy energy level E x are defined as For the effective densities of state we use the threedimensional parabolic band approximation ) 3∕2 (7)   with the reduced Planck constant ℏ, and the effective mass m * n and m * p of electrons and holes, respectively.Poisson equation.Finally, the three continuity equations (Equation (1)) are self-consistently coupled via the electrostatic potential  to the nonlinear Poisson equation In this equation, C is either an effective acceptor or donor background dopant density with dopant charge q C = qz C , described by the background dopant charge number z C ∈ { − 1, 1}.The quasistatic electric permittivity ɛ = ɛ 0 ɛ r is expressed as the product of the electric vacuum permittivity ɛ 0 and the quasi-static relative electric permittivity ɛ r .
Heat equation and carrier mobilities.Generally, large electric fields and high current densities can increase the temperature of the device via Joule heating. [57]Such self-heating, in turn, can influence the electronic properties, e.g., via the temperaturedependent carrier densities (Equation ( 4)) and the mobilities.To consider such effects, we couple the electron and hole continuity equations (Equation ( 1)) to Poisson's equation and the transient heat equation for the temperature T with the given mass density , heat capacity C p at constant pressure, and thermal conductivity .The locally inserted power density causing a potential temperature rise is given by the scalar product of the current density j = j n + j p and the electric field E = − ∇. [58]he influence of the temperature and the electric field on vacancy mobility is then calculated using the ion hopping model of Genreith-Schriever. [59]This model has been thoroughly validated with atomistic simulations and extends the validity of the Mott and Gurney approximation to large electric fields. [59].Similarly, we estimate the effect of velocity saturation for the drift of electrons and holes at large electric fields.The detailed analysis, including the boundary conditions, model setup, and a careful comparison with experimental data, is presented in the Supplementary.The results demonstrate clearly that Joule heating effects and electric field dependencies of the mobilities are negligible for all cases relevant here.Therefore, we can reduce the computation effort by excluding Joule heating and the electric field dependency of the mobilities from the model for the following simulations.

Boundary Conditions at Metal-Semiconductor Contacts
The metal-semiconductor contacts , are assumed to represent physical barriers for mobile defects, and therefore, we impose the zero-flux condition (for x ∈ ) with the outward directed unit normal vector  of the semiconductor domain.For electrons and holes, thermionic emission is described with the flux boundary conditions [60] j n (x, t) Here, v n and v p are the electron and hole recombination velocities given by (with h as Planck constant) with the effective mass m * n and m * p of electrons and holes, respectively.The corresponding quasi-equilibrium carrier densities n n,0 and n p,0 at the contacts are where ϕ 0 = ϕ 0 (x) > 0 is an intrinsic Schottky energy barrier constant in time and For the electrostatic potential at the metal-semiconductor contacts (x ∈ ) we apply the Dirichlet condition with  0 = (ϕ 0 − E c,0 )/q n denoting an intrinsic electrostatic potential barrier, and  a denoting an applied electrostatic potential at the contact.The intrinsic Schottky barrier ϕ 0 is considered as a material and device parameter, i.e., for a given device,  0 and ϕ 0 depend only on the location x.Zero-flux conditions can be used at all non-contact boundaries.

Image-Charge-Induced Schottky Barrier Lowering
The two subsections above describe a system where the electrostatic potential at the boundary is prescribed by an intrinsic Schottky barrier, and only the applied electrostatic potential  a adds a time-dependent contribution.However, any point charge present in the semiconductor induces a charge of opposite polarity (image charge) in the metal electrode.This adds a contribution ϕ i (image-charge potential energy) to the potential energy of the charge with an energy minimum at the contact interface, [60] as schematically illustrated in Figure 2a.The image-charge potential energy ϕ i superposes with the residual potential energy  r ≔ q n  r ( r : residual potential), which corresponds to the potential energy without image charge.As a result, the barrier maximum of the total potential energy at the interface is reduced by a value |Δϕ| compared to the case without image-charge contribution (Figure 2b).While the image-charge potential energy is constant in time, ϕ r depends on the charge distribution within the semiconductor and the applied potential  a ≔ q n  a , which can correspondingly influence Δϕ.Consequently, the contact potentials become functions of the applied voltage and implicitly of the charge carrier densities.
Including image-charge-induced barrier lowering requires modifying the bulk equations as well as the boundary conditions.First, we modify the boundary condition for the electrostatic potential (Equation ( 14)) at the metal-semiconductor contact (x ∈ ) by adding the change Δ ≔ Δ∕q n of the potential barrier induced by the image charge.The new boundary condition reads Within the classical approximation, [60] the change Δ of the potential barrier induced by the image charge is only defined for an upwards-bent conduction band edge (∇ v  r < 0, x ∈ ) and can be expressed as a function of the projected gradient of the residual electrostatic potential In this equation, we refer to  i ≔ Re{ * (f i )} as the image-charge dielectric constant, defined as the real part Re{ɛ*(f i )} of the complex dielectric permittivity ɛ* at the frequency f i = 1/T i , which is connected to the time T i an electron needs to travel from the barrier maximum to the metal-semiconductor interface (see Methods).The residual electrostatic potential  r is defined as the electrostatic potential obtained as a solution to Poisson's equation which is now linear in  r because n  is a function of .At the contact boundaries,  r is obtained using the Schottky contact boundary condition Further, the equations for the equilibrium charge carrier densities at the contacts (x ∈ ) in Equation ( 13) must be altered to To consider the change of the barriers in equilibrium.If the effect of SBL is considered, we consequently solve a system of five coupled partial differential equations (Equation (1), Equation ( 8), and Equation ( 17)) for five solution fields, namely the three quasi-Fermi potentials  n ,  p ,  x , and the two electrostatic potentials  and  r .The equation system is supplemented with the equilibrium solution as initial conditions.

Application to Lateral Memristive Devices
In the following, the equations from Section 2 are applied to the geometry of the lateral devices studied.First, the model geometry and application of boundary conditions are specified in Section 3.1, before the material parameters are discussed in Section 3.2.A summary of all model parameters and values used for the simulations is provided in the Methods.

Model Geometry and Application of Boundary Conditions
An illustration of the geometry of a typical lateral memristive device is shown in Figure 3a, along with its cross-section in the x-y plane in Figure 3b.The device comprises a TMDC flake on a SiO 2 /Si substrate, and the flake is laterally sandwiched by two metal electrodes.Typical TMDC flakes in lateral memristive devices have a thickness D of some nanometers and a width W and length L of up to a few micrometers (e.g., [33,36,42,[61][62][63] ).In this study, we use representative device geometries with two different channel lengths of L ≈ 1 μm and L ≈ 2 μm, thicknesses of D ≈15 nm (≈25 monolayers), and channel widths of W ≈ 10 μm.The flakes consist of exfoliated MoS 2 (details in Section 3.2) and are contacted by two Ti (5 nm)/Au (50 nm) side electrodes (Figure 3ab).Because the devices have side contacts, the charge carriers are expected to migrate primarily along the x-axis.[66] Additionally, conduction through the substrate can be omitted because of the large resistivity of SiO 2 .
These aspects are intrinsically considered by the model geometry shown in Figure 3c.We model the device by a 1D conducting channel defined by the flake length L between the electrodes and the flake cross-section A = DW (Figure 3c).The 1D model comprises the semiconductor (TMDC) domain of length L with two metal-semiconductor contacts at the locations x 1 and x 2 , where we apply the corresponding contact boundary conditions (Sections 2.2 and 2.3), i.e., either Schottky boundary conditions with SBL or without SBL.For the applied potential, we assume  a (x 1 ) = 0 V at the left contact and  a (x 2 ) =  a (x 1 ) + U(t) at the right contact.The voltage U(t) is a piecewise linear interpolation with four support points within the period of length t c and the extrema U max and U min = −U max (Figure 1d).

Material Parameters
We use mechanically exfoliated MoS 2 as a well-known example material for the TMDC domain to gain general insights into the hysteresis mechanism.Because of the high formation energy of sulfur vacancies of up to several eV depending on the conditions, [46,47,67,68] the number of sulfur vacancies remains approximately constant during device operation at sufficiently small electric fields This has been confirmed experimentally by many studies (e.g., [66,69,70] ) since they required electron-beam irradiation [66,69] or plasma treatments [70] to introduce significant vacancy densities in exfoliated MoS 2 flakes.The data we use for validation are taken from devices fabricated with a mechanical shear exfoliation technique followed by a plasma treatment step to induce sulfur vacancies, as described comprehensively in. [23]Auger-electron spectroscopy measurements on these devices confirmed a large density of sulfur vacancies sufficiently mobile to drift continuously under an applied electric field without inducing filament formation. [23]This is consistent with the stable device characteristics measured without current compliance over many voltage cycles under maximum applied voltages of 20 V. [23] Less stable device characteristics and an eventual change from analog to discrete switching was only observed after repeatedly applying voltages between 20 and 50 V, [23] indicating the additional formation of vacancies in this operating regime.All data we use for validation were measured in the stable regime with voltages U ≪ 20 V (Table 4 and Table 5, Methods).Therefore, we can assume that the number of sulfur vacancies remains sufficiently constant during device operation.
While the role of sulfur vacancies for doping and their stable charge states have been debated, [46,66,71] for large concentrations, they are expected to behave as n-type dopants. [71]Here, we assume that the mobile ionic species x represents donor-type sulfur vacancies, single-positively charged with z x = +1 in their ionized state.The maximum density of vacancies N x is estimated from the density of sulfur sites to N x ≈ 4 × 10 28 m −3 (see Methods).Additionally, the negligible vacancy generation and significant ntype doping (i.e., n n ≫ n p ) (see Supplementary) is considered in the model by setting the source term in Equation ( 1) to r  = 0.
The electronic properties of MoS 2 are well investigated and depend significantly on the number of monolayers and other extrinsic and intrinsic factors, such as the stress state and the temperature. [40,72,73]We review and estimate the range of material parameters expected for the most stable phase 2H-MoS 2 and a layer thickness of D = 10 − 15 nm (≈16-25 monolayers [74] ).Note that the parameters are summarized for a temperature of T = 300 K because Joule heating is negligible in the analyzed devices (see Supplementary).
With increasing thickness, the conduction band minimum (CBM) and the valence band maximum (VBM) shift from the K point (in the monolayer) to a midpoint between Γ and K (CBM) and the Γ point (VBM) (in the bulk). [75]A transformation occurs from a direct semiconductor to an indirect one, [61,[76][77][78] with a gradual shift of the electronic bandgap from E g = 1.8 eV (monolayer) to E g = 1.2 eV (bulk).This shift has been explained theoretically by an increasing contribution of interfacial stress to the overall atomic configuration as the material gets thinner [77] and is consistent with ab initio calculations for the influence of stress on the electronic properties of MoS 2 monolayers. [79,80]Assuming no external or internal strain of the MoS 2 layer, e.g., from the fabrication process [81] or defect clusters, [66] simulations predict values of E g ≈ 1.3 eV for a layer thickness of 15 nm. [77]The changing band structure with layer thickness also alters the effective masses m * n and m * p of electrons and holes.For 2H-MoS 2 monolayers, values at the K point of m * n = (0.45 − 0.57) m 0 and m * p = (0.59 − 0.64) m 0 are reported, [75,80,82] with geometric means of m * n ≈ 0.54 m 0 and m * p ≈ 0.61 m 0 [82] considering the slight anisotropy (m 0 : electron rest mass).This is slightly smaller than the effective masses of the bulk with values of m * n ≈ 0.55 m 0 and m * p ≈ 0.71 m 0 . [75]For a layer thickness of 15 nm, we expect effective masses close to the values in the bulk.
Experimental and theoretical studies on the in-plane quasistatic dielectric permittivity of monolayer MoS 2 report a large scatter of values of ɛ r ≈ 3.7, [83][84][85] and ɛ r ≈ 10 − 17. [72] Generally, ɛ r tends to increase as a function of the film thickness and values of ɛ r ≈ 10 were measured for layer thicknesses of D = 10 − 15 nm (close to ɛ r ≈ 12 in the bulk). [86]While these results are in excellent agreement with ab initio simulations, [87] other calculations yield smaller values of ɛ r = 7.1 for the bulk. [85]For the imagecharge dielectric constant, the frequency dependency of the complex electric permittivity becomes relevant, but our estimation shows ɛ i ≈ ɛ r as a reasonable approximation (see Methods).Further, electron affinities of  e ≈ 3.7 eV, [88]  e ≈ 4.0 eV [89,90] and  e ≈ 4.3 eV [91,92] are reported for 2H-MoS 2 .A summary of the parameters and their estimated range is provided in Table 2 (see Methods).
While the abovementioned parameters are obtained directly from the literature and used for all simulations, other parameters are expected to be highly sensitive to the sample's microstructure, defect composition, and interface quality.These sampledependent model parameters include the mobilities   , the intrinsic Schottky barriers ϕ 0 and the defect energy level E x,0 or the corresponding vacancy density n x .For the electron and hole mobilities, a large range of  n,p ≈ 0.1 − 200 cm 2 /(Vs) [40,81,[93][94][95] is reported.This range has been explained by the scattering of charge carriers with interfacial Coulomb impurities. [93]As a result, the mobilities tend to increase with the layer thickness and depend significantly on the gate dielectric and interface quality of the respective device. [81,93]Note that the electric field dependency of the mobilities is negligible in the analyzed devices (see Supplementary).Ab-initio studies find activation energies of diffusion of E a ≈ 0.8 − 2.5 eV for sulfur vacancies in MoS 2 , depending strongly on the spatial configuration of the point defects. [67]his energy barrier is expected to be reduced along grain boundaries or interfaces, [38] and a previous estimation based on Augerelectron spectroscopy measurements resulted in a much smaller value of E a = 0.6 eV [23] .These measurements were performed on lateral memristive devices with vacancy densities of approximately n x = 10 26 − 10 27 m −3 [23] , assuming N x = 4 × 10 28 m −3 (see Methods).With a work function of approximately 4.3 eV for Ti [96] and 4.6 eV [97] for pristine n-type MoS 2 , nonrectifying contacts are anticipated following the Schottky-Mott rule [23] .Experiments demonstrated slightly rectifying contacts with reported Schottky barriers of ≈ 0.05 eV [98] and 0.18 eV [96,99] caused by the effect of Fermi-level pinning.

Results
In the following, we analyze the origin of hysteresis, the influence of vacancy dynamics, and image-charge-induced Schottky barrier lowering in general for lateral memristive devices.We use parameters within the physically realistic range of MoS 2 (Section 3) as representative material and discuss the cases of vanishing Schottky barriers and significant Schottky barrier heights.Finally, we compare additional simulations with various experimental data sets to demonstrate the model's validity.
For the first case, we consider a parameter set S 1 (Table 1, Methods) with omittable Schottky barriers ϕ 0 (x 1 ) = ϕ 0 (x 2 ) ). e-j) Example configurations of the vacancy density n x (x) as a function of space x for various positions in the I-V curve shown in (a), together with the quasi-static equilibrium configuration of n x (gray), i.e., the configuration that would be reached at a constant applied voltage for t → ∞.The measured data are taken from [23] = 1 meV.The experimental data used to obtain the sampledependent model parameters are taken from [23] and are shown in Figure 4a together with the simulation.The current magnitude |I| is a smooth function of the applied voltage U without visible discontinuities and with approximately equal maximum current magnitudes at U = 13 V and U = − 13 V.While |I(U)| is overall symmetric around U = 0 V, the right hysteresis branch (U > 0 V) has a slightly larger area compared to the left hysteresis branch (U < 0 V).As indicated, the direction of the hysteresis is clockwise and counterclockwise in the right and left hysteresis branch, respectively.The simulations of eight consecutive voltage cycles are shown in Figure 4b, revealing a significant difference in the I-V characteristics of the first cycle compared to the following cycles.In particular, the I-V curve of the first cycle is highly asymmetric, with an initial rise in the current (Figure 4b,1), which quickly saturates and almost remains constant until the maximum voltage is reached.In the left hysteresis branch (Figure 4b, U < 0 V) a crossing of the current is visible, resulting in a change of the hysteresis direction from clockwise (Figure 4b, 3-6) to counterclockwise Figures 4b, 4-5).The other voltage cycles (Figure 4b, Cycle 2-8) show a different but highly reproducible I-V behavior and reflect the characteristics of the measured I-V curve mentioned above very well, including the hysteresis direction and the transitions from the low-resistive state (LRS) to the high-resistive state (HRS) and vice versa.
The origin of the two I-V characteristics is the different initial configurations of the system at the beginning of the respective cycle.At the beginning of the first cycle (t = 0 s, U = 0 V), the modeled system is in equilibrium, i.e., all quasi-Fermi levels are constant q n  n = q n  p = q x  x = 0 eV, as shown in the reduced band diagram (Figure 4c).During the first cycle and at the beginning of the second cycle (t = 10.4 s, U = 0 V) the electrons and holes remain close to their quasi-static equilibrium configuration, i.e., the configuration that would be reached at a constant applied voltage for t → ∞.The small mobility of vacancies on the other hand results in a substantial deviation from equilibrium (Figure 4d).This deviation from equilibrium is apparent from the notable drop in the vacancy quasi-Fermi level q x  x at the left metal-semiconductor contact around x = 0.Such a nonequilibrium configuration is reproducibly reached after each cycle following the first one.This phenomenon appears to be superficially similar to electroforming in filamentary devices.There, the formation of a filament is induced during the initial voltage cycles and then modulated in the subsequent cycles. [100]However, in our model, vacancies drift homogeneously distributed over the cross-section of the conducting channel, in line with measurements of the local vacancy density in these devices. [23]We will show in the following that the difference between the first and subsequent hysteresis loops is just a special case of a more general dynamic phenomenon.In the following, we focus the discussion on the second cycle as a more reproducible and representative device characteristic.
In Figure 4e-j, the vacancy density is shown at various voltages in the second hysteresis cycle.Additionally, its quasistatic equilibrium configuration at the respective voltage is shown (gray, faint).At the beginning of the cycle, the drop in q x  x (Figure 4d) results in a local depletion of vacancies at the left electrode (n x ≪ 10 10 m −3 , Figure 4e).This depletion zone initially limits the current, but it vanishes quickly as the voltage increases (Figure 4e-f).Therefore, the device is in a high-conductance (LRS) state at the beginning of the second cycle.As the voltage is further ramped up, the current increases continuously until a depletion zone forms and grows at the right electrode (Figure 4g).It limits the current just before the maximum voltage is reached.The depletion zone continues growing at the right electrode and reduces the current as the voltage is ramped down (Figure 4g-i).Hence, the device is in the HRS.At U = 0 V, the system strives to restore equilibrium where no depletion zone exists.Because of the small vacancy mobility, the depletion zone is not entirely annihilated at the beginning of the left hysteresis branch (Figure 4j, U < 0 V) leading to a vacancy density configuration similar to the one at the beginning of the right cycle.Therefore, the I-V curve is overall symmetric around U = 0 V, and the device starts in a high conductance state (LRS) at the beginning of the left hysteresis branch.The electric fields and electron and hole densities corresponding to the six different points in the hysteresis are provided in the Supplementary.
To analyze the dependency of the I-V curve on the vacancy mobility in detail, we performed simulations using the same parameter set S 1 (Table 4, Methods) as before but with different values for the vacancy mobility, ranging from  x = 1 • 10 −16 m 2 /(Vs) to  x = 1 • 10 −12 m 2 /(Vs) for two different voltage sweep rates f 1 = 5 V/s and f 2 = 2.5 V/s.We define the hysteresis areas A H of the left and right hysteresis branches as the integral over the respective branch and consider the sign of the integral to distinguish the direction of the hysteresis.The hysteresis areas of the second cycle are plotted in Figure 5h as functions of  x for f 1 and f 2 .Example values of A H calculated from the hysteresis curves in As the mobility exceeds small values  x ≤ 1 × 10 −15 m 2 /(Vs), the hysteresis starts to open (A H ≈ 0) and is initially clockwise oriented in both branches (Figure 5a-b).Because the depletion zone is still present when U < 0 V is reached the device is in the low conducting state at the beginning of the left hysteresis branch.A transition of the hysteresis direction occurs in the left hysteresis branch from clockwise (A H < 0) to counterclockwise (A H > 0) at intermediate mobilities around  x ≈ 3 × 10 −14 m 2 /(Vs) (Figure 5c).During the transition, the current in the left hysteresis branch intersects at a voltage U c (Figure 5c, indicated), which continuously moves towards U = 0 V as the mobility increases until the entire hysteresis branch is directed counterclockwise (Figure 5d).Origin of this transition is the shorter time required for vacancies with a larger mobility to reach quasi-static equilibrium, which causes a faster annihilation of the depletion zone before U < 0 V is reached.Further increasing the vacancy mobility to large values (Figure 5d-f) results in the formation of distinct maxima at voltages U < U max (and U > U min ), i.e., a negative differential resistance, an increasingly symmetric I-V curve, and an overall reduced maximum current.This behavior reflects the formation dynamics of the depletion zone.In the high-mobility regime, the depletion zone can form before the maximum positive or negative voltage is reached, which results in a drop in the current, leading to the formation of the maxima.Simultaneously, The measured data are taken from. [23]e vacancy density is close to its quasi-static equilibrium configuration, which results in the pronounced symmetry around U = 0 V and U c ≈ 0 V.The dependency of − U c /U max as a function of the vacancy mobility is detailed in Figure 5g for the first 20 voltage cycles as a measure for the symmetry of the I-V curve.Over the entire range of  x , the voltage U c always reaches U c = 0 after a sufficiently large number of voltage cycles.While > 20 cycles are required to reach the symmetric I-V curve for  x < 2 × 10 −14 m 2 /(Vs), one to two cycles are enough for values close to  x = 1 × 10 −13 m 2 /(Vs).
The results show that the hysteresis and its symmetry are governed by the formation and annihilation dynamics of the vacancy depletion zone, i.e., the deviation of the vacancy density from its quasi-static equilibrium configuration.Hysteresis occurs in a small vacancy mobility range where the vacancies are sufficiently mobile to follow the potential gradient but slow enough to be in nonequilibrium.Therefore, a change in the voltage sweep rate only shifts the A H ( x ) curve on the vacancy mobility axis (Figure 5h, f 2 ).As another consequence, the asymmetry and intersection in the left hysteresis branch are volatile features.They depend on vacancy mobility and the cycle number and eventually vanish after enough sweep cycles.
In previous experimental work [23] a hysteresis reversal similar to the transition shown in Figure 5b-d was achieved by increasing the maximum voltage from 25 V to large values of 50 V.This was accompanied by a significant increase in the maximum current.Additional results presented and discussed in the Supplementary demonstrate that our model can reproduce this behavior very well.It explains the reversal and increase in the maximum current by an increase in the vacancy density induced by the large voltages consistent with the experiments.An increasing maximum voltage alone results in neither a hysteresis reversal nor sufficiently increased maximum currents.
In these first examples, the intrinsic Schottky barriers are so small that SBL is negligible over the entire voltage range.However, the influence of the Schottky barriers and SBL on the I-V curve could generally contribute to the hysteresis behavior and is therefore explored next in the second example.
For the second case, we consider a similar parameter set S 2 (Table 4, Methods) as before but introduce intrinsic Schottky barriers ϕ 0 (x 1 ) = 0.144 eV and ϕ 0 (x 2 ) = 0.11 eV with a slightly smaller value at the right contact. Figure 6a demonstrates the excellent match of the simulated and measured I-V curve.Both have the same clockwise direction (right branch) and counterclockwise direction (left branch) as in Figure 4a, but the asymmetry around U = 0 V is more pronounced.In contrast to the previous I-V curve (Figure 4a), the two maximum current magnitudes are notably different, and the right hysteresis area is smaller than the left one.Because the intrinsic Schottky barriers are small but significant, SBL can occur.In equilibrium at t = 0 (Figure 6b), the SBL notably reduces the contact barriers by ≈25% and 18% to values of ϕ(x 1 ) ≈ 0.11 eV and ϕ(x 2 ) ≈ 0.09 eV.This difference in the reduction of the Schottky barriers is caused by the difference in the intrinsic Schottky barrier ϕ 0 (x 1 ) and ϕ 0 (x 2 ) at the contacts.At a large intrinsic Schottky barrier, the magnitude of the residual potential gradient |∇ v  r | tends to be larger than at a small barrier.This gradient determines the reduction of the Schottky barrier via Equation ( 16), and therefore, SBL tends to reduce the difference in the Schottky barriers compared to the intrinsic barriers.
This affects the band edges E c , E v , and E x in the proximity of the contact interfaces, correspondingly (Figure 6b, E v and E x not shown).In Figure 6d, the applied voltage U, the current magnitude |I|, and the change Δϕ(x 1 ) and Δϕ(x 2 ) of the Schottky barriers are plotted as functions of the time t over the first two voltage cycles.As the voltage increases (U > 0 V) the band edges at the left contact remain bent upwards (∇ v  r < 0), which results in an increased Schottky barrier Δϕ(x 1 ) > 0 (Equation ( 16)).At the right contact, the band edges move to lower energies, and the projected electrostatic potential quickly reaches positive values (∇ v  r ≥ 0) indicating a downwards bent conduction band.As a result, the change in the right Schottky barrier remains zero, i.e., Δϕ(x 2 ) = 0, over most of the right hysteresis cycle (U > 0 V) and does not limit the current.As the voltage approaches negative values, the situation reverses and the change in Schottky barriers is Δϕ(x 1 ) = 0 and Δϕ(x 2 ) > 0 over most of the left hysteresis cycle (U < 0 V).Note that both Δϕ(t) are asymmetric around the time at which the respective maximum positive and negative voltage is applied.This asymmetry is also present in the current magnitude and reflects the hysteresis in the time domain.
To identify the influence of SBL on the I-V characteristics, we performed a simulation with the same parameters as for the curve in Figure 6a but without SBL.The two I-V curves with and without SBL are compared in Figure 6c, and qualitative differences are apparent.Without SBL, the hysteresis direction of the right branch is reversed, and a pronounced asymmetry is visible around U = 0 V, i.e., largely different maximum currents and hysteresis areas of the left and right branches.This asymmetry is caused by a larger difference in the equilibrium Schottky barriers if SBL is not considered.The total hysteresis area without SBL is similar to the case with SBL.An additional I-V simulation was performed without SBL using the equilibrium Schottky barriers ϕ(x) obtained from the simulation with SBL as intrinsic Schottky barriers ϕ 0 (Figure 6c, No SBL, reduced ϕ 0 ).The resulting I-V curve matches the one with SBL.This comparison shows that the hysteresis is not dominated by a change in the interfacial current but a limitation in the bulk conductance via the formation of the depletion zone.Consequently, SBL greatly affects the I-V curve, but it does not contribute significantly to the hysteresis area in this example.
The effect of a change in the Schottky barriers is illustrated with additional simulations in Figure 7, including SBL.We used the parameter set S 2 (Table 4, Methods) as before, but various intrinsic Schottky barriers for the right contact (Figure 7a) and the left contact (Figure 7b) with values from 0.12 eV to 0.18 eV while keeping the respective other Schottky barrier fixed to a value of 0.11 eV.In both cases, we observe a transition from a relatively symmetric I-V curve for the case of similar Schottky barriers to a notably asymmetric curve for significantly different Schottky barriers.This transition is accompanied by a change of the hysteresis direction in the left hysteresis branch if ϕ 0 (x 1 ) < ϕ 0 (x 2 ) and in the right hysteresis branch if ϕ 0 (x 2 ) < ϕ 0 (x 1 ).A comparison of the I-V curve in Figure 7b (right, ϕ 0 (x 2 ) = 0.18 eV) with the hysteresis curve in Figure 6c demonstrates that the behavior observed without SBL is mainly caused by the difference in the two Schottky barriers of the left and right contact.Other features, such as the formation of a maximum at U < U max and the reduction of the maximum current magnitude with increasing barrier difference are also observed (Figure 7).These features are superficially similar to those caused by the increasing vacancy mobility of the device with symmetric contacts, but their origin differs.The asymmetric Schottky barriers result in different restoring forces ).The measured data are taken from, [23] and the simulation parameters are provided in the Methods.compared to the symmetric configuration, which alters the formation and annihilation dynamics of the depletion zone.A detailed explanation based on a comparison of the electric fields, carrier densities, and driving forces of the symmetric and asymmetric case is provided in the Supplementary.
To analyze the I-V dynamics caused by asymmetric Schottky contacts, we simulate the I-V curve of an example configuration over 100 consecutive voltage cycles (Figure 7c).An asymmetry is apparent in cycle one, which reduces slightly until cycle ten is reached.All following voltage cycles are approximately identical to cycle ten and retain a significant asymmetry demonstrated in the example of cycle 100.As a quantitative measure, we extract the total hysteresis area A H,n normalized to its value at cycle 100 and plot it as a function of the number of cycles for simulations with three different vacancy mobilities (Figure 7d).All three functions reach A H,n = 1 and remain constant after enough cycles.A larger vacancy mobility  x reduces the number of cycles required to reach the dynamically stable configuration.Consequently, the asymmetry in the I-V characteristics of different Schottky barriers and hysteresis reversal is not a mere volatile feature, as observed in the case of symmetric contacts.Instead, it remains present after a sufficiently large number of cycles reflecting the asymmetric restoring force in equilibrium.
For a more detailed validation, we fitted the model to additional sets of measured I-V curves with various maximum voltages and different extends of asymmetry.The results in Figure 8 show an overall excellent agreement with the measured data.Only minor deviations are apparent, e.g., in the right hysteresis branch of the asymmetric curves in Figure 8a,c, with a slightly larger current in the simulation compared to the measurement.All parameters are within the expected range estimated in Section 3 (parameter sets S 3 -S 6 , Methods).
For applications such as analog computing and the emulation of synaptic behavior, the devices are operated by applying current pulses to set or reset the resistive state, followed by a lowamplitude voltage pulse to read out the resistance.Running such pulse simulations is challenging and computationally expensive due to the broad range of time scales involved.Each individual pulse lasts just a few microseconds to milliseconds, while the entire sequence of pulses can extend over many seconds to minutes.In the following, we demonstrate that the model can reproduce the expected pulse-update behavior known from experiments.
For the simulations, several set and reset voltage pulse trains are applied with 500 pulse periods, respectively, using the device parameter set S 1 .Each pulse period consists of a set or reset pulse with voltage amplitudes of U set > 0 and U reset < 0 followed by a much shorter, low-amplitude read pulse, as illustrated in Figure 9a.All pulse parameters are summarized in Table 6 in the Methods.
A representative example of five simulated set pulses is shown in Figure 9b.It demonstrates that the current response follows the applied set and read voltage pulses as expected.With each set voltage pulse, the corresponding maximum set and read currents increase slightly, illustrated by the dashed line in Figure 9b in the example of the set current response.This change in the current magnitude becomes apparent over a larger number of pulses.The maximum read current I r is plotted in Figure 9c as a function of the pulse number N for three example pulse parameter sets (see Methods) with various pulse times  6) using the device parameter set S 1 .d) Example measurements taken from [23] performed on two representative lateral devices.
(t set = 2 − 10 ms, t reset = 1.5 − 3 ms) and rest times (t rest = 1 − 2 ms) over a total of 3000 pulse periods (≈ 45 − 60 s).For all three simulations, I r decreases during the reset pulses and increases again during the read pulses, showing analog memristive switching behavior consistent with the continuous I-V characteristics.The detailed behavior depends significantly on the pulse parameters, resulting in different linearity, maximum currents, and symmetry of the set and reset currents.Like the I-V curve simulations, the pulse simulations start from an ideal equilibrium configuration, which is not necessarily ensured in the measurements.Considering this, the simulation reflects the overall behavior of the experimental pulse-update curves of two representative devices from Ref. [23] very well (Figure 9d).

Discussion and Conclusion
We presented a semi-classical charge transport model, which selfconsistently couples the dynamic drift-diffusion equations for electrons, holes, and charged mobile point defects to Poisson's equation for the electrostatic potential.The formulation covers the specific physics of mobile point defects, particularly their nonlinear dynamics and density saturation prescribed by the material structure.Additionally, local self-heating is included via the transient heat equation, allowing us to analyze the mobilities' temperature and electric field dependency.Further, we developed a self-consistent formulation for image-charge-induced Schottky barrier lowering (SBL) to capture the effect of the applied voltage and charge carriers on the Schottky barriers.The system of equations is solved via a physics-preserving finite-volume discretization to analyze the resistive switching mechanism in twodimensional memristive devices.
The results indicate that the hysteresis in all types of I-V curves with and without significant Schottky barriers is caused by the dynamics of mobile charged vacancies.The vacancy dynamics leads to the formation and annihilation of a vacancy depletion region, which locally reduces the conductivity and eventually limits the current through the device.The dynamic response can be modified by many different parameters leading to qualitatively different device characteristics.This switching process explains the I-V curves and pulse characteristics observed experimentally [23] and many distinct features, such as the formation of a maximum at U < U max (i.e., a negative differential resistance), [16,38,101] the reduction of the maximum current magnitude with increasing mobility or decreasing sweep rate, hysteresis crossing, [101] and different hysteresis directions. [16,23,101]or symmetric contacts, the asymmetry of the I-V curve is volatile and reduces with increasing cycle number before it remains constant.Such an initial cycle-to-cycle drift can be avoided by operating the device in the high-vacancy-mobility range by choosing a corresponding voltage sweep rate (Figure 5f) at the expense of the retention time.Large vacancy mobilities reduce the retention time, which is typically desired for transient memristive devices, [102][103][104][105] but for most applications of memristive devices as nonvolatile memories, it is undesired [9,106] Optimizing the cycle-to-cycle drift and retention time implies selecting the narrow range of the sweep rate where the vacancies are Table 1.Overview of results for the activation energy of diffusion and the parameters extracted from the simulations to calculate the activation energy of diffusion from the vacancy mobilities obtained from the fits in Figure 4a and Figure 6a.sufficiently mobile to quickly occupy their dynamically stable configuration but are still sufficiently delayed, causing a large hysteresis.This is the case at the local maximum in Figure 5h, close to the point where the two hysteresis area curves split up.A permanent asymmetry in the I-V curve can be obtained by introducing Schottky barriers.We show how minuscule changes in the Schottky barriers (<0.04 eV) can qualitatively change the I-V characteristic, i.e., reverse the direction and change the hysteresis area and symmetry.Hence, the simulations quantify how deliberate changes can be introduced in the I-V curves for optimizing the device characteristics and explain the often large device-to-device variability. [25,107]Further, we confirmed that SBL contributes to the overall Schottky barrier heights and the I-V characteristics.More specifically, SBL reduces the differences between the intrinsic barriers and contributes significantly to the small Schottky barriers in equilibrium.However, while SBL notably affects the I-V curves, no system configuration was found where SBL dictates the hysteresis area.This result sheds new light on thefrequentlybrought-up hypothesis that SBL could be a major hysteresis mechanism in such devices. [23,48,49]ble 2. Summary of MoS 2 material parameters collected from the literature for a layer thickness of ≈15 nm and comparison with the values used in the simulations.The electron effective mass at rest is denoted by m 0 .Values for the fit parameters  n ,  p and ϕ 0 depend on the respective dataset and are provided in Table 4 and for E a in Table 1 [ 40,81,[93][94][95] 3. Summary of geometry parameters used to simulate the devices from. [23]The two values provided for the channel length L are used together with the parameter sets S 1 -S 6 (Table 4, Table 5).Hysteresis crossing, i.e., the change in thehysteresis direction, is an important feature well explained by the model.Different hysteresis directions have been previously observed experimentally [23,49,101] and interpreted as caused by different switching mechanisms. [23,49]However, the simulations clarify that the vacancy depletion mechanism can explain both hysteresis directions depending on small changes in parameters such as the Schottky barriers, the mobilities, and the number of sweep cycles.As we demonstrated, such a change in the hysteresis could also be induced at large voltages (20-50 V) by forming additional vacancies, consistent with the measurements in. [23]In all cases, the underlying origin of the hysteresis change is a slight modification of the vacancy's dynamic response and not an entirely different mechanism.
Drift-diffusion processes and the electronic properties are generally temperature-dependent.The comprehensive analysis of Joule heating illuminates the relevance of thermal effects on the I-V characteristics considering full coupling to the drift-diffusion equations of electrons and holes and the influence on vacancy mobility.Consistent with measurements, [108] the analysis shows that Joule heating can significantly influence the I-V curve of two-dimensional memristive devices.However, this effect depends strongly on the individual device properties and operating conditions.For all devices analyzed here, thermal effects are negligible, including the dependency of themobilities on the temperature and the electric field.
The excellent match of the simulations with various experimental I-V curves of lateral MoS 2 -based devices supports the validity of our results and model assumption.Furthermore, first pulse simulations reproduce the dynamic behavior of representative devices and show how the pulse update behavior can be tuned by altering the pulse parameters.
The device parameters obtained from the fits match very well within the expected range discussed in Section 3. Equilibrium Table 4. Sample-specific parameter sets S 1 and S 2 were obtained from the fit of the simulation to the experimental data in Figure 4a (S 1 ) and Figure 6a (S 2 ).For all fits, a sweep rate of f = 5V/s was used.8a) S 4 (Figure 8b) S 5 (Figure 8c) S 6 (Figure 8d) Units Left Schottky barrier ϕ 0 (x 1 ) 0.167 1 × 10 vacancy densities in the simulations are around n x ≈ 10 26 m −3 in good agreement with Auger-electron spectroscopy measurements on the modeled devices. [23]][95] A smaller range of values is obtained for the vacancy mobilities of  x ≈ 1 × 10 −14 − 1 × 10 −13 m 2 /(Vs) in the simulations.By using an ion hopping model [59] , these values translate to activation energies of diffusion of E a ≈ 0.50 − 0.53 eV (see Methods), which matches particularly well with previous estimations for the devices of E a = 0.6 eV. [23]Ab-initio studies result in larger activation energies of E a ≈ 0.8 − 2.5 eV, depending on the microscopic environment, such as the configuration of vacancy clusters and lines. [67]However, the activation energy of diffusion is expected to be further reduced along grain boundaries or interfaces. [38]ence, our results support the relevance of such defects for the device performance, consistent with experimental work on MoS 2based lateral devices. [16,38]he equilibrium Schottky barriers obtained from the simulations vary between ≈0 eV and 0.17 eV, close to the experimental values of ≈ 0.05 eV [98] and 0.18 eV [96,99] for Ti/MoS 2 contacts.The device-to-device variations of the Schottky barriers, mobilities, and vacancy densities can be caused by microstructure and interface quality variations, reflecting the fabrication tolerances. [109,110]ble 6.Summary of the pulse parameters used for the simulations in Figure 9c.Such variations can also occur in the geometry of the individual devices and explain minor deviations between measurements and simulations.Additionally, the charging and discharging dynamics of trap-states could contribute to the I-V characteristics depending on the defect and interface structure of the individual sample. [37,111]hile all simulations were performed with MoS 2 as a representative example of TMDC, the model is not limited to the drift and diffusion of positively charged sulfur vacancies in MoS 2. Instead, the model can be applied to other materials and kinds of mobile ionic defects with different charge states by adapting the material and device parameters.
In conclusion, our results corroborate the experimental indications for the relevance of sulfur vacancies in TMDC devices theoretically, and we identified the switching mechanism with its complex dynamics.Deep insights into the physics of the switching mechanism and the origin of hysteresis are gained as a vital step towards tailoring the device functionality for applications.The results bear general implications for the design and analysis of lateral memristive systems and show that the dynamics of mobile ionic defects can dominate their hysteresis.The model and simulation software are made available as powerful tools for the analysis and design of other semiconductor devices in the future.

Experimental Section
Finite volume model: We implemented all equations and numerical examples in the Julia programming language to extend the software package ChargeTransport.jl. [112]It was developed to simulate charge transport in semiconductors and is based on a finite volume method for spatial discretization. [113,114]For the simulations, we use a mesh with an element size of ≈20 nm in the center of the geometry, which refines linearly to 0.5 nm at the contacts.The refinement ensures an efficient solution while still representing the contact potentials accurately.The simulation software we developed is the open-source tool ChargeTransport.jl,version v.0.2.9 (https://github.com/PatricioFarrell/ChargeTransport.jl).
Image-charge dielectric constant: For estimating the image-charge dielectric constant ɛ i , the frequency dependency of the in-plane electric permittivity must be known.Measurements on CVD-grown monolayers [115] and calculations [116][117][118] show a strong dispersive behavior at energies of ≈E > 2 eV [119] and wavelengths of  < 900 nm with real parts of Re{ * r } = 8 − 20 (at  = 900 nm) depending on the study [115][116][117][118] but generally larger than the corresponding quasi-static values.Further, ab initio calculations of the complex permittivity demonstrate that the resonance peaks shift to even higher E and smaller  with an increasing number of MoS 2 layers until the bulk properties are reached. [119]Therefore, we assume an approximately constant electric permittivity  * r for energies E ≪ 2 eV, corresponding to photon frequencies of f ≪ E/h ≈ 5 × 10 14 Hz with values possibly slightly larger than the quasi-static value of ɛ r ≈ 10 r , and therefore, we assume ɛ i ≈ ɛ r .Maximum defect density: We estimate N x considering hexagonal 2H-MoS 2 with an in-plane lattice constant a = 0.316 nm and c = 1.229 nm, [74] which results in a unit cell volume of V = a 2 csin (60°) = 0.1063nm 3 ≈ 10 −28 m 3 .Each unit cell comprises two formula units of MoS 2 , i.e., two Mo atoms and four S atoms, and is part of two stacked monolayers with a respective distance of ≈ c/2.With V and the number of S-sites, we estimate a maximum sulfur vacancy density N x ≈ 4/V ≈ 4 × 10 28 m −3 .
Activation energy of diffusion and vacancy mobility: The activation energy of diffusion E a corresponding to the fitted values of the mobility  x is calculated using the mobility model of Genreith-Schriever. [59]The hopping distance a x is approximated by the in-plane lattice constant and set to a x = 0.316 nm, [74] for the attempt frequency, we use a value of v 0 = 10 12 Hz. [120]We omit the activation entropy of migration S a = 0, and consider the minimum electric field magnitude, i.e., E min = 0, as well as the maximum electric field magnitude of each time step E max averaged over time.With this range of the electric field, we estimate the error in the activation energy that occurs from assuming a constant vacancy mobility  x .The fraction of occupied vacancy sides is in good approximation n i ≔n x ∕N x ≈ 0. The values for E max and the resulting E a are summarized in Table 1 for the two fits from Figure 4a (data set S 1 ) and Figure 6a (data set S 2 ).
Summary of model parameters: In the following, all model parameters and references are summarized.We use a background doping concentration of C = 10 21 m −3 with z c = 1 for all simulations.This concentration value improves numerical stability and is sufficiently small to have no significant influence on the results.Further, a temperature of T = 300 K is used.The mobile ionic species x is assumed to behave as n-type dopants [71] with z x = +1 in their ionized state.All material parameters of MoS 2 obtained from the literature are summarized in Table 2 together with the respective references.If references are provided for a range of values, the value used from within this range is denoted separately.The geometry parameters are listed in Table 3, and the sample-specific parameters, i.e., the parameter sets S 1 and S 2 , in Table 4 and the parameter sets S 3 -S 6 in Table 5.
Pulse parameters.The followingtable summarizes the pulse parameters defined in Figure 9a, which are used for the simulations in Figure 9c.

Figure 1 .
Figure 1.Example of a lateral memristive device structure comprising two electrodes (source and drain) on top of a two-dimensional memristive material.The three illustrated mechanisms have been suggested as the origin of hysteresis in lateral memristive devices based on 2D TMDC materials: a) the charging and discharging dynamics of trap states,[12][13][14][15][16][17][18] b) the migration of charged point defects in the electric field,[23,38] and c) the voltage-induced phase change from a stable, low-conductive semiconducting phase to a metastable, high-conductive metallic phase.[42,44]

Figure 2 .
Figure2.a) Illustration of an electron point charge q n at the semiconductor-metal interface inducing an image charge -q n in the metal electrode, which results in an attractive electrostatic potential energy ϕ i (image-charge potential energy), and b) illustration of the superposition of the residual potential energy ϕ r and the image-charge potential energy ϕ i , which reduces the interfacial potential energy barrier maximum by ϕ relative to the value of ϕ r at the interface.The boundary value of ϕ r at the interface is prescribed by the sum of the intrinsic Schottky energy barrier ϕ 0 and the applied potential energy ϕ a (ϕ a = q n  a , with the applied electrostatic potential  a ).Note that an offset of ≈0.16 eV is added to ϕ i for illustrative purposes and ϕ i = 0eVfar away from the electrode, i.e., at x → −∞ in this example.

Figure 3 .
Figure 3. Illustration of a) the geometry of the modeled device in three dimensions with indicated SiO 2 /Si substrate, TMDC flake, and Ti/Au electrodes.b) 2D view of the x-z cut plane of the memristive device, and c) simplification of the geometry for the one-dimensional (1D) model with indicated channel length L, and locations x 1 and x 2 of the two contact boundaries .d) Illustration of the piecewise linear voltage U(t) as a function of the time t over a period of length t c , and with a maximum voltage of U max and a minimum voltage U min = −U max .

Figure 4 .
Figure 4. Simulations using the parameter set S 1 (Table 4, Methods).a) Comparison of measurements with simulations of the second I-V cycle.b) All eight consecutively simulated current-voltage cycles.c) Simulated reduced band diagram of the initial equilibrium configuration at t = 0s, and an applied voltage U = 0V, and d) of the nonequilibrium configuration at the beginning of the second cycle at t = 10.4s(U = 0V).e-j) Example configurations of the vacancy density n x (x) as a function of space x for various positions in the I-V curve shown in (a), together with the quasi-static equilibrium configuration of n x (gray), i.e., the configuration that would be reached at a constant applied voltage for t → ∞.The measured data are taken from[23]

Figure 5 .
Figure 5. Simulation results illustrating the dependency of the I-V characteristics on the vacancy mobility  x and the voltage sweep rate f using the parameter set S 1 (Table 4, Methods) a-f) Example hysteresis loops of the second cycle for various vacancy mobility values simulated with the voltage sweep rate f 1 = 5V/s.g) Negative cross-over voltage − U c normalized to the maximum voltage U max = 13V for the first twenty voltage cycles as functions of  x .The values from the second-cycle hysteresis curves in (a-e) are indicated with dots.h) Hysteresis areas A H defined as the integrals over the right (bold) and left (dashed) hysteresis branch of the second sweep cycle for two different voltage sweep rates f 1 = 5V/s and f 2 = 2.5V/s.The values from the hysteresis curves in (a-f) are indicated with dots.

Figure
Figure 5a-f are indicated.As expected, A H = 0 for  x → 0 and  x → ∞ because the vacancy density either remains stationary or in quasi-static equilibrium.A significant hysteresis occurs only in the regime between approximately  x = 10 −15 m 2 /(Vs) and  x = 10 −12 m 2 /(Vs).As the mobility exceeds small values  x ≤ 1 × 10 −15 m 2 /(Vs), the hysteresis starts to open (A H ≈ 0) and is initially clockwise oriented in both branches (Figure5a-b).Because the depletion zone is still present when U < 0 V is reached the device is in the low conducting state at the beginning of the left hysteresis branch.A transition of the hysteresis direction occurs in the left hysteresis branch from clockwise (A H < 0) to counterclockwise (A H > 0) at intermediate mobilities around  x ≈ 3 × 10 −14 m 2 /(Vs) (Figure5c).During the transition, the current in the left hysteresis branch intersects at a voltage U c (Figure5c, indicated), which continuously moves towards U = 0 V as the mobility increases until the entire hysteresis branch is directed counterclockwise (Figure5d).Origin of this transition is the shorter time required for vacancies with a larger mobility to reach quasi-static equilibrium, which causes a faster annihilation of the depletion zone before U < 0 V is reached.Further increasing the vacancy mobility to large values (Figure5d-f) results in the formation of distinct maxima at voltages U < U max (and U > U min ), i.e., a negative differential resistance, an increasingly symmetric I-V curve, and an overall reduced maximum current.This behavior reflects the formation dynamics of the depletion zone.In the high-mobility regime, the depletion zone can form before the maximum positive or negative voltage is reached, which results in a drop in the current, leading to the formation of the maxima.Simultaneously,

Figure 6 .
Figure 6.Comparison of simulations with measurements of an asymmetric hysteresis curve using the parameter set S 2 (Table 4, Methods).a) Simulation of the second voltage cycle and measurements.b) Conduction band edges and quasi-Fermi levels in equilibrium at t = 0s with and without SBL.c) Comparison of simulations of the second voltage cycle with and without SBL and without SBL but with intrinsic Schottky barriers equal to the equilibrium Schottky barriers from the case with SBL (No SBL, reduced ϕ 0 ).d) Magnitude of the current, applied voltage U(t), and change of the Schottky barrier ϕ(t) as functions of the time t.The measured data are taken from.[23]

Figure 7 .
Figure 7. Simulations of the second cycle using parameter set S 2 (Table 4, Methods) changing only the intrinsic Schottky barriers.a) Current-voltage (I-V) characteristics for three different values of the left barrier ϕ 0 (x 1 ) and a fixed value of ϕ 0 (x 2 ) = 0.11eV for the right barrier.b) I-V characteristics for three different values of the right barrier ϕ 0 (x 2 ) and a fixed value of the left barrier ϕ 0 (x 1 ) = 0.11eV.c) I-V characteristics of the first 10 voltage cycles and the 100 th cycle, with indicated Schottky barriers and vacancy mobility.d) Total hysteresis area as a function of the cycle number normalized to its respective value at cycle 100 and for three different vacancy mobilities.The parameter set and the Schottky barriers are identical to those used for the calculations in (c).

Figure 8 .
Figure 8.Comparison of simulations with measured data of I-V curves of four devices a-d) with various maximum voltages and extends of asymmetry (data sets S 3 -S 6).The measured data are taken from,[23] and the simulation parameters are provided in the Methods.

Figure 9 .
Figure 9. Simulations of pulse resistance update measurements and qualitative comparison with measurements.a) Illustration of the set and reset pulse shape with subsequent read pulse and definition of the pulse parameters.While the set pulse amplitude U set > 0V, it is U reset < 0V.b) Examples of the simulated voltage function U(t) and the resulting current I(t) in the time domain.c) Three example simulations with various pulse parameters are summarized in the Methods (Table6) using the device parameter set S 1 .d) Example measurements taken from[23] performed on two representative lateral devices.

Table 5 .
Parameters obtained from the model's fit to the experimental data in Figure8a-d.All fits were performed with a sweep rate of f = 5V/s.