Theory of nonvolatile resistive switching in monolayer molybdenum disulfide with passive electrodes

Resistive-memory devices promise to revolutionize modern computer architecture eliminating the data-shuttling bottleneck between the memory and processing unit. Recent years have seen a surge of experimental demonstrations of such devices built upon two-dimensional materials based metal–insulator–metal structures. However, the fundamental mechanism of nonvolatile resistive switching has remained elusive. Here, we conduct reactive molecular dynamics simulations for a sulfur vacancy inhabited monolayer molybdenum disulfide-based device with inert electrode systems to gain insight into such phenomena. We observe that with the application of a suitable electric field, at the vacancy positions, the sulfur atom from the other plane pops and gets arrested in the plane of the molybdenum atoms. Rigorous first principles based calculations surprisingly reveal localized metallic states (virtual filament) and stronger chemical bonding for this new atomic arrangement, explaining the nonvolatile resistive switching. We further observe that localized Joule heating plays a crucial role in restoring the popped sulfur atom to its original position. The proposed theory, which delineates both unipolar and bipolar switching, may provide useful guidelines for designing high-performance resistive-memory-based computing architecture.


INTRODUCTION
Two terminal nonvolatile resistive-memory (NVRM) devices, also known as memristors, are being explored extensively for inmemory computation to overcome the data-shuttling bottleneck between the memory and processing unit in von-Neumann architecture-based modern computers. NVRM devices are conventionally built upon metal-insulator-metal (MIM) structures, usually use TiO 2 as the insulating oxide, and work on the principle of filament formations 1 . With the emergence of various twodimensional (2D) materials, the recent past has seen a surge of experimental demonstrations of memristors, where atomically thin layers of transition metal dichalcogenides (TMDs) [2][3][4][5] or hBN 6,7 are used as insulators. These 2D MIM structures promise to overcome the vertical scaling obstacle of TiO 2 -based devices and offer highly dense, fast, and ultra-low-power technology solutions [8][9][10] for various in-memory computing and communication applications.
Despite their technological importance, the underlying physics behind the ultra-fast resistive switching (RS) in monolayers is not well understood yet, even for an extensively studied material such as 2H-MoS 2 . Although there are proposed mechanisms such as Au/Ag absorption in S vacancy sites 11 , RS in 2H-MoS 2 has also been observed with inert graphene electrodes 3 . The strong crystalline nature of MoS 2 and the ultra-fast switching speed suggest that the traditional process of electrochemical formation and disruption of filament 12 in oxide-based memristors is unlikely to occur in these devices. Even different sample preparation processes did not affect the RS 3 , hinting some intrinsic property of MoS 2 could be the origin of NVRS. Furthermore, a universal theory applicable to both bipolar and unipolar switching devices is still distant.
In this article, we report an atomistic description of RS in vacancy inhabited monolayer MoS 2 with inert electrodes. From reactive molecular dynamics (MD) simulations, we observe that most S atoms residing directly above S vacancies vertically move to the plane of Mo atoms and get stabilized there at a critical applied electric field. Density functional theory (DFT) calculations reveal that this configuration creates localized metallic states around the vacancies (which we term as virtual filaments) and forms stronger Mo-S bonds. A suitable field with reverse or even the same polarity restores the original insulating configuration when localized Joule heating is considered, achieving reversible RS. We show that the proposed theory of virtual filament can explain almost all the observations made in recent experiments 3,[13][14][15] .

RESULTS
The interatomic forces in MD simulations of this study are computed using reactive force field (ReaxFF), which can model chemical reactions such as bond breaking and formation with an accuracy close to that of ab initio simulations with a much lower computational expense 16 . ReaxFF has already been used to simulate the filament formation process in oxide-based memristive devices 12 . Here, we use the ReaxFF potential developed for MoS 2 17 , which predicts the structural and mechanical properties, defect dynamics, and phase change of single-layer MoS 2 with excellent accuracy 17,18 .
Previous theoretical and experimental studies have established S monovacancies to be the most abundant and stable defect in MoS 2 19,20 . Experiments on MoS 2 NVRM also strongly suggest that S vacancy defects play a crucial role in the RS of monolayer MoS 2 11,14,15 . Recently, memristive action caused by phase change from H to H d and T d in 2H-MoTe 2 has been reported 4 . This is, however, extremely unlikely in the case of 2H-MoS 2 because (1) it is a well-studied material and no metastable phase such as H d has been ever reported, (2) The Mo-S bonds of MoS 2 are significantly stronger and tougher to deform compared to the Mo-Te bonds of MoTe 2 (Supplementary Table 1), and (3) the insulating 2H to metallic 1T phase transition is impeded in MoS 2 by a high energy barrier 18,21 .
The popping of the sulfur atom The device structure used for MD simulations is shown in Fig. 1a, c. In most simulations, no explicit electrodes are used due to the unavailability of suitable reactive force fields, but two Lennard-Jones (LJ) walls fitted to emulate the graphene-MoS 2 interactions are placed on top and bottom to provide the MoS 2 layer necessary mechanical support [see "Methods"]. In our extensive test simulations, at the critical field of 2.85 VÅ −1 , a remarkable event is observed at the vacancy points. Due to the force exerted by the applied field, the S atom sitting immediately above a vacancy "pops" to the plane of Mo atoms and gets arrested there (Fig. 1b, d).
To gain insight into the conductivity change because of S atom popping, we calculate the density of states (DOS) using DFT. The DOS of a 10 × 10 supercell harboring a S monovacancy (0.33% or Fig. 1 Details of "popping". a Schematic of the initial atomic configuration of MoS 2 monolayer with a single S vacancy placed between two graphene LJ walls. The purple and yellow balls, respectively, represent the Mo and S atoms. The vacancy site is highlighted with the red dashed circle. Graphene walls are represented using transparent blue sheets on the top and bottom of the MoS 2 monolayer. b Atomic configuration of the system after applying an electric field of 2.85 VÅ −1 for 1 ns. Zoomed-in image of the side view of atoms around the vacancy c before and d after the application of the E-field. The vertical arrow indicates the direction of the electric field (E). DFT calculated density of states of 10 × 10 supercell model of MoS 2 monolayer with a S monovacancy e before and f after the S "popping." g Energyprojected charge density isosurface plot for states very close to the Fermi level. h Cumulative PDOS of Mo and S atoms that fall inside the circular regions of different radii around the popped S atom. ~10 13 cm −2 vacancy concentration) shows a bandgap of about 1.67 eV with zero states at Fermi level and a small defect induced state near the conduction band (Fig. 1e). Surprisingly, S atom popping causes the appearance of new states at and around the Fermi level along with two more defect states in the bandgap region (Fig. 1f), alluding metallic nature of the popped configuration. In Fig. 1g, the energy-projected charge density for states at and around the Fermi level is depicted. In addition to this, projected DOS (PDOS) for all atoms residing in circles with different radii centered around the popped atom are shown (Fig.  1h). Clearly, the major contributions for the metallic states originate from the atoms situated within a 5 Å-radius circle (0.785 nm 2 area) around the popped atom. We term this region as a "virtual filament" since it functionally mimics the conducting filament of oxide-based memristors 1,12 . It is also interesting to note that despite the S atom occupying an interstitial position after popping, the in-plane lattice parameters expand minimally (around 0.26%), making the effect of strain in the electronic structure negligible.
Both DFT and ReaxFF predict the single S popped freestanding MoS 2 structure to be less stable than the parent S vacancy configuration by more than 1.4 eV. Although in freestanding form and at 0 K this energy difference is quite sizeable, things change drastically when the LJ walls, room temperature (300 K), and an electric field are introduced. From a ReaxFF MD run at 300 K with a 2.85 VÅ −1 applied field, the same structure becomes around 0.86 eV more stable after the S popping in terms of average potential energy (excluding the energy gained from the field). When these two (vacancy and popped) structures' snapshots are run at 300 K with no applied field, the popped structure remained more stable than its vacancy-parent structure with an average energy difference of 0.73 eV, indicating the nonvolatile nature of the transition. To gain further insights into the activation energy of this transition, we perform climbing image nudged elastic band (CINEB) 22,23 calculations between these two position-averaged structures using ReaxFF [see "Methods"]. At zero bias, the forward transition barrier from the vacancy to the popped state is around 0.95 eV, while the reverse barrier turned out to be about 1.81 eV, both being relatively high for any diffusion at 300 K. However, the forward barrier, which prohibits S popping, can be monotonically diminished with a positive electric field, as shown in Supplementary Fig. 1a. At a very high field of 3.5 VÅ −1 , the forward barrier becomes only 0.3 eV, whereas at 2.5 VÅ −1 , close to the threshold field for popping, the barrier height is 0.6 eV. To analyze the chemical properties of this popping phenomenon and to test the reliability of the predictions made by ReaxFF, we perform further DFT-based analysis, the result of which are summarized in Table 1. It can be observed that ReaxFF successfully reproduces the DFT predicted trend for charges and bond orders. Although the ReaxFF charge is significantly underestimated compared to the Bader 24 charge, DDEC charge and bond order 25,26 predictions are quite close. Integrated crystal orbital overlap population (ICOOP) values quantify the bond order of specific bonds, whereas integrated crystal orbital Hamilton population (ICOHP) 27,28 values provide an estimation of bond strengths (more negative is stronger). It turns out that after popping, the Mo-S bonds become significantly more robust than the normal Mo-S or vacancy-S-neighboring-Mo bonds, as depicted by DDEC and ReaxFF bond order and ICOHP bond strength analyses. This also reaffirms the relatively higher stability of the popped configuration. Supplementary Fig. 2 illustrates the charge density difference plot of the popped structure respective to the parent-vacancy configuration. This visually confirms that the popped state is indeed an acceptor-type defect that introduces localized metallic states in the system.

Bipolar switching
Because of much stronger bonds and an extremely high reverse energy barrier, significant force is required to bring the popped S atom back to its initial location. Our extensive test simulations reveal that the forces exerted by even high reverse E-fields are not sufficient to bring back the atom. This is consistent with the results of the CINEB calculations with negative fields shown in Supplementary Fig. 1b. The monotonic decrement of the reverse barrier is evident with the increasing absolute value of the reverse field. However, unlike the forward barrier, the electric field's effect on the reverse barrier is hardly relevant compared to the zero-bias height. Nevertheless, a combination of a suitable negative field and a high local temperature, which models the Joule heating due to high current flow in the conductive regions, restores the initial atomic arrangement. The dominant role of Joule heating in the rupture of conductive filaments in traditional oxide-based RS devices is already well established [29][30][31] . In 2D memristive devices, the heating effect due to high current in low-resistance state (LRS) was also observed 32 . Recently, Akinwande et al. reported thermal optimization setup to minimize the temperature rise in hBN switch 9 , and via personal communication, we came to know that significant heating was also observed in their MoS 2 -based devices. Therefore, we may presume Joule heating to be a critical factor in the reset process of 2D NVRMs. In MD, the Joule heating effect is realized by raising the temperature of all atoms situated within a 5 Å radius around the popped atom.
To provide a complete atomistic narrative of the switching mechanism, we have carried out long-time scale MD simulations changing the electric field sequentially, as shown in Fig. 2 3) possesses 12 scattered S monovacancies on the bottom S plane (concentration is about 0.08% or~2 × 10 12 cm −2 , in the same order as mentioned in ref. 3 .). After initial heating for 1 ns and a The order of application of the field is 0→3→0→−2.5→0. Atomistic snapshots zooming-in at the vacancy locations at different points of time displaying various events, c initial vacancy location (4 ns), S atom popping after set (10 ns), popped atom returning after reset (15 ns), d secondary popping (32 ns) during reset and its recovery (35 ns) with the positive field of 2 VÅ −1 in next cycle, and e localized distortion during reset (32 ns) and its recovery (36 ns). Black, red, and green arrows highlight the vacancy locations, secondary popped atoms, and distorted locations.
S. Mitra et al. T = 300 K run for another 1 ns, slowly varying electric field pulses, each 1 ns wide, are applied in sequence. The popping events initiate at the positive field of 2.75 VÅ −1 , where 3 S atoms out of 12 vacancies pop into the Mo plane. As the field is increased to 2.85 VÅ −1 , the pop count increases to 10, while a single S atom moves directly into the vacancy spot. Further increasing the field to 3 VÅ −1 pops the lone remaining atom, elevating the pop count to 11. In Fig. 2a, the variations of the electric field and the number of popped atoms with simulation time are plotted.
The formation of small conductive paths around the vacancy sites should lead to a significant increment of the local temperature owing to Joule heating. Usually, the temperature change is associated with electrical power 31 . As we lack any knowledge of electrical power in atomistic simulations, the temperature is increased in proportion to the applied field. In experimental measurements, a compliance current is set to protect the device from permanent breakdown during the set process, while in the reset process, the current compliance is removed 3,15 . Our test simulations also reveal that increasing the local temperature to a very high value during the set process distorts and damages the structure irreversibly. Thus, we have assigned the temperature of conducting filaments to a well-tested empirical value of 790 K at field 3 VÅ −1 to eliminate the possibility of structural distortion. Therefore, the rate of temperature change in the set process becomes about 164 KV −1 Å. We have considered 19 atoms in the 5 Å-radius surrounding the vacancy centers as the filament zone, and immediately after popping, the local temperature of these zones is increased to assigned temperatures. Even when the field is reduced to zero, the number of popped atoms remains the same, demonstrating a nonvolatile memory feature.
To reset the device, the electric field is applied in the reverse direction, and the temperature of the virtual filaments is simultaneously increased. Again, after performing rigorous test simulations, to ensure a maximum return in the relatively short simulation time without causing any structural distortion, we have empirically adopted the maximum values of the reverse field and corresponding local temperature to be −2.5 VÅ −1 and 1190 K, respectively, which corresponds to a rate of 355 KV −1 Å. Similar settings have been followed previously to successfully simulate carbon-based RRAM devices, where >2200 K localized temperature is applied to initiate a melting-induced conducting-pathway rapture 33 . The variation of the local temperature with E-field is depicted in Supplementary Fig. 4a, b. The reset event starts at −1.75 VÅ −1 . Further increment of the reverse field to −2.25 and −2.5 VÅ −1 with suitably raised local temperature results in a successful return of all 11 popped atoms to their initial locations. The local heating does not significantly alter the overall temperature of the sample, raising it only around 5-7 K. In the three snapshots of Fig. 2c, the initial S monovacancy, popping of S atom during SET, and return during reset are shown. After the maximum reset field, the total pop count drops to 0 from 11, indicating resistance switching. The electric field is further reduced to zero, and the pop count remains the same. The hysteresis plot is shown in Fig. 2b. The videos of set and reset cycles are included in the Supplementary Materials.
To assess the consistency of various processes involved, we have performed another complete set and reset cycle using the same input parameters as the first cycle. After set, the total pop count reaches 11. In the second cycle, some anomalous events, such as secondary pops and a local distortion, occur during the reset process. Secondary pop connotes the popping of a neighboring S atom of the primary popped atom into the Mo plane. In the first snapshot of Fig. 2d, the secondary pop during the reset state is shown. As evident from the figure, the primary popped atom returns to its initial location while another neighboring S atom, not associating with any vacancy, pops into the Mo plane. DFT calculations reveal that the secondary popped atoms are functionally equivalent to the primary popped atoms ( Supplementary Figs 5-7). Hence, the secondary popped atoms are included in the total pop count, and both are treated equally in terms of local temperature. In the second cycle after reset, the total pop count drops to 3. Local distortion during reset is depicted in the first snapshot of Fig. 2e. To examine the effect of secondary pops and local distortions in the consequent cycles, we further sequentially apply shorter (0.5 ns each) positive pulses of 2 and 2.75 VÅ −1 . Interestingly, the secondary pops and local distortions completely recover under the forward fields (Fig. 2d, e).
The drop in popped atoms count for the first cycle is 100%, while in the second cycle it is about 73%. We have simulated another sample with 24 vacancies (0.16% or~4 × 10 12 cm −2 concentration), where the pop count has reduced about 82% after reset ( Supplementary Fig. 8a, b). As resistivity is directly related to the pop count, the presence of substantial vacancies in realistic samples should lead to a significant change in its resistivity.

Unipolar switching
Next, we focus on the mechanism of unipolar switching. This requires the reset field to be applied in the same direction as set. In our MD simulations, for unipolar operations, two types of events are observed. During the set process, the S atoms from the upper S plane pop into the Mo plane, but during the reset process, some S atoms return to their initial location, i.e., the upper S plane (Fig. 3c), and some move into the initial vacancy location, i.e., the lower S plane (Fig. 3d). This occurrence will obviously reset the device, but it would significantly limit the cyclic endurance as the atoms that have moved to the former vacancy locations cannot be popped again with the same polarity of the field. We observed that only with a very high local temperature rate, the number of atoms returning to their initial position becomes sustainable in the long run. In a recent experimental study 15 , it was mentioned that high current and power density is required to reset the device for unipolar switching, indicating Joule heating plays a governing role in retrieving the high-resistance state (HRS). While the bipolar switching is a combined effect of electric field and Joule heating, the unipolar reset mechanism seems only to be dominated by the local heating effect. The high current demonstrated with very low voltage during the unipolar reset process 3 also hints toward this hypothesis. To simulate the experimental observation, we have applied a minimal positive field of 0.02 VÅ −1 and have increased the local temperature of the virtual filament zones to 1300 K. We again simulate two consecutive cycles (Fig. 3a, b) and observe that 80% of total popped atoms have returned to their initial positions in both cases.

Device with explicit graphene electrodes
We extend this study to explore the effects of real graphene electrodes in defect inhabited MoS 2 structure under an external field. For these systems, the C-Mo, C-S, and C-C interactions are described using a reactive force field. Note that the full ReaxFF treatment with explicit graphene electrodes is computationally extremely expensive. Thus, we have only been able to simulate a somewhat "toy system" with a much smaller number of active MoS 2 atoms and a much shorter time scale [see "Methods"]. However, these simulations are essential to verify the graphene electrodes' passive nature and test the consistency of other processes observed in the LJ-walls based systems. Two separate graphene-MoS 2 interfaces are formed with a mean strain of 1.95% and 0.35%, respectively [before relaxation, see "Methods"]. On one side of the MoS 2 flake, a single freestanding graphene layer is placed, whereas on the other side, bilayer graphene is used, and the atoms of the outermost layer are kept fixed to emulate the bulk electrode behavior. This configuration emulates an experimental setup where a bulk substrate on one side supports the device. In addition to the MoS 2 flake, the electric field is applied to the adjacent graphene layers. Like the LJ walls-based system, S atom popping is also observed in the explicit graphene-MoS 2 system. However, the threshold electric field for popping turned out to be strongly dependent on the strain of the interface. In Fig. 4a, we have plotted the number of popped atoms with the magnitude of electric field pulses applied to the graphene-MoS 2 -graphene system with high (mean strain 1.95%) and low (mean strain 0.35 %) interfacial strains. In the case of a highly strained structure, the popping event occurs even at a low field of 1 VÅ −1 , whereas for the low-strain structure, a high field of about 3 VÅ −1 is required to pop S atoms. As evident from the figure, out of ten vacancies, a maximum of six atoms gets popped at electric field pulses of 1.5 and 2 VÅ −1 for the high-strain interface. Gradual increment of the applied field with long-time pulses, as done for the implicit walls based system, would help in enhancing the pop count and reaching its upper limit. In Fig. 4b, we have shown the applied electric field pulses with simulation time and corresponding zoomed-in side view of highly strained (1.95%) graphene-MoS 2 -graphene structure near a S vacancy location. As significant popping occurs at 1.5 VÅ −1 , we have applied an electric field pulse of 1.5 VÅ −1 to set the device and a reverse pulse of −1.5 VÅ −1 to reset it, both being 1 ns wide. During reset, the local temperature is increased to 1200 K. We also performed CINEB calculations for both the low-and high-strain systems to determine the forward and reverse barrier heights for S popping. At zero bias, the forward (reverse) barrier is 1.03 eV (1.59 eV) for the high-strain system (Fig. 4c) and 1.2 eV (1.44 eV) for the low-strain system (Fig. 4d). This is consistent with the variation of the critical electric field for popping with strain (Fig. 4a). We observe a local minimum just after the barrier, which probably occurs because the terminal points are position-averaged structures (before and after popping). However, to be consistent, we have calculated the barrier heights only with respect to the terminal points. Consistent monotonic decrement of the forward barrier with increasing electric field is also observed for this case (Fig. 4c, d). We note that while the forward barrier for the implicit graphene walls MoS 2 system (0.86 eV) is lower than that of the explicit graphene-MoS 2 system (1.03 eV), the critical field for popping is much greater (2.75 and 1 VÅ −1 ) for the former case. There is a general trend of overestimation of the critical field when the electrode-MoS 2 interaction is modeled through a LJ potential instead of a full reactive potential, which is reflected in our other simulations as well [see Discussion and the SI]. The reason for this is not clear yet. However, the intrinsic S atom popping occurs for both systems regardless of the electrode potential used, which qualitatively confirms this mechanism's consistency. Although the graphene implicit walls MoS 2 system overestimates the critical electric field for set and reset, it is an inexpensive system computationally and is suitable to test the cyclic performance of the device, albeit somewhat qualitatively.

DISCUSSION
Compared to the experimental findings 3,11,14 , the applied electric field values might appear to be high in our simulations. One possible reason for this is the underestimation of charge in ReaxFF, as the force experienced by the ions is the product of their charge and the applied electric field. In a neutral state, the formal charge of Mo/S is +4/−2, while in MoS 2 , DFT-based Bader analysis yields this value to be +1.08/−0.54. Even a much higher Bader charge of +1.72/−0.86 has been reported using different DFT settings 34 . However, in our MD simulations, the ionic charge is about +0.4/−0.2, requiring much larger fields to pop the S atoms. Also, large field values help to accelerate the dynamics and achieve rare events in a short amount of time.
Moreover, the simulation using real electrode materials show the threshold electric field varies significantly with interfacial strain. Considering the underestimation of charges, it is found that the threshold electric field for the graphene-MoS 2 system with a higher interfacial strain is comparable with experimental findings. Since in MD simulation, the threshold value of force required to pop an atom is a combined effect of the external field, atomic charge, interfacial strain, and interatomic forces calculated using different potential functions, it is indeed challenging to estimate the exact terminal voltages from the electric field used in the simulations and corroborate it with experimental findings.
There is a significant difference between the rate of local temperature change for set (r set = 164 KV −1 Å) and reset (r reset = 355 KV −1 Å) cycles in the simulations. The precise estimation of local temperature change in such an atomically thin filament is abstruse. Therefore, we simply assign certain test values of localized temperature to identify the consequences of Joule heating in our simulations. It is found that for much lower r reset , S atom returns also occur, but the process becomes extremely slow ( Supplementary Fig. 9). Hence, to accelerate the reset switching, such a high-temperature rate was applied.
The ratio of the number of popped atoms in LRS and HRS obtained from our simulations is much lower than the experimentally found high on-off current ratios. It should be noted that due to computational limitations, it was only possible to conduct MD simulations on MoS 2 flakes with a maximum area of around 500 nm 2 , whereas the experimental studies have used flakes as large as 4 µm 2 . Therefore, in a computational study, it is not possible to attain a reliable statistical distribution for the number of popped atoms during HRS and LRS over the cycles. However, in an equalsized sample with a higher number of vacancies, we do observe the #popped HRS /#popped LRS to go up significantly (Supplementary Fig.  8). It seems that the MD simulations do capture the essential physics of nonvolatile RS, but with a low degree of freedom. There are a total of ten vacancies in both systems. b Plot of the applied electric field with simulation time and corresponding zoomed-in side-view images near a vacancy location of graphene(monolayer)-MoS 2 -graphene(bilayer) system. The mean interfacial strain between graphene and MoS 2 is 1.95%. The outermost graphene layer has been kept fixed during the simulation. A positive field implies a field applied from the top to the bottom electrode, and a negative field signifies a field applied from bottom to top. For popping, an electric field pulse of 1.5 VÅ −1 is applied, and for return, −1.5 VÅ −1 pulse along with 1200 K local temperature is applied. Electric field-dependent activation energy barriers for S atom popping (forward reaction) for the graphene(monolayer)-MoS 2 -graphene(bilayer) c high-strain and d low-strain systems. Note that the energy in the y-axis is relative and only represents the barrier height and relative total energy difference between the two terminal points. Therefore, only the barriers under different electric fields can be compared from these plots, not the absolute values of total energies. The CINEB calculations have been done considering the potential energy contributions from the electric fields.
Apart from the current ratio at on and off states, reported experimental observations corroborate the simulation results well. As evident from the simulation, only the localized regions around popped atoms become conductive after switching. This is supported by experimental findings 15 , where it was observed that during LRS, most of the sample remained resistive while certain highly localized regions with enhanced conductivity were formed. Multilevel switching can also be explained from this mechanism. Multiple resistance states 13,15 may occur due to discrete changes in the number of popped and returned atoms during the set and reset switching, respectively ( Supplementary Fig. 10). It was reported that some devices failed after achieving LRS and could not be reset again 14 . During our simulations, in some cases, multiple unrecoverable secondary pops are generated, causing that particular virtual filament to remain permanently in LRS. Experimental works further report that annealing in a sulfur-rich atmosphere recovers the original HRS for these samples. We have also reproduced this observation, where the extra S atoms push themselves into the secondary popped regions at elevated temperatures, triggering the popped atoms to return to their initial locations ( Supplementary Fig. 11).
In the literature, it is heavily postulated that Au or other electrode atom absorption in S vacancy sites 11 could be a possible mechanism for the RS. Although Au is thought to be inert, migration of Au nanoparticles in CBRAM-based devices has been shown 35 , and the tendency of adsorption of Au at the defects and the edges of MoS 2 have also been reported 36 . Our DFT calculations also confirm induced states at the Fermi level upon Au/Ag adsorption in single/double S vacancy of MoS 2 ( Supplementary Figs 5-7). However, our COHP analyses reveal that there is a high amount of antibonding states at the Fermi level for Au/Ag atom adsorption in S vacancies, while S popping does not generate any such antibonding state, implying that the S popped configuration is chemically much more stable 28 ( Supplementary Fig. 12). Also, the process of S atom popping is independent of the electrode material and can reasonably explain both bipolar and unipolar operations. Although our study extensively probes the devices with inert electrodes, we have also tried to simulate the behavior of devices with implicit and explicit Au electrodes, whose atoms are expected to migrate from the electrodes to the S vacancies. However, the Au-MoS 2 reactive force field has not been developed yet, and thus we are forced to use a Coulomb corrected LJ force field to describe the Au-Au and Au-MoS 2 interactions. The result from these simulations, which ultimately turned out to be inconclusive about the possibility of Au diffusion, is documented in Supplementary Note 1 and Supplementary Fig. 13. Here, S atom popping at the vacancy sites is observed as expected, but it is difficult to comment on the possibility of metal atom migration from electrodes because of the inadequate LJ force field used, which can only describe the nonbonding dispersion and Coulomb interactions, whereas our first principles calculations reveal a strong reactive nature of Au when it diffuses into the vacancy ( Supplementary Figs 5-7). The implicit Au electrodebased system reproduces all events observed in implicit/explicit graphene electrode-based systems ( Supplementary Figs 14-17), but with slightly different values of critical fields (Supplemen- tary Figs 15 and 16). However, this system is clearly unable to explore the possibility of the Au atom migration. Therefore, except for chemically entirely inert electrodes like graphene, the metal ion absorption process may coexist with S atom popping to trigger switching.
There are two important physical aspects not included in our work, which warrant future explorations. The first one is obviously the exclusion of chemically active electrodes such as explicit Au, whose atoms are reported to be diffusing into the S vacancies of MoS 2 to create conducting regions 11,35,36 . Definitely, the development of a rigorous reactive force field is essential to explore this phenomenon through simulations. The second one is our lack of knowledge of the current-voltage-heat relationship at the atomistic level, which prompted us to use empirical temperature values. The large size of the system required to simulate the experimental setup prohibits us from using methodologies like nonequilibrium Green's function (NEGF) to find out the transmission spectra or even the I-V relationship for the various states of the system. Also, convergence issues are well known for ballistic NEGF calculation under high-bias conditions 37 . A drastic improvement in methodologies is required to model the quantum transport phenomenon of these atomically thin devices. Furthermore, a suitable model of the current-voltage-heat relation must be developed for these systems to model the variation of the filament temperatures with varying current and voltage.
In conclusion, we propose a theory of "virtual filament" to explain the nonvolatile RS of monolayer MoS 2 -based memristor. Since this process occurs in a monolayer with a minimal movement of a S atom, it can be exceptionally fast, as observed in experiments 3 . Despite revealing minor recoverable nonidealities, the process is extremely reversible and should result in good cyclic performance. Although we have explicitly investigated 2H-MoS 2 here, the proposed theory might also be applicable for other 2H or even 1T/1T′ TMDs resulting in ultra-fast switching 38 . As the process is electrode independent, a library of TMD materials can be identified with proper computational settings that would show switching performance similar to or even better than MoS 2 . The insights provided in our work may help engineers to design reliable and high-performance 2D memristor-based in-memory computing architecture.

MD simulations setup
The MD simulations of the present study are performed using the LAMMPS 39 package, and the packages OVITO 40 and VESTA 41 are used for visualizations. Especially, the KOKKOS 42 package of LAMMPS has been used extensively to enable GPU acceleration. The interatomic forces are mostly computed using reactive force field (ReaxFF) 16,43 , which can model chemical reactions such as bond breaking and formation with an accuracy close to that of ab initio simulations with much lower computational expense 16 .
For time integration of Newton's equations of motion, the velocity Verlet algorithm with a time step of 0.5 fs is used. Extensive tests are performed to choose the proper time step. The 0.5 fs step reproduces all events that occurred during runs with much smaller steps (0.1-0.2 fs) without any energy drift. The simulations are performed at constant ambient pressure using the Nose-Hoover barostat (NPH ensemble). The temperatures of various zones (the high-resistance bulk and conducting filament zones) are controlled using separate Berendsen thermostats. Periodic boundary conditions are applied in the two in-plane directions, while a >45 Å thick vacuum region is inserted in the out-of-plane direction to avoid spurious interactions between two periodic images. Scattered S monovacancies of different concentrations are created in the top S plane of MoS 2 flakes (Supplementary Fig. 3a). For most simulations, to emulate the C-MoS 2 interactions, two flat walls 3.42 Å away from the MoS 2 surfaces are placed on the top and bottom of the flake. The energy of the wall-MoS 2 interaction is described using LJ 9-3 potential. The binding energies of the graphene-MoS 2 interface with varying interlayer distances are first calculated using DFT [detailed below]. The parameter sets for LJ 9-3 potential are obtained by performing a fitting procedure to this curve of binding energy vs. interlayer spacing (Supplementary Fig. 3b). The employed energy and distance parameters in the potential are ε = 0.1646 eV, σ = 4 Å, and r c = 10.0 Å. After energy minimizations, the samples are slowly heated from T = 10 to 300 K during 1 ns runs. Right after an event, i.e., a significant change of charge of an S atom indicating popping/return, the temperature changes in the filament zones are induced linearly over a 10 ps time period. The external electric field is applied in the out-of-plane direction. As the lower z plane of MoS 2 contains S vacancies, the set process is initiated by applying the field toward positive z direction, while bipolar reset is initiated S. Mitra et al.
using the reverse field. To simulate annealing of the secondary popped S inhabited sample, a one-sided wall is used, and the S atoms are fed from the open side.

Interface design
For the explicit graphene-MoS 2 interface, the equilibrium distance between MoS 2 -graphene and graphene-graphene turned out to be about 3.3 and 3.38 Å, respectively, which matches previous reports 44,45 . For the high-strain graphene-MoS 2 system, the angle between the cell vectors and the amount of rotation between the two surfaces is about 120°and 161°, respectively. The strain is applied to graphene. The individual strain components are ε 11 ¼ À2:93%, ε 22 ¼ À2:93% and mean strain = 1.95%. Small device dimensions are maintained to expedite the simulations with explicit electrodes. Here, it is about 3.8 × 5.4 nm 2 with a total of 3230 atoms, including 720 MoS 2 atoms.
For the low-strain graphene-MoS 2 system, the angle between the vectors and the rotation amount between two surfaces is about 60°and 165°, respectively. The strain is applied to graphene. The individual strain components are ε 11 ¼ 0:53%, ε 22 ¼ 0:53% and mean strain = 0.35%. The device dimension is about 4.1 × 6.0 nm 2 , with a total of 3645 atoms, including 865 MoS 2 atoms. In both the structures, ten S monovacancies are created.
For explicit gold-MoS 2 system, the Au-MoS 2 interaction is described by a Coulomb interaction included LJ 12-6 potential with ε ¼ 22 meV and σ ¼ 0:29 nm, while LJ 12-6 parameter for Au-Au interaction are ε ¼ 412:9 meV and σ ¼ 0:26 nm 46,47 . Au-MoS 2 interlayer equilibrium distance turned out to be about 3.2 Å. The assigned charges of Au atoms are obtained from DFT-based Bader/DDEC data. The device dimension with a higher strained interface is 3.3 × 7.5 nm 2 , including 756 MoS 2 atoms.

DFT calculations
DFT calculations are carried out using the generalized gradient approximation as implemented in the code VASP 48-51 with the PAW 52 method using the Perdew−Burke−Ernzenhof 53 exchange-correlation functional. The following electrons are treated as the valence electrons and are expanded in the plane-wave basis set; Mo: 4d 5 5s 1 , S: 3s 2 3p 4 , Au: 5d 10 6s 1 , Ag: 4d 10 5s 1 . A sufficiently large cutoff energy of 520 eV is used to avoid any Pulay stress. For all structural relaxations, a > 30 1 similar k-mesh is used for all static runs. The BZ is integrated with Gaussian smearing using a well-tested smearing factor of 0.05 eV. Electronic convergence is set to be attained when the difference in energy of successive electronic steps becomes less than 10 −6 eV, whereas the structural geometry is optimized until the maximum Hellmann-Feynman force on every atom falls below 0.05 eVÅ −1 . A large vacuum space of ≥25 Å in the direction of c is applied to avoid any spurious interaction between periodically repeated layers. Considering the high computational cost of DFT, we build a small model of the MD structure with a 10 × 10 supercell of MoS 2 (300 atoms) and a single S mono/divacancy (Fig.  1a). Reducing the cell size for DFT increases the vacancy concentration to about ten times of that of the MD structure. However, the distance of 31.2 Å between two periodic images of a defect is large enough to eliminate any spurious interactions.
For the studies requiring comparison between DFT and ReaxFF, such as formation and activation energy of the popped state and various charge and bond order parameters, the same 10 × 10 supercell of MoS 2 with a single S vacancy used for the DFT studies have been used for ReaxFFbased calculations as well, albeit with necessary implicit or explicit electrodes. The ReaxFF-based static runs, energy minimizations, and CINEB calculations 22,23,54,55 are also performed using LAMMPS. After a gradual heating and equilibration run, the above-mentioned 10 × 10 supercell with a single vacancy and implicit/explicit electrodes is exposed to 300 K temperature, 1 atm pressure, and the threshold electric field for popping. To calculate the average potential energy change due to popping, the energies of 250 snapshots (each snapshot 500 fs apart) before and after popping are averaged and compared. For the CINEB calculations, the average atomic position of the same snapshots is used as the terminal points, while five images between them are used to determine the transition barrier.
The Bader charge analyses are performed using the code developed by Henkelman et al. 24,[56][57][58] , where charge densities generated from DFT static runs are used as inputs. In addition to this, DDEC 25,26,59-62 charge partitioning and bond order analyses, as implemented in the program Chargemol, are also employed on the same charge densities. COHP analyses have been performed using the tool LOBSTER 27,63-66 , where wavefunctions generated by the DFT static runs are used as the inputs. Despite a large number of atoms in the structure, the charge spilling for all COHP projections is less than 1.45%, establishing the reliability of the obtained results.

DATA AVAILABILITY
The authors declare that the main data supporting the findings of this study are available within the paper and its Supplementary Files. Additional supporting videos can be found in the YouTube playlist: https://www.youtube.com/ playlist?list=PLhnED2_HRtqi1748hhfJENVFX_9634o89. Moreover, all highdefinition videos are made available for download at: https://osf.io/ nd2es/?view_only=7beb0d98fa3145c38609735406ade8e3. Other relevant data are available from the corresponding author upon reasonable request.