Nearly hyperuniform, nonhyperuniform, and antihyperuniform density fluctuations in two-dimensional transition metal dichalcogenides with defects
Duyu Chen, Yu Zheng, Chia-Hao Lee, Sangmin Kang, Wenjuan Zhu, Houlong Zhuang, Pinshane Y. Huang, Yang Jiao
NNearly hyperuniform, nonhyperuniform, and antihyperuniform density fluctuations indefected two-dimensional transition metal dichalcogenides
Duyu Chen, ∗ Yu Zheng, Chia-Hao Lee, Sangmin Kang, WenjuanZhu, Houlong Zhuang, Pinshane Y. Huang, and Yang Jiao
7, 2 Tepper School of Business, Carnegie Mellon University, Pittsburgh, PA 15213, United States † Department of Physics, Arizona State University, Tempe, AZ 85287, United States Department of Materials Science and Engineering,University of Illinois Urbana-Champaign, Urbana, IL 61801, United States Semiconductor Research Center, Samsung Electronics, Hwaseong-si, Gyeonggi-do 445701, South Korea Department of Electrical and Computer Engineering,University of Illinois Urbana-Champaign, Urbana, IL 61801, United States Mechanical and Aerospace Engineering, Arizona State University, Tempe, AZ 85287, United States Materials Science and Engineering, Arizona State University, Tempe, AZ 85287, United States (Dated: February 5, 2021)Hyperuniform many-body systems in d -dimensional Euclidean space R d are characterized by com-pletely suppressed (normalized) infinite-wavelength density fluctuations, and appear to be endowedwith novel exotic physical properties. Recently, hyperuniform systems of disordered varieties havebeen observed in the context of various atomic-scale two-dimensional (2D) materials. In this work,we analyze the effects of localized defects on the density fluctuations across length scales and onthe hyperuniformity property of experimental samples of two-dimensional transition metal dichalco-genides. In particular, we extract atomic coordinates from time series annular dark field-scanningtransmission electron microscopy (ADF-STEM) imaging data of 2D tungsten chalcogenides with the2H structure (Te-doped 2H-WSe ) showing continuous development and evolution of electron-beaminduced defects, and construct the corresponding chemical-bonding informed coordination networksbetween the atoms. We then compute a variety of pair statistics and bond-orientational statisticsto characterize the samples. At low defect concentrations, the corresponding materials are nearlyhyperuniform, characterized by significantly suppressed scattering at the zero wave-number limit(omitting forward/ballistic scattering). As more defects are introduced, the (approximate) hyper-uniformity of the materials is gradually destroyed, and the system becomes non-hyperuniform evenwhen the material still contains a significant amount of crystalline regions. At high defect concen-trations, the structures become antihyperuniform with diverging (normalized) large-scale densityfluctuations, mimicking those typically observed at the thermal critical points associated with phasetransitions. Overall, the defected materials possess varying degrees of orientation order, and thereis apparently no intermediate hexatic phase emerging. To understand the observed nearly hyper-uniform density fluctuations in the slightly defected materials, we construct a minimalist structuralmodel and demonstrate that the experimental samples can be essentially viewed as perturbed hon-eycomb crystals with small correlated displacements and double chalcogen vacancies. Moreover,the small correlated displacements alone can significantly degrade hyperuniformity of the perfecthoneycomb structure. Therefore, even a small amount of vacancies, when coupled with correlateddisplacements, can completely destroy hyperuniformity of the system. I. INTRODUCTION
Hyperuniformity is a recently introduced novel conceptthat provides a unified framework to categorize crystals,quasicrystals, and certain unusual disordered systems [1–3]. Hyperuniform many-body systems in d -dimensionalEuclidean space R d are characterized by completely sup-pressed (normalized) density fluctuations at large lengthscales. In particular, the static structure factor S ( k ),which is proportional to the scattering intensity in X -ray, light, or neutron scattering experiments, vanishesin the infinite-wavelength (or zero-wavenumber) limit for ∗ correspondence sent to: [email protected] † Present address: Materials Research Laboratory, University ofCalifornia, Santa Barbara, California 93106, United States hyperuniform systems, i.e., lim k → S ( k ) = 0, where k isthe wavenumber. Here S ( k ) is defined as S ( k ) ≡ ρ ˜ h ( k ) , (1)where ρ is the number density of the system, and ˜ h ( k ) isthe Fourier transform of the total correlation function h ( r ) = g ( r ) − , (2)where g ( r ) is the pair correlation function. Note thatthis definition implies that the forward scattering contri-bution to the diffraction pattern is omitted. Equivalently,the local number variance σ N ( R ) ≡ (cid:104) N ( R ) (cid:105) − (cid:104) N ( R ) (cid:105) (3)of these systems associated with a spherical window ofradius R grows more slowly than the window volume (i.e., a r X i v : . [ c ond - m a t . s o f t ] F e b a scaling of R d in d -dimensional Euclidean space) in thelarge- R limit [1, 3], where N ( R ) is the number of particlesin a spherical window with radius R randomly placed intothe system. Recently, hyperuniformity has been observedin many physical, biological and material systems [4–20].Hyperuniform systems, in particular the disorderedvarieties, appear to be endowed with novel photonic,phononic, transport and mechanical properties, andwave-propagation characteristics [21–27]. For example,disordered hyperuniform dielectric networks were foundto possess complete photonic band gaps comparable insize to photonic crystals, while at the same time main-taining statistical isotropy, enabling waveguide geome-tries not possible with photonic crystals [21, 22]. More-over, disordered hyperuniform patterns can have nearlyoptimal color-sensing capabilities, as evidenced by avianphotoreceptors [28]. In addition, disordered stealthy hy-peruniform two-phase materials and cellular solids wererecently found to possess virtually optimal transportproperties [24, 26, 27]. The reader is referred to Ref.[3] for a thorough overview of hyperuniform systems.Despite the exotic structural characteristics and phys-ical properties that hyperuniform systems possess, inpractice it is very difficult to find perfect hyperuniformsystems of both ordered and disordered varieties due tothe inevitable existence of imperfection [29], except ina few cases where tailored optimization techniques aredesigned and deployed to fabricate perfect hyperuniformmaterials [30–34]. Therefore, it is important to system-atically investigate how various types of imperfectionsaffect hyperuniformity and the associated physical prop-erties of the systems. In the past, the investigation ofthe effect of imperfections on the physical and struc-tural properties of crystals is well documented [35–37].The effects of imperfections/defects on hyperuniformityhave also been extensively investigated theoretically inthe context of perturbed lattices, e.g., see Refs. [29, 38–41] and references therein. Recently, in a seminal study,using various theoretical models, Kim and Torquato [29]demonstrated that while thermal excitation and pointdefects such as vacancies and interstitials destroy hyper-uniformity, uncorrelated random displacements preservehyperuniformity, but could change the class of hyperuni-formity. To quantify the degree of hyperuniformity ofreal systems, the hyperuniformity index [29, 42–44] H ≡ S (0) /S ( k max ) (4)based on the structure factor S ( k ) is often employed,where k max is the position of the largest peak in theFourier space. We note that for most practical pur-poses effective hyperuniform systems with H ≤ − [29, 42–44] behave essentially the same as perfect hy-peruniform systems. However, systematic study of howthe introduction of imperfections affect hyperuniformityof real experimental systems, in particular atomic-scalelow-dimensional materials, is still lacking.Very recently, disordered hyperuniformity (DHU) hasbeen observed in the context of various atomic-scale two- dimensional (2D) materials [45–48]. For example, DHUis found to arise in patterns of electrons emerging from aquantum jamming transition of correlated many-electronstate in a quasi-two-dimensional dichalcogenides, whichleads to enhanced electronic transport [46]. It is alsodiscovered that DHU distribution of localized electronsin 2D amorphous silica results in an insulator-to-metaltransition in the material, which is in contrast to theconventional wisdom that disorder generally diminisheselectronic transport [47]. Moreover, Chen and coworkers[48] have rigorously demonstrated that the introductionof Stone-Wales (SW) topological defects [49] into hon-eycomb network structures preserves hyperuniformity ofthese systems to a large extent, and the resulting amor-phous structural models capture the salient features ofgraphene-like 2D materials at low temperatures.While SW defects are prevalent in graphene-like 2Dmaterials, other types of defects such as chalcogen va-cancies, metal vacancies and trefoil defects are dominantin transition metal dichalcogenide monolayers (TMDCs)such as MoS and WSe [50, 51]. For example, Lin andcoworkers [50] have elaborated how the local structuresevolve as various types of defects are introduced intoMoS .In this work, we conduct a comprehensive characteri-zation of the evolution of global structures as defects aregradually introduced into real experimental samples ofTMDCs, in the context of hyperuniformity. Specifically,we employ deep-learning algorithms to extract the atomicpositions in a sequence of image frames obtained fromADF-STEM, as the scanning electron probe continuouslyintroduces defects into a monolayer crystalline 2D transi-tion metal dichalcogenide alloy, Te-doped 2H-WSe . Wethen employ a multi-step procedure to identify and refinethe chemical-bonding informed coordination networks inthis evolving system. Subsequently, we employ a vari-ety of theoretical and quantitative tools from soft-matterphysics, in particular pair statistics and bond-orientationstatistics to quantify the evolution of global structures.We find that the systems are nearly hyperuniform at lowdefect concentrations, and the (approximate) hyperuni-formity is completely destroyed even when there is stilla significant portion of crystalline sites (specifically, lessthan 20% defects). No intermediate hexatic phase arefound to exist as the defects are gradually introducedinto the systems, which is distinctly different from the2D melting process as temperature increases. Moreover,we generalize an analytical formula to describe the struc-ture factor S ( k ) of crystal with correlated displacementsand vacancies. Using this analytical formula and MonteCarlo simulations, we demonstrate that the experimentalsamples in the early frames can be essentially viewed asperturbed crystals with small correlated displacementsand double chalcogen vacancies. Moreover, our resultsindicate the level of degradation of hyperuniformity thatone should expect due to the finite experimental mea-surement precision in real STEM experiments. In addi-tion, we note that our analysis procedures can be readilyadapted to characterize the structures of other orderedand disordered two-dimensional materials. FIG. 1. A schematic (a) of the 2H-WSe monolayer andits projection (b) when seen from above. Note that the Te-doped 2H-WSe monolayer (at low Te concentration can beviewed as effectively the 2H-WSe monolayer), when pro-jected from above, is mapped into a perfect honeycomb lat-tice (ignoring the local strains introduced by the Te substitu-tions, which are small compared to our measurement preci-sion), with each “particle” in the projected plane possessingthree bonds. Moreover, half of the honeycomb lattice sites areoccupied by the W atoms, and each of the other half sites isoccupied by 2 overlaying Se/Te atoms. Each pair of W sitesare separated by a Se/Te site, and vice versa. The rest of the paper is organized as follows: in Sec.II, we describe the methods that we employ to extractatomic coordinates and determine chemical bonds be-tween atoms from the obtained STEM images, as wellas the statistical descriptors that we use to characterizethese structures. In Sec. III, we employ various statisti-cal descriptors to characterize the evolving global struc-tures of the experimental samples. In Sec. IV, we presenta minimalist structural model of the real experimentalsamples in the early frames, and use analytical formulaand Monte Carlo simulations to validate it. In Sec. V,we provide concluding remarks.
II. METHODSA. Extraction of Atomic Coordinates andIdentification of Chemical-Bonding InformedCoordination Networks
We acquired aberration-corrected scanning transmis-sion electron microscopy images of Te-doped WSe . InADF-STEM, an angstrom-scale electron beam is scannedacross the sample, and scattered electrons are collected as a function of the position of the electron beam. Thematerial studied is WSe − x Te x where x = 0 .
06. Thesynthesis and TEM sample preparation methods werepreviously described in Ref. [51]. Because the Te frac-tion is small and the local lattice distortion induced byTe substitutions (2-4 pm) is below the measurement pre-cision of these frames ( ∼
15 pm), we do not expect theimpact of the Te to be significant (a schematic of the 2H-WSe structure is shown in Fig. 1), and the sample canbe treated as if it were primarily WSe . During imaging,beam-induced defects, primarily vacancies, voids, and lo-cal stripes of phase transformations [52], gradually mod-ify the underlying lattice. In this text, we analyzed a 120-frame atomic-resolution movie to capture the generationand evolution of beam-induced defects. The readers arereferred to the supplementary information of Ref. [51]for more details including sample fabrication, acquisitionparameters, and the movie itself. Next, we extracted2D-projected atomic coordinates from the movie using adeep learning package (AtomSegNet [53]). The extractedatomic coordinates are further processed to remove pos-sible artifacts from the imaging and the deep learningtreatment, particularly for the atomic coordinates thatare too close to each other. For example, we merge to-gether any group of atoms that are within 73.3 pm fromeach other by averaging over their coordinates. Since theaverage shortest projected bond length is around 190 pm,pair of atoms within the 73.3 pm threshold are considerednon-physical and thus merged together.To construct chemical-bonding informed coordinationnetworks from the final extracted atomic positions, wecompute the distance between each pair of atoms, and usea three-step procedure to identify and refine the chemi-cal bonds. We first assign a bond to those pairs with apair distance smaller than the cutoff 241.89 pm, which isfound to identify the bonds in the crystalline region rel-atively well. Next, we add bonds around the defects byallowing a larger cutoff of 322.52 pm for any particle withone or zero bonds identified in the first step. Finally, welimit the maximum number of bonds that any particlecan possess to three and remove the extra bonds by sort-ing bonds according to the bond length and only keepingthe shortest three bonds for each particle. Visual ex-amination indicates that overall our procedure generatesreasonably accurate coordination networks that faithfullyrepresent the actual chemical bonding network. B. Statistical descriptors
We consider a variety of statistical descriptors thatare sensitive in picking up structural information acrosslength scales. A basic quantity of a point configuration isthe aforementioned pair correlation function g ( r ), whichis proportional to the probability density function of find-ing two centers separated by distance r [54]. In practice, FIG. 2. Extracted atomic coordinates and determined chemical-bonding informed coordination networks overlaid with rawimages of different frames obtained using the ADF-STEM technique. g ( r ) is computed via the relation g ( r ) = (cid:104) N ( r ) (cid:105) ρ πr ∆ r , (5)where (cid:104) N ( r ) (cid:105) is the average number of particle centersthat fall into the circular ring at distance r from a centralparticle center (arbitrarily selected and averaged over allparticle centers in the system), 2 πr ∆ r is the area of thecircular ring, and ρ is the number density of the system[54, 55]. The static structure factor S ( k ) is the Fouriercounterpart, and for computational purposes, S ( k ) is theangular-averaged version of S ( k ), which can be obtaineddirectly from the particle positions r j , i.e., S ( k ) = 1 N (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N (cid:88) j =1 exp( i k · r j ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( k (cid:54) = ) , (6)where N is the total number of points in the system [55–57]. The trivial forward scattering contribution ( k =0) in Eq. 6 is omitted, which makes Eq. 6 completelyconsistent with the aforementioned definition of S ( k ) inthe ergodic infinite-system limit.The aforementioned local number variance σ N ( R ) isanother quantity that is often used to characterize den-sity fluctuations of many-body systems. To compute σ N ( R ), we randomly place circular observation windowswith radius R in the system under the constraint that thewindows should fall entirely within the image frame toavoid boundary artifacts [28, 58]. Also, the largest radiusof the window that one can sample must be much smallerthan half of the box length, otherwise density fluctuationsare artificially diminished [11]. We count the numberof particles N ( R ) that fall into the observation window,which is a random variable. The variance associated with N ( R ) is denoted by σ N ( R ) ≡ (cid:104) N ( R ) (cid:105) − (cid:104) N ( R ) (cid:105) , whichmeasure density fluctuations of particles within a windowof radius R .The bond-orientational order metric Q and correla-tion function C ( r ) [59, 60] are often used to study themelting process. Specifically, the order metric Q is de-fined as Q ≡ |(cid:104) Ψ (cid:105)| , (7)where Ψ ( r i ) = 1 n i n i (cid:88) j =1 e θ ij , (8)and (cid:104)· · · (cid:105) denotes ensemble average, n i is the number ofneighbors of vertex i located at r i , and θ ij is the polarangle associated with the vector from vertex i to j -thchemically bonded neighbor of vertex i .The bond-orientational correlation function C ( r ) isdefined as C ( r ) ≡ (cid:104) Ψ ( r )Ψ ∗ (0) (cid:105) , (9)where Ψ ∗ is the complex conjugate of Ψ . We note that Q = 1 and C ( r ) = 1 for a perfect honeycomb network;while for isotropic fluid phase, Q ≈ C ( r ) decayswith an exponential envelop at large r [59, 60]. To avoidartifacts caused by the image boundaries, we exclude thevertices that are within certain distance (439.8 pm in thiswork) from each edge of the bounding box of the image. III. STRUCTURAL CHARACTERIZATION OFEVOLVING GLOBAL STRUCTURE OFDEFECTED TWO-DIMENSIONAL TRANSITIONMETAL DICHALCOGENIDES
We apply the aforementioned scheme to extract atomiccoordinates from different frames of STEM images, asshown in Fig. 2. We then construct chemical-bonding in-formed coordination networks using the aforementionedprocedures and the resulting networks are shown in Fig.2. We note that crystalline Te-doped 2H-WSe , an al-loyed 2D TMDC monolayer, when projected from aboveonto a 2D plane, is mapped into a perfect honeycomblattice and thus hyperuniform, with each “particle” inthe projected plane possessing three bonds. Moreover,half of the honeycomb lattice sites are occupied by theW atoms, and each of the other half sites is occupied by 2overlaying Se/Te atoms. Visual examination of the con-figurations and the associated networks in Fig. 2 suggestthat the samples in the early frames appear to be primar-ily affected by double chalcogen vacancies. The Atom-SegNet [53] identifies any occupied site with significantimage intensity as an atom. As a result, it does not dis-tinguish between metal and chalcogen sites, and it onlyidentifies a site as a “defect” if there is no atom in projec-tion or the intensity is significantly weaker than others.Also, the atomic coordinates in the crystalline region de-viate slightly from the perfect honeycomb crystal, whichmight be due to various experimental factors including 1)detector noise, 2) uncertainties introduced by the deep-learning algorithm, and 3) instrument instabilities (sam-ple drift, scanning errors, mechanical vibrations). Thesefactors limit the measurement precision and accuracy to ∼
15 pm. It is noteworthy that various types of de-fects such as vacancies and interstitials, and correlateddisplacements may destroy or degrade hyperuniformityof the structures, as studied theoretically. In the laterframes, large voids begin to form in the samples. TheTe-doped WSe sample gradually evolves from a nearlyperfect crystal to a highly defective crystal due to thedamage introduced by the electron beam used for imag-ing. The high energy (80keV) electrons induce knock-on damage and radiolysis in the sample, producing vacan-cies, voids, and local lattice reconstruction. The abilityto image and generate atomic-scale defects allows us tosystematically study how hyperuniformity evolves withdefect concentration. A. Pair statistics
To characterize the density fluctuations of the exper-imental samples across different length scales, we com-pute various pair statistics of these samples, and the re-sults are shown in Fig. 3. In particular, as the defectsare gradually introduced into the system from frame 0to frame 20, the magnitudes of the Bragg peaks in S ( k )decrease, which indicates the degradation of the crys-talline order of the system. Moreover, the structure fac-tor S ( k ) of the structures in frames 0-10 appear to de-crease slightly or plateau to relatively small values as k approaches 0, indicating suppressed large-scale fluctu-ations. On the other hand, S ( k ) of frame 20 appearsto converge to value appreciably larger than zero as k decreases at small k , showing the destruction of hype-runiformity at this point. There are significant wigglesin S ( k ) at large k as well in all of these structures, sug-gesting short-scale structures in these materials. Also,the scaling exponent β in | h ( r ) | = | g ( r ) − | ∼ /r β increases from frame 0 to frame 10, i.e., the total cor-relation function h ( r ) decays faster as r increases as de-fects are introduced into the system, which is consistentwith the increasing disorder and loss of large-scale struc-tural correlation in the system. In addition, the normal-ized local number variance σ N ( R ) /R of the structures inframes 0-10 decreases slightly as R increases at large R ,indicating (approximate) hyperuniformity of the struc-tures, while σ N ( R ) scales as R for frame 20, indicatingthe loss of hyperuniformity at this point. To quantifythe degree of hyperuniformtiy, we employ the hyperuni-formity index H for the series of structures in differentframes by extrapolating S ( k ) to k = 0 with a linear fit-ting of S ( k ) within k ∈ [0 . pm − , . pm − ]. Wefind that H (cid:47) − for the first 10 frames and H > − for frame 20, suggesting that the structures in frames0-10 are nearly hyperuniform and the (approximate) hy-peruniformity is essentially destroyed for the structure inframe 20.For structures in frames 20-40, S ( k ) plateaus at small k and converges to finite values as k approaches 0, and σ N ( R ) scales as the volume of the observation window(i.e., R ) as shown in Fig. 4, indicating that the cor-responding systems enter the nonhyperuniform regime.In this regard, these structures are similar to typical liq-uids and glasses [3, 61]. As more defects are introducedand large voids begin to appear beyond frame 60, S ( k )of the corresponding structures start to diverge as k goesto zero, and σ N ( R ) grows faster than the window volume(i.e., σ N ( R ) /R is an increasing function of R ), and suchstructures can be regarded as antihyperuniform with hy-perfluctuations [3], since they are the antithesis of a hy-peruniform system. Note that antihyperuniform systemsare often observed at thermal critical points, all of whichhave fractal structures [3]. To our best knowledge, thisis the first time that antihyperuniformity is observed inreal disordered atomic-scale 2D materials. In addition,the total correlation function h ( r ) decays faster as r in-creases as more defects are introduced into the systembeyond frame 20, which is consistent with the increasingdisorder. B. Bond-orientational statistics
The process of introducing defects into the honey-comb lattice can also be viewed as the “melting” ofhexagonal 2D materials to some extent. It is widelyknown that 2D melting of colloidal systems is a two-stepcrystalline-hexatic-liquid phase transition, and the bond-orientational correlation function changes from oscillat-ing around certain constant to power-law scaling, andthen exponential scaling at large length scale as tem-perature increases [59, 60, 62, 63]. To investigate themelting behavior of our experimental system, we com-pute the bond-orientational order metric Q and bond-orientational correlation function C ( r ), which have beenroutinely used to study 2D melting [59, 60]. The resultsof C ( r ) and Q are shown in Fig. 5. It can be clearlyseen that C ( r ) oscillates around certain constant for allof the investigated structures, i.e., the scaling behaviordoes not change, although that constant decreases as de-fects are gradually introduced, suggesting decreased ori-entational order of the system. Moreover, Q also de-creases relatively smoothly, which is consistent with theresults of C ( r ). We note that Q is small but still muchlarger than 0 even for the structure in frame 119, whichmeans that there is remaining orientatinal order in thesystem, a reflection of the presence of the remaining crys-talline regions in the system. These results indicate thatthere is no intermediate hexatic emerging in the “melt-ing” of 2D TMDC monolayer, which is consistent withthe fact that there is no mechanism in the system lead-ing to the formation of grain boundary and the loss oflong-range orientational order. This behavior is differ-ent from the 2D melting of colloidal systems and similarto the observation in structural models of graphene-likematerials where disorder is introduced through the SWtopological defects [48].To quantify the “perfectness” of the honeycomb lat-tice during the damaging process, we compute the defectconcentration p d defined as: p d = 1 − N c /N , (10)where N c is the number of crystalline sites in a struc-ture (excluding the region within 439.8 pm of the edges),and N is the number of particles in a reference perfecthoneycomb lattice (if it were to occupy the same region) where the side length of a hexagon in the perfect crys-tal is set as the same as the average bond length of allthe crystalline sites in the experimental structure underconsideration. Here we consider a particle to reside ina crystalline site if the following conditions are met: 1)this particle has three bonds; 2) the bond length dif-ference of its longest bond and shortest bond is withincertain threshold (set as 36.65 pm here); 3) all of its bondangles are in the vicinity of 120 degrees (set as 100 ∼ p d is muchsmaller than 20%, i.e., when the material still contains asignificant amount of crystalline sites, and the structuresenter the antihyperuniform regime when p d exceeds 40%. IV. STRUCTURAL MODEL OF DEFECTEDTWO-DIMENSIONAL TRANSITION METALDICHALCOGENIDES
As mentioned above, the experimental samples in theearly frames appear to be primarily affected by doublechalcogen vacancies and correlated displacements. Tofully understand the effect of imperfections on the struc-tures, we construct a simplistic structural model of theexperimental samples that enables us to tune differentfactors independently and see the outcomes. Here weconsider a simple honeycomb lattice where each parti-cle is connected to its chemically-bonded neighbors bysprings of spring constant K . Specifically, the energy E of the system is given by E = (cid:88) i 500 particles atdifferent σ . Moreover, to introduce chalcogen vacanciesinto the displaced crystal, we randomly remove pN/ p , where p is defined as the ratioof double chalcogen vacancies over the total number ofchalcogen sites in the lattice. Representative configura-tions at σ = 0 . p are shown in Fig. 6.We note that chalcogen sites comprise half of the totalsites in the honeycomb lattice, and every pair of chalco-gen sites is separated by a metal site. We then compute S ( k ) of these stochastically displaced configurations withor without defects. To better compare the results withthe experimental systems, we normalize the distance inour simulated systems so that the characteristic lengthscale d = 2 π/k is the same for the initial experimentalframe and the vacancy-free displaced crystal at a giventemperature, where k is the position of the first peak.We find that S ( k ) of vacancy-free displaced honeycombcrystal at σ = 0 . S ( k ) at small k further in-creases as p increases, as shown in Fig. 7(a), indicatingthe loss of hyperuniformity. In addition, the scaling ex-ponent β in | h ( r ) | = | g ( r ) − | ∼ /r β slightly increasesas p increases, i.e., the total correlation function h ( r )decays faster as r increases as vacancies are introducedinto the system, as shown in Fig. 7(b). These trends areconsistent with those observed in experimental systems,and demonstrate that the series of experimental samplesstudied in the early frames can be essentially viewed as honeycomb crystals with small correlated displacementsand varying concentration of double chalcogen vacancies.To further validate our simulation results, we gener-alize a previously known analytical formula for S ( k ) ofdefective point process in any space dimension with spa-tially uncorrelated point vacancies [29] to displaced hon-eycomb crystals with double chalcogen vacancies. Wenote that this formula was first known for defective crys-tals [65, 66]. Specifically, following similar procedures inRef. [29] , we derive the generalized formula for S ( k ) as S ( k ) = 1 + (1 − p (cid:48) )[ S ( k ) − , = 1 + (1 − p S ( k ) − , = p − p S ( k ) , (14)where S ( k ) is the structure factor of the vacancy-freedisplaced crystal, and p (cid:48) = p/ S ( k ) of displaced honeycomb crystal with p concentra-tion of double chalcogen vacancies using this analyticalformula by plugging in the simulated S ( k ) of vacancy-free displaced crystals. We find that the analytical resultsmatch the directly simulated S ( k ) of defected displacedcrystals well, as shown in Fig. 8, which demonstratesthat our assumption is indeed valid. V. CONCLUSION AND DISCUSSION In this work, we investigated structural features oftimelapse STEM images of 2D TMDCs across lengthscales as an electron probe was used to gradually intro-duce various types of defects into the 2D materials. In (a)(b)(c) FIG. 6. Representative configurations of honeycomb crystalwith correlated Gaussian displacements with variance σ =0 . p . (a) p = 0. (b) p = 0 . 02. (c) p = 0 . particular, we quantified density fluctuations and the as-sociated hyperuniformity/antihyperuniformity propertyof these defected 2D materials. We find that in the earlyframes the defects are dominantly double chalcogen va-cancies, and at very low defect concentrations the cor-responding materials are nearly hyperuniform, as man-ifested in various pair statistics and quantified by thehyperuniformity index H . However, as additional de-fects are introduced, the (approximate) hyperuniformityof the materials is completely destroyed, when there issignificant amount of crystalline regions in the system.No intermediate hexatic phase emerged, which is differ-ent 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 struc- FIG. 7. Pair statistics of simulated honeycomb lattice withcorrelated Gaussian displacements with variance σ = 0 . S ( k )of the simulated frames and representative reference experi-mental frames. (b) Log-log plot of | g ( r ) − | with a linearfitting (dashed line). tures. By constructing a minimalist structural modelfor the samples in the early frames, we were able todemonstrate that the experimental samples can be es-sentially viewed as perturbed honeycomb crystals withsmall correlated displacements and double chalcogen va-cancies. Moreover, the correlated displacements alone,which is the result of various uncontrollable experimentalnoises and perturbations that one should usually expectwhen acquiring STEM images of atomic-scale 2D mate-rials, largely degrade hyperuniformity of the system, andlow concentration of vacancies, when coupled with corre-lated displacements, basically destroy (approximate) hy-peruniformity.We note that here we primarily studied the effect ofdouble chalcogen vacancies and correlated displacementson the density fluctuations of defected TMDCs. How-ever, if other more complex types of defects such as trifolddefects can be introduced into experimental samples ofTMDCs in a controllable manner, in principle we can em-ploy similar procedures to investigate the effect of those0 FIG. 8. Structure factor S ( k ) of honeycomb crystal with cor-related Gaussian displacements with variance σ = 0 . p = 0 . 02. (b) p = 0 . defects. Also, in the future it would be interesting tolook at how SW defects, when coupled with other factors,affect structural features and physical properties of realgraphene-like materials in experiments. Previously, it hasbeen demonstrated [48] that the introduction of SW de-fects and local structural relaxation alone preserves hy-peruniformity of honeycomb lattice while the disorderof the structure increases, and are accompanied by theemergence of exotic physical properties. More generally,to our best knowledge we for the first time introduce avariety of theoretical machinery from soft-matter physicsto study the structural evolution of experimental sam-ples of atomic-scale 2D materials, and these techniquescan be readily adapted and applied to other ordered anddisordered 2D materials. The structural studies of 2Dmaterials in this paper and related works will strengthenour fundamental understanding of the physics underlyingthese materials, and serve as the basis for future func-tional materials design. ACKNOWLEDGMENTS We are grateful for Dr. JaeUk Kim and Dr. RunzeTang for very helpful discussion. H.Z. thank the start-up funds from ASU. Y. J. thanks ASU for support andPeking University for hospitality during his sabbaticalleave. C.-H. L. and P. Y. H. acknowledge funding sup-port from the U.S. Department of Energy, Office of Ba-sic Energy Sciences, Division of Materials Sciences andEngineering under award [1] S. Torquato and F. H. Stillinger, Phys. Rev. E , 041113(2003).[2] S. Torquato, J. Phys. Condens. Matter , 414012(2016).[3] S. Torquato, Phys. Rep. , 1 (2018).[4] A. Donev, F. H. Stillinger, and S. Torquato, Phys. Rev.Lett. , 090604 (2005).[5] C. E. Zachary, Y. Jiao, and S. Torquato, Phys. Rev. Lett. , 178001 (2011).[6] C. E. Zachary and S. Torquato, Phys. Rev. E , 051133(2011).[7] R. D. Batten, F. H. Stillinger, and S. Torquato, J. Appl.Phys. , 033504 (2008).[8] R. D. Batten, F. H. Stillinger, and S. Torquato, Phys.Rev. Lett. , 050602 (2009).[9] Y. Jiao and S. Torquato, Phys. Rev. E , 041309 (2011).[10] J. L. Lebowitz, Phys. Rev. A , 1491 (1983).[11] R. Dreyfus, Y. Xu, T. Still, L. A. Hough, A. G. Yodh,and S. Torquato, Phys. Rev. E , 012302 (2015).[12] D. Hexner and D. Levine, Phys. Rev. Lett. , 110602 (2015).[13] M. Hejna, P. J. Steinhardt, and S. Torquato, Phys. Rev.B , 245204 (2013).[14] A. Gabrielli, M. Joyce, and F. S. Labini, Phys. Rev. D , 083523 (2002).[15] D. Hexner and D. Levine, Phys. Rev. Lett. , 020601(2017).[16] D. Hexner, P. M. Chaikin, and D. Levine, Proc. Natl.Acad. Sci. U.S.A. , 4294 (2017).[17] J. H. Weijs and D. Bartolo, Phys. Rev. Lett. , 048002(2017).[18] Z. Ding, Y. Zheng, Y. Xu, Y. Jiao, and W. Li, Phys. Rev.E , 063101 (2018).[19] Q.-L. Lei, M. P. Ciamarra, and R. Ni, Sci. Adv. ,eaau7423 (2019).[20] Q.-L. Lei and R. Ni, Proc. Natl. Acad. Sci. U.S.A. ,22983 (2019).[21] M. Florescu, S. Torquato, and P. J. Steinhardt, Proc.Natl. Acad. Sci. U.S.A. , 20658 (2009).[22] W. Man, M. Florescu, E. P. Williamson, Y. He, S. R. Hashemizad, B. Y. C. Leung, D. R. Liner, S. Torquato,P. M. Chaikin, and P. J. Steinhardt, Proc. Natl. Acad.Sci. U.S.A. , 15886 (2013).[23] G. Gkantzounis, T. Amoah, and M. Florescu, Phys. Rev.B , 094120 (2017).[24] G. Zhang, F. H. Stillinger, and S. Torquato, J. Chem.Phys. , 244109 (2016).[25] Y. Xu, S. Chen, P. Chen, W. Xu, and Y. Jiao, Phys. Rev.E , 043301 (2017).[26] D. Chen and S. Torquato, Acta Mater. , 152 (2018).[27] S. Torquato and D. Chen, Multifunct. Mater. , 015001(2018).[28] Y. Jiao, T. Lau, H. Hatzikirou, M. Meyer-Hermann, J. C.Corbo, and S. Torquato, Phys. Rev. E , 022721 (2014).[29] J. Kim and S. Torquato, Phys. Rev. B , 054105 (2018).[30] S. Torquato, G. Zhang, and F. H. Stillinger, Phys. Rev.X , 021020 (2015).[31] G. Zhang, F. Stillinger, and S. Torquato, Phys. Rev. E , 022119 (2015).[32] R. A. DiStasio Jr, G. Zhang, F. H. Stillinger, andS. Torquato, Phys. Rev. E , 023311 (2018).[33] J. Kim and S. Torquato, Acta Mater. , 143 (2019).[34] J. Kim and S. Torquato, Phys. Rev. E , 052141 (2019).[35] N. W. Ashcroft and N. D. Mermin, Solid State Physics(Brooks/Cole, Cengage Learning, Belmont, 1976).[36] C. Kittel and P. McEuen,Introduction to solid state physics, vol. 8 (Wiley,New York, 1976).[37] P. M. Chaikin and T. C. Lubensky,Principles of condensed matter physics (Cambridgeuniversity press, Cambridge, 2000).[38] T. R. Welberry, G. H. Miller, and C. E. Carroll, ActaCrystallogr. Sect. A: Found. Crystallogr. , 921 (1980).[39] A. Gabrielli, Phys. Rev. E , 066131 (2004).[40] A. Gabrielli and S. Torquato, Phys. Rev. E , 041105(2004).[41] M. A. Klatt, J. Kim, and S. Torquato, Phys. Rev. E ,032118 (2020).[42] F. Martelli, S. Torquato, N. Giovambattista, and R. Car,Phys. Rev. Lett. , 136002 (2017).[43] M. A. Klatt and S. Torquato, Phys. Rev. E , 012118(2018).[44] D. Chen, E. Lomba, and S. Torquato, Phys. Chem.Chem. Phys. , 17557 (2018).[45] W. Zhou, X. Zou, S. Najmaei, Z. Liu, Y. Shi, J. Kong,J. Lou, P. M. Ajayan, B. I. Yakobson, and J.-C. Idrobo, Nano Lett. , 2615 (2013).[46] Y. A. Gerasimenko, I. Vaskivskyi, M. Litskevich,J. Ravnik, J. Vodeb, M. Diego, V. Kabanov, and D. Mi-hailovic, Nat. Mater. , 1078 (2019).[47] Y. Zheng, L. Liu, H. Nan, Z.-X. Shen, G. Zhang, D. Chen,L. He, W. Xu, M. Chen, Y. Jiao, et al., Sci. Adv. ,eaba0826 (2020).[48] D. Chen, Y. Zheng, L. Liu, G. Zhang, M. Chen, Y. Jiao,and H. Zhuang, Proc. Natl. Acad. Sci. U.S.A. ,e2016862118 (2021).[49] A. J. Stone and D. J. Wales, Chem. Phys. Lett. , 501(1986).[50] Y.-C. Lin, T. Bj¨orkman, H.-P. Komsa, P.-Y. Teng, C.-H.Yeh, F.-S. Huang, K.-H. Lin, J. Jadczak, Y.-S. Huang,P.-W. Chiu, et al., Nat. Commun. , 6736 (2015).[51] C.-H. Lee, A. Khan, D. Luo, T. P. Santos, C. Shi, B. E.Janicek, S. Kang, W. Zhu, N. A. Sobh, A. Schleife, et al.,Nano Lett. , 3369 (2020).[52] Y.-C. Lin, D. O. Dumcenco, Y.-S. Huang, and K. Sue-naga, Nat. Nanotechnol. , 391 (2014).[53] R. Lin, R. Zhang, C. Wang, X.-Q. Yang, and H. Xin,arXiv preprint arXiv:2012.09093 (2020).[54] S. Torquato, Random Heterogeneous Materials(Springer-Verlag, New York, 2002).[55] S. Atkinson, F. H. Stillinger, and S. Torquato, Phys. Rev.E , 062208 (2013).[56] C. E. Zachary and S. Torquato, J. Stat. Mech. Theor.Exp. , P12015 (2009).[57] D. Chen, Y. Jiao, and S. Torquato, J. Phys. Chem. B , 7981 (2014).[58] D. Chen, W. Y. Aw, D. Devenport, and S. Torquato,Biophys. J. , 2534 (2016).[59] K. Zahn, R. Lenke, and G. Maret, Phys. Rev. Lett. ,2721 (1999).[60] K. Wierschem and E. Manousakis, Physical Rev. B ,214108 (2011).[61] J.-P. Hansen and I. R. McDonald,Theory of simple liquids (Academic Press: London,1986).[62] J. M. Kosterlitz and D. J. Thouless, J. Phys. C , 1181(1973).[63] X. Xu and S. A. Rice, Phys. Rev. E , 011602 (2008).[64] K. Hun, Proc. R.Soc. London, Ser. A , 102 (1947).[65] P. H. Dederichs, J. Phys. F: Met. Phys. , 471 (1973).[66] H. Peisl, J. Phys. Colloq.37