Amorphous chalcogenides as random octahedrally bonded solids: I. Implications for the first sharp diffraction peak, photodarkening, and Boson peak Free

Generating trustworthy structures for glass-forming liquids near the conventional laboratory glass transition on the seconds-to-hours time scale remains out of reach for the present-day computer technology. According to the random first order transition (RFOT) theory,^1–4 structural relaxation near laboratory glass transition occurs via relatively compact, cooperative events involving several hundred rigid molecular units, implying a cooperativity size approaching a thousand atoms or so. Conversely, simulated systems of smaller size do not sample the appropriate portion of the phase space. This can be reflected in the amorphous samples being overly strained or, for instance, exhibiting an excess of homonuclear bonds in a binary alloy.⁵ State-of-the-art simulations of molecular glassformers⁶ appear to have reached microsecond relaxation times. It is expected that a glass transition on the microsecond scale still exhibits much of the physics relevant near the typical laboratory T_g,^1,7 though significantly modified by mode-coupling effects.⁸ At these shorter time scales, structural relaxations can be pictured as a compact core containing a relatively modest number of rigid molecular units, 50 or so,⁸ but decorated by string-like excitations that each measure a few such molecular units.^9,10 Again, one must simulate systems containing thousands of atoms to sample the appropriate portion of the phase space. These notions are consistent with severe finite size effects seen in simulations of model binary mixtures.¹¹ 

The present work aims to generate high quality structures for covalently bonded amorphous materials exhibiting distorted octahedral coordination. We hereby develop a novel structure-generation methodology and apply it to a family of alloys As_xSe_1−x, with a particular emphasis on the stoichiometric alloy As_0.4Se_0.6 (As₄₀Se₆₀). These alloys belong to an important class of intermetallic compounds often called the chalcogenides, which contain elements S, Se, and Te, and various combinations of elements from groups 14 and 15, such as Ge₂Se₂Te₅. The chalcogenides have enjoyed renewed interest in the past decade because of their potential applications in non-volatile computer memory and smart optics, among others.^12–14 The mechanism of current-driven “switching” between the amorphous and crystalline states of these fascinating compounds¹⁵ is still poorly understood, despite a half-century worth of effort.^16,17 Additionally, basic challenges posed by the amorphous chalcogenides include the puzzling midgap electronic states, which reveal themselves as light-induced ESR signal and midgap optical absorption,^18,19 and anomalous photo-luminescence.^20–22 Many other interesting phenomena occurring in chalcogenides²³ await detailed explanation, including: photo-darkening,²⁴ field-induced changes in ultrasonic attenuation,²⁵ light-induced circular dichroism,²⁶ reversible nano-distortion,²⁷ photomelting at cryogenic temperatures,²⁸ and the exponential (“Urbach”) tail in light absorption.²⁹ 

Similar to other classes of glassy compounds, the chalcogenides exhibit the so-called first sharp diffraction peak (FSDP),³⁰ whose detailed structural origin has been debated for decades, see Refs. 31–34 and references therein. The FSDP is engendered by structural motifs on a supra-atomic length scale that is often referred to as medium-range order (MRO) to distinguish it from the short-range order (SRO), which corresponds to the atomic length scale, and the longer-range “persistence” of oscillations in the pairwise correlation function g(r), which may extend for several nanometers in network materials.³⁵ Despite its large magnitude, the persistence length is essentially an integer multiple of the atomic length scale.³⁵ In contrast, the FSDP is not strictly tied to the atomic spacing.^32,36 One of the most interesting properties of the FSDP is that its strength in a glass modestly increases with temperature,³⁷ contrary to the expectation that a static structural feature should be progressively smeared as the vibrational amplitude increases. In addition, the FSDP is suppressed at increased pressures and is negatively correlated with photodarkening.³⁸ 

Glassy materials also display dynamic phenomena on lengthscales exceeding the atomic length. Arguably the best known example of such dynamic phenomena is referred to as the “Boson peak.” The latter is an umbrella term referring to degrees of freedom seen by vibrational spectroscopies at frequencies near 1 THz and has been also associated with the “bump” in the heat capacity and excess phonon scattering observed around the so-called plateau in the heat conductivity.^39–41 The term “Boson” refers to the linear temperature dependence of the peak’s strength. Novikov and Sokolov⁴² have reported an intriguing correlation between the position of the peak and FSDP.

When in crystalline form, the stoichiometric compound As₂Se₃ consists of covalently bonded double layers.⁴³ Inter-layer bonding is significantly weaker and can be formally classified as a closed-shell interaction. Except for inherent malcoordinated motifs,⁴⁴ the bonding in the amorphous state is expected to be similar in that each arsenic atom will typically form three covalent bonds and selenium two bonds with its immediate neighbors;⁴³ most atoms are surrounded by a perfect Lewis octet. Similar notions underlied early approaches to the structure of amorphous chalcogenides as network glasses.^45–48 

Accordingly, the stoichiometric compound As_0.4Se_0.6 can be thought of as consisting of corner-sharing pyramidal units, such as AsSe_3/2. For other stoichiometries, motifs other than the heteronuclear chain patterns Se–As–Se and As–Se–As must be present. This structural model has been supported by many studies.^49–52 Importantly in the present context, Zhugayevych and Lubchenko (ZL) have argued^43,44 that similar to their crystalline counterparts, glassy chalcogenides can be thought of as symmetry broken, distorted versions of higher symmetry parent structures defined locally on the simple cubic lattice. In this way, the aforementioned pyramidal AsSe_3/2 units can be thought of as distorted versions of a molecular unit consisting of an arsenic atom coordinated strictly octahedrally to three selenium atoms, see illustrations in the supplementary material.

The connection with octahedrality is important for two related reasons: On the one hand, it rationalizes the notion that the stability of the compound is largely due to $ppσ$ bonding while sp-mixing is a perturbation;⁴³ this is consistent with electronic structure calculations.⁵³ On the other hand, it suggests that the appropriate parent structures for chalcogenides should be coordinated octahedrally, not tetrahedrally. These conclusions are not obvious based on geometric grounds alone because bond-bond angles in the chalcogenides are numerically close to $100°$⁠, which is just midway between the $90°$ and $110°$ angles symptomatic of strictly octahedral and tetrahedral coordination, respectively.

Covalently bonded glass-formers, such as the chalcogenides, pose additional challenges to the already difficult task of simulating glassy melts because the effective forces between the atoms, as mediated by the electrons, are many-body and much more complicated than the relatively simple force fields employed when simulating common model glass-formers. Arguably the most significant realization of such electron-mediated many-body effects in octahedrally bonded compounds is the synergic relation between two adjacent, nearly co-linear covalent and secondary bonds and its interplay with sp-mixing.^43,54 Much like in the classic $π$-back bonding,⁵⁵ adjacent covalent and secondary, close-shell interactions can share electrons so that their respective strengths and the corresponding bond lengths change in accord, as has been documented in thousands of compounds.⁵⁶ 

Current approaches^53,57–60 to generating the structure of amorphous chalcogenides seem to converge on the so-called first principles molecular dynamics⁶¹ (FPMD) and related methodologies. Hereby one effectively infers the force fields governing atomic motions based on the ground state electronic energy evaluated using, for instance, the density functional theory (DFT). Aside from finite size effects and limitations on the simulation length, as discussed above, the accuracy of structure prediction by such methodologies is limited only by the accuracy of the DFT-based evaluation of the electronic ground state. Because of the aforementioned computational costs, present-day FPMD simulations must be performed at temperatures significantly exceeding typical laboratory glass transition temperatures and for modestly sized systems. Despite these potential limitations, FPMD has been remarkably successful in reproducing many features of the pair correlation functions in chalcogenide melts,⁵³ likely helped by the apparently weak temperature dependence of such pair correlation functions.^37,62 Notwithstanding these successes, establishing the structural signature of the FSDP using direct molecular modeling has been difficult. The systematic trend of the FSDP with the partial quantity of arsenic in As_xSe_1−x⁶³ has not been consistently reproduced in theoretical studies to date; the peak is found for some stoichiometries but not others.⁵³ 

The present approach deliberately gives up the objectivity of first principles molecular dynamics toward any given specific structural motif. Instead, here we essentially engineer amorphous structures so as to control the amount of octahedrality in the coordination and, to an extent, the number of homo- and hetero-nuclear contacts and other properties of the resulting bond network. We build our structures using a random close-packed lattice of monodisperse hard spheres as a scaffold, add a desired number of lattice points that are maximally octahedrally coordinated, assign specific atoms to the sites of the lattice, introduce vacancies as needed for stoichiometry, and then geometrically optimize the so-obtained parent structure.

This structure-building algorithm is motivated by the following notions: Random close-packed structures of hard spheres are easily generated even for relatively large systems. At the same time, the simple cubic structure (in which every atom is exactly octahedrally coordinated!) can be generated by filling the octahedral cavities in a particular periodic close-packed structure, viz., the fcc lattice, see Fig. 1(a). Likewise, here we explore the possibility that high density amorphous structures with distorted octahedral coordination can be generated by filling a subset of cavities in a random close-packed structure, as in Fig. 1(b). Within the spectral range of its spatial profile corresponding to the random close-packed scaffold, these structures have the highest value of coordination attainable by an aperiodic lattice.

We presently generate multiple structures under a variety of conditions. When optimized using relatively high quality force fields, the samples exhibit the distorted octahedral coordination and a spatial coexistence of covalent and secondary bonding expected of the chalcogenides.^43,54 At the same time, no extended double layers characteristic of crystalline As₂Se₃ are found. We consistently recover the FSDP at several stoichiometries and the variation of its position with arsenic content. We subsequently argue the FSDP is a spectral signature of the (aperiodic) pattern formed by regions of weak, secondary bonding. The complementary regions of stronger bonding correspond to covalent interactions. The resulting, spatially modulated density pattern of valence electrons represents a particular type of charge density wave^54,64,65 (CDW), whose crests and troughs—corresponding to regions of stronger and weaker bonding—are located off the atomic sites, and is often called a bond-order CDW. (An archetypal instance of such a bond-order wave is represented by the bond-alternation pattern in trans-polyacetylene.⁶⁶) Furthermore, the hereby established connection between the CDW and the FSDP will be used to discuss effects of pressure and temperature on the latter peak, and its connection with photodarkening. A connection can be made with much earlier work of Elliott,^32,36,67 who has ascribed the FSDP to ordering of interstices around cation-centered clusters. In contrast with relatively recent, FPMD-based studies by Bauchy, Kachmar, and Micoulaut⁵³ (BKM), we find the secondary bonding and the peak itself are mostly due to homo-nuclear (secondary) contacts. Also in contrast with the BKM study, the covalent bonds in the present structures are almost exclusively heteronuclear, consistent with experiment.⁶⁸ 

The notion of relatively compact, strongly bonded regions separated by regions of relatively weak bonding further allows one to connect the present results to the random first order transition (RFOT) theory.^3,4 The effective particle of the latter theory is a rigid molecular unit,^1,8 often called the “bead,” that is not significantly deformed during activated transport in a glassy liquid. While easily identifiable in some cases—think, e.g., of monomers in a polymeric substance⁶⁹—the identity of the bead in a covalently bonded substance is less clear. The bead represents the effective particulate degree of freedom relevant to configurational equilibration in glassy liquids as opposed, for instance, to mixing effects⁷⁰ or vibrations of individual bonds. Accordingly, the bead size naturally emerges in the Bevzenko-Lubchenko formalism^71,72 as the ultra-violet cut-off for the configurational degrees of freedom. Lubchenko and Wolynes⁸ have proposed and tested a prescription by which the effective bead size is determined by calibrating the fusion entropy of a substance by that of a weakly interacting liquid with a relatively rigid core, such as argon. This prescription is unambiguous and has enabled dozens of quantitative predictions,^3,4 consistent with the notion that the configurational entropy of a liquid, which pertains to the translational degrees of freedom, is numerically close to the excess liquid entropy relative to the corresponding crystal.⁷³ Yet we do not have a direct way to determine the bead size that does not explicitly refer to the crystal phase, except for some model liquids amenable to replica-symmetry breaking calculations.⁷⁴ The present results thus advance us toward resolving this long-standing challenge. Furthermore, combined with Lubchenko and Wolynes’ findings^40,41 on the relation between the Boson peak frequency and the size of the cooperativity size in glassy melts—which is directly related to the bead size—this notion sheds light on the aforementioned correlation between the FSDP and the Boson peak.

The article is organized as follows: The present algorithm to generate the octahedrally bonded amorphous structure is described in Sec. II. In Sec. III, we discuss structures optimized using DFT and argue on the detailed origin of the first sharp diffraction peak. In Sec. IV, we delve deeper into the structural analysis, by inducing artificial excess of tetrahedral coordination and comparing the resulting structures with those obtained in Sec. III. Section V summarizes the findings and argues that the presently proposed microscopic picture accounts for the observed pressure and temperature dependence of the FSDP, as well as its empirical correlation with photodarkening and the Boson peak.

The intuition behind the present algorithm is graphically illustrated in Fig. 1. First, we note that the simple cubic structure is composed of two equivalent, inter-penetrating fcc structures. The two can be visualized, respectively, as the anions and cations of the rock salt structure. Furthermore, each vertex of one fcc lattice is located in the center of an octahedral cavity of the other fcc lattice, Fig. 1(a). (There are two tetrahedral cavities per each octahedral cavity in the fcc lattice.) By analogy, we anticipate that by identifying non-tetrahedral cavities in the random close-packed structure and placing atoms in the centers of such cavities, one will obtain an aperiodic octahedrally bonded solid, as alluded to in Fig. 1(b). The algorithm is summarized in the block diagram in Fig. 2.

In the first step, we generate a densely packed, disordered lattice with both uniform density and coordination characteristic of a glassy liquid. As a practical matter, we generate a maximally coordinated disordered lattice with a large repeat unit. (The periodic boundary conditions—as opposed to a sample with a free surface–will allow us to obtain the electronic ground state energy for the ultimate solid using the computationally efficient plane-wave DFT.) To this end, we generate our disordered, periodically continued sample using the Lubachevsky-Stillinger^75,76 (LS) algorithm. The latter algorithm produces a jammed collection of monodisperse hard spheres at a filling fraction approaching that of the putative random closed-packed (rcp) solid, which exhibits the maximum coordination attainable by a bulk-disordered assembly of monodisperse hard spheres.^77,78 Specifically, our structures have the filling fraction 0.638. From here on, we will designate the centers of the spheres as the sites of the scaffold RCP lattice, as exemplified by the red spheres in Fig. 1(b).

Monodisperse spheres crystallize readily below the melting temperature,⁷⁸ the filling fraction of the liquid in equilibrium with the corresponding crystal being around 0.49.⁷⁹ The LS algorithm employs a combination of particle moves and particle swelling to jam the liquid into an aperiodic solid before it could crystallize; the resulting lattice does not correspond to an equilibrated liquid. Alternatively, one may envision generating an aperiodic scaffold structure by equilibrating a liquid that does not readily crystallize, such as the venerable Kob-Andersen binary mixture.⁸⁰ In this case, the eventual octahedrally bonded system could be regarded as equilibrated in a limited sense, viz., with regard to density fluctuations at half-wavelengths exceeding the particle spacing in the scaffold structure. A connection can be made with earlier methodologies, in the work of Torquato, Stillinger et al.,^81,82 in which the spectral content of density fluctuations is controlled selectively during the generation of aperiodic structures.

There are two sources of ambiguity in identifying cavities in a random lattice, as we illustrate in Figs. 3 and 4 using a two-dimensional lattice for clarity: On the one hand, one must decide on the rank of the cavity. Clearly the top-left situation in Fig. 3 corresponds to two triangular cavities, while the top-right shows a square cavity. To classify an intermediate situation, with respect to the compactness of the cavities involved, one must adopt a specific prescription. We will adhere to the following convention: We consider the cavity shown in red in Fig. 3 to be quadrangular if the area A of the shared face of the two Voronoi cells in contact is less than a certain threshold value. This quadrangular cavity can be thought of as resulting from merging the two triangular cavities separated by the line connecting the host particles of those two Voronoi cells; the line is perpendicular to the shared face. For brevity, we shall call such connecting lines Delaunay connections. The other source of ambiguity in cavity identification is illustrated in Fig. 4. The quadrangles ABCE and ACDE both represent acceptable choices for quadrangular cavities, if the lengths xy and yz are both less than some preset threshold value. Yet these two choices are mutually exclusive and lead to distinct tilings of the sample by the cavities. The ensuing ambiguity is inconsequential only if the distinct ways of tiling are statistically equivalent.

Hence we adopt the following two-step strategy to identify cavities and tile the sample: In the first step, we perform Delaunay triangulation of the scaffold structure (the red lines in Fig. 4), which is complementary to the Voronoi partition; both are generally unique for disordered lattices. In the pertinent case of three spatial dimensions, Delaunay triangulation produces a unique set of the most compact tetrahedra that tile space. In the second step, we single out a subset of face-sharing tetrahedra that are relatively non-compact and can be merged to create higher order cavities that are relatively compact. To do so, we remove each Delaunay connection corresponding with a face shared by two Voronoi cells, if the area of the face is less than a certain threshold value denoted with A. We will designate the latter condition as the area rule. As a result, no fewer than three tetrahedral cavities will be merged upon the removal of a single connection, resulting in the appearance of a hexahedral or higher order cavity. At the same time, the number of tetrahedral cavities decreases. If the merged cavities are exclusively tetrahedral, the result of the merger is a higher-order cavity that is convex. Conversely, the removal of edges in the presence of non-tetrahedral cavities may result in the appearance of non-convex cavities, and, furthermore, in a percolation of cavities. To prevent the latter, we take the edges of non-tetrahedral cavities off the list of the removable Delaunay connections. We will designate the latter restriction as the percolation rule.

As already mentioned, the tiling pattern is generally not unique and depends on the precise order in which the pertinent Delaunay connections are removed. The simplest strategy, which will be employed below, is to remove the connections randomly, subject to the area and percolation rule. The statistics of the resulting tiling are summarized in Fig. 5, where we show the partial quantity of the cavities of distinct order as functions of the threshold area A. Incidentally, we observe the results are not very sensitive to the system size. One may imagine an alternative, entirely deterministic strategy, in which the Voronoi faces are first ranked by area and then the corresponding Delaunay connections are removed in ascending order of the areas (subject to the percolation rule). This procedure produces smoother but qualitatively similar curves, see Fig. S3 in the supplementary material. Presumably the most systematic way to assign the priority for removal of Delaunay connections would be via a smooth cost function such as the Boltzmann distribution with judiciously chosen parameters.

We finish by pointing out that one could conceive many alternatives to the area rule. For instance, one may choose to remove Delaunay connections that are longer than a certain threshold value L, see Fig. 3. Although this alternative strategy will not be pursued here, we do provide the pertinent analogs of Fig. 5 and Fig. S3 in the supplementary material as Figs. S4 and S5, for the reader’s reference.

To create an octahedrally coordinated parent structure, we now place additional–or secondary-lattice–sites inside non-tetrahedral cavities in the scaffold structure, at most one secondary site per cavity, as in Fig. 1(b). For concreteness, we place each secondary site in the respective cavity’s center of mass, the latter determined with respect to the corners of the cavity. As part of the procedure, we decide on the chemical identity of the atoms that will be assigned to each site. If necessary for maintaining the stoichiometry, vacancies in either the scaffold or secondary lattice will be created. After both lattices are finalized, an atom is assigned to each lattice site, and the lattice spacing is adjusted to achieve the desired density, the resulting structure is called the parent structure. While possessing chemically inspired structural motifs, this structure is not optimized with respect to the electronic energy; the latter step will be detailed in Subsection II D.

Zhugayevych and Lubchenko⁴³ have noted that the contacts in the parent structures for the crystalline structure of the stoichiometric compound As₂Se₃ that would give rise to covalent bonds in the eventual, distorted structure are exclusively hetero-nuclear. In keeping with this notion—which is consistent with Pauling’s rules,⁸³ of course—we adopt the following two conventions: In one, the scaffold lattice is populated exclusively with chalcogens and the secondary lattice with pnictogens. The resulting structure will be referred to as C-type. In the other convention, the scaffold and secondary lattices are occupied with pnictogens and chalcogens, respectively; the corresponding structure is called P-type.

Consider a C-type structure and the stoichiometric compound As₄₀Se₆₀ for concreteness. To build such a structure, we first identify non-tetrahedral cavities according to Subsection II B. If the number of cavities is less than two-thirds the number of scaffold lattice sites, the secondary lattice will not have the requisite number of sites to achieve the As_0.4Se_0.6 stoichiometry. Thus we create vacancies in the scaffold structure, so as to remove extra chalcogens as needed to achieve the stoichiometry. In the opposite situation when the number of non-tetrahedral cavities exceeds the number of sites needed to accommodate the pnictogens, we create vacancies in the secondary lattice.

Every time the algorithm calls for the creation of vacancies, one must make a decision as to their spatial distribution. Here we adopt the simplest strategy by choosing vacancy sites randomly. Additional optimization could be performed, such as with respect to spatial uniformity of the density. Such optimization, exemplified in the supplementary material, leads to a modest increase in hetero-nuclear covalent contacts.

The resulting statistics of sites filled with chalcogens, pnictogens, and the number of created vacancies are shown in Fig. 6(a), all as functions of the threshold area A. Here we observe that for A numerically less than 0.1 or so, there are indeed not enough sites to accommodate the pnictogens and so vacancies in the scaffold lattice must be made. Beyond that threshold value, there are too many potential sites for pnictogens, which effectively means one must introduce vacancies in the secondary lattice. Note that near the latter threshold value of A, the partial quantity of octahedral cavities reaches its maximum value, see Fig. 5. Most of the data presented here pertain to A = 0.1. Other values of A yield very similar structures, at least at the level of pairwise correlation functions, see the supplementary material. There, we also show that the parent structures are indeed octahedrally coordinated.

The vacancy statistics for the P-type structures are shown in Fig. 6(b). Here we observe that for all values of A, vacancies must be made in the scaffold structure. Despite this difference between the C- and P-type structures, we find that the pairwise correlation functions are not overly sensitive to the structure type, see Figs. S18 and S19 of the supplementary material. Henceforth, we shall focus exclusively on C-type structures.

To perform geometric optimization, we must assign chemically appropriate coordinate values to the lattice sites. We do so by matching the mass density of each sample to the experimentally available mass density^84–89 of each compound near its glass transition. Note that the reported density varies. For instance, the figures pertaining to the stoichiometric composition As₄₀Se₆₀ are distributed between 4.455 and 4.62 g/cc.^85–89 Two types of geometric optimization will be adopted: In one, labeled “const-P,” we allow the size of the unit cell to optimize at constant pressure. In the other, labeled “const-V,” the volume is fixed.

For geometric optimization, we employed two quantum-chemical approximations: One is plane-wave DFT, as implemented in packages VASP (Vienna Ab Initio Simulation Package)^90–93 and CP2K/QUICKSTEP.^94,95 VASP utilizes projector augmented-wave potentials with plain wave basis set. The Perdew-Wang (PW91^96–100) generalized gradient approximation was used for exchange-correlation functions. Geometric optimization was performed using the conjugate-gradient algorithm. VASP allows for unit cell optimization. The CP2K/QUICKSTEP⁹⁴ utilizes a Gaussian basis set¹⁰¹ and Goedecker, Teter, and Hutter (GTH) pseudopotentials.^102–104 The exchange-correlation functions are computed using Padé approximants. Geometric optimization is performed using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) procedure.^105–108 CP2K does not support unit-cell optimization, therefore CP2K simulations were performed at several, fixed values of the density. We have found that VASP yields secondary bonding that is more consistent with experiment than CP2K, see the supplementary material. Most of our DFT-based optimization has been performed using VASP.

A separate set of optimizations were performed using the semiempirical quantum chemistry package MOPAC with PM6 parametrization,^109,110 which employs Dewar and Thiel’s NDDO approximation.¹¹¹ Like CP2K, it utilizes the BFGS procedure for structure optimization. MOPAC’s built-in parameters tend to significantly overestimate the length of the secondary bond, according to our tests on crystalline As₂Se₃. With the default parameters, MOPAC-optimized structures tend to exhibit excessive tetrahedral coordination within the double layers and an artificially large spacing between the layers. We have obtained an alternative, re-trained set of parameters, using numerical optimization¹¹² that produces significantly improved structures, crystalline and amorphous alike; there are some structural artifacts in the amorphous samples however, see Sec. IV. Re-parametrization of MOPAC is detailed in the supplementary material.

Figures 7 and 8 show the evolution of the pairwise correlation function g(r) and structure factor S(q) with the progress of geometric optimization, for an individual realization of the structure. Individual realizations all use the same scaffold structure, while the cavities and atom placement vary from realization to realization. [Note that to evaluate the quantities g(r) and S(q), we properly smear the atoms’ positions to account for vibrations,^114,115 see the supplementary material for details.] The quantity N gives the number of steps the optimization algorithm has performed. The optimization was performed at a constant pressure and resulted in some reduction in density.

The most noticeable feature in the time-series of g(r) is the progressive decrease in the nearest and next nearest neighbor spacings toward lower values, and significantly more so for the former than the latter. This process is straightforwardly understood from the intrinsic instability of approximately linear motifs in octahedral bonded $ppσ$ networks toward dimerization,^43,44,116 see the illustration in Fig. 9. This instability, in turn, can be understood at the simplest level as arising from the Peierls instability of linear chains at half filling.^66,117–119 We note that the position of the first peak matches very well with the length of the covalent As-Se bond, viz., 2.4 Å, while the position of the second peak matches the length of the secondary bond, both homo- and heteronuclear, in crystalline As₂Se₃. Indeed, the latter bonds range between 3.3 and 4 Å in length, see Table S2 in the supplementary material. There, we also exemplify a pyramidal unit in the distorted structure.

In contrast with g(r), whose gross features appear already in the parent structure, the structure factor S(q) of the parent structure almost appears structureless. Still, it is seen to recover the experimental curve near quantitatively once optimization is completed.

Next we compare our fully optimized structures with the experiment and, also, with results obtained relatively recently using first principles molecular dynamics (FPMD), from the work of Bauchy, Kachmar, and Micoulaut⁵³ (BKM). Figures 10 and 11 pertain to g(r) and S(q), respectively. The present results are shown for one specific realization of a parent structure and data averaged over ten realizations. For the reader’s reference, we provide the curves for individual realizations in the supplementary material. We observe that, overall, the present algorithm and FPMD alike fare well in comparison with the experimental data obtained near the laboratory glass transition. (For As₄₀Se₆₀, $Tg≃466$ K.)

In the BKM study,⁵³ the sample is equilibrated but at a high temperature of 1200 K—which much exceeds the laboratory T_g—and then quenched to a temperature of interest. The present data are obtained using an effectively hybrid protocol, of which some steps correspond to effectively zero temperature (generation of the scaffold structure and the geometrical optimization) and some to infinite temperature (cavity identification and atom placement). Despite their substantial size, individual samples show many fine features that are not present in the actual bulk material. Averaging over already a few realizations produces rather smooth curves.

The most noteworthy difference between the present results and those of BKM⁵³ is in the vicinity of the so-called first sharp diffraction peak (FSDP) in the structure factor S(q). The latter peak is observed near q = 1.3 Å in As₄₀Se₆₀ (it is indicated by the vertical line in Figs. 8 and 11); the position varies somewhat with stoichiometry. The FSDP corresponds to the smallest non-trivial supra-atomic lengthscale; one’s ability to reproduce the peak thus serves as a rather stringent criterion for the ability of a microscopic model to reproduce bonding beyond the immediate coordination shell. Renninger and Averbach’s⁶³ experiment revealed that the strength of the FSDP positively correlates with the partial amount of arsenic, within the 10%-50% range of arsenic content. BKM have simulated several stoichiometries within the latter range. They have recovered the FSDP in the total structure factor at 35% As content (As₃₅Se₆₅), but not other concentrations, where, instead, they saw a soft shoulder or no discernible features at all. Nevertheless, the latter authors consistently observe a peak in the partial structure factor $SAsSe(q)$ at all probed values of arsenic concentration, as in Fig. 15. BKM’s As₄₀Se₆₀ data are shown in Fig. 11.

We have tested the present structure-building algorithm for several stoichiometries, viz., corresponding to 20%, 35%, 40%, and 50% arsenic content, respectively. The results of averaging over 5 realizations, superimposed on experimental data, are shown in Fig. 12. (The As₄₀Se₆₀ data are averaged over 10 realizations.) We observe that the FSDP is reproduced robustly for all tested stoichiometries. Its position, as a function of arsenic content, is shown in Fig. 13; it is seen to follow the experimental trend semi-quantitatively.

With the exception of the most arsenic-rich sample, As₅₀Se₅₀, the FSDP is quantitatively similar to its experimental counterpart not only in position but also in strength. In particular, the strength of the peak does positively correlate with arsenic content, except for the As₅₀Se₅₀ stoichiometry. To avoid confusion, we reiterate that individual realizations exhibit significant variation in both position and strength of the peak, see the supplementary material.

Since the BKM simulations were performed at constant volume, it seems appropriate to test whether some of the difference between the predictions of that and the present study stems from our allowing the sample size to relax. The results of geometric optimization at constant volume, shown in Fig. 14, indicate that the FSDP is not recovered nearly as robustly when the unit cell is not allowed to relax. In addition, the latter figure shows that the FSDP is consistently shifted toward higher qs relative to its experimental value. We thus observe that on the one hand, the FSDP is surprisingly sensitive to the density variations. On the other hand, the results in Fig. 14 suggest that the relative weakness of the FSDP in the first principles MD could be at least in part due to the latter simulations being performed at constant volume.

Yet there seem to be more factors at play than the density alone. In Fig. 15, we compare the partial structure factors for the present structures, as produced both at constant pressure and volume, and the corresponding structure factors obtained by BKM. To set the stage, we point out that the three partial structure factors contribute unequally to the total structure factor.^5,122 Specifically at 40% arsenic content, $S(q)=0.416⋅SSe, Se(q)+0.458⋅SAs, Se(q)+0.126⋅SAs, As(q)$⁠, implying the two main contributions to the total structure factor are from Se–Se and As–Se pairwise correlations. First of all, the data in Fig. 15 indicate that the structures obtained by geometric optimization at constant volume and pressure, respectively, are different quantitatively, not qualitatively. There is some divergence in the strength and position of the peaks but not their identity. In contrast, there is a significant difference between the present and FPMD-generated samples in that the former exhibit prominent Se–Se and As–As correlations in the spectral region around the FSDP. Furthermore, the present As–Se partial structure factor exhibits a relatively broad peak with a dip at the wavevector of the FSDP. In contrast, the FPMD-produced As–Se partial structure factor is the only one exhibiting a peak at that wavevector.

To elucidate the apparent As–As and Se–Se correlations on supra atomic scales, we next examine the partial pairwise correlation functions in direct space. These are shown in Fig. 16 alongside the corresponding data from Ref. 53. The presently generated samples clearly exhibit significantly fewer homonuclear covalent bonds than FPMD and a commensurately greater amount of hetero-nuclear covalent bonds. Conversely, the present structures exhibit substantially more secondary bonding between like atoms than the FPMD data. The amounts of hetero-nuclear secondary bonds are comparable between the two methodologies; however, the present structures exhibit heteronuclear secondary bonds that are shorter and thus stronger than those produced in FPMD. Both methodologies exhibit an excess of bonds at 2.7-2.8 Å, relative to experimental structures. This excess is apparently due to hetero- and homo-nuclear contacts for the present and FPMD methodology, respectively. It seems likely that this artifact is due to unrelaxed regions of excessive density that are structural relics of the high temperature portions of the methodologies. We remind the reader that the secondary bond decreases in length at higher density.^43,54,56,116

The relative preponderance of heteronuclear covalent bonding in our structures is not surprising: The parent structures are constructed so as to have as few homonuclear contacts as possible, subject to the stoichiometry. Importantly, the optimized structures exhibit even fewer homonuclear bonds than the parent structures, as can be seen in Fig. 17(a). The latter trend is more pronounced when extra effort is made to make the secondary lattice more spatially uniform by placing the vacancies more homogeneously in space, see Fig. 17(b). This apparent propensity for making heteronuclear covalent bonds becomes somewhat weaker for large values of the threshold area A, which correspond to relatively large numbers of vacancies in the parent structure and relative facility for particle transport. Data on particle displacement during optimization can be found in the supplementary material. There, we also demonstrate the converse notion that the movement is partially suppressed, if geometric optimization is performed at constant volume.

Now, having exclusively heteronuclear covalent bonds is unfavorable entropically at finite temperatures, and more so, the higher the temperature; the relatively large number of homonuclear bonds in FPMD-produced samples is consistent with the latter samples obtained by quenching a melt equilibrated at very high temperatures. According to the experiment, effects due to the mixing entropy are essentially negligible by the glass transition temperature. Indeed, although the presence of some As–As homopolar bonds in the stoichiometric compound As_0.4Se_0.6 has been inferred using high-resolution XPS studies,⁴⁹ their quantity is very small according to recent NMR experiments.⁶⁸ 

Let us now compile a partial summary of FSDP-related notions. First, the partial structure factor in the present study exhibits a peak at the FSDP wave-vector in the homonuclear partial structure factors but not the heteronuclear one. Second, the FSDP is at a length $≃7.7/1.3≈5.9$⁠, which clearly exceeds the covalent bond length and there are few homonuclear covalent bonds in the first place. [The numerical constant 7.7 comes about because the first $q≠0$ peak of the function $sin(qr)/q$ is at $q≈7.7/r$⁠, see Eq. (S3) of the supplementary material and Ref. 31.] Third, the time progression in Figs. 7 and 8 indicates that the FSDP reaches its steady-state shape and position only by the time the number of secondary bonds reaches its steady state value; the comparison of snapshots at N = 10 and N = 50 is particularly revealing. Note that the secondary bonds emerge in part at the expense of the covalent bonds, as can be seen in two ways: The bonds comprising the right flank of the first peak in the g(r) move toward the second peak, while the “progenitor” peak of the eventual FSDP moves toward shorter wavevectors. Thus we conclude that the FSDP is mostly due to a pattern formed by homonuclear secondary bonds. The characteristic lengthscale of the pattern matches well the sum of the covalent and secondary bonds. Thus the pattern is formed by the regions housing regions of relatively weak, secondary bonding that are separated by regions of relatively strong, covalent bonding. As mentioned in Sec. I, the pattern thus represents a charge density wave (CDW). The latter can be visualized directly by graphing the density of valence electrons as a function of the spatial coordinate. In Fig. 18(a), we show a 2D slice of the CDW; panel (b) shows the charge distribution for crystalline As₂Se₃, for the reader’s reference. In the supplementary material, we show an alternative way to visualize spatial localization of valence electrons, viz., using the so-called electron localization function (ELF),^123,124 which has been used to locate metallic regions in solids. The ELF-produced patterns echo those in Fig. 18.

The microscopic identification above the pattern underlying the FSDP is consonant with several aspects of the early analysis by Elliott,^32,36,67 to be discussed in some detail in Sec. V. It is also consistent, in a sense, with the BKM data. In the latter data, the FSDP is due to the heteronuclear component of S(q) (Fig. 15), consistent with the relatively large number of heteronuclear secondary bonds (Fig. 16). The latter come about because of the relatively large number of homonuclear covalent bonds. It is remarkable that the total structure factors produced in that and the present study are so similar despite the significant difference in local stoichiometry. In any event, both studies point at the pattern formed by weak, secondary bonding as causing the FSDP.

Here we present a complementary perspective on local bonding as afforded by analyzing the bond angles. The latter perspective is particularly useful in that it allows one to connect relatively extended structural motifs to hybridization of electronic wavefunctions with distinct symmetries on individual atoms. The semi-empirical package MOPAC is known to overestimate the extent of sp-mixing in the chalcogenides, which results in a significant departure of MOPAC-optimized structures from distorted octahedral bonding, toward tetrahedral coordination.⁴³ As a byproduct of the bias toward sp³ hybridization, there is a decrease in the number and strength of secondary bonds because of steric hindrance.

Here we take advantage of this (otherwise unwanted) complication to investigate effects of tetrahedral coordination on the pairwise correlation functions, by analyzing structures optimized using MOPAC. To minimize direct effects of excessive sp³ hybridization, if any, on the length of the secondary bond, we first modify some of MOPAC’s parameters to maximally reproduce the structure of crystalline As₂Se₃ and, in particular, the bond lengths for a representative subset of the secondary bonds. (The default MOPAC parameters do lead to secondary bonds that are too long, see the supplementary material for details.) After the re-parametrization, MOPAC produces crystal structures that are hardly distinguishable from those determined by X-ray diffraction. Nonetheless, MOPAC-optimized disordered parent structures still exhibit secondary bonding that is on average weaker than what is suggested by the experimental g(r). Indeed, Fig. 19(a) indicates that MOPAC overestimates the length of secondary bonds by 0.2 Å on average. In addition to the latter feature, g(r) in Fig. 19(a) also exhibits a small peak at bond lengths less than that expected for the nearest neighbor bond, either homo- or heteronuclear. This small peak is an artifact stemming from our use of modified parameters. Structures optimized using default parameters do not exhibit this artefactual feature; at the same time, they exhibit a covalent bond that is significantly longer than that found in the experiment. In addition, the discrepancy in the secondary bond length exceeds that found for the re-trained set of parameters from Fig. 19, see the supplementary material.

Consistent with the discussion in the beginning of this section, the weakening of secondary bonding apparent from Fig. 19(a) is positively correlated with increased tetrahedrality. Indeed, according to Fig. 20, the bond angles in the MOPAC-optimized structures center around $110°$⁠, a value symptomatic of tetrahedral coordination. In contrast, the bond-angles in structures produced using the DFT calculations are consistently found to have smaller values that are appropriate for distorted-octahedral bonding. As a concrete reference for assessing the degree of octahedrality, we will use the bond angles in crystalline As₂Se₃. There are two inequivalent arsenics in the crystal, each of which makes three covalent bonds with its neighbors. The corresponding sets of bond angles, sorted in ascending order, are¹²⁵ 

(The label “X” signifies crystal or Xtal.) There are also three inequivalent seleniums, each making two covalent bonds with their neighbors where the bond angles are

We observe that VASP-produced bond angles in disordered samples are consistent with the angle values listed in Eqs. (1) and (2). Incidentally, we point out that VASP-optimized structures for crystalline As₂Se₃ are very close to those determined using X-ray crystallography; see the supplementary material for further detail. In any event, the position of the second peak in g(r) for VASP and MOPAC optimized structures, Figs. 10 and 19(a), suggests the VASP-derived lengths for secondary bonds match experiment significantly better.

To further quantify the above discussion, we note that each triplet of Se–As–Se bond-angles in Eq. (1) can be thought of as a three-dimensional vector, call them $θ→1,X$ and $θ→2,X$⁠, respectively. We next identify each pyramidal unit in the disordered structure whose apex hosts an arsenic atom covalently bonded to three atoms, irrespective of their identity; calculate the three respective bond angles; and sort them in ascending order. We call the resulting “vector” $θ→As$ and compute the lengths of the deviations $θ→As−θ→i,XAs$ (i = 1, 2). The smallest of the two,

is then used to measure deviation of local coordination in the disordered structure from that in the crystal. The analogous quantity for the bond angle on selenium atoms is computed according to the following equation:

where i = 1, 2, 3. The quantities from Eqs. (3) and (4) are histogrammed in Fig. 21. Again, we observe that DFT-optimized structures exhibit angles that are tightly distributed in values around their analogs in the crystal structure.

One does not fail to notice a small but discernible amount of $60°$ bond angles in Fig. 20 in DFT-optimized structures, also seen as a feature near $20°$ in Fig. 21. Since the latter figures report angles exclusively between covalent bonds, the $60°$ bond angle automatically implies a presence of triangular units, each unit contributing three incidences to the histograms in Figs. 20 and 21. Such triangular units also have been reported previously in first principles MD performed at high temperatures.⁵⁹ Such triangular units are expected to be scarce in actual materials near the laboratory glass transition because they necessarily imply the presence of homonuclear contacts; the latter are indeed rare according to recent experimental studies.⁶⁸ Note that triangular units are virtually absent in structures obtained using MOPAC, again consistent with its propensity for open structures.

Perhaps the most contrasting view of the DFT- and MOPAC-optimized samples is afforded by the statistics of closed rings.¹²⁶ Here we analyze the latter statistics using the computational tool R.I.N.G.S.,¹²⁷ by counting the number of closed rings of size 15 or less, where each bond distance is 2.8 Å or less. The most conspicuous distinction is that MOPAC-produced structures exhibit many more rings and thus are much more connected than the DFT-optimized samples, see Fig. 22. In addition, the former are dominated by even-numbered rings, while in the latter, even- and odd-numbered rings occur with comparable likelihood. Note that the presence of odd-numbered rings in our structures necessarily implies some presence of homonuclear bonds, consistent with findings from Sec. III.

The relative dearth of connectivity in DFT-optimized structures indicates that the latter do not exhibit significant long-range ordering such as the double layering so characteristic of crystalline As₂Se₃, see also Fig. 18. This is despite the coordination not deviating too much from that found in a crystalline material, as witnessed by the statistics of nearest neighbors displayed in Fig. 23. According to the latter figure, the majority of arsenics are coordinated with three seleniums, while most seleniums are coordinated with two arsenics. Our calculations show that, on average, arsenics have 3.0 neighbors and seleniums 2.0, consistent with existing structural models.^49–52 

Conversely, it might appear at a first glance that the extensive connectivity of the MOPAC-optimized samples is analogous to that exhibited by tetrahedrally bonded materials such as diamond. Consistent with this expectation, a substantial number of atoms have four covalently bonded neighbors, including seleniums that should normally have only two, see Fig. 24. On average, arsenics have 3.5 neighbors and seleniums 2.5. Yet the consequences of these vestiges of tetrahedrality turn out to be only limited. On one hand, there are plenty of 4-member rings in the MOPAC-generated structures. On the other hand, the six member rings do not display the bond angles characteristic of traditional tetrahedrally bonded networks, as can be seen in Fig. 25 where we show representative 6-member rings for both DFT and MOPAC-optimized samples. Despite exhibiting relatively open angles that are compatible with sp³ hybridization, the MOPAC-produced 6-rings do not show the boat and chair conformation symptomatic of a sp³ hybridized ring-molecule such as cyclohexane.

Finally, we address the effects of tetrahedrality on the structure factor S(q). The latter quantity obtained for the MOPAC-optimized sample is shown in Fig. 19(b). The structure factor exhibits an across-the-board shift toward shorter wavevectors relative to the experimental data even though the density does not change on average following optimization. The leftward shift of the second and higher order peaks seems to be fully accounted for by the increase in the secondary bond length, relative to actual structures. Note that the latter length is distributed around 3.8 Å or so, which corresponds to the spectral region around $q≈7.7/3.8≈2.0$ Å⁻¹. It is interesting that MOPAC-generated structures exhibit a double peak in the spectral region corresponding to the FSDP that is not unlike that in the under-optimized structure from Fig. 8 at N = 50, with an important distinction that the r.h.s. peak is more of a shoulder shifted to higher qs relative to the experimental position of the FSDP. These observations may well be explained by a decrease in the differentiation between covalent and secondary bonds, as in Fig. 9(b), that would be expected owing to the increased tetrahedrality. In any event, the above discussion is consistent with the earlier, general notion that increased tetrahedrality reduces the amount of secondary bonding.

We have designed and implemented a methodology for generating structures of amorphous solids characterized by distorted octahedral coordination. The present framework allows one to partially decouple the establishment of local density, coordination, network topology, and chemical ordering in such solids. This is in contrast with ab initio based methods, in which all of these aspects of the structure are generated simultaneously—and self-consistently—by equilibrating a collection of atoms subject to quantum-chemical force fields and the Boltzmann distribution. Our motivations are twofold: On one hand, present day computer technology only allows one to equilibrate very modestly sized sample at a temperature that is very much higher than the laboratory glass transition temperature T_g. The resulting samples, when quenched to T_g, likely contain bonding motifs that are uncharacteristic of samples equilibrated near T_g. For instance, one may expect the so-obtained samples of even stoichiometric compounds, such as As_0.4Se_0.6, to host an excess of homonuclear covalent bonds. The latter are entropically favored at high temperatures but should be quite rare at low temperatures, consistent with experimental evidence.⁶⁸ On the other hand, the ability to separately manipulate those distinct aspects of structure formation may help one directly test microscopic proposals for the origin of specific observables such as the first sharp diffraction peak (FSDP).

Each step of the algorithm, Fig. 2, has a stochastic aspect that can be managed probabilistically using a cost function. We anticipate that a cost function can be designed, for each step, that corresponds to the Boltzmann distribution at a specified temperature. In the current implementation, some steps are performed at either effectively zero or infinite temperature. The former set covers the generation of the close-packed scaffold structure, homo-vs. hetero-nuclear assignment of nearest neighbor bonds in the parent structure, and the geometric optimization. The “infinite temperature” set covers the identification of compact non-tetrahedral cavities and the assignment of vacancies in the lattice, when it is needed to satisfy the stoichiometry. We have found that increasing selectivity in these infinite T steps does not qualitatively change the outcome of the overall procedure.

When applied to a technologically important family of chalcogenide alloys As_xSe_1−x, the present methodology produces structures that are similar to those produced⁵³ using first-principles molecular dynamics (FPMD) in some respects, but significantly different in others. At the level of the total pairwise correlation function g(r), the two methods fare comparably well against the experiment; the respective functions are remarkably similar, in fact. The same thing can be said about the statistics of bond angles and shortest rings. The distinctions between the outcomes of the two methodologies thus appear to be less with respect to the overall topology of the network than to the chemical identity of the nodes of the network. (The two aspects are often referred in the literature as the topological and chemical ordering, respectively.³⁵) The presently generated structures have fewer homonuclear covalent bonds and commensurately more homonuclear secondary bonds. Both approaches indicate that the puzzling first sharp diffraction peak (FSDP) is due to the aperiodic modulation pattern formed by regions of weak, secondary bonding that are separated by regions housing covalent interactions. In significant distinction, the present studies indicate the spatial correlation underlying the FSDP is largely that between like atoms, even in stoichiometric compounds. In addition, the FSDP is recovered consistently at all tested values of arsenic content in the As_xSe_1−x compounds and so is the apparent tendency of the peak’s position to decrease with As content. The FPMD studies seem to lack such consistency by often producing either a shoulder or no discernible features in the structure factor at the wave lengths in question.

The present method of, essentially, engineering structures so as to minimize the number of homonuclear covalent bonds could be perceived as unconventional but it achieves the intended goal of partially circumventing the intrinsic limitation of ab initio methods to exceedingly high temperatures, where the mixing free energy is much too large compared with its value near the laboratory T_g. Likewise, we use an admittedly unconventional device of using the known bias of the semi-empirical quantum-chemistry package MOPAC toward tetrahedrality to artificially increase the extent of tetrahedral bonding and monitor the resulting effects. The so generated network is significantly more connected while the atoms, as expected, exhibit significantly more sp³ hybridization. The resulting weakening of the secondary bonding leads to partial suppression of the FSDP, consistent with preceding observations. The overall picture emerges in which the chalcogenides structurally interpolate between octahedrally and tetrahedrally bonded network but in a rather non-obvious way so as to give rise to an aperiodic modulation pattern of relatively strongly and weakly bonded regions, as in Fig. 18. An essential feature of this interpolation is that it exhibits a spatially non-uniform bond saturation. In contrast, the extreme cases of strictly octahedral and tetrahedral bonding—exemplified by simple-cubic (pressurised) arsenic and diamond—exhibit strictly uniform bond strength throughout.

The appearance of such a modulation pattern of bonding strength was anticipated by Zhugayevych and Lubchenko,^43,44,116 who argued that crystalline and amorphous $ppσ$-bonded networks can be thought of as symmetry-broken versions of relatively symmetric structures locally defined on the simple cubic lattice. These structures consist of intersecting chains at half-filling and thus are Peierls unstable toward dimerization. The resulting alternation pattern of covalent and secondary bonds, as in Fig. 9, gives rise to a disordered checker-board-like pattern of regions housing relatively strong and weak bonding from Fig. 18. This pattern—as formed by regions of, respectively, relatively high and low density of valence electrons—is often referred to as “bond-order” or “off-site” charge density wave (CDW).⁵⁴ In terms of the electronic Hamiltonian, it corresponds, respectively, to high and low values of the hopping matrix element and the corresponding entry of the electronic density matrix. Thus we argue that the first sharp diffraction peak in the chalcogenides is due to a bond-order, aperiodic CDW. We have seen that the results from the work of Bauchy et al.⁵³ could be interpreted as suggesting the FSDP arises not from a CDW in which those crests and troughs reside, instead, on anions and cations. This, “on-site” type of CDW could result, for instance, from a spatial variation of electronegativity as in ionic insulators. Both off-site and on-site CDW are present in the chalcogenides, see further discussion in Ref. 54. The CDW observed in the presently generated samples apparently does not have the “stripe” morphology seen in crystalline chalcogenides. Note that aside of the precise mechanism of the instability leading to the formation of the CDW, the eventual geometry is affected by the precise amount of sp-mixing.¹²⁸ 

The CDW-based microscopic picture above is consistent with the temperature and pressure dependence of the FSDP in the chalcogenide glasses. First note that the semiconducting gap is due to both the on-site and off-site components of the charge density wave. The simplest way to appreciate this is to write the following Hückel Hamiltonian for a one-dimensional chain, one orbital per site, the on-site energy, and transfer integral both exhibiting a perfect alternation pattern:

where $cn,s†$ (c_n,s) creates (annihilates) an electron on site n with spin s, and

is the hopping matrix element between sites n and n + 1, and $2ε$ is equal to the electronegativity variation. It is straightforward to show¹²⁹ that at half filling, the chain in Eq. (5) is an insulator with a gap

The hopping matrix element t_n,n+1 decreases roughly exponentially with the distance between sites n and (n + 1), and so does the corresponding entry of the electronic density matrix,

where $ψkm$ are the expansion coefficients of the normalized kth solution of the Schrödinger equation corresponding to Eq. (5). The strength of the off-site CDW is reflected by the differential of the amount of valence electrons between even and odd sites,

Now suppose for concreteness that the even n correspond to the weaker, secondary bond: $t2<t1$⁠. Heating will preferentially lead to the expansion of the more-weakly bonded, softer regions; this will lead to a further decrease in t₂ and, hence, a decrease in the density matrix entry P_2n,2n−1, but less so in P_2n+1,2n. The strength $Δρ(off)$ of the CDW will thus increase, consistent with the apparent increase of the FSDP with temperature. Conversely, application of pressure will preferentially lead to the compression of the more weakly bonded regions thus leading to a decrease in the strength of the CDW and, hence, the FSDP. These notions are buttressed by the present results: On one hand, we have seen that the FSDP is partially suppressed when geometric optimization is performed at constant volume. On the other hand, such optimization also results in secondary bonds that are somewhat shorter than those obtained at constant pressure, see the supplementary material.

Elliott^32,36,67 has proposed a general mechanism by which the FSDP stems from chemical ordering of interstices around cation-centered clusters that was built on an earlier proposal by Blètry¹³⁰ that certain amorphous materials can be thought of as mixtures of spherical particles and voids. Although it appears that there is generally no strict correspondence between chemical ordering³³—in the Bhatia-Thornton sense¹³¹—and the FSDP, the central aspect of Elliott’s microscopic proposal can be connected with the present results insofar as one is willing to associate Elliott’s interstices with regions housing the weaker, secondary bonds.

The phenomenon of reversible photodarkening²⁴ is observed in well annealed chalcogenide films, whereby soaking the samples for hours with intense light above the gap frequency leads to substantial changes in how the material absorbs light. The apparent changes are lasting and can be reproducibly “erased” by thermal annealing. The changes are two-fold: On one hand, the material begins to absorb at frequencies as low as half-gap. On the other hand, the gross absorption edge becomes red-shifted, by 10% or so for macroscopic amounts of absorbed photons. The shift can be seen clearly after the midgap absorption is saturated. It is this overall shift that is usually referred to as “photodarkening.” Zhugayevych and Lubchenko^44,116 (ZL) have argued that the midgap absorption is due to special midgap states that are an intrinsic consequence of the structural degeneracy of glass-forming liquids and are tied to the structural motifs corresponding to transition state configurations for activated transitions between the distinct free energy minima. These special midgap states are of topological origin; they are filled and absorb at supra-gap frequencies in pristine samples but become singly occupied and absorb at frequencies much below the gross optical gap, upon light soaking. The origin of photodarkening, on the other hand, is less clear.

Photodarkening is interesting in the present context because it is known to be negatively correlated with the strength FSDP. This is puzzling because photodarkening is accompanied by a modest increase in volume, which seems to contradict the logic from the preceding discussion: On one hand, one would seem to be warranted in associating an increase in volume with weakening of the secondary bonds and an increase of the optical gap, according to Eq. (7). On the other hand, heating, which is accompanied by some expansion, leads to an increase in the FSDP intensity, as just discussed.

This intriguing behavior can be rationalized within the present microscopic picture. First we note that the T-dependent FSDP studies and the light-soaking experiments are distinct in a key way. In the former, one changes the ambient temperature and, ordinarily, maintains it only for a time necessary to perform X-ray diffractometry. Some aging will undoubtedly occur—especially closer to T_g—but not very much; the sample remains largely arrested in the same state it was residing in following vitrification. (In any event, aging and its effects on the FSDP, if any, could be tested for reversibility with regard to volume changes by temperature cycling and by increasing the wait time at the target temperature of the experiment.) In contrast, exposing a sample to intense light so as to have it absorb macroscopic quantities of photons at a gap frequency or above¹³² clearly leads to structural changes—including the aforementioned volume increase—that last for apparently long times unless the sample is annealed. Such long-lasting changes directly indicate the sample has indeed aged toward a different structure as a result of irradiation.

To gain qualitative insight as to those structural changes we first recall that according to the RFOT theory,^3,4 glass ages via local relaxation of a quenched structure toward the state that would be equilibrium at the ambient temperature/pressure; the driving force has a contribution due to the enthalpy difference between the frozen-in and equilibrated state and a contribution due to the configurational entropy.¹³⁷ The latter contribution vanishes if the ambient temperature is below the (putative) Kauzmann temperature.¹³⁸ The RFOT-based analysis can be accommodated for the presence of steady state irradiation. Such irradiation shifts equilibrium toward a greater population of the antibonding orbitals, at the expense of the bonding orbitals, see the schematic in Fig. 26. This creates a driving force toward structural states whose LUMOs and adjacent orbitals are stabilized compared with the original state. As aging toward a LUMO-stabilized structure progresses, the optical gap will shrink. Thus, photodarkening can be argued to be of a very general nature, not limited just to chalcogenides. The only requirement is that the material be structurally degenerate. This is, of course, consistent with crystalline materials not displaying photodarkening.²⁴ Furthermore, according to the molecular orbital (MO) theory, LUMO is stabilized when the length of the pertinent subset of bonds is increased. Thus, one can expect quite generally that the material should expand, even if modestly, consistent with observation. Conversely, pressurization would destabilize such LUMO-stabilized states and partially suppress photodarkening, as is borne out in experiment. Golden et al.⁵⁴ have made a general, quantitative case that the band gap is positively correlated with the strength of the charge density wave, whether the latter is on-site or off-site or a combination of the two. Thus we conclude that the lowering of the gap corresponds to a decreased strength of the CDW and, hence, a suppression of the FSDP, again consistent with experiment. Likewise, one can use the RFOT-based notions of aging to explain on general grounds why samples aged using polarized light exhibit slightly less photodarkening for light polarized in the same direction than in the perpendicular direction.²⁴ In the spirit of Le Chatelier’s principle, there will be some photo-bleaching for transitions induced by the original polarized light but not in the perpendicular direction. (This partial bleaching is still a much weaker effect than the overall photodarkening.)

Deciding on detailed, system-specific structural changes accompanying the lowering of LUMO is however less straightforward; the reader is referred to the review of existing microscopic proposals in Ref. 24. Clearly, the overall expansion of the lattice cannot lead to uniform expansion of every bond since the resulting destabilization of the occupied orbitals in the valence band would more than offset the stabilization of the conduction-band orbitals occupied by optically promoted electrons. The only logical possibility, then, is that the overall width of the valence band is increased in the aged sample; this is reflected in Fig. 26. This way, the destabilization of the HOMO and adjacent orbitals, if any, due to the overall lengthening of the bonds, is compensated by stabilization and, hence, shortening, of some of the bonds. These conclusions are consistent with structural studies by Yang, Paesler, and Sayers,¹³² who find that irradiated samples exhibit a significantly increased Debye-Waller factor and some enhancement in the secondary bond length.

Finally, we address the apparent correlation between the Boson peak and the FSDP observed by Novikov and Sokolov.⁴² First we recall that Lubchenko and Wolynes^2,40,41 have used the RFOT theory to show that a small subset of the reconfigurations in a glassy liquid near T_g will have very small barrier. While fully dampened near T_g, such transitions become facile at sufficiently low temperatures where they occur by tunneling. The density of states for resulting structural resonances can be estimated and quantitatively account for the so-called two-level systems³⁹ (TLS) and two signature phenomena associated with the Boson peak (BP), viz., the bump in the heat conductivity and the plateau in thermal conductivity. Each TLS represents the two lowest energy levels of such a structural resonance, while the BP is caused by a sub-set of higher energy states that corresponds to the vibrations of the boundary of the resonance. We have called these modes “ripplons.”^40,41 Although exhibiting the density of states typical of a harmonic mode, the ripplons do not correspond to a static structural inhomogeneity but, instead, are of dynamical origin and can be thought of intrinsically stemming from the quantum-mechanical uncertainty in the position of the boundary of the structural resonance. According to Refs. 2, 40, and 41, the basic length scale behind the Boson peak is given by the size $ξ$ of the cooperatively reconfiguring region. The latter size, in terms of the volumetric size a of the rigid molecular unit (or “bead”), is largely system-independent and depends only logarithmically on the time scale of the glass transition so that reconfiguring regions encompass $(ξ/a)3≈200$ beads near the typical laboratory glass transition.^1,2 The resulting (angular) frequency spectrum of the ripplons is given by^40,41

where $ωD$ is the Debye frequency defined as $(6π2)1/3cs/a$⁠, c_s standing for the speed of sound. Mode l is (2l + 1)-degenerate. [To avoid confusion we note that the Debye frequency is not strictly defined in amorphous materials, see discussion in Ref. 139, and is understood here as the ultra-violet cut-off of the theory. The odd-looking factor $(6π2)1/3$ stems from our defining the bead size volumetrically.] The quantum number l is greater or equal than 2, by symmetry, and is not expected to exceed ten, because the wave length of the ripplon cannot be less than 2a. Sokolov and Novikov (SN)⁴² specifically referred to the BP frequency $ωBP$ as determined by the Raman shift in inelastic neutron scattering, for which the pertinent frequency scale would be set—in the present framework—by frequency differences $ωl1−ωl2$⁠. According to the SN study (see their Fig. 4), $qFSDP≈30ωBP/2πcs$⁠. This, combined with Eq. (10) and $ωD=(6π2)1/3cs/a$ yields

where the angular brackets denote averaging over the ripplon harmonics, properly weighted by their multiplicity and Boltzmann factors. In view of the high anharmonicity of the resonance and its being embedded in a solid matrix, there are likely no selection rules, consistent with the results for phonon scattering obtained in Ref. 40 The overall scaling of $qFSDP$ with the inverse bead size a is certainly consistent with our earlier identification of the FSDP as reflecting the length scale describing spatial modulation of bonding strength. Numerically, the expression in the angular brackets in Eq. (11) varies between ca. 0.6 and 4.2 (up to a sign), the two extremes corresponding to $(ω3−ω2)$ and $(ω10−ω2)$⁠, respectively. The average in Eq. (11) is thus not too different, numerically, from 2.5. (A more accurate estimate seems hardly warranted because of a lack of detailed knowledge of the anharmonicity and the effects of damping, some of which have been estimated in Ref. 40.) The resulting expression on the r.h.s. of Eq. (11), $≃7.5a−1$ is obviously consistent with the lengthscale corresponding to FSDP estimated earlier to be near $7.7/qFSDP$⁠. Thus, on one hand, the above notions help rationalize the empirical correlation between the position of the FSDP and the Boson peak frequency⁴² and, on the other hand, support our earlier conclusions on the intrinsic connection between the FSDP and the rigid structural unit in amorphous chalcogenides.

See supplementary material for information on (I) cavity statistics when mergeable Delaunay tetrahedra are sorted in a deterministic fashion, (II) illustration of an alternative criterion for cavity merger, (III) quantification of octahedrality of parent structures, (IV) details of evaluation of g(r) and S(q), (V) re-parametrization of MOPAC, (IV) performance of CP2K vs. VASP, (VII) g(r) and S(q) for individual realizations of optimized structures, (VIII) dependence of g(r) and S(q) for different values on A, (IX) comparison of C-type and P-type structures at the level of pairwise correlation functions, (X) details of pre-optimization with respect to local density fluctuations, (XI) particle mobility during geometric optimization, (XII) alternative way to visualize the CDW, (XIII) additional ring statistics, (XIV) effects of pressure on secondary bonds.

We thank Dr. Jon C. Golden for his expert advice. This work has been supported by the National Science Foundation Grant Nos. CHE-0956127 and CHE-1465125, and the Welch Foundation Grant No. E-1765. We gratefully acknowledge the use of the Maxwell/Opuntia Cluster and the untiring support from the Center for Advanced Computing and Data Systems at the University of Houston. Partial support for this work was provided by resources of the uHPC cluster managed by the University of Houston and acquired through NSF Award No. ACI-1531814.