On the switching mechanism and optimisation of ion irradiation enabled 2D MoS₂ memristors

Abstract

Memristors are prominent passive circuit elements with promising futures for energy-efficient in-memory processing and revolutionary neuromorphic computation. State-of-the-art memristors based on two-dimensional (2D) materials exhibit enhanced tunability, scalability and electrical reliability. However, the fundamental of the switching is yet to be clarified before they can meet industrial standards in terms of endurance, variability, resistance ratio, and scalability. This new physical simulator based on the kinetic Monte Carlo (kMC) algorithm reproduces the defect migration process in 2D materials and sheds light on the operation of 2D memristors. The present work employs the simulator to study a two-dimensional 2H-MoS₂ planar resistive switching (RS) device with an asymmetric defect concentration introduced by ion irradiation. The simulations unveil the non-filamentary RS process and propose routes to optimize the device's performance. For instance, the resistance ratio can be increased by 53% by controlling the concentration and distribution of defects, while the variability can be reduced by 55% by increasing 5-fold the device size from 10 to 50 nm. Our simulator also explains the trade-offs between the resistance ratio and variability, resistance ratio and scalability, and variability and scalability. Overall, the simulator may enable an understanding and optimization of devices to expedite cutting-edge applications.

---

1. Introduction

The development of Information and Communication Technology (ICT), such as 5G circuits and the internet of things, demands breakthroughs in non-volatile memory¹ since the state-of-the-art flash technology is hitting its physical limits in terms of power consumption and device scalability.² Among the emerging memory technologies, memristors have attracted much research interest and become the herald of next-generation computational architecture, revolutionizing ICT.^3–6 The potential of memristors is rooted in their superior performance (e.g., ultrafast switching, low power consumption, data retention and endurance), intrinsically high scalability, stackability, compatibility with complementary metal oxide semiconductor (CMOS) technology and flexibility for wearable applications.^1,2,7 However, device variability, including cell-to-cell variability and cycle-to-cycle variability, hinders the industrial deployment of memristor technology. The cell-to-cell variability is the inhomogeneity between devices via the same fabrication process and cycle-to-cycle variability is related to the operation of individual devices,^1,7 which emerges from the stochastic processes during the resistive switching (RS). For example, in filamentary RS devices the location and morphology of conducting filaments may vary between cycles and cells.^8,9 Device variability is inevitable in RS devices which relies on material defects,^10–13 while the forming process exacerbates the problem by varying defect concentration and distribution.¹⁴ Although the cycle-to-cycle variability can be insignificant in some applications exemplified by neuromorphic imaging recognition,¹ it imposes major limitations on memristor applications, rendering a range of high-end applications impractical (e.g., multilevel information processing and long-term storage⁷). For example, device variability causes signal degradation in crossbar arrays¹⁵ and hampers projections of device lifetime, demanding excessive budget in testing.¹⁶ Verification and iterative approaches may mitigate the resistance variability in multilevel information¹⁷ and radiofrequency applications respectively.¹⁸ However, a consensus on a figure of merit for the variability issue is yet to be achieved and a lack of comparable statistics on device variability remains a main obstacle to the technology.⁷

Recently, memristive behaviour has been observed in a range of two-dimensional (2D) materials.^19–25 This may expedite the industry deployment of the memristor technology since 2D memristors exhibit superior tunability, scalability and electrical reliability.^4,19,26–34 These characteristics arise from the ultrathin nature and unique mechanical, electric and optoelectronic properties of 2D materials.^26,27 The 2D memristor landscape shows diverse device architectures and switching mechanisms. For example, the switching of 2D vertical memristors^20,35–38 depends on the formation and rupture of conductive filaments. 2D planar memristors may rely on phase transition,^39,40 charge trapping/de-trapping,^41,42 electron tunnelling modulated by polarization,⁴³ electrochemical processes⁴⁴ and Schottky barrier modulation.²⁹ Defect migration plays a crucial role in these processes.^32,45–47 Planar 2D memristors show great promise in the future of energy-efficient neuromorphic computations. The planar architecture allows effective gate tuning and multi-terminal operation, enhancing device controllability and enabling complex neural functionalities.³² Their low intrinsic capacitance supports fast switching and low power consumption. We note the planar 2D memristor is still in its infancy. The prototype devices operate at much higher voltages (∼10s V) than the state-of-the-art vertical 3D memristors (∼1–2 V). It is imperative to improve the device performance. Compared with their 3D counterparts, the device variability of 2D memristors has rarely been explored, while it is of utmost importance to gain in-depth knowledge of the switching process and mitigate the variability issues in 2D memristors.

Physical simulators are indispensable for understanding the resistive switching process and they can greatly facilitate investigations on device viability.^7,48,49Ab initio calculations can accurately relate the resistive switching to the defect creation,⁵⁰ the electronic structure and transport properties in nanostructures.⁵¹ However, they are limited to relatively small volumes (a few nm as maximum) and short times (shorter than ns) and can hardly reproduce RS processes at the device level.⁴⁸ Continuous models are apt to describe the average behaviours of individual devices or even circuits. However, they overlook the microscopic characteristics of the system, such as particle migration,^52–55 and are hence not suitable to investigate the stochasticity that emerges from the evolution of microscopic configuration of the system during the switching process. The kinetic Monte Carlo (kMC) algorithm is an established technique to study the microscopic evolution and its impacts on device performance.^{11,12,56–59} For example, the kMC algorithm can reproduce particle diffusion and the formation and rupture of percolation paths involving several RS cycles of memristors.^11,56–58 Although the kMC algorithm and continuous models have been employed to explore defect accumulation and electrical conduction in 2D materials,^29,38 the simulation of 2D memristors at the device level is scarce and the variability of 2D memristors has been rarely explored. In this work, we build a new physical kMC simulator for the defect migration process in 2D materials and shed light on the operation of 2D memristors. We investigate the effects of initial vacancy distribution, the scaling limits and the factors that regulate device endurance and variability. We use experimental data from the MoS₂ memristor enabled by site-specific defects to collate and verify the model.¹⁹ The simulator helps understand the physics of resistive switching in 2D materials, offering practical guidance to optimize device performance.

2. Experimental results and simulation approach

Fig. 1a shows a device schematic. The device exhibits a planar metal–semiconductor–metal structure and the semiconductor channel is a 2D 2H-MoS₂ with an asymmetric defect concentration introduced by focused helium ion irradiation (see more details about the device fabrication in ESI, section 1†). We note that the irradiation-induced defects enable the resistive switching, while devices of pristine MoS₂ do not switch.¹⁹ The defect region, referred to as the fissure, bisects the channel with an asymmetric concentration along the horizontal direction, i.e., across the electrodes. For example, the right tail of the peak (tens of nanometres) in Fig. 1a is much wider than the left (<1 nm).

Fig. 1 Defect enabled memristive device. (a) Schematic device structure shows the fissure of defects in the MoS₂ channel. The 50 × 50 nm simulation domain represents the defect profile of a skewed Gaussian distribution, while the channel is on the micrometre scale. The arrows in the lattice structure indicate the possible routes of vacancy migration, i.e., a S atom (yellow) exchanges its position with a nearest-neighbour vacancy site (grey). (b) 14 representative experimental I–V hysteresis loops sampled from continuous 1200 cycles. (c) The evolution of the resistance ratio (calculated at −10 V) over cycles (time in the top x-axis). The dashed line is an exponential fitting, and the shaded indicates the 95% prediction interval.

In our simulation, we have assessed three initial defect distributions with different asymmetries (see ESI, section 2†), i.e., a skewed Gaussian distribution (Fig. S1a†), a triangle function (Fig. S1b†) and a step function (Fig. S1c†). The skewed Gaussian distribution emulates the asymmetric distribution of defects observed experimentally,¹⁹ and the triangle and step functions of different asymmetries are to explore possible routes for device optimization. The source electrode as convention is grounded and placed at the left side of the asymmetric peak. The polarity of applied voltage regulates the direction of the fields with respect to the asymmetric defect distribution, i.e., a positive voltage (V_d > 0) indicates the electric field points from the abrupt edge of the fissure to the long tail side.

Fig. 1b shows 14 experimental quasi-static pinched hysteresis loops by selecting every 15th loop from 1200 consecutive cycles measured at a triangular voltage ramp between ±35 V with a rate of 6 V s⁻¹. The device shows bipolar operation, and the SET (RESET) process occurs at a positive (negative) voltage. The current level increases with the cycling. Fig. 1c shows the resistance ratio over cycle, which shows that cycling the device leads to an exponential decay (shown as a dashed line) in resistance ratio (21% after 1000 cycles). These degradation processes will be further discussed alongside the simulation work.

To investigate the device switching process, we implemented the kinetic Monte Carlo (kMC) algorithm⁶⁰ using MATLAB in a simulation domain of 50 × 50 nm. The vacancy distribution within the fissure evolves with the external electric field applied via the electrodes,¹⁹ while no defects escape the simulation domain during the switching. Since the fissure region dominates the device resistance, the size of the simulation domain is sufficient to include the main physical processes involved in the resistive switching. The defects are doubly charged sulfur vacancies because the helium-ion irradiation preferentially removes sulfur atoms from the MoS₂ lattice.³³ The generation of new defects under the applied electric field is negligible due to the high activation energy (5.85 eV (ref. 61) for vacancy and >5 eV for antisite defect⁶²). The sulfur vacancy migrates via the exchange of the vacancy position with one of the adjacent sulfur atoms,^29,63 as shown in the Fig. 1a. In the quasi-static switching, the vacancy migration events occur at a much longer time scale (∼s) than the lattice vibration (10⁻¹³ s), so the system is in thermodynamic equilibrium for any vacancy distributions.⁶⁰ Furthermore, the field-driven migration renders the reverse migration negligible, validating the kMC approach. In our simulation, we combine the electric response with the thermal effect since the lattice temperature modulates the transition rate, which is given by Maxwell–Boltzmann statistics and transition state theory (TST), i.e. Γ = ν exp(−E_A/K_BT), where ν = 7 × 10¹³ s⁻¹ is the vibration constant of the particle, K_B the Boltzmann constant, T the temperature and E_A the activation energy of the migration. The activation energy is modulated by the local electric field as follows: E_A = E⁰_A − b·F(x,y),¹² where F(x,y) is the electric field, b the polarization factor and E⁰_A = 2.297 eV the activation energy for migration in the zero-field condition.⁴ The kMC algorithm weighs all the possible migrations and chooses the evolvement path. It should be noted that the larger the transition rate, the smaller the time step t = −ln(rand)/∑Γ, where rand stands for a random number between 0 and 1 and ∑Γ is the summation of the transition rates for all possible migrations. For each defect distribution during the switching, the local vacancy concentration is averaged over 6 × 6 grid points (3.2 nm²), which determines the local resistance via the empirical relationship ,³³ (n is a parameter extracted from the experimental results, see ESI, section 4†). The electric field screening is assumed to be a linear function of the local resistance since the dielectric constant in MoS₂ strongly depends on the distribution and number of sulfur vacancies.⁶⁴ For a given applied voltage, the local electric field decreases with the defect density, indicating the defect migration is a self-limiting process. Further details about the simulation can be found in the ESI, section 3.†

3. Simulation results and discussion

Fig. 2a shows a typical simulated pinched hysteresis loop from a defect distribution of skewed Gaussian profile under a voltage ramp rate of 0.71 V s⁻¹ between 35 V and −35 V. Prior to the switching, the fissure region is 8 nm wide with a peak density ρ = 5.64 V_s nm⁻² (see the probability distributions in ESI, section 2†). The I–V sweep follows the directions indicated by the arrows. The loop shows a SET process with a positive voltage and a RESET process with a negative voltage, suggesting the same bipolar switching behaviour observed experimentally. The simulation corroborates that the operation of the device does not need a forming process, facilitating circuit simplicity.^65,66 The switching is progressive, in contrast to an abrupt resistance change, suggesting the absence of filamentary conduction. This explains the observed low power consumption.⁷ The resistance ratio (r = R_off/R_on) is 1.44 calculated at −4 V during the RESET process and the maximum current level is 3 μA. It is interesting to note the simulated loop reproduces the asymmetry found experimentally between the SET and RESET processes (see Fig. 1b).

Fig. 2 Switching mechanism. The initial defect profile of the device exhibits a skewed Gaussian distribution with a peak density of 5.64 V_s nm⁻² and a width of 8 nm. The voltage ramp starts from positive polarity with a rate of 0.71 V s⁻¹ between 35 V and −35 V. (a) A typical simulated I–V pinched hysteresis loop. Numerals on the loop mark four representative states of the switching process. The SET process takes place from I to II, and for the RESET from III to IV. The colour maps shown in the right panel correspond to the microscopic configuration of defects in the LRS (top) and HRS (bottom). (b), (c) and (d) are the defect density, local resistance and electric field profiles along the x-axis during the SET, respectively. (e), (f) and (g) are the correspondent profiles during the RESET processes.

The switching is due to the reconfiguration of the defect distribution within the fissure by the electric field. The two colour-maps attached to Fig. 2a show the microscopic distributions of defects in the high resistance state (HRS) and the low resistance state (LRS), corresponding to the state labelled by (I) and (II) on the loop respectively (more microscopic details about the states I–IV can be found in ESI, section 5†). The defects accumulate within a small region in the HRS, leading to a much higher density than the LRS. Fig. 2b details the defect evolution during the SET process. The defect profiles are extracted at a series of sequential voltages from state I to II (see Fig. 2a) and the colours of the curves shifting to blue with increasing time. The initial vacancy concentration exhibits a skewed Gaussian distribution (the green curve), mimicking the defect profile by the ion irradiation. The external field of positive polarity gradually drives the defects away from the fissure region, lowering the peak and leading to a plateau of the distribution across a 10 nm wide region. The local resistance varies with the defect profile (see Fig. 2c). The spatial-varying resistance regulates the electric potential distribution in the channel when an external voltage is applied (see Fig. 2d). The larger the resistance of a region, the larger the electric field, and the more significant the vacancy drifting. Therefore, the self-adaptive electric field reduces the vacancy concentration and hence the resistance. From state I to II, the peak resistance reduces by an order of magnitude, leading to an overall reduction in the channel resistance and hence the SET process.

Fig. 2e reveals the evolution of the defect profile during the RESET process from state III to IV (see Fig. 2a). The defects move towards and accumulate at the abrupt edge of the fissure (x = 22 nm in Fig. 2e), which increases the local resistance (see Fig. 2f) and hence the local electric field (see Fig. 2g). On the long tail of the peak (i.e., x > 22 nm), the field pushes defects from a wide region (22 nm < x < 30 nm) towards the peak, while on the left side of the peak the field drops drastically within a 1 nm region. This asymmetric distribution of the field limits the escape of the defects from the fissure into the left side of the channel and causes the defect accumulation at the peak, recovering the initial defect configuration. The drift of defects to the left side of the fissure may become important if the field is sufficiently high where the RESET process will fail, leading to device failure (see Fig. S4†). The simulator can predict the operational range of the voltage for the device. For example, a device with a peak density of 5.64 V_s nm⁻² and a width of 12 nm can be stressed with −30 V for 17 s or −40 V for 2 s before the device fails (see Fig. S4a†).

The simulator allows us to investigate the device performance. Fig. 3a shows the temporal evolution of the device resistance to a triangular voltage ramp (2.1 V s⁻¹) between 25.2 and −25.2 V. The initial defect profile has a skewed Gaussian distribution with a peak density of 5.64 V_s nm⁻² and a width of 9 nm. The voltage ramp starts from a positive polarity +V_d. The device switches over 45 cycles (2161 s). Both the R_on and R_off drop 57% after the first ten cycles and reaches a steady state where R_on and R_off drop slower (17% in 35 cycles). The resistance drop occurs due to a progressive reduction in the peak density of the defect profile, which dominates the device resistance. This is evident in Fig. 3b and c, which show the evolution of the density profile with consecutive RESET and SET processes. Here, we can see that cycling gradually relocates the defects into the originally low-density (right tails) regions (below 3.5 V_s nm⁻²) at the expense of the peak density. This reduces the defects population of the peak region from 31% in cycle 1 to 22% in cycle 45, as can be seen in Fig. S5a.† The defect accumulation in the tail region (>25 nm onward in the x-axis) stems from the low electric field in the region, which is 16% of the electric field found in the peak (see Fig. S5b†). The low field modulates the migration barrier by 1% (see Fig. S5c†), leading to negligible field-driven migration in the tail region. The further the defects migrate away from the peak into the tail during the SET process, the harder for them to move back to the peak during the RESET process. Consequently, the defect distribution flattens across the fissure during the cycling leading to a gradual reduction in the device resistance. This flattening of the defect profile progressively reduces the difference between the OFF and ON states, reducing the ON/OFF ratio and eventually leading to device failure (see also the resistive switching loops in Fig. S6†). As observed experimentally, the power consumption increases (see the rise from 20 μW to more than 50 μW in 45 cycles in Fig. S7†) with the resistance reduction. Nevertheless, the simulation qualitatively explains the fatigue behaviour (see Fig. 1c) and allows quantitative prediction. For devices with long cycling tolerance, we can fit the behaviours of R_on and R_off (red and blue dashed lines in Fig. 3a respectively) and project the evolvement of ON/OFF ratio for longer times (Fig. 3d).

Fig. 3 Device fatigue. The initial defect profile of the device exhibits a skewed Gaussian distribution with a peak density of 5.64 V_s nm⁻² and a width of 9 nm. The voltage ramps from a positive polarity with a rate of 2.1 V s⁻¹ between 25.2 V and −25.2 V for 45 cycles (2161 s). (a) Resistance evolution over time with two exponential fittings for the LRS (in red) and the HRS (in blue) of the cycles. The shaded regions correspond to the 95% prediction interval of the fitting. (b) and (c) correspond to the density profile along the x-axis (averaged over the y-axis) for successive RESETs and SETs respectively. (d) Resistance ratio projection based on the simulation.

The most significant advantage of the simulator is its capability of exploring device variability since the kMC algorithm is apt to investigate stochastic processes. Fig. 4a shows the cycle-to-cycle variability of the resistance ratio in one 45-cycle simulation, starting with a skewed Gaussian distribution of vacancies, negative polarity, with a peak density ρ₁ = 5.64 V_s nm⁻² and a width of 9 nm. The value of resistance ratio distributes uniformly in the range of 1.2 to 1.4 with the mean of 1.29. We define the cycle-to-cycle variability of two consecutive cycles as Δr = |r_i − r_i+1|. The standard deviation of the cycle-to-cycle variability is 0.03. In Fig. 4b, we investigate the cell-to-cell variability by running 26 independent simulations for a given macroscopic distribution (see section 3 in ESI† for more details about running different simulations). Each simulation starts with a unique initial microscopic vacancy configuration in the lattice, while the macroscopic defect distribution remains the same. This scenario mimics the device fabrication process, where the sputtering process and the creation of vacancies are stochastic at the nanometre scale,^19,33 introducing cell-to-cell variability. We note the device fabrication is limited by many other parameters (e.g., ion beam stability, focusing, sample cleanliness, etc.) and the cell-to-cell variability here represents the upper limit of the ideal situation. For each such simulation, we average the resistance over 15 cycles and use the standard deviation of the cycle-to-cycle variability as the error. Here, we can observe a uniform distribution of resistance ratios around the mean of 1.31 with a standard deviation of 0.02.

Fig. 4 Device variability. The initial defect profile of the device exhibits a skewed Gaussian distribution with a peak density of 5.64 V_s nm⁻² and a width of 9 nm. The voltage ramps from a negative polarity with a rate of 2.1 V s⁻¹ between 25.2 V and −25.2 V. (a) Shows simulated cycle-to-cycle variability of the resistance ratio with an average of 1.29 (the dashed line). (b) Corresponds to the resistance ratio variability of 26 independent 15-cycle simulations initiated with the same parameters but varying microscopic defect configurations. For each simulation, the resistance ratio is the average over the 15 cycles and the error bar is the standard deviation of the cycle-to-cycle variability.

We investigate the potential of device scalability by evaluating the effects of the device dimensions (i.e., height and width) on the device resistance and the resistance ratio. For the simulations shown in Fig. 5, the initial defect profile has a skewed Gaussian distribution with a peak density of 5.64 V_s nm⁻². The device is under a voltage ramp rate of 2.1 V s⁻¹ between 25.2 and −25.2 V, which starts from the positive polarity. Fig. 5a displays experimental and simulated pinched I–V hysteresis loops with varying device height, while Fig. 5b shows the simulated temporal evolution of resistance. Both the experiment and simulation show the decrease in device height results in an overall resistance increase. This is due to the reduced conducting channels along the vertical direction. This corroborates the non-filamentary conduction. Fig. 5c demonstrates the limit of scaling down the vertical dimension. Both the resistance ratio and the device variability deteriorate as the device height reduces. A trade-off needs to be identified between the device conductance and the resistance ratio (variability), indicating a limit on the scalability of the y-dimension. Scaling along the horizontal dimension is regulated by the fissure width and the range of defect drift. Fig. 5d shows typical experimental and simulated I–V hysteresis loops from devices with varying fissure width. The temporal resistance evolution is simulated in Fig. 5e. The device resistance decreases with decreasing the fissure width, but the resistance ratio and device variability deteriorate as the fissure width is scaled down (Fig. 5f). Our simulations reveal the limits of device scaling.

Fig. 5 Device scalability. (a), (b) and (c) The effects of the device height on the switching, (d), (e) and (f) the effects of the fissure width on the switching. In (a) and (d), the experimental (simulated) data are plotted against bottom (top) and left (right) axes. For the simulation, the initial defect profile of the device exhibits a skewed Gaussian distribution with a peak density of 5.64 V_s nm⁻². The voltage ramps from a positive polarity with a rate of 2.1 V s⁻¹ between 25.2 V and −25.2 V for 15 cycles. The resistance ratio is averaged over the 15 cycles and the error is the standard deviation of the cycle-to-cycle variability.

The asymmetric nature of the initial defect distribution is crucial to the resistive switching. This indicates that device optimization may be possible by tuning the initial defect distributions. In Fig. 6 we show the resistance ratio for an initial defect distribution with a triangle shape (Fig. 6a) and a step function shape (Fig. 6b) with varying peak densities and fissure region widths. The skewed Gaussian distribution case is in Fig. S8.† In the triangle case (Fig. 6a), when the peak density is higher than 3.95 V_s nm⁻², the resistance ratio appears to exhibit a maximum when the fissure width is varied from 4 nm to 32 nm. The value of the maximum ratio decreases and occurs at a larger fissure width when the peak density decreases. Devices with the step function distribution (Fig. 6b) exhibits a similar maximum ratio. In contrast to the triangle case, the maximum appears in all the peak densities simulated and on the right side of the peak the resistance ratio falls more rapidly with increase in the fissure width compared with the triangle distribution. However, the triangle distribution enables a higher resistance ratio at a narrower fissure than both the Gaussian and step function cases, so it may offer a better option for device scaling and further performance optimization. We also note that in all the distributions simulated, the resistance ratio and variability can be enhanced by increasing the initial peak density. This is because a larger defect population enables more significant differences between the HRS and LRS.

Fig. 6 Device optimization. The dependence of the resistance ratio on the initial peak density (5.64 V_s nm⁻², 5.08 V_s nm⁻², 4.52 V_s nm⁻² and 3.95 V_s nm⁻²) and device width for two density profiles: (a) a triangle and (b) step distributions. The devices are stressed under consecutive voltage ramps with a rate of 2.1 V s⁻¹ between 25.2 to −25.2 V. The resistance ratio is averaged over 15 cycles and the error is the standard deviation of the cycle-to-cycle variability.

Finally, we note that although the voltage is high in this MoS₂-based device, the current is low, so the device energy consumption is reasonable. Besides, as the voltage drops mainly in the fissure region, we cannot address the scaling of the switching voltage by reducing the length between the drain and the source. In this sense, the minimum electric field needed to move the vacancies determine the voltage scale. Hence, it might be possible to lower the operating voltage and further reduce the power consumption by increasing the density peaks (see in Fig. S4† how the electric field has a stronger influence in higher densities), selecting materials with suitable activation energies⁶³ or by defect engineering (defects migrate easier through grain boundaries⁴).

4. Conclusions

We have developed a kMC simulator for 2D planar memristors based on defect migration using the case of a MoS₂-based device enabled by a helium ion beam microscope. The simulator reproduces the asymmetric resistive switching cycle with a performance close to observed experimentally. Besides, this approach is helpful for insights into the switching mechanism, in addition to study the device variability, the device endurance and the resistance ratio. We studied the device degradation with cycling, which produces defect relocations into low-density regions, causing a drop in the resistance ratio and reducing the switching window. Some trade-offs are proposed to tailor the device features by controlling the device size, the number of defects introduced in the channel and their distribution. For example, the device can be miniaturized at the expense of reducing the resistance ratio, increasing the variability and the resistance, or higher peak densities can be used to increase the resistance ratio and variability. We note that to reduce the switching voltage of a device based on defect migration, we must enhance the migration of defects by employing defect engineering, using higher peak densities or other transition metal dichalcogenides with lower activation energy for defect migration. Besides, the conduction mechanism during the ON state is a distributed one, which explains the drop of resistance when the device size in the y-axis is increased. We have also assessed different distributions of defect density in the channel to find some routes for device optimization. Comparing the skewed Gaussian distribution, the step function and the triangle distribution, we can conclude the triangle shape enables larger resistance ratios for narrower fissure regions. Besides, the fissure width of these distributions also affects the resistance ratio.

Author contributions

J. J. fabricated the devices and collected the experimental data. S. A. developed the simulator and conducted the simulation. S. A. and H. Z. analysed the data and wrote the manuscript. H. Z. conceived the study and supervised the project. All authors have given approval to the final version of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors gratefully acknowledge financial support by the Science Foundation Ireland under 20/FFP-P/8727.

References

1. V. Gupta, S. Kapur, S. Saurabh and A. Grover, IETE Tech. Rev., 2020, 37, 377–390.
2. S. Munjal and N. Khare, J. Phys. D: Appl. Phys., 2019, 52, 24.
3. Z. R. Wang, C. Li, W. H. Song, M. Y. Rao, D. Belkin, Y. N. Li, P. Yan, H. Jiang, P. Lin, M. Hu, J. P. Strachan, N. Ge, M. Barnell, Q. Wu, A. G. Bartos, Q. R. Qiu, R. S. Williams, Q. F. Xia and J. J. Yang, Nat. Electron., 2019, 2, 115–124.
4. L. Wang, W. G. Liao, E. E. H. Wong, Z. G. Yu, S. F. Li, Y. E. F. Lim, X. W. Feng, E. E. C. Tan, X. Huang, L. Chen, L. Liu, J. S. Chen, X. Gong, C. X. Zhu, X. K. Liu, Y. W. Zhang, D. Z. Chi and K. W. Ang, Adv. Funct. Mater., 2019, 29, 10.
5. P. Yao, H. Q. Wu, B. Gao, J. S. Tang, Q. T. Zhang, W. Q. Zhang, J. J. Yang and H. Qian, Nature, 2020, 577, 641–646.
6. M. Lanza, A. Sebastian, W. D. Lu, M. Le Gallo, M. F. Chang, D. Akinwande, F. M. Puglisi, H. N. Alshareef, M. Liu and J. B. Roldan, Science, 2022, 376, eabj9979.
7. M. Lanza, H. S. P. Wong, E. Pop, D. Ielmini, D. Strukov, B. C. Regan, L. Larcher, M. A. Villena, J. J. Yang, L. Goux, A. Belmonte, Y. C. Yang, F. M. Puglisi, J. F. Kang, B. Magyari-Kope, E. Yalon, A. Kenyon, M. Buckwell, A. Mehonic, A. Shluger, H. T. Li, T. H. Hou, B. Hudec, D. Akinwande, R. J. Ge, S. Ambrogio, J. B. Roldan, E. Miranda, J. Sune, K. L. Pey, X. Wu, N. Raghavan, E. Wu, W. D. Lu, G. Navarro, W. D. Zhang, H. Q. Wu, R. W. Li, A. Holleitner, U. Wurstbauer, M. C. Lemme, M. Liu, S. B. Long, Q. Liu, H. B. Lv, A. Padovani, P. Pavan, I. Valov, X. Jing, T. T. Han, K. C. Zhu, S. C. Chen, F. Hui and Y. Y. Shi, Adv. Electron. Mater., 2019, 5, 28.
8. A. Chen, Solid-State Electron., 2016, 125, 25–38.
9. A. Fantini, L. Goux, R. Degraeve, D. J. Wouters, N. Raghavan, G. Kar, A. Belmonte, Y. Y. Chen, B. Govoreanu and M. Jurczak, Intrinsic Switching Variability in HfO₂ RRAM, IEEE, Monterey, CA, 2013.
10. G. S. Park, X. S. Li, D. C. Kim, R. J. Jung, M. J. Lee and S. Seo, Appl. Phys. Lett., 2007, 91, 3.
11. S. Aldana, P. Garcia-Fernandez, R. Romero-Zaliz, M. B. Gonzalez, F. Jimenez-Molinos, F. Gomez-Campos, F. Campabadal and J. B. Roldan, J. Phys. D: Appl. Phys., 2020, 53, 11.
12. A. Padovani, L. Larcher, O. Pirrotta, L. Vandelli and G. Bersuker, IEEE Trans. Electron Devices, 2015, 62, 1998–2006.
13. G. Bersuker, D. C. Gilmer, D. Veksler, P. Kirsch, L. Vandelli, A. Padovani, L. Larcher, K. McKenna, A. Shluger, V. Iglesias, M. Porti and M. Nafria, J. Appl. Phys., 2011, 110, 12.
14. H. D. Kim, H. M. An, S. M. Hong and T. G. Kim, Phys. Status Solidi A, 2013, 210, 1822–1827.
15. A. Chen and M. R. Lin, Variability of resistive switching memories and its impact on crossbar array performance, IEEE, Monterey, CA, 2011.
16. Process Integration, Devices, and Structures section, International Technology Roadmap for Semiconductors, 6th edn, 2013, https://www.semiconductors.org/resources/2013-international-technology-roadmap-for-semiconductors-itrs/.
17. D. M. Nminibapiel, D. Veksler, P. R. Shrestha, J. P. Campbell, J. T. Ryan, H. Baumgart and K. P. Cheung, IEEE Electron Device Lett., 2017, 38, 736–739.
18. S. Lopez-Soriano, J. M. Purushothama, A. Vena and E. Perret, Sci. Rep., 2022, 12, 13.
19. J. Jadwiszczak, D. Keane, P. Maguire, C. P. Cullen, Y. B. Zhou, H. D. Song, C. Downing, D. Fox, N. McEvoy, R. Zhu, J. Xu, G. S. Duesberg, Z. M. Liao, J. J. Boland and H. Z. Zhang, ACS Nano, 2019, 13, 14262–14273.
20. C. B. Pan, Y. F. Ji, N. Xiao, F. Hui, K. C. Tang, Y. Z. Guo, X. M. Xie, F. M. Puglisi, L. Larcher, E. Miranda, L. L. Jiang, Y. Y. Shi, I. Valov, P. C. McIntyre, R. Waser and M. Lanza, Adv. Funct. Mater., 2017, 27, 10.
21. Y. Li, S. B. Long, Q. Liu, H. B. Lv and M. Liu, Small, 2017, 13, 35.
22. K. Qian, R. Y. Tay, V. C. Nguyen, J. X. Wang, G. F. Cai, T. P. Chen, E. H. T. Teo and P. S. Lee, Adv. Funct. Mater., 2016, 26, 2176–2184.
23. C. L. Tan, Z. D. Liu, W. Huang and H. Zhang, Chem. Soc. Rev., 2015, 44, 2615–2628.
24. J. Q. Liu, Z. Y. Zeng, X. H. Cao, G. Lu, L. H. Wang, Q. L. Fan, W. Huang and H. Zhang, Small, 2012, 8, 3517–3522.
25. Y. B. Li, A. Sinitskii and J. M. Tour, Nat. Mater., 2008, 7, 966–971.
26. S. Bertolazzi, P. Bondavalli, S. Roche, T. San, S. Y. Choi, L. Colombo, F. Bonaccorso and P. Samori, Adv. Mater., 2019, 31, 35.
27. L. Zhang, T. Gong, H. D. Wang, Z. N. Guo and H. Zhang, Nanoscale, 2019, 11, 12413–12435.
28. V. K. Sangwan, H. S. Lee, H. Bergeron, I. Balla, M. E. Beck, K. S. Chen and M. C. Hersam, Nature, 2018, 554, 500–504.
29. D. Li, B. Wu, X. J. Zhu, J. T. Wang, B. Ryu, W. D. Lu, W. Lu and X. G. Liang, ACS Nano, 2018, 12, 9240–9252.
30. J. Jadwiszczak, C. O'Callaghan, Y. B. Zhou, D. S. Fox, E. Weitz, D. Keane, C. P. Cullen, I. O'Reilly, C. Downing, A. Shmeliov, P. Maguire, J. J. Gough, C. McGuinness, M. S. Ferreira, A. L. Bradley, J. J. Boland, G. S. Duesberg, V. Nicolosi and H. Z. Zhang, Sci. Adv., 2018, 4, 10.
31. M. Wang, S. H. Cai, C. Pan, C. Y. Wang, X. J. Lian, Y. Zhuo, K. Xu, T. J. Cao, X. Q. Pan, B. G. Wang, S. J. Liang, J. J. Yang, P. Wang and F. Miao, Nat. Electron., 2018, 1, 130–136.
32. V. K. Sangwan, D. Jariwala, I. S. Kim, K. S. Chen, T. J. Marks, L. J. Lauhon and M. C. Hersam, Nat. Nanotechnol., 2015, 10, 403–406.
33. D. S. Fox, Y. B. Zhou, P. Maguire, A. O'Neill, C. O'Coileain, R. Gatensby, A. M. Glushenkov, T. Tao, G. S. Duesberg, I. V. Shvets, M. Abid, M. Abid, H. C. Wu, Y. Chen, J. N. Coleman, J. F. Donegan and H. Z. Zhang, Nano Lett., 2015, 15, 5307–5313.
34. B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti and A. Kis, Nat. Nanotechnol., 2011, 6, 147–150.
35. M. Kim, R. J. Ge, X. H. Wu, X. Lan, J. Tice, J. C. Lee and D. Akinwande, Nat. Commun., 2018, 9, 7.
36. R. J. Ge, X. H. Wu, M. Kim, J. P. Shi, S. Sonde, L. Tao, Y. F. Zhang, J. C. Lee and D. Akinwande, Nano Lett., 2018, 18, 434–441.
37. M. M. Rehman, G. U. Siddiqui, Y. H. Doh and K. H. Choi, Semicond. Sci. Technol., 2017, 32, 10.
38. C. B. Pan, E. Miranda, M. A. Villena, N. Xiao, X. Jing, X. M. Xie, T. R. Wu, F. Hui, Y. Y. Shi and M. Lanza, 2D Mater., 2017, 4, 9.
39. P. Zhang, C. X. Gao, B. H. Xu, L. Qi, C. J. Jiang, M. Z. Gao and D. S. Xue, Small, 2016, 12, 2077–2084.
40. M. Yoshida, R. Suzuki, Y. J. Zhang, M. Nakano and Y. Iwasa, Sci. Adv., 2015, 1, 6.
41. G. H. Shin, C. K. Kim, G. S. Bang, J. Y. Kim, B. C. Jang, B. J. Koo, M. H. Woo, Y. K. Choi and S. Y. Choi, 2D Mater., 2016, 3, 9.
42. M. M. Rehman, G. U. Siddiqui, J. Z. Gul, S. W. Kim, J. H. Lim and K. H. Choi, Sci. Rep., 2016, 6, 10.
43. X. Y. Xu, Z. Y. Yin, C. X. Xu, J. Dai and J. G. Hu, Appl. Phys. Lett., 2014, 104, 5.
44. T. J. Echtermeyer, M. C. Lemme, M. Baus, B. N. Szafranek, A. K. Geim and H. Kurz, IEEE Electron Devices Lett., 2008, 29, 952–954.
45. H. Zhang, W. Z. Bao, Z. Zhao, J. W. Huang, B. Standley, G. Liu, F. L. Wang, P. Kratz, L. Jing, M. Bockrath and C. N. Lau, Nano Lett., 2012, 12, 1772–1775.
46. G. N. Panin, O. O. Kapitanova, S. W. Lee, A. N. Baranov and T. W. Kang, Jpn. J. Appl. Phys., 2011, 50, 6.
47. B. Standley, W. Z. Bao, H. Zhang, J. Bruck, C. N. Lau and M. Bockrath, Nano Lett., 2008, 8, 3345–3349.
48. L. L. Gao, Q. Y. Ren, J. W. Sun, S. T. Han and Y. Zhou, J. Mater. Chem. C, 2021, 9, 16859–16884.
49. D. Ielmini and V. Milo, J. Comput. Electron., 2017, 16, 1121–1143.
50. H. P. Komsa, J. Kotakoski, S. Kurasch, O. Lehtinen, U. Kaiser and A. V. Krasheninnikov, Phys. Rev. Lett., 2012, 109, 5.
51. M. Ghorbani-Asl, A. N. Enyashin, A. Kuc, G. Seifert and T. Heine, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 7.
52. G. Gonzalez-Cordero, M. Pedro, J. Martin-Martinez, M. B. Gonzalez, F. Jimenez-Molinos, F. Campabadal, N. Nafria and J. B. Roldan, Solid-State Electron., 2019, 157, 25–33.
53. F. Jimenez-Molinos, G. Gonzalez-Cordero, P. Cartujo-Cassinello and J. B. Roldan, SPICE modeling of RRAM thermal reset transitions for circuit simulation purposes, IEEE, Barcelona, SPAIN, 2017.
54. G. Gonzalez-Cordero, M. B. Gonzalez, H. Garcia, F. Campabadal, S. Duenas, H. Castan, F. Jimenez-Molinos and J. B. Roldan, Microelectron. Eng., 2017, 178, 26–29.
55. F. Jimenez-Molinos, M. A. Villena, J. B. Roldan and A. M. Roldan, IEEE Trans. Electron Devices, 2015, 62, 955–962.
56. S. Aldana, J. B. Roldan, P. Garcia-Fernandez, J. Sune, R. Romero-Zaliz, F. Jimenez-Molinos, S. Long, F. Gomez-Campos and M. Liu, J. Appl. Phys., 2018, 123, 8.
57. S. Dirkmann, J. Kaiser, C. Wenger and T. Mussenbrock, ACS Appl. Mater. Interfaces, 2018, 10, 14857–14868.
58. S. Aldana, P. Garcia-Fernandez, A. Rodriguez-Fernandez, R. Romero-Zaliz, M. B. Gonzalez, F. Jimenez-Molinos, F. Campabadal, F. Gomez-Campos and J. B. Roldan, J. Phys. D: Appl. Phys., 2017, 50, 10.
59. S. Menzel, P. Kaupmann and R. Waser, Nanoscale, 2015, 7, 12673–12681.
60. A. F. Voter, Erice, ITALY, 2004.
61. D. Le, T. B. Rawal and T. S. Rahman, J. Phys. Chem. C, 2014, 118, 5346–5351.
62. J. H. Hong, Z. X. Hu, M. Probert, K. Li, D. H. Lv, X. N. Yang, L. Gu, N. N. Mao, Q. L. Feng, L. M. Xie, J. Zhang, D. Z. Wu, Z. Y. Zhang, C. H. Jin, W. Ji, X. X. Zhang, J. Yuan and Z. Zhang, Nat. Commun., 2015, 6, 8.
63. A. Yoshimura, M. Lamparski, N. Kharche and V. Meunier, Nanoscale, 2018, 10, 2388–2397.
64. A. Singh, S. Lee, H. Bae, J. Koo, L. Yang and H. Lee, RSC Adv., 2019, 9, 40309–40315.
65. Y. S. Chen, H. Y. Lee, P. S. Chen, T. Y. Wu, C. C. Wang, P. J. Tzeng, F. Chen, M. J. Tsai and C. Lien, IEEE Electron Device Lett., 2010, 31, 1473–1475.
66. F. Kreupl , Carbon Memory Assessment, white paper for the ITRS meeting on emerging research devices (ERD), August 25–26, 2014.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2nr06810a

---