Atomistic dynamics of sulfur-deficient high-symmetry grain boundaries in molybdenum disulfide†

As a common type of structural defect, grain boundaries (GBs) play an important role in tailoring the physical and chemical properties of bulk crystals and their two-dimensional (2D) counterparts such as graphene and molybdenum disulfide (MoS2). In this study, we explore the atomic structures and dynamics of three kinds of high-symmetry GBs (α, β and γ) in monolayer MoS2. Atomic-resolution transmission electron microscopy (TEM) is used to characterize their formation and evolutionary dynamics, and atomistic simulation based analysis explains the size distribution of α-type GBs observed under TEM and the inter-GB interaction, revealing the stabilization mechanism of GBs by pre-existing sulfur vacancies. The results elucidate the correlation between the observed GB dynamics and the migration of sulfur atoms across GBs via a vacancy-mediated mechanism, offering a new perspective for GB engineering in monolayer MoS2, which may be generalized to other transition metal dichalcogenides.


Introduction
Grain boundaries (GBs) in general, including phase boundaries (PBs), [1][2][3][4][5][6][7] are local structural imperfections in crystalline materials. They can profoundly influence the mechanical, electronic, optical, magnetic and chemical properties of materials. [8][9][10][11][12][13] Precisely imaging the structures and/or the dynamics of GBs is essential to reveal the structure-property correlation in materials, including the development of two-dimensional (2D) materials such as graphene and MoS 2 . Aberration corrected scanning transmission electron microscopy (AC-STEM) 14 enables direct atomic-resolution imaging of structural defects (including GBs) in 2D materials such as graphene, transition metal dichalcogenides (TMDs), 15 and their atomistic dynamics [16][17][18] under electron beam irradiation and/or at elevated temperatures. In monolayer MoS 2 , Lin et al. 19 reported the nucleation and growth of GBs (identified as α-type GB in this work) around Re/Au dopant centers, as activated by the electron beam irradiation at elevated temperatures. The energetics and kinetics of the relevant 2H-1T phase transition were recently elucidated by Zhao and Ding based on first-principles calculations. 20 In monolayer MoSe 2 , several groups have confirmed the formation of mirror twin-boundaries [21][22][23] (also termed β-type GBs) due to the Se deficiency, which could be induced by the electron beam irradiation or grown intrinsically during the molecular beam epitaxy (MBE) process. For the effect of GBs on the electronic properties of 2D materials, Barja et al. 24 characterized onedimensional charge density waves at the β-type GB in MoSe 2 via low-temperature scanning tunneling microscopy. Huang et al. 25 showed that GBs could cause a local modulation of band-gaps that depend on the mis-orientation between neighboring grains in the monolayer MoS 2 , introducing in-gap states that may benefit the catalytic function in hydrodesulfurization, 21,24,26 and modifying the magnetic properties. Ly et al. 27 confirmed experimentally that the electronic transport across GBs is highly sensitive to the inter-domain misorientation. All these previous results have unambiguously demonstrated that GBs have a remarkable impact on the physical and chemical properties of 2D TMDs.
Despite these previous studies on GB dynamics, the S-vacancy induced atomic sliding-migration mechanism of GB evolution in the MoS 2 system has not been reported due to the lack of statistically sufficient data on GBs and analysis of the inter-GB interaction and S-vacancy-assisted kinetic processes.
Hence, an in-depth view of the formation and evolution of GBs in TMDs is still necessary. To address this point, we present direct characterization of the atomic processes during the formation and annihilation of high-symmetry GBs in monolayer 2H-MoS 2 . We focus on these GBs because of their frequent appearance under electron beam irradiation, as well as wellresolved atomic structures that can be used to construct models for atomistic simulations. The simulation results show that the GB structures with a three-Mo-column width and a length of 6-7 Mo-units are most energetically favorable for the α-type GB, which are in good agreement with our experimental evidence and the findings reported by Lin et al. 19 Our TEM observations show that these GBs are always accompanied by the presence of S-vacancies near the GB, and the migration of S-atoms in the pristine lattice participates in the process of GB evolution. Further calculations demonstrate that an increase in the amount of S-vacancies and reduction in the vacancy-GB distance could dramatically reduce the formation energy of GBs, indicating that the S-vacancies induced by irradiation damage under the experimental conditions promote the nucleation and growth of GBs. By calculating the energy barriers involved in the mass transport of S-atoms across the GBs, we show that the S-vacancies aid the migration of S-atoms in the pristine lattice, leading to the observed GB dynamics, where the vacancy-mediated atomic sliding-migration process and the beam-atom interaction are responsible. This mechanism is different from the doping effect of Re impurities reported previously. 19 Our investigation of GBs in monolayer 2H-MoS 2 not only deepens our understanding of the energetics and dynamic evolution of structural defects in 2D materials, but also offers insights into the design of GB networks in TMDs for their catalysis function, for example in the hydrogen evolution reaction.

Sample preparation
Monolayer MoS 2 samples were prepared through standard micromechanical exfoliation of mineral-form bulk MoS 2 (SPI supplies), followed by deposition onto a silicon substrate using scotch tape. Under an optical microscope, identified monolayers were transferred onto lacey carbon TEM grids via a wet-chemistry lift-off process.

TEM characterization and image simulation
All the annular dark field (ADF)-STEM experiments were conducted with a spherical-aberration corrected TEM (FEI Titan ChemiSTEM, probe corrected) at an acceleration voltage of 80 kV. The spherical aberration C s has been optimized down to 2 μm. A probe current of ∼60 pA was chosen for ADF imaging to minimize the beam damage. The convergence angle α of the incident electron probe was set to 21 mrad, and the detection angle β of the ADF detector was set to 50-200 mrad. Experimental ADF images were processed through the well-known Wiener filtering to enhance the signal-to-noise ratio. ADF image simulation was carried out by QSTEM 27 under the same parameter settings as the experiments such as C s , α, and β besides the probe size ∼1.2 Å.

First-principles calculations
We explored the structural and electronic properties of monolayer 2H-MoS 2 with GBs identified in the experiments by performing first-principles calculations using plane-wave basis sets based density functional theory (DFT) methods. Ultrasoft (US) pseudopotentials 28 were used for the core-valence electron interaction and the Perdew-Burke-Ernzerhof (PBE) parameterization 29,30 of the generalized gradient approximation (GGA) was used for the exchange-correlation functional. We used the Vienna Ab initio Simulation Package (VASP) 31 for DFT calculations, with the atomic models illustrated in Fig. 1a-c. For all results presented, an energy cutoff of 280 eV was used for plane-wave basis sets. For the study of GBs, we considered a 10 × 1 supercell, and a vacuum layer of 20 Å was used to isolate the single layers. A 1 × 6 × 1 Monkhorst-Pack mesh grid for sampling k points was set up for the Brillouin zone integration. These settings were verified by achieving a total energy convergence of less than 1 meV per atom. Geometrical relaxation was carried out before calculating the structural properties and electronic structures, 32 until the force on atom was converged below a threshold of 0.01 eV Å −1 . The formation energy of a GB, E f = E GB − E 0 , is calculated from the relaxed structure of monolayer 2H-MoS 2 with GBs and compared to the native monolayer 2H-MoS 2 . Here, E GB and E 0 are the total energy of the GB-containing and pristine structures, respectively. The DFT calculation for α-GB has also been independently checked using CASTEP 43 using equivalent parameters, with a similar qualitative result. The potential energy barrier E b for the migration of a S atom across the GBs was calculated using the nudged elastic band (NEB) technique. 33,34

Empirical potential based calculations
DFT methods offer first-principles and an accurate description of interatomic interactions in materials, which is especially valuable for materials with complex structures and defects. However, because of their prohibitively high computational cost in the exploration of large systems with several hundreds of atoms, we used empirical potential based molecular dynamics (MD) simulations instead to study the structures for discussing the width/length distribution of GBs, and the interaction between mis-aligned GBs, where the model consists of hundreds to thousands atoms (see the ESI † for details). MD simulations are performed using the large-scale atomic/molecular massively parallel simulator (LAMMPS). 35 The Stillinger-Weber (SW) potential recently developed for 2H-MoS 2 was used in the calculations, which was validated by the calculations of its mechanical and thermal properties. 36 Periodic boundary conditions were applied to a 2D supercell of the 2H-MoS 2 . In the MD simulations, the atomic structures of GBs were firstly optimized, and then equilibrated under ambient conditions (temperature T = 300 K) under a Nosé-Hoover thermostat 37 for 50 ps. A time step of 0.5 fs was used to integrate the equations of motion. The atomic structures were further relaxed for 50 ps to evaluate the averaged potential energy for the formation energy E f of 2H-MoS 2 structures with GBs. The formation energies calculated from MD simulations are comparable to those obtained from DFT calculations.

Atomic structures and electronic properties of GBs in single-layer 2H-MoS 2
Three major types of high-symmetry GBs (labeled as α, β and γ-types) are frequently observed along the zigzag direction of a monolayer 2H-MoS 2 matrix, as illustrated in Fig. 1, where the experimental ADF images are in good agreement with the relaxed structural models obtained from DFT calculations. The cyan and green triangles in Fig. 1a mark the orientation of both sides of these GBs. Among them, the α-type GB can be regarded as a 0°GB, or named the strained T phase in the matrix of the 2H-phase as defined by Lin et al. 19 For consistency, here we have followed their terminology in naming this defect "α-type GB" where GB normally refers to the interface between grains, although no new grains are observed during the structural evolution. β-type and γ-type GBs shown in Bs, respectively. 6,22 The β-type GB is parallel to the zigzag direction (named 4|4P), while the γ-type GB is characterized by the shared zigzag edges (4|4E). It should be noted that the β-type GB is also known as the inversion domain boundary (IDB) 21 or mirror twin boundary (MTB) 23 in the literature.
DFT calculations were employed to determine the formation energies of these GBs. The results are summarized in Fig. 1d, which shows that the β-type GB has the lowest formation energy of 0.21 eV Å −1 , compared to 0.32 eV Å −1 for the α-type GB and 1.22 eV Å −1 for the γ-type GB, consistent with our experimental findings that the dominant GBs are β-type rather than other types of high-symmetry GBs in the MBE-grown MoS 2 and MoSe 2 monolayers. 21 Fig. 1e shows the electronic density-of-states (DOS) of 2H-MoS 2 with and without the GBs. The pristine monolayer 2H-MoS 2 sheet is found to be semiconducting with a band gap of 1.62 eV, aligning with other reports. [38][39][40] To characterize the spatial distribution of electrons near the Fermi level, the partial charge densities with energy levels ranging from E F −0.5 eV to E F +0.5 eV are plotted in Fig. S1, † where E F is the Fermi energy. This charge analysis shows that electronic states near the Fermi level are localized at the GB. These GB-induced states located within the native gap could lead to a semiconductor-to-metal transition at the GB, where one-dimensional (1D) metallic states are reported. 6,21 These in-gap states induced by GBs could also, in principle, benefit the catalytic function of TMDs in hydrodesulfurization or hydrogen evolution reaction where GBs behave as the active sites for catalysis.

Size statistics of α-type GBs in 2H-MoS 2
In this section, we present a quantitative analysis of the formation and structural evolution of α-type GBs in the monolayer 2H-MoS 2 . Fig. 2a-e show a time sequence of ADF images for the birth and subsequent growth of α-type GBs at room temperature. The α-type GB (marked by yellow polygons in Fig. 2b), or the strained-T-phase, is formed coherently within the lattice of the 2H-MoS 2 matrix, along the zigzag 1,10 direction. Fig. 2f shows the atomic structures of α-type GBs, where one can see significant lattice distortion near the GB. Specifically, equilateral triangles of Mo sub-lattices have turned into isosceles. To quantify the in-plane strain distribution, we employed geometric phase analysis (GPA) on the TEM images. The results are summarized in Fig. 2g, from which we conclude that a typical α-type GB leads to a shear strain ε xy up to +20% (Fig. 2h-i). The nucleation process of α-type GBs at room temperature without dopants involved here is different from the thermally activated mechanism near dopant atoms as reported in the literature. 19 From this experimental evidence and our ab initio energetic calculation, we conclude that the sulfur vacancies nearby, marked by red arrows in Fig. 2a-e (also see Fig. S2 † for the image simulation), could serve as important stimuli to the α-type GB formation.
TEM observation and statistics in the size distribution of over 60 α-type GBs (Fig. 3) show that these α-type GBs in the shape of nano-strips grow laterally only at the very early stage of GB evolution ( Fig. 2 and 3b), while the widths are all three-Mo-column wide (Fig. 3a). This specific width preference corresponds to the width for the lowest GB formation energy ΔE in our empirical potential based calculations (Fig. 3c). According to the Boltzmann distribution factor p B = exp(−ΔE/ k B T ) calculated for thermal equilibrium at T = 300 K, the α-type GBs with this optimal width has ∼10 3 and ∼10 15 higher probability to be observed than the GBs that are two-or morethan-three-Mo-columns in width (Fig. S3 †). On the other hand, the length distribution of GB shown in Fig. 3b is centralized with a peak value of 6-7 Mo-sublattice units, which could be explained as follows. The formation energy of GBs increases almost linearly with the length of the embedded α-type GBs (inset of Fig. 3d), as indicated by our empirical potential based MD simulation results. However, the lateral growth/decay of GBs is a kinetic process. The nucleation of a single Mo unit with two terminations costs an energy of 2E e = 0.021 eV per atom as the first step of α-type GB formation, where E e is the formation energy of one Mo unit and the termination. Additional cost for increasing the GB length by one Mo unit, E u = 0.00485 eV per atom, is however much lower. Consequently, the formation energy can be written as E f = N u E u + 2E e where the N u (>0) is the length measured in the Mo unit. During the kinetic evolution process, short α-type GBs with very few (e.g. two) Mo unit lengths are very unstable due to the high formation energy of termination, and thus are not easy to be captured within the experimental time window. The longer α-type GBs, however, has a finite lifetime that could be captured in the time window of TEM observation, as thermal fluctuation can activate lowenergy excitation. The abovementioned physical mechanisms can be quantified by applying a kinetic model 41 and considering the Boltzmann factor of equilibrium distribution (see details in the ESI †). From the calculation results (Fig. 3d), we conclude that the probability of α-type GBs characterized in experiments is a function of the GB length, and the most probable length of GB is N u = ∼7, which aligns with the fact that 6-7 Mo-sub-lattices are most common in the experiments.
Once several α-type GBs are nucleated, the interaction between two α-type GBs also becomes important in generating high-density GB networks in monolayer 2H-MoS 2 . From our DFT calculations, we find that the average distance between neighboring Mo-Mo atoms, d Mo-Mo , is 2.98 Å within the α-type GB, 3.18 Å in pristine matrix 2H-MoS 2 , and 3.31 Å in the GB/ matrix interface. These results suggest that the GB is under compression with a strain of ∼−6.0%, while the region between the GB and pristine 2H-MoS 2 is under tension with a strain of ∼4.1%. This indicates that two parallel α-type GBs are subjected to repulsion due to the tension strain in between when they are close to each other, while two non-parallel GBs could release the strain energy by approaching each other, where the compression and tension strain fields counteract.
Specifically, our empirical potential based calculations show that the formation energy of two parallel α-type GBs decreases with the inter-GB distance, indicating the repulsive nature of interactions, while that of two α-type GBs mis-aligned by 60°a ngle increases with the inter-GB distance, which indicates an effective attraction between them ( Fig. 3e and f ). These results explain why (from the TEM images) two parallel GBs are all very well separated while the GBs mis-aligned by 60°angle prefer to merge by joining the ends (Fig. 2e).

Atomistic dynamics of inversion domains enclosed by β-type and γ-type GBs
It has been found that βand γ-type GBs could be formed within the monolayer MoSe 2 due to the deficiency of chalcogen atoms caused by the electron beam irradiation, 22 which were also frequently observed in MoS 2 under our experimental conditions (Fig. S4 †). Interestingly, in addition to the GB growth, we also sometimes captured the inversed processthe annihilation of GBs, which has not been well documented for TMDs in the literature. Fig. 4 demonstrates that the inversed domain is enclosed by a triangle composed of one β-type GB and two γ-type GBs. We have resolved structural models corresponding to the experimental series (I-IV), highlighting the atomic processes involved. It should be noted that the MoS 2 lattice inside the domain is still in the 2H phase, but in a 180°inversion with respect to the outer matrix. These GBs are high-symmetry 60°GBs along the zigzag direction of 2H-MoS 2 . From the TEM images of the domain evolution, we find that the γ-type GBs migrate towards the domain center and lengths of both the β-type and γ-type GBs decrease, through the migration of the S-atom from the γ-type GB to the domain center. As a result, the GBs are shortened, the inversed domain shrinks and eventually reconstructed to the 2H-matrix. The calculated energies along this evolution pathway, using the empirical potential method (Fig. 4), verify this mechanism and show that the total energy decreases dramatically during the annihilation of the IDB by relieving the lattice strain. As both β-type and γ-type GBs are sulfur deficient, sulfur atoms in the 2H phase outside the triangles must participate in the beam-driven conveyance of S-atoms to heal the GBs, leaving S-vacancies in the 2H-MoS 2 lattices accompanied by the annihilation of GBs. This indicates that both experimentally observed beam-induced growth 22 and annihilation of GB-bound inversion domains are complex dynamical processes involving different kinetic pathways. The former domain growth process is initialized at the beam-induced vacancy sites (line defects) in normal 2H lattices under beam irradiation, while the latter annihilation process starts from S-deficient GBs of the inverse domains and turns into normal 2H lattices with the S supply nearby the domains. The roles played by interaction with either the chalcogen vacancy and/or mobile chalcogen adatoms under beam irradiation will be further investigated.
A common feature between α-, βand γ-type GB nanostructures is the remarkable lattice distortion (Fig. 2 and 4). Compared to the α-type GBs, strain distributions at the β-type and γ-type GBs exhibit local dipoles, indicating the presence of a dislocation pair in the atomic structures (Fig. S5 †). The lattice discontinuity and deformation revealed by our GPA analysis exactly reflect the intrinsic dislocation characteristics of β-type and γ-type GBs, aligning with the fact that GB is often regarded as assembled dislocations in the literature. 42

The effect of S-vacancies on the GB dynamics
For the structural evolution of GBs characterized in Fig. 2 and 4, we always observed that there were sulfur vacancies nearby the GBs, which may account for the observed GB dynamics. To clarify this point, DFT calculations were carried out to reveal the dependence of the S-vacancy presence on the relative formation energy, E f , of the GBs and the results are summarized in Fig. 5. For the sake of simplicity, we plot the dependence of E f on the vacancy-GB distance d with the presence of only one vacancy nearby (Fig. 5a). We find that E f is reduced as d decreases, indicating that the incorporation of S vacancy nearby stabilizes the α-type and γ-type GBs by releasing strain in the distorted lattice containing GBs. We then consider multiple S-vacancies in the 2H phase with the same distance of three Mo-units to the center of GB and the calculated values of E f behave as a function of the number of vacancies (Fig. 5b). Fig. 4 Atomic processes of the GB evolution. Experimentally observed annihilation processes (I-IV) of the inversion domain consisting of one β-type GB (the red rectangle) and two γ-type GBs (green rectangles). Scale bar: 0.5 nm. In the resolved atomic models corresponding to the experimental image series, green rectangles indicate the γ-type GBs and the red ones correspond to β-type GBs. Red arrows mark the migration path of S-atoms in the lattice, from within the γ-type boundary to the central inversion domain in the 2H phase. S-atoms in the green γ-type GBs migrate onto the hollow sites inside as indicated by the red arrows, leading to the shift of γ-type GBs. Empirical potential based MD simulation results show that the total energy decreases dramatically with the evolution of GBs. The formation energy of GBs decreases with the number of S-vacancies, indicating that a high concentration of S-vacancies can promote the formation of GBs.
As highlighted by the red arrows in Fig. 2 and 4, the S-atom migration via a vacancy-mediated mechanism occurs during the formation of α-type GBs and in the dynamics of γ-type GBs, for example. Migration of S atoms mediated via the vacancy mechanism actually dominates the atomic processes (Fig. S6 †) of α-type (Fig. 2) and γ-type (Fig. 4) GBs. To visualize this process, we highlight the slipped S atoms in green color ( Fig. 4 and S6 †), underneath which are other S-atoms (yellow) in their sub-lattices. Atomic migration and lattice collapse induced by these vacancies will then lead to α-type GB formation. Similarly, S atoms within the γ-type GBs migrate to the neighboring hollow sites along the red arrows in Fig. 4, leading to a collective shift of the γ-type GB. As the inversed domain shrinks, the migration of S-atoms compensates the S-deficiency in β-type and γ-type GBs and gradually turns the shortening β-type and γ-type GBs into the normal 2H phase.
To explore the energetics of S-vacancy migrating from α-type and γ-type GBs to the neighboring 2H phase of MoS 2 , we perform DFT calculations for the potential barriers of such processes across different types of GBs using the NEB technique. From the values of energy barriers, E b , summarized in Fig. 6, we conclude that E b is 0.86 and 1.97 eV per atom for S-atom migration across α-type and γ-type GBs, respectively. In our experiments, the GB dynamics observed at room temperature only occurs in the presence of an electron beam and thus the beam-specimen interaction must be a part of the driving force as well. The energy barriers for the S-atoms to migration across GBs are much higher than that available by thermal activation. Although the overall process is exothermal, beam irradiation is still necessary to accelerate the dynamics by transferring energies to GB atoms to overcome the barriers. The energy transfer due to elastic electron-nuclei collision may cause the slippage of S-layers during the GB evolution dynamics. As shown in Fig. S7, † the energy transfer onto the lattice S atom is about 5.87 eV, which is much higher than the energy barriers calculated (<2 eV) and sufficient to drive the S-atom migration via the vacancy mechanism. The beam irradiation thus assists sulfur vacancy to act as an active conveyance of S-atoms that dominates the GB evolution.
With these results, one may wonder why the experimentally observed α-type GBs are the most common GBs characterized, as our DFT calculation shows that the β-type GBs have the lowest formation energy among all the three types. During the normal growth process of TMDs in thermal equilibrium, the β-type GBs were also observed experimentally to be the most abundant in STM and TEM for MBE-grown TMDs, which supports the conclusion from our DFT calculations. However, under the electron beam irradiation, α-GBs appear extrinsically and more frequently within the normal 2H domain. Compared to the simple S-atom migration involved in the α-type GB formation, both Mo and S atoms will have to slide (driven by beam radiation) to form the β-type and γ-type GBs. Consequently, among these extrinsic GBs induced by beam radiation, the formation of α-type GBs is more probable than the other types in in situ TEM observation.

Conclusions
In summary, we present a joint experiment-theory investigation of high-symmetry GBs (α-, βand γ-types) in monolayer 2H-MoS 2 . The study shows that α-type GBs are always three-Mo-wide and six-or seven-Mo-long, consistent with the calculated formation energies and evolutionary kinetics of the GBs. Moreover, the interaction between parallel α-type GBs is repulsive, but attractive for those mis-aligned by 60°angle. The role of S-vacancy in GB formation is studied and we find that the formation energies of αand γ-type GBs decrease as the density of S-vacancy increases and the vacancy-GB distance decreases. The S-vacancy helps to reduce the formation energy of GBs and thereby promotes the GB formation. The vacancyinduced sliding-migration processes of S-atoms are analyzed to understand the dynamic evolution of GBs. It is believed that electron irradiation induced S-vacancy serves as an active conveyance of S-atoms in promoting the GB evolution. In principle, Fig. 6 Migration of S-atoms across a GB via a vacancy-mediated mechanism, which constitutes the basic step of GB dynamics observed in the experiments. (a) S-migration from a GB to the S-vacancy nearby for α-type and γ-type GBs. (b) The energy profiles calculated from DFT for a S-atom to migrate across a GB to the S-vacancy nearby for α-type and γ-type GBs.
our investigation on high-symmetry GBs can be extended to GBs with ordinary inter-domain mis-orientation angles in other chemical vapor deposited TMDs. The findings will also advise the control of GBs in 2D materials toward a specific GB network for catalytic applications such as the hydrogen evolution reaction.