Nearly hyperuniform, nonhyperuniform, and antihyperuniform density fluctuations in two-dimensional transition metal dichalcogenides with defects

Hyperuniformity is a recently introduced novel concept that provides a unified framework to categorize crystals, quasicrystals, and certain unusual disordered systems Torquato and Stillinger (2003); Torquato (2016, 2018). Hyperuniform many-body systems in $d$-dimensional Euclidean space $\mathrm{R}^{d}$ are characterized by completely suppressed (normalized) density fluctuations at large length scales. In particular, the static structure factor $S(k)$, which is proportional to the scattering intensity in $X$-ray, light, or neutron scattering experiments, vanishes in the infinite-wavelength (or zero-wavenumber) limit for hyperuniform systems, i.e., $\lim_{k\rightarrow 0}S(k)=0$, where $k$ is the wavenumber. Here $S(k)$ is defined as

where $\rho$ is the number density of the system, and $\tilde{h}(k)$ is the Fourier transform of the total correlation function

where $g_{2}(r)$ is the pair correlation function. Note that this definition implies that the forward scattering contribution to the diffraction pattern is omitted. Equivalently, the local number variance

of these systems associated with a spherical window of radius $R$ grows more slowly than the window volume (i.e., a scaling of $R^{d}$ in $d$-dimensional Euclidean space) in the large-$R$ limit Torquato and Stillinger (2003); Torquato (2018), where $N(R)$ is the number of particles in a spherical window with radius $R$ randomly placed into the system. Recently, hyperuniformity has been observed in many physical, biological and material systems Donev et al. (2005); Zachary et al. (2011); Zachary and Torquato (2011); Batten et al. (2008, 2009); Jiao and Torquato (2011); Lebowitz (1983); Dreyfus et al. (2015); Hexner and Levine (2015); Hejna et al. (2013); Gabrielli et al. (2002); Hexner and Levine (2017); Hexner et al. (2017); Weijs and Bartolo (2017); Ding et al. (2018); Lei et al. (2019); Lei and Ni (2019); Ru19 ; San19 ; San20 .

Hyperuniform systems, in particular the disordered varieties, appear to be endowed with novel photonic, phononic, transport and mechanical properties, and wave-propagation characteristics Florescu et al. (2009); Man et al. (2013); Gkantzounis et al. (2017); Zhang et al. (2016); Xu et al. (2017); Chen and Torquato (2018); Torquato and Chen (2018). For example, disordered hyperuniform dielectric networks were found to possess complete photonic band gaps comparable in size to photonic crystals, while at the same time maintaining statistical isotropy, enabling waveguide geometries not possible with photonic crystals Florescu et al. (2009); Man et al. (2013). Moreover, disordered hyperuniform patterns can have nearly optimal color-sensing capabilities, as evidenced by avian photoreceptors Jiao et al. (2014). In addition, disordered stealthy hyperuniform two-phase materials and cellular solids were recently found to possess virtually optimal transport properties Zhang et al. (2016); Chen and Torquato (2018); Torquato and Chen (2018). The reader is referred to Ref. Torquato (2018) for a thorough overview of hyperuniform systems.

Despite the exotic structural characteristics and physical properties that hyperuniform systems possess, in practice it is very difficult to find perfect hyperuniform systems of both ordered and disordered varieties due to the inevitable existence of imperfection Kim and Torquato (2018), except in a few cases where tailored optimization techniques are designed and deployed to fabricate perfect hyperuniform materials Torquato et al. (2015); Zhang et al. (2015); DiStasio Jr et al. (2018); Kim and Torquato (2019a, b). Therefore, it is important to systematically investigate how various types of imperfections affect hyperuniformity and the associated physical properties of the systems. In the past, the investigation of the effect of imperfections on the physical and structural properties of crystals is well documented Ashcroft and Mermin (1976); Kittel and McEuen (1976); Chaikin and Lubensky (2000). The effects of imperfections/defects on hyperuniformity have also been extensively investigated theoretically in the context of perturbed lattices, e.g., see Refs. Welberry et al. (1980); Gabrielli (2004); Gabrielli and Torquato (2004); Kim and Torquato (2018); Klatt et al. (2020) and references therein. Recently, in a seminal study, using various theoretical models, Kim and Torquato Kim and Torquato (2018) demonstrated that while thermal excitation and point defects such as vacancies and interstitials destroy hyperuniformity, uncorrelated random displacements preserve hyperuniformity, but could change the class of hyperuniformity. To quantify the degree of hyperuniformity of real systems, the hyperuniformity index Martelli et al. (2017); Kim and Torquato (2018); Klatt and Torquato (2018); Chen et al. (2018)

based on the structure factor $S(k)$ is often employed, where $k_{max}$ is the position of the largest peak in the Fourier space. We note that for most practical purposes effective hyperuniform systems with $H\leq 10^{-4}$ Martelli et al. (2017); Kim and Torquato (2018); Klatt and Torquato (2018); Chen et al. (2018) behave essentially the same as perfect hyperuniform systems. However, systematic study of how the introduction of imperfections affect hyperuniformity of real experimental systems, in particular atomic-scale low-dimensional materials, is still lacking.

Very recently, disordered hyperuniformity (DHU) has been observed in the context of various atomic-scale two-dimensional (2D) materials Zhou et al. (2013); Gerasimenko et al. (2019); Zheng et al. (2020); Chen et al. (2021). For example, DHU is found to arise in patterns of electrons emerging from a quantum jamming transition of correlated many-electron state in a quasi-two-dimensional dichalcogenides, which leads to enhanced electronic transport Gerasimenko et al. (2019). It is also discovered that DHU distribution of localized electrons in 2D amorphous silica results in an insulator-to-metal transition in the material, which is in contrast to the conventional wisdom that disorder generally diminishes electronic transport Zheng et al. (2020). Moreover, Chen and coworkers Chen et al. (2021) have rigorously demonstrated that the introduction of Stone-Wales (SW) topological defects Stone and Wales (1986) into honeycomb network structures preserves hyperuniformity of these systems to a large extent, and the resulting amorphous structural models capture the salient features of graphene-like 2D materials at low temperatures.

While SW defects are prevalent in graphene-like 2D materials, other types of defects such as chalcogen vacancies, metal vacancies and trefoil defects are dominant in monolayer transition metal dichalcogenides (TMDCs) such as MoS₂ and WSe₂ Lin et al. (2015); Lee et al. (2020). For example, Lin and coworkers Lin et al. (2015) have elaborated how the local structures evolve as various types of defects are introduced into MoS₂. There are also a few preliminary experimental studies Lehnert et al. (2019); Leiter et al. (2020); Komsa et al. (2013) examining the evolution of structures as various types of defects are introduced into samples of TMDCs.

In this work, we conduct a comprehensive characterization of the evolution of global structures as defects are gradually introduced into real experimental samples of TMDCs, in the context of hyperuniformity. Specifically, we employ deep-learning algorithms to extract the atomic positions in a sequence of image frames obtained from ADF-STEM, as the scanning electron probe continuously introduces defects into a monolayer crystalline 2D transition metal dichalcogenide alloy, Te-doped 2H-WSe₂. We then employ a multi-step procedure to identify and refine the chemical-bonding informed coordination networks in this evolving system. Subsequently, we employ a variety of theoretical and quantitative tools from soft-matter physics, in particular pair statistics and bond-orientation statistics to quantify the evolution of global structures. We find that the systems are nearly hyperuniform at low defect concentrations, and the (approximate) hyperuniformity is completely destroyed even when there is still a significant portion of crystalline sites (specifically, less than $20\%$ defects). No intermediate hexatic phase are found to exist as the defects are gradually introduced into the systems, which is distinctly different from the 2D melting process as temperature increases. Moreover, we generalize an analytical formula to describe the structure factor $S(k)$ of crystal with correlated displacements and vacancies. Using this analytical formula and Monte Carlo simulations, we demonstrate that the experimental samples in the early frames can be essentially viewed as perturbed crystals with small correlated displacements and double chalcogen vacancies. Moreover, our results indicate the level of degradation of hyperuniformity that one should expect due to the finite experimental measurement precision in real STEM experiments. In addition, we note that our analysis procedures can be readily adapted to characterize the structures of other ordered and disordered two-dimensional materials.

The rest of the paper is organized as follows: in Sec. II, we describe the methods that we employ to extract atomic coordinates and determine chemical bonds between atoms from the obtained STEM images, as well as the statistical descriptors that we use to characterize these structures. In Sec. III, we employ various statistical descriptors to characterize the evolving global structures of the experimental samples. In Sec. IV, we present a minimalist structural model of the real experimental samples in the early frames, and use analytical formula and Monte Carlo simulations to validate it. In Sec. V, we provide concluding remarks.

We apply the aforementioned scheme to extract atomic coordinates from different frames of STEM images, as shown in Fig. 3. We then construct chemical-bonding informed coordination networks using the aforementioned procedures and the resulting networks are shown in Fig. 3. We note that crystalline Te-doped 2H-WSe₂, an alloyed 2D TMDC monolayer, when projected from above onto a 2D plane, is mapped into a perfect honeycomb lattice and thus hyperuniform, with each “particle” in the projected plane possessing three bonds. Moreover, half of the honeycomb lattice sites are occupied by the W atoms, and each of the other half sites is occupied by 2 overlaying Se/Te atoms. Visual examination of the configurations and the associated networks in Fig. 3 suggest that the samples in the early frames appear to be primarily affected by double chalcogen vacancies. The AtomSegNet Li20 identifies any occupied site with significant image intensity as an atom. As a result, it does not distinguish between metal and chalcogen sites, and it only identifies a site as a “defect” if there is no atom in projection or the intensity is significantly weaker than others. For example, although single chalcogen vacancies are also prevalent Lee et al. (2020) in these Te-doped 2H-WSe₂ samples, these vacancies still contain single chalcogen atoms and will not be identified as defects in the projected structures. Also, the atomic coordinates in the crystalline region deviate slightly from the perfect honeycomb crystal, which might be due to various experimental factors including 1) detector noise, 2) uncertainties introduced by the deep-learning algorithm, and 3) instrument instabilities (sample drift, scanning errors, mechanical vibrations). These factors limit the measurement precision and accuracy to $\sim 15$ pm. It is noteworthy that various types of defects such as vacancies and interstitials, and correlated displacements may destroy or degrade hyperuniformity of the structures, as studied theoretically. In the later frames, large voids begin to form in the samples. The Te-doped WSe₂ sample gradually evolves from a nearly perfect crystal to a highly defective crystal due to the damage introduced by the electron beam used for imaging. The high energy (80keV) electrons induce knock-on damage and radiolysis in the sample, producing vacancies, voids, and local lattice reconstruction. The ability to image and generate atomic-scale defects allows us to systematically study how hyperuniformity evolves with defect concentration. In addition, we note that when analyzing any real experimental samples from images, one is almost inevitably limited by the finite precision for the determination of atomic coordinates, which effectively adds random uncorrelated displacements to the particle positions. However, as previously proven in a theoretical study, random uncorrelated displacements preserve hyperuniformity, so finite measurement precision should not affect our ability of determining the particular hyperuniformity property of a given experimental system.

As mentioned above, the experimental samples in the early frames appear to be primarily affected by double chalcogen vacancies and correlated displacements. To fully understand the effect of imperfections on the structures, we construct a simplistic structural model of the experimental samples that enables us to tune different factors independently and see the outcomes. Here we consider a simple honeycomb lattice where each particle is connected to its chemically-bonded neighbors by springs of spring constant $K$. Specifically, the energy $E$ of the system is given by

where $H_{ij}=1$ indicates that vertices $i$ and $j$ are connected by a chemical bond, and $\mathbf{u}_{i}$ is the vector displacement of vertex $i$ from its corresponding reference honeycomb lattice site. We introduce correlated displacements to the particles in the honeycomb lattice according to a multivariate Gaussian distribution:

where $\sigma^{2}=k_{B}T/(Ka^{2})$ is an effective variance associated with the Gaussian distribution and can be viewed as “fictitious” dimensionless temperature. Here $k_{B}$ is the Boltzmann’s factor, $T$ is the “fictitious” temperature, and $a$ is the side length of the hexagon in the original honeycomb lattice. By varying $T$, one can effectively vary the variance $\sigma^{2}$. It is noteworthy that uncorrelated stochastic displacements of crystals preserve perfect hyperuniformity, and for this reason we employ correlated displacements in our model to account for the loss of perfect hyperuniformity observed in the relatively defect-free experimental samples in the first few frames. Experimentally, correlated displacements can result from as scan distortions, sample drift, mechanical vibration, noise, and local phase transitions. However, thermal motion can be excluded from the potential source of these displacements since the time scale associated with thermal motion is much shorter than the STEM image acquisition time, and thus the obtained atomic positions are time-averaged and the thermal motion information filtered out.

Note that introducing correlated displacements in the above fashion is mathematically equivalent to generating equilibrium structures at finite positive temperature, which would allow us to utilize standard Monte Carlo simulations to obtain perturbed honeycomb lattices with correlated displacements according to Gaussian distribution with different variance $\sigma^{2}$. Specifically, at each trial move, each vertex is allowed to randomly move within a prescribed maximal distance from its old position in each dimension and the trial move is accepted with the probability

where $E_{old}$ and $E_{new}$ are the energies of the system before and after the trial move as defined in Eq. 13. To fully equilibrate the system, we perform $1000N$ trial moves, where $N$ is the number of particles/vertices in the system.

Using this scheme, we generate configurations of displaced honeycomb crystals with $N=2,500$ particles at different $\sigma^{2}$. Moreover, to introduce chalcogen vacancies into the displaced crystal, we randomly remove $pN/2$ particles at chalcogen sites according to prescribed vacancy concentration $p$, where $p$ is defined as the ratio of double chalcogen vacancies over the total number of chalcogen sites in the lattice. Representative configurations at $\sigma^{2}=0.005$ with varying $p$ are shown in Fig. 7. We note that chalcogen sites comprise half of the total sites in the honeycomb lattice, and every pair of chalcogen sites is separated by a metal site. We then compute $S(k)$ of these stochastically displaced configurations with or without defects. To better compare the results with the experimental systems, we normalize the distance in our simulated systems so that the characteristic length scale $d_{0}=2\pi/k_{0}$ is the same for the initial experimental frame and the vacancy-free displaced crystal at a given temperature, where $k_{0}$ is the position of the first peak.

We find that $S(k)$ of vacancy-free displaced honeycomb crystal at $\sigma^{2}=0.005$ shown in Fig. 8(a) captures the salient features of the experimental sample in the initial frame, which contains relatively few vacancies. Also note here that correlated displacements already largely degrades hyperuniformity even in the absence of vacancies, similar to the effect of thermal motions Hun (1947); Dederichs (1973) . Moreover, as more chalcogen vacancies are introduced into the simulated system, $S(k)$ at small $k$ further increases as $p$ increases, as shown in Fig. 8(a), indicating the loss of hyperuniformity. In addition, the scaling exponent $\beta$ in $|h(r)|=|g_{2}(r)-1|\sim 1/r^{\beta}$ slightly increases as $p$ increases, i.e., the total correlation function $h(r)$ decays faster as $r$ increases as vacancies are introduced into the system, as shown in Fig. 8(b). These trends are consistent with those observed in experimental systems, and demonstrate that the series of experimental samples studied in the early frames can be essentially viewed as honeycomb crystals with small correlated displacements and varying concentration of double chalcogen vacancies. In Fig. 8(a) we also show $S(k)$ of the reference perfect honeycomb crystal with $N=2,500$ particles in Fig. 7(a), and $S(k)$ is strictly zero in the intervals between the sharp Bragg peaks.

To further validate our simulation results, we generalize a previously known analytical formula for $S(k)$ of defective point process in any space dimension with spatially uncorrelated point vacancies Kim and Torquato (2018) to displaced honeycomb crystals with double chalcogen vacancies. We note that this formula was first known for defective crystals Peisl (1976); Dederichs (1973). Specifically, following similar procedures in Ref. Kim and Torquato (2018) , we derive the generalized formula for $S(k)$ as

where $S_{0}(k)$ is the structure factor of the vacancy-free displaced crystal, and $p^{\prime}=p/2$ is the effective vacancy concentration in the system, i.e., the ratio of chalcogen vacancies over the total number of lattice sites. Since in our model vacancies can only be introduced randomly at chalcogen sites, strictly speaking the vacancies are not completely spatially uncorrelated, but we assume that this spatial correlation can be neglected. We compute $S(k)$ of displaced honeycomb crystal with $p$ concentration of double chalcogen vacancies using this analytical formula by plugging in the simulated $S_{0}(k)$ of vacancy-free displaced crystals. We find that the analytical results match the directly simulated $S(k)$ of defected displaced crystals well, as shown in Fig. 9, which demonstrates that our assumption is indeed valid.

Because certain correlated displacements of atomic positions may be due to various experimental measurement errors rather than intrinsic features of the experimental samples, we also use computer simulations to investigate the large-scale structural features of honeycomb crystals with double chalcogen vacancies and no correlated displacements. Specifically, we randomly remove particles at chalcogen sites and compute the structure factor $S(k)$ of the defected honeycomb crystals at different double chalcogen vacancy concentrations without correlated displacements, as shown in Fig. 10. It is noteworthy that $S(k)$ is relatively flat away from the Bragg peaks associated with the honeycomb crystals, and the values of $S(k)$ in these intervals between Bragg peaks are roughly proportional to the effective vacancy concentration, which is consistent with previous theoretical predictions Kim and Torquato (2018).

In this work, we investigated structural features of timelapse STEM images of 2D TMDCs across length scales as an electron probe was used to gradually introduce various types of defects into the 2D materials. In particular, we quantified density fluctuations and the associated hyperuniformity/antihyperuniformity property of these defected 2D materials. We find that in the early frames the chemical bonding-informed coordination network is mainly influenced by double chalcogen vacancies, and at very low defect concentrations the corresponding materials are nearly hyperuniform, as manifested in various pair statistics and quantified by the hyperuniformity index $H$. However, as additional defects are introduced, the (approximate) hyperuniformity of the materials is completely destroyed, when there is significant amount of crystalline regions in the system. No intermediate hexatic phase emerged, which is different from the 2D melting process for colloidal systems. In later frames large voids begin to form in the samples, leading to the rise of antihyperuniformity of the structures. By constructing a minimalist structural model for the samples in the early frames, we were able to demonstrate that the experimental samples can be essentially viewed as perturbed honeycomb crystals with small correlated displacements and double chalcogen vacancies. Moreover, the correlated displacements alone, which is the result of various uncontrollable experimental noises and perturbations that one should usually expect when acquiring STEM images of atomic-scale 2D materials, largely degrade hyperuniformity of the system, and low concentration of vacancies, when coupled with correlated displacements, basically destroy (approximate) hyperuniformity.

We note that here we primarily studied the effect of double chalcogen vacancies and correlated displacements on the density fluctuations of defected TMDCs. However, if other more complex types of defects such as trefoil defects can be introduced into experimental samples of TMDCs in a controllable manner, in principle we can employ similar procedures to investigate the effect of those defects. Also, it is noteworthy that while chalcogen vacancies are dominant in the projected structures of certain samples of TMDCs such as those in the early frames of the current work, SW defects are prevalent in graphene-like 2D materials. In the future it would be interesting to look at how SW defects, when coupled with other factors, affect structural features and physical properties of real graphene-like materials in experiments. Previously, it has been demonstrated Chen et al. (2021) that the introduction of SW defects and local structural relaxation alone preserves hyperuniformity of honeycomb lattice while the disorder of the structure increases, and is accompanied by the emergence of exotic physical properties. It is important to note that while a single experimental realization of amorphous graphene was previously found to be hyperuniform Chen et al. (2021), even the samples of TMDCs in the very early frames in the current work are found to be only nearly (or approximately) hyperuniform. More generally, to our best knowledge we for the first time introduce a variety of theoretical machinery from soft-matter physics to study the structural evolution of experimental samples of atomic-scale 2D materials, and these techniques can be readily adapted and applied to other ordered and disordered 2D materials. The structural studies of 2D materials in this paper and related works will strengthen our fundamental understanding of the physics underlying these materials, and serve as the basis for future functional materials design.