## Abstract

Disordered hyperuniformity (DHU) is a recently proposed new state of matter, which has been observed in a variety of classical and quantum many-body systems. DHU systems are characterized by vanishing infinite-wavelength normalized density fluctuations and are endowed with unique novel physical properties. Here, we report the discovery of disordered hyperuniformity in atomic-scale two-dimensional materials, i.e., amorphous silica composed of a single layer of atoms, based on spectral-density analysis of high-resolution transmission electron microscopy images. Moreover, we show via large-scale density functional theory calculations that DHU leads to almost complete closure of the electronic bandgap compared to the crystalline counterpart, making the material effectively a metal. This is in contrast to the conventional wisdom that disorder generally diminishes electronic transport and is due to the unique electron wave localization induced by the topological defects in the DHU state.

## INTRODUCTION

Disordered hyperuniform (DHU) systems are a unique class of disordered systems that suppress large-scale density fluctuations such as crystals and yet have no Bragg peaks (*1*, *2*). For a point configuration (e.g., a collection of particle centers of a many-body system), hyperuniformity is manifested as the vanishing structure factor in the infinite-wavelength (or zero–wave number) limit, i.e., *k* = 2π/λ is the wave number. In this case of a random field, the hyperuniform condition is given by *2*). It has been suggested that hyperuniformity can be considered as a new state of matter (*1*), which has a hidden order in between of that of a perfect crystal and a totally disordered system (e.g., a Poisson distribution of points).

Recently, a wide spectrum of physical and biological systems have been identified to have the remarkable property of hyperuniformity, which include the density fluctuations in early universe (*3*), disordered jammed packing of hard particles (*4*, *5*), certain exotic classical ground states of many-particle systems (*6*, *7*), jammed colloidal systems (*8*, *9*), driven nonequilibrium systems (*10*, *11*, *12*, *13*), certain quantum ground states (*14*), avian photoreceptor patterns (*15*), organization of adapted immune systems (*16*), amorphous silicon (*17*), a wide class of disordered cellular materials (*18*), dynamic random organizing systems (*19*–*21*), and even the distribution of primes on the number axis (*22*). In addition, it has been shown that hyperuniform materials can be designed to have superior physical properties including large isotropic photonic bandgaps (*23*).

In this paper, we report the discovery of hyperuniformity in amorphous two-dimensional (2D) silica [conventionally modeled as “continuous random networks” (*24*)], based on the analysis of aberration-corrected transmission electron microscopy (TEM) images of the material. To the best of our knowledge, this is the first discovery of disordered hyperuniformity in atomic-scale 2D materials (i.e., those composed of a single layer of atoms), which can have unique novel electronic, magnetic, and optical properties compared to their bulk counterparts.

We show via simulations that the observed DHU in amorphous silica is closely related to the strong topological and geometrical constraints induced by the local chemical order in the system. In addition, our density functional theory (DFT) calculations show that DHU significantly reduces the electronic bandgap in 2D amorphous silica, leading to almost complete closure of the bandgap compared to the crystalline counterpart. This is in contrast to the conventional wisdom that disorder generally diminishes electronic transport and is due to the unique electron wave localization induced by DHU.

## RESULTS AND DISCUSSION

### Hyperuniformity in 2D amorphous silica

We first analyze the high-resolution TEM images of 2D amorphous silica, see Fig. 1A and fig. S1 for a large-field image. The material samples were fabricated using chemical vapor deposition, and the procedure for obtaining the imaging dataset was reported in detail in (*25*) and briefly described in the Supplementary Materials. As shown in Fig. 1A, the black spots (with diffusive boundaries) represent the silicon atoms. The micrographs are processed to retain the distribution information of the silicon atoms by thresholding and fitting the grayness intensity distribution associated with each silicon atom using a Gaussian function, i.e., *G*(*x*) = *I*_{0}*e*^{−∣x − xi∣2/σ2}, where *I*_{0} is the maximal intensity, *x _{i}* is the center of the silicon atom, and σ is an effective radius.

The associated spectral density *26*) and shown in the inset of Fig. 1C. The angularly averaged *k* =∣**k**∣) is shown in Fig. 1C. We note that the spectral density analysis (instead of structure factor) is used here to efficiently use all of the information on Si atom distributions contained in the TEM micrographs and minimize the possible systematic errors induced by converting the intensity map into center distributions. In addition, we note that spectral density would be a more natural indicator of hyperuniformity in quantum systems. The TEM images can be considered as a mapping of electron density of the system, which achieves peak values surrounding the Si atoms, allowing us to also extract the well-defined center distribution of Si atoms. However, for a general system, the electron density map (e.g., associated with extended states) might be better represented by a random field, for which the spectral density would be a better metric for hyperuniformity. It can be seen that *k* values, which indicates that the 2D amorphous silica samples analyzed are hyperuniform (see the Supplementary Materials for detailed scaling analysis).

### Numerical models of DHU 2D amorphous silica networks

Next, we computationally generate disordered hyperuniform 2D silica networks. We note that unlike 3D amorphous materials such as metallic glasses, which can be simulated by numerically quenching a high temperature liquid state, even very rapid quenching of a 2D material will lead to a highly crystalline material with a small number of local defects. These defected crystalline materials cannot represent the experimentally obtained DHU amorphous silica network. Together with the observation that the 2D silica systems are disordered hyperuniform, this motivates us to use a structure-based method, which is a two-step approach.

In the first step, a three-coordinated DHU network is generated using a modified “collective coordinate” approach (*7*), in which a random initial configuration of points is gradually evolved to match a prescribed targeted structure factor *S**(*k*) while simultaneously satisfying mutual exclusion volume constraints. The targeted *S**(*k*) = 0 for *k* ≤ *K** drives the system to a hyperuniform state, and the exclusion volume constraints ensure that the final configuration can be feasibly mapped to an amorphous silica network. In the second step, the obtained point configuration is mapped to a three-coordinated network by connecting a point with its three nearest neighbors. This network is further converted to a silica network, by placing a silicon atom centered at each point and placing an oxygen atom at the midpoint of the two connected silicon atoms. We then perform molecular statics simulations using the Si─O potential based on the Tersoff parameterization (*27*) as implemented in the LAMMPS program (*28*) to optimize the constructed silica network to physically metastable states by minimizing their total potential energy.

The spectral density of the simulated amorphous silica network (see Fig. 1B) is computed by placing a Gaussian kernel function at the center of each silicon atom and is shown in Fig. 1D. It can be clearly seen that *k* limit, i.e., *k* values (see the Supplementary Materials for detailed scaling analysis and comparison to experimental results). We also compute and analyze the structure factors associated with the center distribution of Si atoms in both the experimental and simulated 2D amorphous silica, which have the same small-*k* behavior, i.e., *S*(*k*) ~ *k* and lim_{k → 0}*S*(*k*) = 0. In addition, we directly analyze the Si atom center number variance σ^{2}(*R*) as a function of the radius of circular observation window *R*, which also confirms the same scaling behavior, i.e., σ^{2}(*R*) ∼ *R* ln (*R*) for large *R* (see the Supplementary Materials for details).

We note that in the second step of the simulation (i.e., the potential energy minimization), hyperuniformity [i.e., the small *k* values of *S*(*k*)] is not constrained anymore. The evolution of the system is dominated by the interactions of atoms and constrained by the topology, and the positions of the atoms have been significantly perturbed compared to the final configuration obtained in the first step. Yet, the resulting system is still hyperuniform. This result suggests the strong geometrical and topological constraints, i.e., the bond length and angle associated with the Si-O bonds, as well as three-coordinated configurations, induced by the local chemical order could lead to the observed hyperuniformity in the system. We note that the emergence of disordered hyperuniformity from geometrical constraints has also been observed in jammed particle packings (*4*, *5*), avian photo receptor patterns (*15*), and amorphous laser speckle patterns (*29*).

We also obtain and compare the number of *n*-fold rings (where *n* = 3,4,…) formed by silicon atoms in both the experimental and simulated networks, see Fig. 1 (E and F). The ring statistics for the two systems again agree very well with one another. In addition, the comparison of local structural statistics of the experimental and numerical systems, including the pair correlation function and nearest neighbor distribution of the Si atoms, also shows excellent agreement (see the Supplementary Materials for details). These results indicate that our numerical network model can provide a statistically accurate structural representation of the 2D amorphous silica system by capturing key correlations on both large- and small-length scales, as well as the local topological order. Therefore, we expect that the physical properties computed based on the numerical network model should also be representative of those of the experimental system.

### Stone-Wales defects preserve hyperuniformity in 2D materials

What is the origin of DHU in 2D amorphous silica? We note that as a first approximation, the 2D silica glass can be considered as obtained from a 2D crystalline silica network by continuously introducing the Stone-Wales (SW) defects (*30*), which change the topology of the network. However, the SW defects do not affect the number of particles within a large observation window, since the SW transformation is localized and only affects a pair of atoms on the single-bond scale. Nonetheless, SW transformations of atom pairs on the boundary of the observation window might lead to a bounded fluctuation of particle numbers, which are scaled with the surface area of the window. Therefore, the SW defects should preserve hyperuniformity in the system. We provide numerical evidence for this speculation in the Supplementary Materials. Although the actual 2D amorphous silica has a structure that deviates from the ideal SW transferred crystalline network, the above argument could provide a possible explanation of the observed DHU in the system.

### Disordered hyperuniformity significantly reduces electronic bandgap

We use the simulated DHU SiO_{2} networks to calculate its density of states (DOS) at the DFT-Perdew-Burke-Ernzerhof (*31*) level of theory. Specifically, the DHU structure consists of three sublayers (1800 atoms; 600 Si atoms and 1200 O atoms). For comparison, we also create a supercell of 2D crystalline SiO_{2} with the same number of Si and O atoms. Several models of 2D crystalline SiO_{2} have been studied using DFT calculations in the literature (*32*). Here, we refer 2D crystalline SiO_{2} to the hexagonal bilayer crystalline network observed in the experiment (*33*).

We use numerical atomic orbitals (*34*) as implemented in the Atomic-orbital Based Ab-initio Computation at UStc (ABACUS) package (*35*) for calculating the electronic structure. The simulation methods and parameters can be found in Materials and Methods. For comparison, we also compute the DOS for 2D hexagonal crystalline SiO_{2}. The energy difference between the 2D DHU and hexagonal crystalline SiO_{2} calculated from the Tersoff potential and DFT is both positive and comparable (0.074 and 0.134 eV/atom, respectively). The positive energy differences indicate that the crystalline structure is more energetically stable than the DHU structure. Nevertheless, the TEM image (see Fig. 1A) shows that the experimental atomic structure of 2D amorphous SiO_{2} is drastically different from the crystalline model. By contrast, the high similarity between the DHU model and experimentally observed atomic structure reveals the metastable nature of DHU systems.

Figure 2A shows that the 2D crystalline SiO_{2} is essentially an insulator with a predicted bandgap of 5.31 eV, consistent with the previously reported bandgap of 5.48 eV calculated for a unit cell at the same level of theory (*36*). By contrast, we observe from Fig. 2B that a small but finite number of states occupy the Fermi level of the DHU structure, showing metallic behavior of the electrons with a typical bandgap of ~50 meV. This is comparable to the thermal fluctuations at room temperature of ~25 meV. In other words, the disordered hyperuniformity fundamentally changes the electrical transport behavior of 2D SiO_{2}, from an effective insulator at room temperature (as in the crystalline form) to an effective metal (as in the DHU form).

From Fig. 2B, we estimate the density ρ of the electrons that contribute to the electrical conductivity of the DHU structure at room temperature. By integrating the number of states in the energies ranging from 25 meV (corresponding to the thermal energy) to the Fermi level in Fig. 2B, we determine ρ as 2.33× 10^{12} cm^{−2}. This magnitude belongs to the category of “high doping” (e.g., 6.0 × 10^{11} and 9.2 × 10^{12} cm^{−2}) applied to common 2D semiconductors such as WS_{2} and MoS_{2} (*37*). To put it another way, considering the electron density alone, the conductivity associated with the DHU structure could be comparable to those of the 2D materials.

To better understand this metallic behavior of DHU 2D SiO_{2}, we compute the charge densities within an energy window of 0.5 eV below the highest occupied molecular orbital (HOMO) level for the crystalline structure and below the Fermi level for the DHU structure. For the former structure, Fig. 3A (a complete version of Fig. 3 can be found in the Supplementary Materials) shows that the electrons are distributed around Si and O atoms in the entire structure, i.e., fully occupying the valence bands. The amount of these electrons is significant, as can be seen from the large DOS below the HOMO level. But these electrons cannot be thermally excited at room temperature to the conduction bands due to the large bandgap, leading to zero electrical conductivity for pure 2D crystalline SiO_{2}.

On the other hand, for the DHU structure, the number of valence electrons in the same energy window is much less, resulting in a low carrier density but, nevertheless, a nonzero conductivity. A closer look at the slightly wider energy window associated with lowest DOS (e.g., −2 to 0 eV) reveals that the distribution of states still forms an almost continuous spectrum of peaks, see the inset in Fig. 2B. Figure 3B reveals that the valence electrons contributing to the conductivity originate from a small portion of the Si and O atoms in the DHU system.

### Topological defects in DHU SiO_{2} lead to electron localization with high energies

To understand why the electrons are localized around these atoms, we show in Fig. 3 (C and D) the distributions of the computed potential energies using the Tersoff potential for the crystalline and DHU systems, respectively. Consistent with Fig. 3A, Fig. 3C shows that the potential energies are homogeneously distributed in the crystalline system. On the contrary, for the DHU system, we observe from Fig. 3D that the potential energies have a highly heterogeneous distribution, with significantly higher energies localized in regions where the atomic arrangements substantially deviate from the sixfold hexagonal configuration. These are topological defects, which are necessary to achieve DHU in amorphous 2D silica. The high-energy localizations perfectly coincide with the electron localizations (see Fig. 3B), indicating that the local potential energies due to the topological defects induced by DHU are sufficiently high to activate the electrons to the energy levels near the Fermi energy.

Electron localization in a disordered system is generally associated with a number of exotic physical phenomena (*38*) such as metal-to-insulator transitions suggested by Anderson (*39*). The electron localization in the DHU silica system appears to be phenomenologically different from the Anderson localization, as the DHU localization gives rise to the opposite (i.e., insulator-to-metal) transition. A mathematical model based on model Hamiltonian that considers the disordered potential needs to be devised to investigate the electrical transport property of the localized electrons. Furthermore, whether the insulator-to-metal transition also occurs in other 2D DHU semiconductors/insulators remains unknown and is certainly worth exploring.

In summary, we have discovered, for the first time, disordered hyperuniformity in 2D amorphous materials and showed that DHU fundamentally changes the electronic transport behaviors in the material, making 2D DHU silica metallic. This interesting prediction awaits experimental confirmation. We also showed that the metallic behavior (i.e., with a virtually continuous spectrum of DOS without notable bandgaps in DHU silica) is resulted from the localization of high-energy states due to the topological defects induced by DHU. Since the observed DHU in 2D silica is associated with the local topological and geometrical constraints common in many 2D materials, we would expect to observe DHU and, thus, the metallic behavior in the amorphous states of other 2D materials, should these states be metastable at least. With the increasing interest in 2D amorphous materials, we expect our methods of building realistic structural model of amorphous 2D material systems along with large-scale DFT calculations to be applicable to a wide range of other 2D materials such as graphene (*40*) and molybdenum disulfide (*41*) in the amorphous form.

Last, we emphasize that the conclusions about electron localization, closure of bandgaps, and enhanced transport behaviors are mainly based on the numerical simulations. In future work, we will develop theories and carry out experiments that prove the predicted electron localization and enhanced electron transport behaviors in DHU 2D silica. On the theory side, one possible approach is to analytically examine the propagation of electron waves via homogenization of Schrödinger’s equation with model potentials having HU distributions (e.g., mimicking the effects of Si atoms in the system), following a similar contrast expansion procedure for obtaining optimal bounds on photonic bandgaps described in (*42*). On the experimental side, we will carry out surface nanoscale conductivity measurements and tunneling spectroscopy, which have been successfully used to identify emergent hyperuniformity in many-electron glasses formed via a quantum jamming transition (*43*).

## MATERIALS AND METHODS

### Numerical realization of hyperuniform configurations of Si atoms

As the first step in the generation of the DHU silica network models, we use a generalized collective coordinate (*7*) method to generate DHU realizations of Si atoms. Specifically, we first generate a random initial configuration of nonoverlapping circular disks with diameter σ, which is taken to be the average nearest neighbor distance between a pair of Si atoms in the SiO_{2} network. Next, the positions of the Si atoms are randomly perturbed, subject to the nonoverlapping constraints, to generate new configurations of Si atoms. For each new configuration, an effective energy *E*, defined as the sum of the absolute value of the structure factor *S*(*k*) associated the configuration for *k* < *K**, i.e., *E* = ∑_{k < K*}∣*S*(*k*)∣. This effective energy is minimized using the simulated annealing method: A new configuration (generated by random perturbation of the Si atom positions) will be accepted and replace the old configuration with the probability *p* = min {1, exp [(*E*_{old} − *E*_{new})/*T*]}, where *T* is the annealing parameter that is gradually decreased during the simulation. The final “ground state” configuration has a *S*(*k*) = 0 for *k* < *K**, which is then mapped to a three-coordinated network, with an O atom inserted in the middle of every Si─Si bond, to generate a model SiO_{2} network.

### Molecular statics simulations

We use the Tersoff potential (a representative force field for SiO_{2}) to describe the energy dependence of the SiO_{2} network on the atomic coordinates. We then apply the conjugate gradient algorithm as implemented in LAMMPS to minimize the energy of the system by optimizing the in-plane lattice constants and atomic locations until the energy convergence criterion of 1.0 × 10^{−16} eV is satisfied. This energy minimization process is called the molecular statics simulations, which are distinct from molecular dynamics simulations, where atoms are given initial velocities and adjust their positions following Newton’s second law.

### DFT calculations

The DFT-based electronic structure calculations are computed by ABACUS with numerical atomic orbitals. We adopt the PBE functional and pseudopotential and use double-zeta plus polarization basis sets (Si: 2*s*2*p*1*d*, O:2*s*2*p*1*d*), the cutoffs for Si and O orbitals are set to 8 and 7 arbitrary units, respectively. The energy cutoff is set to 100 Ry. A single *k*-point (Γ point) is used for the two large supercells each containing 1800 atoms.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/16/eaba0826/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank P. Y. Huang for helpful discussion and sharing TEM images of amorphous silica. Y.Z. and H.N. thank the Arizona State University (ASU) for the University Graduate Fellowship. The molecular statics simulations were performed on the ASU Agave computer cluster. The DFT calculations were performed on the USTC HPC facilities. We also thank the anonymous reviewers for extremely helpful and constructively comments and suggestions for improving the paper.

**Funding:**This research was supported by the start-up funds from ASU.

**Author contributions:**W.X., M.C., Y.J., and H.Z. conceived the project and supervised the research. Y.Z., H.N., G.Z., and D.C. performed the DHU characterizations. L.L., Z.-X.S., and L.H. performed molecular statics and DFT calculations. All authors analyzed the results and wrote the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).