High-throughput numerical modeling of the tunable synaptic behavior in 2D MoS2 memristive devices

Memristive devices are emerging as neuromorphic candidates for next-generation computing systems because of their compact size, high switching speed, multi-bit memory states, and energy-efficient analog computing operation^1,2,3. By exploiting pinched hysteresis loops in current-voltage ((I)–(V)) characteristics, these devices can store thousands of different resistive states that can be updated and read by applying simple voltage or current pulses^4,5,6. Recently, two-dimensional (2D) materials have expanded the library of memristive systems, enabling a variety of additional features such as ultimate scalability down to atomic dimensions, atomically sharp interfaces, anomalous memory in moiré heterostructures, tunable electronic structure, reduced electrostatic screening, and strong electric field effects^{7,8,9,10,11,12}. Moreover, 2D materials have been exploited for multi-terminal memristive devices^{13,14,15,16,17,18,19,20,21,22}, typically in lateral geometries with additional gate electrodes that permit tuning of the device characteristics to achieve desirable performance features such as high linearity and symmetry in the synaptic response, tunable synaptic learning rate, heterosynaptic plasticity, and an increased number of resistive states^{13,16,23,24,25}. These multi-terminal memristive devices are also referred to as memtransistors because they unify nonvolatile memristive properties with the volatile gate tunability of field-effect transistors¹⁷. While this device technology is still in its infancy, recent work has achieved rapid improvements in the performance metrics of memtransistors based on 2D transition metal dichalcogenides (TMDCs)^26,27. For example, memtransistors have been scaled to sub-micron channel lengths to achieve low switching energies down to 20 fJ/bit, low switching voltages below 1 V, and minimal sneak path currents of 0.1 nA in crossbar architectures^14,23,28,29.

Despite these early empirical achievements, further advances in memtransistor performance and integration can be expected by developing and exploiting a thorough understanding of the physical origins of memristive switching mechanisms⁵. Specifically, optimizing 2D memristive devices relies on understanding and controlling pulse-controlled synaptic responses through device design parameters including channel electronic properties, stoichiometry, defect density, grain boundaries, contact materials, and electrode geometries^5,30. A critical challenge is to develop a model that concurrently captures all nuances of charge transport, electrostatic screening, and memristive switching in order to establish correlations between device design parameters and performance. For example, most 2D memristors and memtransistors based on TMDCs (e.g., MoS₂) as a channel material^{27,31,32,33,34} are believed to involve migration of mobile defects (e.g., sulfur vacancies) as the origin of resistive switching in these devices^13,35,36,37. Since the field effects from the gate electrode further complicate the switching mechanism in memtransistors, the first step is to develop a physical model of resistive switching in lateral 2D memristors without the gate electrode. Indeed, previous theoretical literature reproduced (I)–(V) characteristics of such 2D memristive systems by the formation and annihilation dynamics of a local vacancy depletion zone in the conduction channel^22,38. However, a larger number of material and device parameters complicate the physical model of the charge transport, resulting in a lack of insight into the dominant drift-diffusion processes. This incomplete understanding further hinders the establishment of correlations between materials parameters and characteristics of 2D memristive systems as synapses in artificial neural networks (ANNs). For example, hardware implementation of ANNs with backpropagation benefits from linear and symmetric weight-update rules^{2,39,40,41,42,43}. Several schemes have been proposed to achieve linear and symmetric responses, such as pulse shaping and filament tuning in conventional oxide memristors^{2,39,40,41,42,43}. However, no predictive model directly connects materials parameters to the tunable synaptic learning behavior in 2D memristive systems.

Computational, compact, and Monte Carlo models present a first step toward providing the physical insight to build robust empirical correlations between materials parameters and device performance. To date, relatively few computational models have been presented to study memristive switching in memtransistors and other 2D memristive devices^22,44,45,46. Compact models use an equivalent circuit formulation comprising simplified physics-inspired descriptions of the device to reproduce its (I)–(V) characteristics⁴⁴. These models are efficient for integrating the device characteristics on a circuit/system level but require many fitting parameters and often do not directly reflect the microscopic physics of the switching mechanism. In contrast, kinetic Monte Carlo methods^45,46 consider the stochasticity of defect migration over nanoscale volumes, but fall short of directly incorporating materials parameters and semiconductor-metal interfaces. Recently, quasi-static²² and fully transient³⁸ mesoscale models have been presented to describe the switching processes in 2D memristive devices based on TMDCs. These models solve systems of partial differential equations for the charge and mass transport in the entire device, thus providing some connections between microscopic switching mechanisms and device performance^47,48. Given the complexity of the described physics in the spatial and temporal scales, these simulations can be computationally demanding, which makes detailed analyses of large parameter spaces prohibitively time-consuming⁴⁸.

Here, we propose device metrics and combine them with high-throughput charge-transport simulations to screen the vast material and device parameter space for 2D memristive systems. Specifically, we employ a transient and fully coupled semiconductor charge-transport model, which considers electron and hole dynamics in addition to the drift and diffusion of point defects that are coupled through Poisson’s equation for the electrostatic potential. High-throughput simulations are performed for 2D memristive systems over a wide range of parameters to establish a robust data set of pulsed learning behavior that serves as the basis for synaptic weight-update protocols in practical neuromorphic circuits. Synaptic learning metrics, such as linearity and symmetry, are extracted from the simulations to identify correlations with the design and operating parameters. In this manner, this model can not only be used as a tool for designing 2D memristive synaptic elements but also serve as a basis for comprehensive memtransistor models involving multiple terminals.

Quasi-2D MoS₂ is used as a representative material to investigate resistive switching processes through vacancy migration in lateral memristive devices. A charge-transport model is used to describe the resistive switching process by solving the fully coupled transient drift-diffusion equations of electrons, holes, and positively charged vacancies together with Poisson’s equations for the electrostatic potential (psi). The equations consider the nonlinear dynamics of the charge carriers and the different physical nature of the vacancies compared to the free electrons and holes (e.g., their limited maximum density prescribed by the material structure and distinct nonlinear dynamics). The highly anisotropic mobilities of the charge carriers are simulated by confining them in the plane and allowing carrier migration only along the channel direction. Figure 1a illustrates the geometry of a 2D lateral memristor comprising the semiconducting conduction channel and two metal-semiconductor interfaces at the metal contact electrodes. We assume Schottky barriers at the contacts and consider their dependency on the local electric field via image-charge-induced Schottky barrier lowering (SBL). While the left contact is grounded (source), a voltage (V) can be applied to the right contact (drain). In contrast to other models, we self-consistently include SBL by solving for an additional auxiliary potential together with the fully coupled time-dependent semiconductor equations and the drift-diffusion of vacancies. The model equations consider both drift and nonlinear diffusion using the quasi-Fermi potentials as unknowns, which has several numerical advantages⁴⁹. The model is implemented in the Julia programming language utilizing the VoronoiFVM solver with a finite volume discretization for the spatial domain^50,51,52. The finite volume method is particularly suitable for such conservative flow problems because it locally conserves the flux and is constrained not to produce unphysical negative densities⁴⁹. The specialized code allows for efficient parallelization of the simulations leading to a high throughput of approximately 6550 simulations per day, distributed among three central processing units. Further details on the model and simulations are provided in the “Methods” section. All simulations start by finding the initial equilibrium configuration of all charge carriers and the electrostatic potential for (V=0,{rm{V}}) before applying a time-dependent voltage impetus.

a An illustration of a lateral memristor geometry and its cross-section along the conduction channel. b Comparison of experimental current–voltage (I–V) curves from ref. ³⁵ with simulated curves obtained by applying a triangular voltage function shown in the inset. c Logarithmic representation of the current magnitude ({|I|}) of the data shown in (b). d Example results of the electrostatic potential (psi), the vacancy density ({n}_{{rm{x}}}), and the electron density ({n}_{{rm{n}}}), at three different points (I–III) in the current-voltage characteristic indicated in (b) and (c).

Full size image

Current–voltage and pulse characteristics

We first apply triangular voltage sweeps (Fig. 1b, inset) with a maximum voltage ({V}_{max }=10,{rm{V}}) and a voltage sweep rate ({rm{d}}V/{rm{d}}t=5,{rm{V}}/{rm{s}}) to simulate the I–V characteristics. Example results of the second cycle are shown in Fig. 1b (linear scale) and Fig. 1c (logarithmic scale), using the material parameters of a 15-nm-thick MoS₂ layer. A pinched hysteresis loop is visible in the (I)–(V) curves, resulting from non-filamentary bipolar resistive switching. The switching is caused by the formation and annihilation dynamics of a local vacancy depletion zone (Fig. 1d). The vacancies migrate homogeneously along the conduction channel, driven by the applied voltage. Because the vacancies are positively charged, their local depletion is accompanied by a locally reduced electron density in the channel, which decreases the local conductivity and the net current. The (I)–(V) simulations match well with measurements on MoS₂-based lateral memristors without the gate electrode in addition to previous calculations^35,38, demonstrating that the model captures the essential mechanism of resistive switching in these devices. The image-charge-induced Schottky barrier lowering does not contribute significantly to the hysteresis because it is largely dominated at the contacts by the completely reversible electric field from the applied voltage and not by the shift in the charge carrier densities³⁸.

To further validate the memristor model, we compare experimental measurements of the synaptic behavior from ref. ³⁵ with simulations from the present charge-transport model. Specifically, multi-bit memory in the reset and set process is used to realize long-term potentiation (LTP) and long-term depression (LTD) in synaptic devices, respectively, that can be used to mimic ANNs in crossbar architectures^2,3,6. The non-filamentary soft switching in lateral devices is particularly well-suited for multi-bit memory operations rather than filamentary hard switching in vertical devices that produce digital-like on and off states. For the simulations, we apply a sequence of 1000 voltage pulses to set the device into a defined resistance state (set pulses) followed by 1000 voltage pulses to reset the resistance to its initial state (reset pulses). After each set and reset pulse, a short read pulse is applied to extract the read current ({I}_{{rm{r}}}) and the resistance state. Hence, each pulse period comprises a set or reset pulse and a read pulse as illustrated in Fig. 2a. For the simulations, we select geometry and pulse voltages matching those in the experiments, and the same material parameters are used both in (I)–(V) curve simulations (Fig. 1b) and synaptic behavior simulations (Fig. 2b). Only the charge carrier mobilities are slightly modified to accommodate for expected device-to-device variability³⁸ and the pulse durations to maintain numerical stability. The results shown in Fig. 2b qualitatively match the experimental nonlinear pulse update behavior, suggesting that the model captures the essential nonlinear dynamics of the memristive system.

a Illustrations of a single set and reset pulse period used in simulations in b providing definitions of the pulse parameters for the case of pulse periods with negative set pulses SET(-­) and positive reset pulses RES(+). b Comparative plot of experimental long-term potentiation (reset) and long-term depression (set) curves from ref. ³⁵ and simulation over four sequences of set and reset pulse trains with 1000 pulse periods each.

Full size image

We expect these essential aspects of the model to remain valid for longer and shorter pulse durations as long as the memristive mechanism is dominated by vacancy dynamics. However, for very long pulse durations (and/or large voltages), additional processes may become relevant that cause permanent changes in the materials and are not included in the model, for example, Joule heating and electromigration at the semiconductor–metal interfaces^53,54,55. Based on previous calculations, such thermal effects are omitted in the operation regime investigated here³⁸. Electromigration and other effects causing permanent changes to the device were ruled out in the experimental study of the devices for small operation voltages and pulse times³⁵.

High-throughput pulse simulations

A large number of pulse simulations are performed to develop empirical relationships between material properties and synaptic plasticity. First, we analyze the pulse-dependent read current of the system for applications in neuromorphic computing in an idealized case of equal Schottky barrier heights at the two contact electrodes.

Definition of pulse schemes and polarities

During a single pulse simulation, we apply an alternating set and reset pulse sequences with 50 pulse periods each, leading to a total of 10 pulse sequences and 500 pulse periods per simulation (see Fig. 3a–c). We refer to all pulse sequences with an odd number as set pulse sequences, i.e., the first, third, fifth, seventh, and ninths, and all pulse sequences with an even number as reset pulse sequences. For the simulations, we select set and reset pulses with variable maximum voltage magnitudes ({V}_{{rm{set}}}) and ({V}_{{rm{res}}}) (between −5 and 5 V), respectively, and equal time durations of ({t}_{{rm{set}}}={t}_{{rm{res}}}=5,{rm{ms}}). Each set and reset pulse is followed by a short read pulse with a voltage of ({V}_{{rm{read}}}=0.2,{rm{V}}) and duration ({t}_{{rm{read}}}=0.14,{rm{ms}}) to extract the read current ({I}_{{rm{r}}}). The time duration between the set/reset pulse and the read pulse is ({t}_{{rm{rest}}}=1,{rm{ms}}).

a Symmetric rectangular pulse scheme (scheme 1) with an alternating series of set and reset pulse sequences using constant pulse voltages. b Illustration of the asymmetric pulse scheme (scheme 2) with a constant set voltage, but linearly increasing reset voltage. c Symmetric triangular pulse scheme (scheme 3) with linearly increasing set and reset pulse voltages. All pulse schemes are illustrated in the example of pulse case II over the first six sequences and define the set pulse sequence as the pulse sequences with odd numbers, i.e., the first, third, fifth sequence. d Definition of the pulse case I, comprising positive set pulses SET(+) and positive reset pulses RES(+). e Definition of the pulse case II, with negative set pulses SET(-­) and positive reset pulses RES(+). f Definition of pulse case III, with negative set and reset pulses, SET(-­) and SET(+). g Definition of pulse case IV with positive set pulses SET(+) and negative reset pulses RES(-­). All pulse simulations use five alternating set and reset pulse trains with 50 pulse periods each, resulting in a total of 500 pulse periods.

Full size image

Since we vary the magnitude and the sign of the voltage pulses, we refer to set pulse sequences with ({V}_{{rm{set}}} > 0,{rm{V}}) as SET(+) and those with ({V}_{{rm{set}}} < 0,{rm{V}}) as SET(-­). Correspondingly, reset pulse sequences with ({V}_{{rm{res}}} > 0,{rm{V}}) are indicated as RES(+) and those with ({V}_{{rm{res}}} < 0,{rm{V}}) as RES(-­). This variation in the pulse amplitude leads to four different cases (cases I–IV) with different combinations of positive and negative pulse sequences within the two-dimensional parameter space spanned by ({V}_{{rm{set}}}) and ({V}_{{rm{res}}}) (Fig. 3d–g).

Each pulse case is analyzed for three types of pulse schemes, illustrated in Fig. 3a–c (schemes 1–3) for case II pulses. These pulse schemes differ in how the pulse amplitudes change within each pulse sequence. For the symmetric rectangular pulse scheme (scheme 1), ({V}_{{rm{set}}}) and ({V}_{{rm{res}}}) are kept constant (Fig. 3a). For the asymmetric pulse scheme (scheme 2), ({V}_{{rm{set}}}) is kept constant, while the reset voltage amplitude is linearly ramped up from 0 V at the first pulse period to the final value of ({V}_{{rm{res}}}) at the last pulse period (Fig. 3b). In the symmetric triangular pulse scheme (scheme 3), both pulse voltages are linearly ramped up from 0 V to their final values of ({V}_{{rm{set}}}) and ({V}_{{rm{res}}}), resulting in an alternating triangular-shaped pulse train (Fig. 3c). The three pulse schemes in Fig. 3a–c are illustrated for ({V}_{{rm{set}}} < 0,{rm{V}}) and ({V}_{{rm{res}}} > 0{,rm{V}}) (Case II). In particular, the second pulse scheme is included as an intermediate scheme between the fully rectangular and the fully triangular pulse schemes to analyze the influence of individual triangular or rectangular pulse sequences on the read current. We perform 2601 pulse simulations for each pulse scheme by varying ({V}_{{rm{set}}}) and ({V}_{{rm{res}}}) between (-5) and (5) V with steps of 0.2 V to capture all four pulse cases, resulting in a total of 7803 pulse curves.

Influence of pulse polarities and vacancy dynamics

Figure 4a shows two example pulse curves that were simulated with the first pulse scheme the pulse case II and IV, respectively. For the first simulation, ({V}_{{rm{set}}}=-5) V, ({V}_{{rm{res}}}=4,{rm{V}}) (case II) and for the second simulation, the polarity of the voltages was reversed (i.e., ({V}_{{rm{set}}}=5) V, ({V}_{{rm{res}}}=-4,{rm{V}}), case IV). Note that the read current of both curves is approximately identical. It decreases during the set process and increases during the reset process. Over the five sequences of set and reset processes, this behavior repeats periodically. The origin for the independence of the read current on the voltage polarity can be explained by the vacancy dynamics in the system.

a Read current ({I}_{{rm{r}}}) as a function of the pulse number (N) for two example simulations using the pulse scheme 1, and opposite pulse polarities, namely case II with SET(-) and RES(+) pulses and case IV with SET(+) and RES(-); ({I}_{{rm{r}}}) at (N=25) and (N=50) is indicated. b Vacancy densities ({n}_{{rm{x}}}) as a function of space from the case II simulation in (a) in the initial equilibrium ((N=0)) and at (N=25) and (N=50), and c the same for the case IV simulation in (a).

Full size image

Figure 4b, c shows the corresponding vacancy density ({n}_{{rm{x}}}) along the center of the conduction channels in equilibrium and after (N=25) and (N=50) pulses in the first set process. In equilibrium (for both simulations), the vacancy density (({n}_{{rm{x}}})) is approximately constant at ({n}_{{rm{x}}}approx 6,cdot, {10}^{23},{{rm{m}}}^{-3}) in the middle of the channel and is depleted in the narrow space charge region around 20 nm to the electrode contacts resulting in ({n}_{{rm{x}}}approx 5,cdot, {10}^{20},{{rm{m}}}^{-3}) at the interfaces. Because the two Schottky barriers are identical in equilibrium, the vacancy density is symmetric around the center of the channel. During the set process and with increasing pulse number, the vacancy density depletes at one end of the channel (and increases at the other end) reaching a minimum of ({n}_{{rm{x}}}approx 8,cdot, {10}^{16},{{rm{m}}}^{-3}) at (N=50). This process results in the reduction of the conductivity consistent with the hysteresis in the I–V curves in Fig. 1b, c (I–III). Hence, under these idealized conditions, pulses of both polarities induce an equal change in the read current owing to the symmetry of drift-diffusion of mobile vacancies with equal Schottky contacts. Depending on the sign of the pulse voltages, the depletion zone forms either at the left or right side of the conduction channel and leads to an almost identical reduction of conductivity.

To analyze the memristive process over several set and reset sequences, we plot the read current in Fig. 5a, together with the vacancy densities ({n}_{{rm{x}},{rm{L}}}) at the left contact and ({n}_{{rm{x}},{rm{R}}}) at the right contact in Fig. 5b as functions of the pulse number for the case IV (V_set = 5 V, V_res = −4 V). Consistent with Fig. 4c, the vacancy density ({n}_{{rm{x}},{rm{R}}}) decreases during the set process by several orders of magnitude, while ({n}_{{rm{x}},{rm{L}}}) increases only slightly. Because ({n}_{{rm{x}},{rm{R}}}ll {n}_{{rm{x}},{rm{L}}}) for most of this process, the read current is mainly limited by the vacancy depletion at the right contact. These changes in the vacancy density are restored during the following reset process where ({n}_{{rm{x}},{rm{L}}}) decreases and ({n}_{{rm{x}},{rm{R}}}) increases. Halfway through the reset process ({n}_{{rm{x}},{rm{L}}}le {n}_{{rm{x}},{rm{R}}}) (i.e., a shallow depletion zone forms at the left contact), which limits the further increase of the read current. This set/reset behavior repeats periodically over the following pulse sequences. As a result, the read current in this particular case is dominated by the formation and subsequent annihilation of the vacancy depletion at the right contact, while the shallow depletion at the left contact only has a small effect on the slope of the read current.

a Read current ({I}_{{rm{r}}}) simulated with V_set = 5 V, V_res = −4 V as a function of the pulse number, and b corresponding vacancy density at the left and right contact. c Read current simulated with V_set = 1 V, V_res = −1.5 V as a function of the pulse number, and d corresponding vacancy density at the left and right contacts.

Full size image

Even with constant memristor model parameters, the vacancy dynamics strongly depends on biasing conditions, resulting in different synaptic plasticity characteristics. Depending on the pulse voltages, not only can the incremental change in the read current be reversed (LTP versus LTD) but different pulse trains can also produce complex and irregular features. For example, Fig. 5c shows reversed set and reset processes obtained using the same pulse scheme (scheme 1, case IV) as in Fig. 5a, b but smaller pulse magnitudes (({V}_{{rm{set}}}=1) V, and ({V}_{{rm{res}}}=-1.5,{rm{V}})). The read current changes aperiodically initially and only becomes periodic over the last six pulse sequences. The final periodic changes are in opposite directions compared to the simulations in Fig. 5a (i.e., the read current increases during the set process and decreases during the reset process).

To understand this behavior, we correlate the read current in Fig. 5c with the vacancy density at the left and right contact in Fig. 5d. The vacancy density behaves similarly to Fig. 5b (i.e., decreasing at the right contact and increasing at the left contact during the set process and vice versa during the reset process). However, in contrast to the simulation in Fig. 5a, b, the vacancy density oscillates at both contacts over only one order of magnitude (shallow depletion). Additionally, the periodic oscillation of the vacancy density is superposed by a drift of ({n}_{{rm{x}},{rm{R}}}) and ({n}_{{rm{x}},{rm{L}}}) to higher and lower values, respectively. Hence, ({n}_{{rm{x}},{rm{R}}}) and ({n}_{{rm{x}},{rm{L}}}) diverge over the series of pulse sequences until the minimum of ({n}_{{rm{x}},{rm{R}}}) is close to the maximum of ({n}_{{rm{x}},{rm{L}}}) in the last set sequence. This drift is the consequence of the larger pulse voltage ratio (left|{V}_{{rm{res}}}/{V}_{{rm{set}}}right|=1.5) compared to (left|{V}_{{rm{res}}}/{V}_{{rm{set}}}right|=0.8) for the simulation in Fig. 5a, b. Because of the large voltage ratio, the change of the vacancy density created during one pulse sequence is barely compensated during the following sequence, leading to the observed drift. This drift, together with the shallow and symmetric oscillation of ({n}_{{rm{x}},{rm{R}}}) and ({n}_{{rm{x}},{rm{L}}}) can explain the observed read current behavior.

Because of the symmetric change of the vacancy density, the read current is largely dominated by vacancy depletion during the first set process (({n}_{{rm{x}},{rm{R}}} < {n}_{{rm{x}},{rm{L}}})) and the following reset process (({n}_{{rm{x}},{rm{L}}} < {n}_{{rm{x}},{rm{R}}})). This changes gradually over the following pulse sequences. Owing to the shallow oscillation of the vacancy density, the drift has a significant influence on the ratio of ({n}_{{rm{x}},{rm{L}}}) to ({n}_{{rm{x}},{rm{R}}}). As a consequence, the set process becomes increasingly dominated by the annihilation of the depletion at the left contact. The net effect is an irregular read current during the first few pulse sequences until the set process is entirely dominated by the annihilation of the depletion zone at the left contact. Hence, compared to the simulation in Fig. 5a, b, the read current is eventually dominated by ({n}_{{rm{x}},{rm{L}}}), which causes the reversed potentiation-depression behavior. This example shows that the formation and annihilation mechanism can lead to a complex non-trivial dynamic behavior even if only the bias conditions are modified. Consequently, comparing different pulse schemes requires screening a larger parameter space instead of analyzing individual pulse curves.

Pulse metrics

To quantify the relation between materials parameters and operating conditions with the synaptic behavior of 2D memristive devices, we define seven different metrics to characterize pulsed LTP and LTD curves and study their dependence on different operating conditions. These metrics are: (1) drift (D), (2) symmetry (S), (3) average current change (Delta {I}_{{rm{r}}}) per set pulse, (4,5) linearity of the set and reset process ({L}_{{rm{set}}}) and ({L}_{{rm{res}}}), and (6,7) current pulse overshoot ({P}_{{rm{set}}}) and ({P}_{{rm{res}}}) to characterize the operation regime.

Drift metric D

First, we define the drift metric (D) as the slope obtained from a linear least mean squared fit to the pulse curve, normalized to the mean value of ({I}_{{rm{r}}}) to reduce the dependency of the drift on the absolute value of ({I}_{{rm{r}}}) (see Fig. 6a). The magnitude ({|D|}) of the extracted drift is plotted in Fig. 6b–d for all three pulse schemes (scheme 1–3) as a function of ({V}_{{rm{set}}}) and ({V}_{{rm{res}}}). Qualitatively, the behavior of (D) and ({|D|}) for three pulse schemes is very similar. Over the entire parameter space, the drift is negative, and its magnitude varies over several orders from ({|D|}={10}^{-4}) to ({|D|}approx 1). The largest values in the range of ({10}^{-2}) and 1 are obtained in the first and third quadrants of the plots, i.e., in regions where ({V}_{{rm{set}}}) and ({V}_{{rm{res}}}) have identical signs. This behavior is intuitively obvious because for identical signs, both (set and reset sequences) drive the vacancy density in the same direction and thus a reset of the device does not occur. The smallest values of ({|D|}) approximately between ({10}^{-4}) and ({10}^{-3}) are achieved in a narrow parameter band (a “low-drift valley”) that extends from the second quadrant through the middle of the plot down to the fourth quadrant (Fig. 6, white dashed lines). In this low-drift valley, the effect of the set pulse sequence on ({I}_{{rm{r}}}) is mostly compensated by the subsequent reset sequence with a voltage of opposite sign, and the positive read pulse leading to a splitting into two sub bands. Because of this compensation, the low-drift valley lies approximately on the diagonal where ({V}_{{rm{set}}}approx {V}_{{rm{res}}}) for the two symmetric pulse schemes (scheme 1 and 3). A slight deviation from this diagonal is apparent as ({V}_{{rm{set}}}to -5,{rm{V}}) and ({V}_{{rm{res}}}to +5,{rm{V}}) (i.e., the low-drift valley does not end in the corner of the second quadrant but at a slightly larger ({V}_{{rm{set}}})). This asymmetry occurs because a positive pulse voltage does not have the same effect on the vacancy density as a negative one owing to the hysteresis and nonlinearity in the system, and the minor influence of the read voltage pulse. Hence, it is relevant whether the first pulse sequence has a positive or negative pulse voltage. In the asymmetric pulse scheme (scheme 2), the low-drift valley is rotated further clockwise compared to schemes 1 and 3. Hence, in scheme 2, for a set pulse sequence with a given (|{V}_{{rm{set}}}|), a corresponding reset sequence with a larger (|{V}_{{rm{res}}}|) and opposite sign is required to compensate the effect of the set process in ({I}_{{rm{r}}}) to obtain low drift. This overall smaller influence of the reset sequence on ({I}_{{rm{r}}}) results from its triangular character, where all pulses (except the last one) have a voltage (< |{V}_{{rm{res}}}|). Consequently, increasing (|{V}_{{rm{set}}}|) and/or (|{V}_{{rm{res}}}|) of a triangular pulse sequence (as in scheme 2 and 3) has a smaller influence on the drift than a corresponding increase in the rectangular pulse scheme (scheme 1). Therefore, the low-drift valley is narrowest in scheme 1, slightly wider in scheme 2, and significantly wider in scheme 3.

a Illustration of the drift pulse metric, and b plot of drift as a function of the set and reset voltage for scheme 1, c for scheme 2, and d for scheme 3. e Illustration of symmetry metric, and f plot of symmetry as a function of set and reset voltage for scheme 1, g for scheme 2, and h for scheme 3. i Illustration of the read current increment per pulse (Delta {I}_{{rm{r}}}) within the set pulse train, and j plot of (Delta {I}_{{rm{r}}}) for scheme 1, k for scheme 2, and l for scheme 3. The red dots in f–h indicate the maximum symmetry values, and the arrows point in the direction of increasing symmetry.

Full size image

Symmetry metric S

As a second metric, we define the symmetry (S) via the root mean square deviation of the normalized read current of the last set process to the reversed read current of the last reset process, as illustrated in Fig. 6e. This metric is defined within (0le {S}le 1), where (S=1) indicates perfect symmetry. Hence, for (S=1), the read current in the reset process returns the same way as it was changed during the set process. The results of the symmetry evaluation for the three pulse schemes (scheme 1–3) are shown in Fig. 6f–h. In contrast to the drift, significant qualitative differences are visible in (S) for the three pulse schemes. For all three schemes, the largest symmetry is found within the low-drift valley (Fig. 6f–h), which is conceivable because drift always introduces asymmetry to some extent. Therefore, low symmetry tends to occur in the 1st and 3rd quadrants. This behavior is particularly apparent for the symmetric rectangular pulse scheme (scheme 1), where two high-symmetry regions with maximum values of (S > ,0.9) in the 2nd and 4th quadrants (cases II and IV) are separated from each other by a low-symmetry diagonal band with minimum values of (Sapprox 0.35) that extends from the 3rd quadrant to the 1st quadrant. Along the low-drift valley, the symmetry increases with increasing magnitudes of the pulse voltages (Fig. 6f, arrows) to a maximum of (S=0.943), at ({V}_{{rm{set}}}=-4,{rm{V}}) and ({V}_{{rm{reset}}}=3.5,{rm{V}}).

In the symmetric triangular pulse scheme (scheme 3), the two high-symmetry regions are rotated clockwise relative to those of scheme 1 with an overall slightly lower symmetry of mostly around (Sapprox 0.8). The low-symmetry band is concentrated around the origin and reaches minimum symmetries down to (Sapprox 0.2). In contrast to the rectangular pulse scheme (scheme 1), the symmetry is maximum at small voltages ((S=0.93) at ({V}_{{rm{set}}}=-0.5,{rm{V}}) and ({V}_{{rm{res}}}=0.5,{rm{V}})) and decreases along the low-drift valley with increasing voltages (Fig. 6h, arrows). In contrast to schemes 1 and 3, the asymmetric pulse scheme (scheme 2) offers two additional regions with large symmetry mostly around S = 0.84 (left and right side of the plot, Fig. 6g). Compared to the symmetric pulse schemes (schemes 1 and 3), the second pulse scheme also has significantly larger values in the low-symmetry ((S > 0.6)) band and a larger minimum value of (S > 0.4). The largest symmetry is obtained in the 4th quadrant at the left edge of the low-drift valley with a maximum of (S=0.97) at ({V}_{{rm{set}}}=1,{rm{V}}) and ({V}_{{rm{res}}}=-2.5,{rm{V}}) (Fig. 6g, red dot).

In all three schemes, the transition between low and high symmetry occurs mostly smoothly. However, in schemes 2 and 3, abrupt changes from low to high symmetry are visible within the low-drift valley. It will be shown later that this phenomenon can be attributed to qualitative changes in the character of the pulse curves. Overall, the largest region of high symmetry is observed in the rectangular symmetric pulse scheme at large voltages (scheme 1, 2nd quadrant), and the lowest symmetry in the symmetric triangular scheme at small voltages (scheme 3, around the center). Using an asymmetric pulse scheme (scheme 2) extends the parameter range where symmetry is high while increasing the minimum symmetry and offering a narrow band of extremely high symmetry close to (S=1) at intermediate voltage values.

Read current increment metric (Delta {{boldsymbol{I}}}_{{bf{r}}})

Based on these results, a high symmetry and low drift can be obtained by selecting voltages either in the second quadrant (scheme 1), at the edge of the low-drift valley (scheme 2), or close to zero pulse voltages (scheme 3). These three options are not equivalent because they lie at different pulse voltage ranges, which influence the number of potentially available resistance states, as discussed later. Whether two configurations of the memristive system can be distinguished as distinct resistance states depends on their difference in the read current magnitude relative to the noise density of the memristive system⁶. Generally, other factors, such as random fluctuations in the operating voltage, can also cause noise and thereby limit the number of distinguishable conductance states⁶. Assuming that the noise density is independent of the pulse voltages, the incremental change of the read current can be used as a qualitative measure for the number of potentially possible resistive states. We define the read current increment metric (Delta {I}_{{rm{r}}}) as the average change of ({I}_{{rm{r}}}) per pulse number over the last set voltage sequence of the respective simulation (Fig. 6i). The results of (Delta {I}_{{rm{r}}}) are shown in Fig. 6j–l for the three pulse schemes.

In all three schemes, (Delta {I}_{{rm{r}}}) is minimum around the origin and tends to increase with the pulse voltage magnitudes. For all pulse schemes, it is small in the parameter space where the effect of the set sequence and the reset sequence on ({I}_{{rm{r}}}) compensate. This compensation leads to a low-(Delta {I}_{{rm{r}}}) band, extending along the low-drift valley in the 2nd and 4th quadrants. A second low-(Delta {I}_{{rm{r}}}) band extends along the 1st and 3rd quadrants where both pulse voltages have the same sign and cause a large drift with low symmetry. In the symmetric rectangular scheme 1 pulse scheme, the two bands lie at (|{V}_{{rm{set}}}|approx |{V}_{{rm{res}}}|) (Fig. 6j). In the asymmetric pulse scheme, the reset sequence is triangular with voltages pulses (le |{V}_{{rm{res}}}|) (Fig. 6k). Therefore, both low (Delta {I}_{{rm{r}}}) bands are rotated inwards, off the diagonals. For the symmetric triangular scheme (scheme 3), the low-(Delta {I}_{{rm{r}}}) band in the 1st and 3rd quadrant is more localized around the zero-voltage point and approximately follows the low symmetry region (Fig. 6l). The results show that there is an apparent trade-off between small drift, maximum symmetry, and large (Delta {I}_{{rm{r}}}) when selecting the pulse voltages.

Overshoot ({{boldsymbol{P}}}_{{bf{set}}}) and ({{boldsymbol{P}}}_{{bf{res}}})

To analyze the non-linear and non-monotonic behavior of the pulse curves, we define the overshoot metrics ({P}_{{rm{set}}}) and ({P}_{{rm{res}}}). The set overshoot metric ({P}_{{rm{set}}}) indicates the location of the maximum ({I}_{{rm{r}}}) in the set pulse sequence relative to the first set pulse on a normalized scale (0le {P}_{{rm{set}}}le 1) (Fig. 7a). Hence, for ({P}_{{rm{set}}}=0), the largest ({I}_{{rm{r}}}) in the respective set sequence occurs after the first set pulse, and for ({P}_{{rm{set}}}=1), the largest ({I}_{{rm{r}}}) is after the last set pulse (i.e., for (N=50)). The reset overshoot metric ({P}_{{rm{res}}}) is defined equivalently for the reset pulse sequence, but relative to the last pulse (Fig. 7b), which ensures that ({P}_{{rm{set}}}={P}_{{rm{res}}}=0) for a set process where ({I}_{{rm{r}}}) steadily decreases followed by a reset sequence, where ({I}_{{rm{r}}}) increases again.

a Illustration of the set overshoot metric P_set and b reset overshoot metric P_res. c plot of P_set as a function of the set voltage V_set and reset voltage V_res for the symmetric rectangular pulse scheme (scheme 1), d the asymmetric pulse scheme (scheme 2), and e the symmetric triangular pulse scheme (scheme 3). f Plot of P_res as a function of V_set and V_res for scheme 1, g scheme 2, and h scheme 3. i Illustration of the set linearity metric L_set for the set process and j the reset linearity metric L_res for the reset process. k Plot of L_set as a function of V_set and V_res for scheme 1, l scheme 2, and m scheme 3. n Plot of L_res as a function of V_set and V_res for scheme 1, o scheme 2, and p scheme 3.

Full size image

The overshoot metrics ({P}_{{rm{set}}}) and ({P}_{{rm{res}}}) are calculated for the last set and reset sequences of each simulation and shown in Fig. 7c–e and Fig. 7f–h over the entire voltage range. In all three schemes, two large triangular-shaped regions with ({P}_{{rm{set}}}=1) are visible in the top and bottom half of the plots. These large-overshoot regions are separated by a narrow transition edge from the rest of the parameter space where ({P}_{{rm{set}}}=0). Only for the third pulse scheme, this transition regime is smeared out over the entire parameter space. There, ({P}_{{rm{set}}}) is mostly around 0.2 and only reaches ({P}_{{rm{set}}}approx 0) in a small region around the origin. The smallest region with ({P}_{{rm{set}}}=1) is obtained in the asymmetric pulse scheme (scheme 2). Here, the triangular high-overshoot regions are limited to smaller set voltage magnitudes of (left|{V}_{{rm{set}}}right|approx 2.5,{rm{V}}) because the rectangular set sequence has a larger influence on ({I}_{{rm{r}}}) than the triangular reset sequence. The behavior of ({P}_{{rm{res}}}) is similar to ({P}_{{rm{set}}}). In all three pulse schemes, there are two triangular regions in the top and bottom half of the plots with ({P}_{{rm{res}}}ne 0). In contrast to the set process, ({P}_{{rm{res}}}) changes gradually in all three schemes with only one narrow transition edge within the low-drift valley (indicated) where ({P}_{{rm{res}}}=1). These sharp transition regions coincide approximately with the locations of the transition regions in the set process.

Linearity ({{boldsymbol{L}}}_{{bf{set}}}) and ({{boldsymbol{L}}}_{{bf{res}}})

We quantify the linearity of the set and reset process via the root mean squared deviation of the normalized ({I}_{{rm{r}}}) during the set or reset process to an ideal linear ({I}_{{rm{r}}}) connecting the first and last pulse of the respective pulse sequence as illustrated in Fig. 7i, j. The linearity metrics ({L}_{{rm{set}}}) and ({L}_{{rm{res}}}) of the set and reset process are scaled to ensure that a large value is associated with a small deviation (i.e., a high linearity). The linearity metrics are calculated for the last set and reset pulse sequences of each simulation and plotted in Fig. 7k–m and Fig. 7n–p.

In all three pulse schemes, ({L}_{{rm{set}}}) is large with values of ({L}_{{rm{set}}}ge 0.8) for pulse voltages where the maximum read current is located at the beginning or the end of the pulse sequence (i.e., where ({P}_{{rm{set}}}=1) or ({P}_{{rm{set}}}=0)). Therefore, for the first and second pulse schemes, small set linearity only occurs in narrow cross-arranged parameter bands that correspond to the locations of the transition edges of ({P}_{{rm{set}}}) in Fig. 7c–e. These low-linearity bands lie directly within the low-drift valley. In the third pulse scheme, the smeared-out set overshoot transition leads to only two triangular-shaped regions at the top and bottom of the plot where ({L}_{{rm{set}}}ge 0.8) and close to 1. Outside these two regions, the set linearity is mostly around ({L}_{{rm{set}}}approx 0.7) with minimum values around ({L}_{{rm{set}}}approx 0.3).

The reset linearity behaves equivalent to the set linearity. However, because the reset overshoot changes gradually over a larger parameter regime and has only two regions with ({P}_{{rm{res}}}=0), there are only two corresponding triangular-shaped regions of high linearity with ({L}_{{rm{res}}}approx 1). While this holds for schemes 2 and 3, the reset linearity of the symmetric rectangular scheme (scheme 1) behaves qualitatively like the set linearity (i.e., mostly ({L}_{{rm{res}}}ge 0.8)) but with low ({L}_{{rm{res}}}) bands at locations corresponding to the edges of the reset overshoot.

Selecting operating parameters

Computational results from a simple 2D memristor model reveal empirical tradeoff relationships between linearity, symmetry, and drift such that optimized operating conditions of 2D memristive synaptic devices depend on the target applications. For all three pulse schemes, low drift is only found in the low-drift valley, where the overshoot has a narrow transition region causing nonlinearities. The linearity and symmetry are therefore mostly large at the edge of the valley or outside of it. Simultaneously, large (Delta {I}_{{rm{r}}}) are obtained outside the low drift valley and far away from the origin. However, depending on the applied pulse scheme, different combinations of pulse metrics are possible.

To demonstrate the implications of these results for selecting suitable pulse voltages, we select six voltage points, three on the left edge of the low drift valley (points 1–3) and three on the right edge (points 4–6). They are indicated in Fig. 8a–c on the set overshoot map of the respective pulse scheme. The corresponding pulse curves are shown in Fig. 8d–f (points 1–3) and Fig. 8g–i (points 4–6). All three pulse schemes behave similarly in their fundamental pulse update behavior. For the first three points on the left edge (points 1–3), the overshoot is ({P}_{{rm{set}}}=1) and ({P}_{{rm{res}}}approx 1) and correspondingly, the read current changes irregularly during the first two to four pulse sequences. After these initial sequences, the positive set voltage causes an increase in ({I}_{{rm{r}}}), while the negative reset voltage reduces ({I}_{{rm{r}}}) (Fig. 8d–f). The other three points on the right edge (4–6) lie outside the overshoot regime (({P}_{{rm{set}}}={P}_{{rm{res}}}=0)). Consequently, their pulse update behavior is reversed. Here, a positive set voltage reduces ({I}_{{rm{r}}}), and a negative reset voltage increases ({I}_{{rm{r}}}) (Fig. 8g–i). This reversal of the read-current update behavior follows the mechanism explained in Section “High-throughput pulse simulations”. Differences between the three pulse schemes become apparent when considering linearity and symmetry.

a–c Overshoot metric maps for the three pulse operation schemes with symmetric rectangular pulse trains (scheme 1), asymmetric pulse trains (scheme 2), and symmetric triangular pulse trains (scheme 3). Six points (1–6) are indicated in each map to denote the set voltages V_set and reset voltages V_res that correspond to the pulse curves in (d–f) and (g–i). d–f Read current I_r as a function of the pulse number N for the three different pulse operation schemes (scheme 1–3) and three different V_set and V_res combinations at the right edge of the overshoot transition regime as indicated in the metric maps in (a–c) (white, dashed line). g–i Read current I_r as a function of the pulse number N for the three different pulse operation schemes (scheme 1–3) and three different V_set and V_res combinations at the left edge of the overshoot transition regime as indicated in the metric maps in (a–c).

Full size image

In the symmetric rectangular pulse scheme (scheme 1), points 1–3 on the left edge of the low-drift valley offer high symmetry (0.75–0.86) as well as moderate to high set and reset linearities (0.77–0.93). The set/reset linearities increase towards the origin from 0.77/0.81 (point 1) to 0.85/0.87 (point 2) and 0.93/0.93 (point 3), which is reflected in the pulse curves in Fig. 8d by an increasingly triangular shape of the set–reset sequences from point 1 to point 3. On the right edge of the low-drift valley (points 4-6), the symmetry is slightly higher (0.85–0.9) but the set/reset linearities are smaller, with values of 0.78/0.72 (point 4), 0.81/0.77 (point 5), and 0.84/0.82 (point 6). Therefore, a triangular-shaped pulse curve as for the point 3 is not obtained (Fig. 8g). Instead, the reset-set sequences are asymmetric and concave-shaped.

In the asymmetric pulse scheme (scheme 2), the points 1–3 on the left edge of the low-drift valley offer the highest symmetry among the analyzed points, with values between 0.9 and 0.94. The set/reset linearities are notably smaller compared to those in scheme 1, with values of 0.74/0.67 (point 1), 0.82/0.71 (point 2), and 0.89/0.74 (point 3). As a result, the pulse curves show mostly symmetric set-reset sequences with nonlinear, concave shapes. On the right edge of the low-drift valley, the symmetry is significantly smaller with values between 0.71 (point 4) and 0.65 (point 6). Simultaneously, the set/reset linearities are the highest observed among the analyzed points, with values of 0.8/0.98 (point 4), 0.82/0.97 (point 5), and 0.85/0.94 (point 6). As a result, all three pulse curves (4–6) show triangular-shaped and mostly piecewise linear reset–set sequences.

In the symmetric triangular pulse scheme (scheme 3), the points 1–3 on the left side of the low-drift valley offer intermediate symmetries with values between 0.76 (point 1) and 0.81 (point 3). While the set linearity has values close to one with 0.93 (point 1), 0.98 (point 2), and 0.96 (point 3), the reset linearity is small with values of 0.66 (point 1), 0.70 (point 2), and 0.74 (point 3). On the right edge of the low-drift valley, the symmetry and set linearity are both the smallest among the analyzed points with values of 0.45–0.54, and 0.57–0.61, respectively. The linearity of the reset process is large between 0.89 (point 3) and 0.95 (point 1). Hence, for both sets of points, clear differences are apparent between the set and reset linearity, which leads to a reduced symmetry metric and to reset-set sequences that are neither triangular nor symmetric and concave, but instead asymmetrically tilted and nonlinear (Fig. 8f, i).

This comparison demonstrates that the three pulse schemes offer different features and combinations of these features for tuning the synaptic behavior of memristive devices. It is shown that the third pulse scheme results in the most nonlinear and asymmetric behavior, while the first scheme offers mostly nonlinear but more symmetric pulse curves. Depending on the bias conditions, the asymmetric (second) pulse scheme allows for symmetric but nonlinear pulse curves or highly linear triangular pulse curves with a small asymmetry. In all three schemes, changing the voltage regime can result in a reversal of the read current update behavior from LTD to LTP. Hence, for this case of symmetric Schottky barriers and fixed material and device parameters, a large range of synaptic tuning can be achieved by modifying the voltage conditions. This analysis also shows that the complexity of memristive device behavior can arise from the dynamics of a seemingly simple but fundamental mechanism. Generally, the dependency of linearity and symmetry on the pulse schemes seems to be qualitatively consistent with experiments on MoS₂ memtransistors barring the additional complexity of gate-dependent synaptic plasticity¹⁹, where improved linearity and symmetry were observed for pulse scheme 2, and relatively larger asymmetry for scheme 3, albeit with smaller nonlinearity. However, our results also imply that the linearity-symmetry relations among the three pulse schemes discussed here are not intended as general rules for the entire parameter space. Instead, they depend significantly on the pulse operation regime and other material and device parameters, such as the amplitude of Schottky barrier heights and the relative symmetry of Schottky barrier heights at two contacts (see Section 2.5). Understanding this behavior represents the first step towards a comprehensive understanding of more complex memristive systems and memtransistors that involve multiple memristive mechanisms such as trap states^{56,57,58,59,60}, phase transitions^59,61,62, and the application of additional gate bias^21,23,37,63.

Schottky barriers and nonidealities

In this work, we analyzed the behavior of the read current for the idealized case of small and symmetric Schottky barriers at the contacts and demonstrated that varying V_set and V_res can result in a complex and highly nonlinear behavior of the read current that changes qualitatively depending on the voltage regime. In a more realistic scenario, variations in the material properties can significantly alter the read-current behavior⁶⁴. For example, the charges introduced by additional vacancies or grain boundaries are expected to influence the electron and vacancy mobilities⁶⁵, alter the band gap and trap density of states^66,67,68, and increase the background doping by introducing stationary charges. Further, defects and inhomogeneities at the interfaces can alter the interface potential⁶⁴, which can be considered in the model by varying the effective Schottky barrier height. Previously, it was shown that asymmetry in Schottky barrier heights at two contacts can qualitatively change the ({I})–({V}) characteristics, such as the hysteresis direction and symmetry of the memristive loop³⁸. Therefore, we expect that these parameters also have a considerable influence on the pulse characteristics.

To assess the influence of semiconductor-metal contacts, we perform a series of simulations with variations in the left and right Schottky barrier heights ({phi }_{0,{rm{L}}}) and ({phi }_{0,rm{R}}). All other model parameters are the same as in Fig. 4a (scheme 1, case II). First, we consider three different Schottky barrier configurations with small variations in the barrier height, namely, a symmetric case with ({phi }_{0,{rm{L}}}={phi }_{0,{rm{R}}}=0.001,{rm{eV}}), a symmetric case with ({phi }_{0,{rm{L}}}={phi }_{0,{rm{R}}}=0.05,{rm{eV}}), and an asymmetric case with ({phi }_{0,{rm{L}}}=0.001,{rm{eV}}) and ({phi }_{0,{rm{R}}}=0.05,{rm{eV}}). The resulting read currents are shown in Fig. 9a as functions of the pulse number. All three configurations result in qualitatively identical curves with only minor differences. Secondly, we consider a Schottky barrier ({phi }_{0,{rm{L}}}=0.1,{rm{eV}}) at the left contact and calculate the read current curves for several values of the Schottky barrier at the right contact, namely ({phi }_{0,{rm{R}}}=0.22,{rm{eV}}), (0.25,{rm{eV}}), (0.3,{rm{eV}}), and (0.35,{rm{eV}}) (Fig. 9b). All four read current curves differ significantly and have smaller read currents for larger ({phi }_{0,{rm{R}}}). Additionally, the first three curves differ also qualitatively. The first configuration (({phi }_{0,{rm{R}}}=0.22,{rm{eV}})) yields a highly nonlinear and irregular read current curve owing to the same dynamic phenomenon as elaborated in Fig. 5c earlier. With increasing ({phi }_{0,{rm{R}}}), the read current curves start to recover the typical potentiation-depression behavior.

a Read current as a function of the pulse number for small variations in the left and right Schottky barriers ({phi }_{0,{rm{L}}}) and ({phi }_{0,{rm{R}}}), and b for larger variations in ({phi }_{0,{rm{R}}}), while ({phi }_{0,{rm{L}}}=0.1,{rm{eV}}) is kept constant. All simulations use pulse amplitudes of V_set = −5 V and V_res = 4 V.

Full size image

These results demonstrate that small variations and asymmetries (< 0.05 eV) in the Schottky barriers are expected to cause only minor quantitative changes in the read current curves for initially symmetric and small Schottky barriers. However, the influence of variations in the Schottky barriers depends significantly on the configuration of the Schottky barriers and can lead to qualitative changes for slightly larger variations and asymmetries.

We used high-throughput charge-transport simulations to develop comprehensive correlations between voltage pulse parameters and quantitative metrics of synaptic learning behavior in a lateral 2D memristive system. The pulse dynamics for three different voltage pulse schemes describe a vacancy-dominated hysteresis mechanism in these devices that agrees with experimental data³⁵. This mechanism leads to non-filamentary soft resistive switching that could be particularly suited for multi-bit memory operations and neuromorphic learning^2,3,6. The synaptic learning behavior was defined by the linearity of the set and reset process, the symmetry, drift, and overshoot of the set and reset process, and the current increment per pulse. The results provide valuable insight in selecting operation voltages. Voltage conditions that simultaneously yield high linearity, high symmetry, low drift, and a large current increment per pulse require making compromises. For all three pulse schemes, low drift is mostly found in a narrow voltage regime, where the current increment and the linearity are small. The linearity and symmetry are in most cases either large at the edge of this voltage regime or outside of it. Despite these similarities, we find that the three pulse schemes lead to different combinations of high or low symmetry, and linearity. In particular, the asymmetric (second) pulse scheme allows for either symmetric but nonlinear pulse curves or highly linear triangular pulse curves with a small asymmetry, depending on the bias conditions. The other two pulse schemes are mostly nonlinear and asymmetric (scheme 3) or offer mixed scenarios (scheme 2). However, the linearity and symmetry of potentiation and depression curves depend significantly on several parameters such as voltage regime and (a)symmetry of the Schottky barrier heights at two contacts. In all three schemes, changing the voltage regime can result in high nonlinearity and, eventually, in reversing the read current update behavior from LTD to LTP, which is explained in detail by the vacancy dynamics.

The analysis shows that a large range of synaptic tuning can be achieved by modifying the voltage conditions even for the simplified case of symmetric Schottky barriers and fixed material and device parameters. It is demonstrated that complexity can arise from a seemingly simple but fundamental hysteresis mechanism that considers the agglomeration and depletion dynamics of ionic vacancies as the primary source of hysteresis. The migration of such ionic point defects is believed to be the origin of resistive switching in most 2D memristors and memtransistors based on TMDCs^13,35,36,37. Hence, the presented analysis marks the first step towards a comprehensive theoretical picture of more complex memristive systems and memtransistors that involve multiple mechanisms and multi-terminal geometries^21,23,37,63. Understanding the memristive behavior of such systems could open a gateway to tailoring and controlling the bio-realistic synaptic behavior that is uniquely enabled by the physical properties of 2D materials and van der Waals heterojunctions¹². Overall, this model not only explains the experimental results from previously reported 2D memristive devices but can also serve as a valuable tool in the design of future neuromorphic systems enabled by 2D analog memories and synaptic devices.

Definition of pulse metrics

For the definition of the linearity metrics ({L}_{{rm{set}}}) and ({L}_{{rm{res}}}), the read current of the selected set and reset sequence are scaled to obtain scaled read current values (0le {I}_{{rm{r}},{rm{n}},beta }(N)le 1), following

with the read currents ({I}_{{rm{r}},beta }) of the set and reset sequence, the minimum value ({I}_{{rm{r}},beta ,min }) and the maximum value ({I}_{{rm{r}},beta ,max }) of the respective ({I}_{{rm{r}},beta }). Similarly, we also scale the pulse axis of each pulse sequence to values between (0le Nle 1). Next, we interpolate ({I}_{{rm{r}},{rm{n}},beta }) linearly on the new pulse axis between its initial value at (N=0) and the last value, which is now at (N=1) on the scaled pulse axis. The linear interpolated functions for the selected set and reset sequences are then given by

with the offset (b) on the scaled read current axis. Using the two equations above, we obtain the root-mean-squared deviation (Delta {I}_{{rm{r}},{rm{n}},beta }) between the linear interpolation ({I}_{mathrm{int},beta }) and the scaled read current ({I}_{{rm{r}},{rm{n}},beta }), as

with the number of ({N}_{beta }=50) pulses in the set/reset sequence. Because the read current axis and the pulse number axis are scaled to values between zero and one, it is (Delta {I}_{{rm{r}},{rm{n}},beta }le 1). This allows us to define the linearity metrics as

For calculating the overshoot metrics ({P}_{{rm{set}}}) and ({P}_{{rm{res}}}), we numerically determine the locations ({N}_{max ,beta }) (with (beta in left{{rm{set}},{rm{res}}right})) of the local maximum in the scaled read current (Eq. 1) on the scaled pulse axis ((0le Nle 1)) of the selected set or reset sequence. The set and reset overshoot metrics are then defined as

For calculating the symmetry metric, we use the scaled read current of the set and the following reset process ({I}_{{rm{r}},{rm{n}},beta }(N)) (with (beta in left{{rm{set}},{rm{res}}right})) (Eq. 1) on the scaled pulse axis. From these, we obtain the root-mean-squared deviation between the read current in the set sequence ({I}_{{rm{r}},{rm{n}},{rm{set}}}(N)) and the reversed reset process ({I}_{{rm{r}},{rm{n}},{rm{res}}}(1-N)), as

which allows us to define the symmetry metric with (0le Sle 1) as

The read current increment metric (Delta {I}_{{rm{r}}}) is defined as the average change of the read current per pulse within the selected set sequence

For calculating the drift metric (D), we first normalize the read current ({I}_{{rm{r}}}) to its mean value considering all 500 pulses. A linear least squared fit is performed on the normalized read current. The slope of the linear equation obtained from the fit is defined as the drift (D).

Numerical charge-transport model

The numerical charge-transport model is based on a finite volume discretization in the spatial domain and is implemented in the Julia programming language utilizing the VoronoiFVM solver^51,69. It solves the drift-diffusion equations for the quasi-Fermi potentials ({varphi }_{alpha }) with (alpha in left{{rm{n}},{rm{p}},{rm{x}}right}), of electrons n, holes p, and mobile, single-positively charged sulfur vacancies x

with the respective carrier charges ({q}_{alpha }), carrier densities ({n}_{alpha }), and the current densities

Here, the carrier mobilities are denoted as ({mu }_{alpha }), the positive elementary charge is (q), and the charge numbers are ({z}_{{rm{n}}}=-1), ({z}_{{rm{p}}}=+1), and ({z}_{{rm{x}}}=+1). The charge-carrier densities are expressed as (with (alpha in {{rm{n}},{rm{p}},{rm{x}}}))

with the temperature (T), the Boltzmann constant ({k}_{{rm{B}}}), and the effective densities of state of the conduction band and valence band ({N}_{{rm{c}}}={N}_{{rm{n}}}) and ({N}_{{rm{v}}}={N}_{{rm{p}}}) in parabolic approximation. The intrinsic band edge energies of the conduction and valence bands for electrons and holes are denoted as ({E}_{{rm{c}},0}={E}_{{rm{n}},0}), ({E}_{{rm{v}},0}={E}_{{rm{p}},0}), respectively, and the vacancy energy level is denoted as ({E}_{{rm{x}},0}). Different statistics functions ({{mathscr{F}}}_{alpha }) distinguish the physical nature of the charge carriers, namely the Fermi-Dirac integral of order ½ (({{mathscr{F}}}_{{rm{n}}}) and ({{mathscr{F}}}_{{rm{p}}})) for electrons and holes, and the Fermi-Dirac integral of order -1 (({{mathscr{F}}}_{{rm{x}}})) for vacancies⁷⁰. The above equations are coupled to Poisson’s equation via the electrostatic potential (psi)

with the effective background donor doping density (C) with charge ({q}_{{rm{C}}}=-q), the electric vacuum permittivity ({{rm{varepsilon }}}_{0}) and the quasi-static relative electric permittivity ({{rm{varepsilon }}}_{{rm{r}}}). At the electrode–semiconductor contacts at the locations ({x}_{1}) and ({x}_{2}), Schottky-contact boundary conditions are applied with the intrinsic Schottky barriers ({phi }_{0,{rm{L}}}={phi }_{0}({x}_{1})) and ({phi }_{0,{rm{R}}}={phi }_{0}({x}_{2})). Image-charge-induced Schottky barrier lowering (SBL) is included by solving for an additional auxiliary electrostatic potential. Details about the equation system, the formulation of boundary conditions, and SBL can be found elsewhere³⁸.

Model throughput

To evaluate the efficiency, we performed numerical experiments measuring the simulation time over 100 simulations with random V_set and V_res values on a single-CPU system (Intel Core i9-14900KF). We achieved a throughput of approximately 91 ± 1 simulations per hour. For the high-throughput simulations, we distributed the simulations on three CPU systems to achieve approximately 273 simulations/h translating to 6550 simulations/day.

Model parameters

The material parameters used for all high-throughput pulse simulations and the (I)–(V) curve simulations are listed in Table 1 (Set 1). The (I)–(V) curve simulations are performed with the voltage parameters in Table 2, and the pulse voltage parameters for the high-throughput simulations are listed in Table 3. The device geometry parameters used for all simulations are listed in Table 4 To improve numerical stability, we ramp up and down the pulse voltage linearly for each set, reset, and read pulses, as described to in Fig. 2a. For the simulation in Fig. 2a, we increased the read pulse duration compared to the measurement from 100 µs to 400 µs to maintain numerical stability at reasonably large time steps. The simulations in Figs. 4, 5, and 9 use the simulation parameters of the high-throughput simulations. For the high-throughput pulse simulations, the set/reset pulse time of 5 ms includes 1.5 ms in the beginning of the pulse for ramping it up, and 1 ms at the end of the pulse for ramping it down. The read pulse of length 0.14 ms includes 0.02 ms at the beginning and end of the pulse for the ramp. The read current is always extracted from the simulations as the current occurring during the last time step of the respective read pulse before the ramp down.