Boson peak in covalent network glasses: Isostaticity and marginal stability

TL;DR

This study demonstrates that the boson peak in covalent network glasses like silica is governed by near-isostatic constraints and marginal stability, unifying understanding across different glassy systems through molecular dynamics simulations and experimental data.

Contribution

It provides the first direct evidence linking the boson peak in covalent glasses to isostaticity and marginal stability, extending concepts from packing-based glasses.

Findings

Boson peak in silica is governed by near-isostatic constraints.

The boson peak manifests as a wavenumber-independent band in the dynamical structure factor.

Results are consistent with inelastic X-ray scattering data.

Abstract

The boson peak (BP) stands as a key feature in understanding glasses and amorphous materials. It directly underlies their anomalous material properties, including thermal behaviors such as excess specific heat and low thermal conductivity, as well as mechanical characteristics such as nonaffine elasticity and fragile plasticity. Despite its importance, understanding of the BP remains limited in covalent network glasses. The most promising concepts are isostaticity and marginal stability, which have been established in theories of rigidity percolation and the jamming transition. While these concepts, supported by comprehensive data, account for the BP in packing-based glasses, comparable explanations have not yet been demonstrated for covalent network glasses. Here we study silica glass, a prototypical covalent network glass, using molecular dynamics simulations. We show that the BP in silica glass is governed by near-isostatic constraints and marginal stability, supporting their universality across diverse glassy systems. Furthermore, we reveal that these principles manifest as a wavenumber-independent band in the dynamical structure factor, and we demonstrate consistency with inelastic X-ray scattering data on silica glass. Our results provide a unified, experimentally testable framework for deciphering the BP and for refining the interpretation of scattering data in amorphous materials.

Tables

Table 1. Table 1: Parameters for the potential of silica glass. q Si = 1.7755 ​ e q_{\text{Si}}=1.7755\,e , q O = − q Si / 2 q_{\text{O}}=-q_{\text{Si}}/2 , r c ​ S = 8.0 r_{cS}=8.0  Å, r c ​ L = 10.0 r_{cL}=10.0  Å, and γ S = γ L = 0.2 \gamma_{S}=\gamma_{L}=0.2  Å.

$\alpha$–$\beta$	$A_{\alpha \beta }$ (eV)	$B_{\alpha \beta }$ (Å-1)	$C_{\alpha \beta }$ (eVÅ6)	$D_{\alpha \beta }$ (eVÅ24)
Si–O	$23107.8$	$5.098$	$139.7$	$66.0$
O–O	$1120.5$	$2.893$	$26.1$	$16800.0$
Si–Si	$2797.9$	$4.407$	$0.0$	$3423204.0$

Table 2. Table 2: Physical quantities including elastic moduli and Debye values. For silica glass, the quantities are measured as mass density ρ \rho  (g/cm 3 ), elastic moduli K , G K,G  (GPa), sound speeds c α c_{\alpha}  (m/s), wavenumber q q  (Å -1 ), Debye level A D A_{D}  (THz -3 ), and frequency ω \omega  (THz). See also Table S1 of SI for values of the overconstrained-network system and the isostatic-network system.

Glass	$\rho$	$K$	$K_A$	$K_N$	$\frac{K_N}{K_A}(\%)$	$G$	$G_A$	$G_N$	$\frac{G_N}{G_A}(\%)$	$\nu$	$c_L$	$c_T$	$\frac{c_L}{c_T}$	$q_D$	$A_D$	${\omega }_D$	${\omega }_{\text{BP}}$	${\omega }_0$
Silica	$2.20$	$40.9$	$172$	$131$	$76.3$	$30.5$	$104$	$73.6$	$73.7$	$0.202$	$6090$	$3720$	$1.64$	$1.58$	$0.00257$	$10.5$	$1.21$	$0.46$
HSH	$1.40$	$0.544$	$0.674$	$0.131$	$19.4$	$0.122$	$0.344$	$0.222$	$64.4$	$0.395$	$0.712$	$0.296$	$2.40$	$4.36$	$0.965$	$1.46$	$0.0970$	$0.033$
HSL	$1.25$	$0.332$	$0.478$	$0.145$	$30.4$	$0.0406$	$0.281$	$0.240$	$85.5$	$0.441$	$0.556$	$0.180$	$3.09$	$4.20$	$4.70$	$0.81$	$0.0379$	$-$
LJ	$1.015$	$61.2$	$61.7$	$0.530$	$0.859$	$13.6$	$35.8$	$22.2$	$61.9$	$0.396$	$8.84$	$3.67$	$2.41$	$3.92$	$0.000699$	$16.3$	$1.05$	$0.38$
SS	$1.015$	$40.8$	$40.8$	$0.00$	$0.00$	$6.21$	$14.7$	$8.50$	$57.8$	$0.428$	$6.96$	$2.47$	$2.81$	$3.92$	$0.00225$	$11.0$	$0.798$	$0.35$

Table 3. Table 3: Characteristics of Si–O bonds and O–O contacts in silica glass. “Separation” lists the typical nearest-neighbor separations. “Spring constant” lists ϕ S ′′ \phi^{\prime\prime}_{S} and ϕ L ′′ \phi^{\prime\prime}_{L} , the second derivatives of the short-range and long-range Coulomb pair potentials, respectively. Neighbor (pair) definitions use distance cutoffs taken at r ≤ 1.8 r\leq 1.8  Å for Si–O bonds and r ≤ 2.9 r\leq 2.9  Å for O–O contacts. See also Fig.  1 .

Pair	Separation (Å)	Spring constant (eVÅ-2)	Number of pairs
Si–O	$1.6$	${\varphi }_S^{\prime \prime }=35.8$	$4.00N_{\text{Si}}$
		${\varphi }_L^{\prime \prime }=-11.1$
O–O	$2.6$	${\varphi }_S^{\prime \prime }=4.55$	$3.00N_{\text{O}}$
		${\varphi }_L^{\prime \prime }=1.29$

Table 4. Table S1: Physical quantities including elastic moduli and Debye values. For silica glass, the quantities are measured as mass density ρ \rho (g/cm 3 ), elastic moduli K , G K,G (GPa), sound speeds c α c_{\alpha} (m/s), wavenumber q q (Å -1 ), Debye level A D A_{D} (THz -3 ), and frequency ω \omega (THz). In addition to values of the original system (reported in Table 3 of the main text), we list results for the overconstrained-network system and the isostatic-network system.

		$\rho$	$K$	$K_A$	$K_N$	$\frac{K_N}{K_A}(\%)$	$G$	$G_A$	$G_N$	$\frac{G_N}{G_A}(\%)$	$\nu$	$c_L$	$c_T$	$\frac{c_L}{c_T}$	$q_D$	$A_D$	${\omega }_D$	${\omega }_{\text{BP}}$
Silica	Original	$2.20$	$40.9$	$172$	$131$	$76.3$	$30.5$	$104$	$73.6$	$73.7$	$0.202$	$6090$	$3720$	$1.64$	$1.58$	$0.00257$	$10.5$	$1.21$
	Overcon	$2.20$	$53.6$	$203$	$149$	$73.6$	$29.4$	$121$	$92.4$	$75.8$	$0.268$	$6500$	$3660$	$1.78$	$1.58$	$0.00265$	$10.4$	$1.93$
	Isostatic	$2.20$	$\approx 0$	$196$	$196$	$100$	$\approx 0$	$118$	$118$	$100$	$-$	$\approx 0$	$\approx 0$	$-$	$1.58$	$+\infty$	$\approx 0$	$\approx 0$
HSH	Original	$1.40$	$0.544$	$0.674$	$0.131$	$19.4$	$0.122$	$0.344$	$0.222$	$64.4$	$0.395$	$0.712$	$0.296$	$2.40$	$4.36$	$0.965$	$1.46$	$0.0970$
	Overcon	$1.40$	$0.548$	$0.641$	$0.0932$	$14.4$	$0.246$	$0.385$	$0.138$	$36.0$	$0.305$	$0.793$	$0.420$	$1.89$	$4.36$	$0.351$	$2.04$	$0.469$
HSL	Original	$1.25$	$0.332$	$0.478$	$0.145$	$30.4$	$0.0406$	$0.281$	$0.240$	$85.5$	$0.441$	$0.556$	$0.180$	$3.09$	$4.20$	$4.70$	$0.81$	$0.0379$
	Overcon	$1.25$	$0.347$	$0.475$	$0.128$	$26.9$	$0.0757$	$0.285$	$0.209$	$73.4$	$0.398$	$0.598$	$0.246$	$2.43$	$4.20$	$1.88$	$1.17$	$0.146$
LJ	Original	$1.015$	$61.2$	$61.7$	$0.530$	$0.859$	$13.6$	$35.8$	$22.2$	$61.9$	$0.396$	$8.84$	$3.67$	$2.41$	$3.92$	$0.000699$	$16.3$	$1.05$
	Overcon	$1.015$	$60.6$	$61.0$	$0.415$	$0.680$	$21.6$	$36.6$	$15.0$	$41.1$	$0.341$	$9.39$	$4.61$	$2.04$	$3.92$	$0.000360$	$20.3$	$4.54$
SS	Original	$1.015$	$40.8$	$40.8$	$0.00$	$0.00$	$6.21$	$14.7$	$8.50$	$57.8$	$0.428$	$6.96$	$2.47$	$2.81$	$3.92$	$0.00225$	$11.0$	$0.798$
	Overcon	$1.015$	$35.4$	$35.4$	$0.00$	$0.00$	$17.6$	$21.2$	$3.66$	$17.2$	$0.287$	$7.61$	$4.16$	$1.83$	$3.92$	$0.000499$	$18.2$	$6.51$

Equations

$$v_S(r)=A_{\alpha \beta }\exp \left(-B_{\alpha \beta }r\right)-\frac{C_{\alpha \beta }}{r^6}+\frac{D_{\alpha \beta }}{r^{24}},$$

$$v_L(r)=\frac{q_{\alpha }{q}_{\beta }}{r},$$

$${\varphi }_S(r)$$

$${\varphi }_L(r)$$

$$G_{cS}(r)$$

$$G_{cL}(r)$$

$$\varphi (r)=\frac{ϵ}{2}{\left(1-\frac{r}{\sigma }\right)}^{2}H(\sigma -r),$$

$$v(r)=4ϵ\left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right],$$

$$\varphi (r)=v(r)-v(r_c)-(r-r_c)v^{\prime }(r_c),$$

$$v(r)=ϵ{\left(\frac{{\sigma }_{\alpha \beta }}{r}\right)}^{12},$$

$$g(\omega )=\frac{1}{3N}\sum _{k=1}^{3N}\delta (\omega -{\omega }_k),$$

$$g_D(\omega )=A_D{\omega }^2=\frac{3}{{\omega }_D^3}{\omega }^2\propto {\omega }^2,$$

$${\omega }_D={\left(\frac{c_L^{-3}+2c_T^{-3}}{3}\right)}^{-1/3}{q}_D={\left(\frac{18{\pi }^2\hat{\rho }}{c_L^{-3}+2c_T^{-3}}\right)}^{1/3},$$

$$\delta e_{k,\text{Si–O}}^2=\frac{m_{\text{O}}}{N_{\text{Si–O}}}\sum _{⟨ij⟩\in {\mathcal{B}}_{\text{Si–O}}}{\left|\left(\frac{{\vec{e}}_{k,i}}{\sqrt{m_{\text{Si}}}}-\frac{{\vec{e}}_{k,j}}{\sqrt{m_{\text{O}}}}\right)\cdot {\vec{n}}_{ij}\right|}^2,$$

$${\mathcal{P}}_k=\frac{1}{N}{\left[\sum _{i=1}^N{\left({\vec{e}}_{k,i}\cdot {\vec{e}}_{k,i}\right)}^2\right]}^{-1}.$$

$$S_{\alpha }(q,\omega )=\frac{k_{B}T}{2N}\frac{q^2}{{\omega }^2}\sum _{k=1}^{3N}{F}_{k,\alpha }(q)\delta \left(\omega -{\omega }_k\right),$$

$$F_{k,T}(q)$$

$$F_{k,L}(q)$$

$${\int }_0^{q_D}\left\{\frac{S_T(q,\omega )}{k_{B}T}+\frac{S_L(q,\omega )}{k_{B}T}\right\}dq$$

$$=\frac{1}{2N{\omega }^2}\sum _{k=1}^{3N}\delta \left(\omega -{\omega }_k\right)\left(\sum _{i,j=1}^N\frac{{\vec{e}}_{k,i}\cdot {\vec{e}}_{k,j}}{\sqrt{m_i{m}_j}}{\int }_0^{q_D}{e}^{\mathrm{i}\vec{q}\cdot ({\vec{r}}_i-{\vec{r}}_j)}{q}^{2}dq\right).$$

$${\int }_0^{q_D}{e}^{\mathrm{i}\vec{q}\cdot ({\vec{r}}_i-{\vec{r}}_j)}{q}^{2}dq=\frac{1}{4\pi }{\int }_{0\le |\vec{q}|\le q_D}{e}^{\mathrm{i}\vec{q}\cdot ({\vec{r}}_i-{\vec{r}}_j)}{d}^3\vec{q}$$

$$=\frac{1}{4\pi }{\delta }_{i,j}\left({\int }_{0\le |\vec{q}|\le q_D}{d}^3\vec{q}\right)=\frac{q_D^3}{3}{\delta }_{i,j}.$$

$${\int }_0^{q_D}\left\{\frac{S_T(q,\omega )}{k_{B}T}+\frac{S_L(q,\omega )}{k_{B}T}\right\}dq$$

$$=\frac{q_D^3}{2{\omega }^2}\frac{1}{3N}\sum _{k=1}^{3N}\delta \left(\omega -{\omega }_k\right)\left(\sum _{i=1}^N\frac{|{\vec{e}}_{k,i}{|}^2}{m_i}\right)=\frac{q_D^3}{2M(\omega )}\frac{g(\omega )}{{\omega }^2},$$

$$N_{\text{const}}=\sum _b{N}_b\left[\frac{b}{2}+(2b-3)\right],$$

$$N_{\text{const}}=N_{\text{O}}\cdot 1+N_{\text{Si}}\cdot (2+5)=9N_{\text{Si}}=N_{\text{dof}},$$

$$\frac{g(\omega )}{{\omega }^2}=\frac{2M(\omega )}{q_D^3}{\int }_0^{q_D}\left\{\frac{S_T(q,\omega )}{k_{B}T}+\frac{S_L(q,\omega )}{k_{B}T}\right\}dq,$$

$$g(\omega )$$

$$S(q,\omega )$$

$$\frac{g(\omega )}{{\omega }^2}=\frac{2m}{q_D^3}{\int }_0^{q_D}3\frac{S(q,\omega )}{k_{B}T}dq,$$

Full text

I Introduction

The vibrational states of atoms and molecules in crystals are well characterized as phonons, providing a clear understanding of their material properties [1]. In contrast, the vibrational states in glasses deviate markedly from phonons, producing anomalous properties relative to crystals [2]. Whereas the vibrational density of states (vDOS) in crystals is accurately captured by Debye theory based on phonons, glasses exhibit an excess over the Debye prediction, an anomaly known as the boson peak (BP) [2, 3]. The BP has been observed in the dynamical structure factor by light, inelastic X-ray, and neutron scattering, establishing it as an experimental hallmark of glasses [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. The BP is directly linked to the excess specific heat of glasses and strongly influences transport properties including sound propagation and heat conduction [17, 18, 19, 20, 21, 22]. It also affects nonaffine elastic responses and plastic deformation in glasses [23, 24, 25, 26]. Thus, the BP is central to understanding glasses and amorphous materials.

Understanding of the BP has advanced over decades. Several theoretical frameworks have been proposed, including elastic heterogeneities [27, 28, 29, 30, 31], soft-potential models [32, 33, 34, 35], and random-matrix approaches [36, 37]. Among these, early works [38, 39, 40, 41] introduced the concept of isostaticity based on the topology of random networks. They accounted for the rigidity of covalent network glasses by counting constraints (bonds) against degrees of freedom. When the number of constraints falls below the number of degrees of freedom, the system can deform continuously at no energy cost, and zero-energy vibrational modes, known as floppy modes, appear. Conversely, when constraints exceed degrees of freedom, excess constraints render energy-free deformations costly, shifting floppy modes to finite frequencies. The boundary where constraints equal degrees of freedom defines the isostatic state, the rigidity threshold known as Maxwell’s criterion [42]. Thus, by tracking how constraints surpass degrees of freedom, one can assess rigidity and the character of low-frequency vibrational states. This concept was developed within the framework of rigidity percolation [43, 44, 45, 46] and has also been employed to explain low-energy dynamics and two-level tunneling states in glasses [47, 48].

Subsequently, isostaticity was incorporated into the physics of jammed materials, providing insights into the jamming and glass transitions [49, 50, 51]. In this context, the topological picture above was unified with energetic considerations through the notion of frustration [52, 53, 54]. Dense and disordered packings generate strong repulsive forces, giving rise to frustration within the system. It has been proposed that glasses are solids that marginally counteract this frustration. From the vibrational perspective, frustration shifts modes that would be floppy at isostaticity but are pushed to finite frequencies by excess constraints back toward low frequencies. This produces a gapless vDOS that extends down to zero frequency and signals marginal stability. These modes contribute to the excess low-frequency vDOS and form the BP. Consequently, the BP can be viewed as a manifestation of isostaticity and marginal stability. Building upon this scenario, a mean-field theory was constructed using random spring networks, which successfully explains the BP [55, 56].

Numerical simulations have tested this scenario and its mean-field predictions. In particular, the packing model of harmonic spheres (HS) has been extensively studied in the context of the jamming transition [49, 50, 51]. Jamming scaling laws for the vDOS and the BP have been established [57, 58, 59, 60, 61], consistent with isostaticity and marginal stability [52, 53, 54] and with mean-field theory [55, 56]. Isostaticity and marginal stability have also been confirmed for the BP of Lennard–Jones (LJ) glass, a representative model of van der Waals glasses [62, 63, 64], and hard-sphere glass, a representative model of colloidal glasses [65, 66]. Across these simulations, quasi-localized vibrations (QLVs) emerge alongside the BP as a consequence of marginal stability driven by frustration [58, 60, 67, 68]. Thus, the BP is now well rationalized in packing-based glasses, including particle packings, van der Waals glasses, and colloidal glasses, in terms of isostaticity and marginal stability.

In contrast, the BP remains less well understood in network-forming glasses, that is, covalent network glasses such as oxide glasses, especially when compared with packing-based systems including van der Waals and colloidal glasses. Extensive light, inelastic X-ray, and neutron scattering experiments have measured the dynamical structure factor of covalent network glasses and consistently observed the BP [4, 6, 7, 8, 10, 12, 13, 14, 16]. The concept of isostaticity was originally proposed to interpret these experimental observations in covalent network glasses [38, 39, 40, 41]. However, a detailed validation of isostaticity, analogous to that achieved for packing, van der Waals, and colloidal glasses, has not yet been performed. Moreover, marginal stability driven by frustration remains largely unexplored in this class of materials. In addition, how isostaticity and marginal stability manifest in the experimentally measured dynamical structure factor has not been clarified. Given the breadth of scattering data on covalent network glasses, establishing the BP mechanism in these materials and its expression in the dynamical structure factor is a crucial step toward bridging theoretical frameworks and experiments.

In this work, we focus on silica ($\text{SiO}_{2}$) glass, a representative oxide glass, to investigate its BP. Previous numerical studies have examined the vDOS and QLVs in silica glass [69, 70, 71, 72]. In contrast to these studies, our primary goal is to understand the BP in terms of isostaticity and marginal stability. We further analyze the dynamical structure factor to clarify how isostaticity and marginal stability are expressed in this experimentally accessible observable. We find that the BP in silica glass follows from near-isostatic constraints and marginal stability, as in HS and LJ glasses, providing strong evidence for the universality of these principles across diverse glassy materials. We further show that these principles manifest as a wavenumber-independent broad band in the dynamical structure factor, and we demonstrate that this band is consistent with inelastic X-ray scattering data for silica glass. These results advance understanding of the BP in network-forming glasses and offer practical criteria for interpreting scattering data.

II Methods

II.1 Silica glass

We perform molecular dynamics (MD) simulations of silica ($\text{SiO}_{2}$) glass. Silica glass consists of $N_{\text{Si}}$ silicon (Si) atoms and $N_{\text{O}}$ oxygen (O) atoms in three-dimensional space, where $N_{\text{O}}=2N_{\text{Si}}$ and the total number of atoms is $N=N_{\text{O}}+N_{\text{Si}}=3N_{\text{Si}}$. The masses of the Si atom and the O atom are denoted as $m_{\text{Si}}$ and $m_{\text{O}}$, respectively, with the ratio $m_{\text{Si}}/m_{\text{O}}=1.755$. The interatomic forces are modeled using the BKS family of potentials, originally proposed in Ref. [73] and subsequently modified and extended in several aspects [74, 75, 76, 77, 78]. BKS-type models are pairwise-additive and comprise a short-range two-body term and a long-range Coulomb term, and therefore they do not include explicit many-body angular contributions that encode covalent directionality. For a more explicit treatment of covalent bonding, many-body potentials with angular terms have been developed, notably the Vashishta potential [79], which adds three-body interactions associated with O–Si–O and Si–O–Si angles. In this study we adopt the SHIK parameterization of the BKS family [78], which has been rigorously tested to reproduce experimental data and first-principles simulations, including thermodynamic quantities, structural properties (radial distribution function, static structure factor, bond-angle distribution), and elastic moduli [74, 75, 76, 77, 78]. Previous numerical studies have investigated mechanical responses and vibrational states of silica glasses using BKS-type models [80, 81, 82], and a recent simulation explored the effects of densification on the BP [83].

The potential comprises two contributions: the short-range term and the long-range Coulomb term. The short-range term is defined as

[TABLE]

where $r$ is the interparticle distance, and $\alpha$ and $\beta$ denote either O or Si. The parameters $A_{\alpha\beta}$, $B_{\alpha\beta}$, $C_{\alpha\beta}$, and $D_{\alpha\beta}$ are specified in Table 1. The long-range Coulomb term is given by

[TABLE]

where the charge of silicon is set as $q_{\text{Si}}=1.7755\,e$ (with $e$ the elementary charge) and, to ensure charge neutrality, the charge of oxygen is set as $q_{\text{O}}=-q_{\text{Si}}/2$. The potentials $v_{S}(r)$ and $v_{L}(r)$ are truncated at $r=r_{cS}=8$ Å and $r=r_{cL}=10$ Å, respectively. For $v_{L}(r)$, this truncation is referred to as Wolf truncation [74, 75]. To prevent discontinuities in the potential and its force (the first derivative of the potential) at the cutoff distances, $v_{S}(r)$ and $v_{L}(r)$ are smoothed as follows:

[TABLE]

with the smoothing functions

[TABLE]

where $\gamma_{S}=\gamma_{L}=0.2$ Å.

We performed three independent MD simulations using LAMMPS [84]. We fixed the mass density of silica glass at $\rho=2.20$ g/cm3 and controlled the temperature with a Nose-Hoover thermostat [85, 86]. Initial Si and O positions were randomized and the system was equilibrated at $T=3500$ K for $100$ ps to obtain a homogeneous liquid. The system was then quenched to $T=300$ K at a rate of $1$ K/ps, followed by equilibration at $T=300$ K for $100$ ps. Subsequently, all atomic velocities were set to zero and the configuration was relaxed (energy minimization) to obtain the inherent structure $\vec{r}=[\vec{r}_{1},\vec{r}_{2},\ldots,\vec{r}_{N}]$.

To access a broad range of wavenumbers $q$ and frequencies $\omega$, we used several system sizes ranging from $N=1.5\times 10^{4}$ to $1.2\times 10^{5}$ atoms. All reported quantities are averages over the three independent configurations.

In addition to silica glass, we analyze HS, LJ, and soft-sphere (SS) glasses in three dimensions. The HS model, extensively studied in the context of the jamming transition [49, 50, 51], is a packing glass known to exhibit isostaticity and marginal stability together with the BP and QLVs [52, 53, 54, 57, 58, 59, 60, 61]. We consider two HS preparations, a high-connectivity sample (HSH) generated at high pressure and a low-connectivity sample (HSL) generated at low pressure. Furthermore, LJ and SS glasses, representative van der Waals systems, have been investigated for their BP and QLVs in prior work [62, 63, 64, 20, 28, 21, 29].

II.2 HS glass

HS glass is a one-component system in three dimensions that has been extensively studied in the context of the jamming transition [49, 50, 51, 57, 58, 59, 60, 61]. Two particles $i$ and $j$ interact via the harmonic potential,

[TABLE]

where $r=r_{ij}$ is the interparticle distance, $\sigma$ is the particle diameter, and $H$ is the Heaviside step function. All particles have the same mass $m$. Physical quantities are measured in units of length $\sigma$, mass $m$, and energy $\epsilon$; accordingly, $q$ and $\omega$ are reported in units of $\sigma^{-1}$ and $\sqrt{\epsilon/(m\sigma^{2})}$, respectively. We prepared HS glasses at pressures $P=0.05$ and $0.005$, referred to as HSH and HSL in this paper. At the higher pressure $P=0.05$ (HSH), particles are more densely packed and exhibit higher connectivity, whereas at the lower pressure $P=0.005$ (HSL), particles are less dense and have lower connectivity. As with silica glass, we simulated multiple system sizes, ranging from $N=1.6\times 10^{4}$ to $1.024\times 10^{6}$ particles, to access broad ranges of $q$ and $\omega$.

II.3 LJ glass

LJ glass is a one-component system in three dimensions, studied in previous works [20, 63, 87]. Two particles $i$ and $j$ interact via the Lennard–Jones (LJ) potential,

[TABLE]

where $r=r_{ij}$ is the interparticle distance and $\sigma$ is the particle diameter. The potential is truncated at $r=r_{c}=2.5\sigma$. To avoid artificial anharmonicities caused by the discontinuity at $r=r_{c}$ [88], we employ the smoothed form

[TABLE]

which ensures that both the potential and the force (the derivative of the potential) vanish smoothly at $r=r_{c}$. All particles have the same mass $m$. Physical quantities are measured in units of length $\sigma$, mass $m$, and energy $\epsilon$. The number density is set to $\hat{\rho}=N/V=1.015$, where $V$ is the system volume. We simulated several system sizes, ranging from $N=1.6\times 10^{4}$ to $1.024\times 10^{6}$.

II.4 SS glass

SS glass is a binary mixture of large ($L$) and small ($S$) particles in three dimensions, which we have previously studied in Refs. [29, 21]. Particles $i$ and $j$ of types $\alpha$ and $\beta$ ($\alpha,\beta\in\{L,S\}$) interact via a 12-inverse-power-law potential,

[TABLE]

where $r=r_{ij}$, $\sigma_{\alpha\beta}=(\sigma_{\alpha}+\sigma_{\beta})/2$, and $\sigma_{\alpha}$ is the diameter of the large or small species ($\sigma_{L}$ or $\sigma_{S}$). The size ratio is $\sigma_{S}/\sigma_{L}=0.7$, and the composition is equimolar with $x_{L,S}=N_{L,S}/N=1/2$ and $N=N_{L}+N_{S}$. As in the LJ case, we employ the smoothed form $\phi(r)$ defined in Eq. (7). All particles have the same mass $m$. Physical quantities are measured in units of length $\sigma=\left(\sum_{\alpha,\beta=L,S}x_{\alpha}x_{\beta}\sigma_{\alpha\beta}^{3}\right)^{1/3}$, mass $m$, and energy $\epsilon$. The number density is set to $\hat{\rho}=N/V=1.015$. We simulated several system sizes, ranging from $N=1.6\times 10^{4}$ to $1.024\times 10^{6}$.

II.5 Vibrational modes

Using the inherent structure $\vec{r}=[\vec{r}_{1},\vec{r}_{2},\cdots,\vec{r}_{N}]$, we perform a standard normal-mode analysis by solving the eigenvalue problem of the dynamical matrix (a $3N\times 3N$ matrix) to obtain the eigenvalues $\lambda_{k}$ and eigenvectors $\vec{e}_{k}=[\vec{e}_{k,1},\vec{e}_{k,2},\cdots,\vec{e}_{k,N}]$ for modes $k=1,2,\cdots,3N$ [89]. We remove the three zero-frequency translational modes from the analysis. The eigenfrequencies are $\omega_{k}=\sqrt{\lambda_{k}}$, and the eigenvectors are orthonormalized such that $\vec{e}_{k}\cdot\vec{e}_{\ell}=\sum_{i=1}^{N}\vec{e}_{k,i}\cdot\vec{e}_{\ell,i}=\delta_{k\ell}$, where $\delta_{k,l}$ is the Kronecker delta.

In the present work, we analyze several system sizes, ranging from $N\sim 10^{4}$ to $10^{6}$, as in Refs. [60, 63]. We compute all vibrational modes for the smallest system and only the low-frequency modes for larger systems. We then merge the mode datasets from different sizes as a function of $\omega_{k}$; results from different sizes connect smoothly, thereby extending access to the lower-frequency regime. $g(\omega)$ in Eq. (9) and $S_{\alpha}(q,\omega)$ in Eq. (14) are computed from these combined datasets, and results from different sizes are presented together in the figures.

II.5.1 vDOS

From the set of eigenfrequencies $\{\omega_{k}\}$ for vibrational modes $k=1,2,\ldots,3N$, we compute the vDOS $g(\omega)$ as

[TABLE]

where $\delta$ denotes the Dirac delta function.

II.5.2 Debye vDOS

In an isotropic elastic medium, vibrational modes are phonons [1]. Continuum elasticity gives the linear dispersions $\omega=c_{T}q$ and $\omega=c_{L}q$, where $c_{T}$ and $c_{L}$ are the transverse and longitudinal sound speeds, respectively. Counting phonon states yields the Debye vDOS,

[TABLE]

where $A_{D}=3/\omega_{D}^{3}$ is the Debye level, and the Debye frequency $\omega_{D}$ is

[TABLE]

with the Debye wavenumber $q_{D}=(6\pi^{2}\hat{\rho})^{1/3}$ and the number density $\hat{\rho}=N/V$.

II.5.3 Stretching and compression of Si–O bonds in silica glass

For silica glass, we quantify the stretching and compression of Si–O bonds in vibrational mode $k$ via the dimensionless measure $\delta e_{k,\text{Si--O}}^{2}$, defined as

[TABLE]

where $\mathcal{B}_{\text{Si--O}}$ is the set of nearest-neighbor Si–O pairs, $N_{\text{Si--O}}=|\mathcal{B}_{\text{Si--O}}|=4N_{\text{Si}}$ is the total number of Si–O bonds, and ${\vec{n}}_{ij}=(\vec{r}_{i}-\vec{r}_{j})/|\vec{r}_{i}-\vec{r}_{j}|$ is the unit vector along the Si–O bond, with $i$ the Si atom and $j$ the O atom. This quantity captures the mean-squared relative displacement of a Si atom and its O neighbor projected along the bond direction; the mass factors render it dimensionless.

II.5.4 Participation ratio

We measure the fraction of particles that participate in vibrational mode $k$ using the participation ratio $\mathcal{P}_{k}$ [89]:

[TABLE]

For reference, $\mathcal{P}_{k}=1$ corresponds to a perfectly extended mode with equal amplitude on all particles, whereas $\mathcal{P}_{k}=1/N$ corresponds to a mode localized on a single particle.

II.6 Dynamical structure factors

From the data of $\omega_{k}$ and $\vec{e}_{k}=[\vec{e}_{k,1},\vec{e}_{k,2},\cdots,\vec{e}_{k,N}]$, we compute the dynamical structure factors $S_{\alpha}(q,\omega)$, where $\alpha\in\{T,L\}$ denotes the transverse and longitudinal polarizations, respectively, following [90, 89]:

[TABLE]

with

[TABLE]

where $k_{B}$ is the Boltzmann constant; $q=|\vec{q}|$ and $\omega$ are the wavenumber and frequency; and $\hat{\vec{q}}=\vec{q}/q$. It is important to note that scattering experiments typically access the longitudinal component $S_{L}(q,\omega)$ [91].

Using $S_{T}(q,\omega)$ and $S_{L}(q,\omega)$ in Eqs. (14) and (15), we obtain

[TABLE]

Assuming an isotropic elastic medium, we have

[TABLE]

Substituting Eq. (17) into Eq. (16) and using $g(\omega)$ from Eq. (9), we obtain

[TABLE]

where we define the effective mass $M(\omega)$ via $M^{-1}=\sum_{i=1}^{N}|\vec{e}_{k,i}|^{2}/m_{i}$, evaluated at $\omega=\omega_{k}$. Thus $M(\omega)$ depends on the mode $k$ and hence on the frequency. Equation (18) corresponds to Eq. (21).

In theoretical analyses [27, 28, 30, 61, 55, 56, 14], Eq. (21) is applied to compute the vDOS from the dynamical structure factor (or from Green’s functions). As shown above, Eq. (21) rests on the isotropic-medium assumption in Eq. (17), which is generally a good approximation in the low-frequency regime.

II.7 Elastic moduli and Debye values

We calculate the elastic moduli, the bulk modulus $K$ and the shear modulus $G$, using the harmonic formulation based on the linear response theory [92, 93, 89]. From these moduli, we compute Debye quantities such as the Debye frequency $\omega_{D}$ and the Debye level $A_{D}$.

Table 2 summarizes the elastic moduli and related quantities. The elasticity of glasses involves not only affine but also nonaffine deformation [92, 94, 87, 93]. Consequently, an elastic modulus $M$ decomposes as $M=M_{A}-M_{N}$, where $M_{A}$ and $M_{N}$ are the affine and nonaffine contributions. HS (HSH and HSL), LJ, and SS glasses exhibit strong nonaffine contributions to their shear moduli, whereas the nonaffine contribution to the bulk modulus is comparatively small; notably, in LJ and SS glasses the bulk modulus is almost entirely determined by the affine component.

By contrast, silica glass exhibits significant nonaffine contributions to both bulk and shear moduli: the nonaffine components exceed $70\%$ of the affine components. Such strong nonaffinity in both moduli arises from the lack of centrosymmetry in the tetrahedral network structure [95, 96], which has also been observed in amorphous silicon [97] and in noncentrosymmetric crystals such as $\alpha$-quartz [98]. As a result, the bulk and shear moduli become closer in magnitude. The Poisson’s ratio $\nu=(3K-2G)/(6K+2G)$ is approximately $0.2$, about half the value $\nu\approx 0.4$ observed in HS, LJ, and SS glasses. Values near $\nu=0.2$ categorize silica glass as a strong glass, whereas $\nu\approx 0.4$ indicates that HS, LJ, and SS glasses are fragile glasses, as discussed in Refs. [99, 100].

III Results

Figure 1(a) shows the radial distribution function $g(r)$ of the inherent structure. The first peaks mark the typical nearest-neighbor separations for Si–O and O–O, at $r\approx 1.6$ Å and $r\approx 2.6$ Å, respectively (purple and green arrows). Throughout, we refer to Si–O nearest-neighbor pairs as Si–O bonds because they represent covalent Si–O linkages, whereas we refer to O–O nearest-neighbor pairs as O–O contacts because there is no covalent O–O bond; these contacts arise from the geometry of O–Si–O linkages, with the directionality of Si–O covalency represented implicitly by the pairwise potential [78, 73]. Specifically, we define Si–O bonds as pairs with $r\leq 1.8$ Å and O–O contacts as pairs with $r\leq 2.9$ Å (black arrows). Panels (b) and (c) of Fig. 1 visualize the atomic network using these definitions. These visualizations confirm a tetrahedral network in which each Si atom sits at the center and four O atoms occupy the vertices, yielding four Si–O bonds and six O–O contacts per tetrahedral unit.

In what follows, we present results for silica glass, HSH glass, and LJ glass in the main text. Complementary data for HSL and SS glasses are provided in the Supporting Information (SI) to corroborate and strengthen our conclusions.

III.1 Isostaticity

Figure 2(a) displays the vDOS $g(\omega)$ for silica glass (purple curve), which agrees with previously reported results [76, 77, 78, 80, 81, 82, 83]. For reference, the Debye vDOS $A_{D}\,\omega^{2}$ is also shown; in the low-frequency regime $g(\omega)$ clearly exceeds this Debye vDOS, revealing excess modes. We further resolve the low-frequency behavior in Fig. 3(a,b). In Fig. 3(b), the reduced vDOS $g(\omega)/\omega^{2}$ lies above the Debye level $A_{D}$, thereby confirming the presence of the BP. We define the BP frequency $\omega_{\text{BP}}$ as the frequency at which the reduced vDOS $g(\omega)/\omega^{2}$ attains its maximum. For silica glass, this yields $\omega_{\text{BP}}\approx 1.21$ THz, consistent with scattering experiments [4, 6, 7, 8, 10, 12, 13, 14, 16].

To understand the vDOS and the BP, we begin by examining isostaticity in silica glass. As shown in Fig. 1, the Si and O atoms form a tetrahedral network with typical nearest-neighbor separations of $1.6$ Å for Si–O bonds and $2.6$ Å for O–O contacts. Accordingly, we construct a network system by connecting pairs with unstressed (neither pre-stretched nor pre-compressed) springs: Si–O pairs within the cutoff $1.8$ Å are linked, and O–O pairs within $2.9$ Å are likewise connected. Table 3 lists typical values of the spring constants $\phi^{\prime\prime}_{S}$ and $\phi^{\prime\prime}_{L}$, i.e., the second derivatives of the short-range and Coulomb potentials, respectively, evaluated at $r=1.6$ Å for Si–O and $r=2.6$ Å for O–O, together with the counts of springs corresponding to Si–O bonds and O–O contacts. We first focus on the short-range interaction $\phi^{\prime\prime}_{S}$, neglecting the Coulomb contribution $\phi^{\prime\prime}_{L}$; that is, we analyze a network system in which these pairs are connected by springs with stiffness $\phi^{\prime\prime}_{S}$.

We can demonstrate that the network system constructed above is isostatic, meaning that the number of constraints $N_{\text{const}}$ equals the number of degrees of freedom $N_{\text{dof}}$. First, the total number of degrees of freedom is $N_{\text{dof}}=3(N_{\text{Si}}+N_{\text{O}})=9N_{\text{Si}}$. Next, following Ref. [41], we count constraints as follows. To match the situation in Ref. [41] within our spring representation, we regard the O–O springs along tetrahedral edges as effective angular constraints associated with the $\angle\text{O--Si--O}$ bending [78, 73], and we treat only the Si–O springs as central-force bonds (this situation is visualized in Fig. 1(b)). The total number of constraints is then

[TABLE]

where $b$ is the coordination number and $N_{b}$ is the number of atoms with coordination $b$; the first term counts bond-stretching (two-body) constraints, and the second term counts bond-bending (three-body) constraints centered on the atom. Oxygen has $b=2$ and contributes only the bond-stretching part $b/2=1$ per O (we do not assign O-centered angular constraints), whereas silicon has $b=4$ and contributes $b/2=2$ bond-stretching constraints and $(2b-3)=5$ angular constraints associated with $\angle\text{O--Si--O}$. Substituting into Eq. (19) gives

[TABLE]

which confirms that the network representation is isostatic.

We present the vDOS for this isostatic-network system in Fig. 2(a) and Fig. 3(a) (cyan curves). The vDOS remains finite as $\omega\to 0$, indicating a gapless spectrum at isostaticity that comprises both strictly zero-frequency floppy modes and many additional soft, low-frequency modes of isostatic origin [52, 53, 54].

A notable point is that Si–O bonds are significantly stiffer than O–O contacts: $\phi^{\prime\prime}_{S}=35.8$ eV Å*-2* for Si–O versus $\phi^{\prime\prime}_{S}=4.55$ eV Å*-2* for O–O, as listed in Table 3. Consequently, vibrations that involve stretching and compression of Si–O bonds appear at high frequencies [82]. Since each tetrahedral unit contains four Si–O bonds, the total number of such modes can reach $4N_{\text{Si}}$. To probe these modes, we quantify Si–O bond stretching and compression using the dimensionless measure $\delta e_{k,\text{Si--O}}^{2}$ (see Materials and Methods). Figure 2(b) plots $N\,\delta e_{k,\text{Si--O}}^{2}$ versus $\omega$. $N\,\delta e_{k,\text{Si--O}}^{2}$ is nearly zero in the low-$\omega$ regime, including the BP region, indicating minimal Si–O bond deformation there. Above $\omega\approx 12.5$ THz, however, $N\,\delta e_{k,\text{Si--O}}^{2}$ increases, signaling the onset of bond-stretching and bond-compression character. Thus $\omega\approx 12.5$ THz marks a boundary: modes above this frequency predominantly occupy the $4N_{\text{Si}}$ sector associated with Si–O bond stretching and compression, whereas modes below it comprise the remaining $5N_{\text{Si}}$ modes.

Because the short-range interaction $\phi^{\prime\prime}_{S}$ couples not only atoms within a tetrahedral unit but also atoms belonging to different units, it is natural to consider spring-network variants that deviate from isostaticity. These interunit couplings correspond to van der Waals interactions and are often referred to as weak interactions [40, 41, 45]. Specifically, we examine a network in which all O–O contacts are connected by unstressed springs. In this case the system is overconstrained, $N_{\text{const}}>N_{\text{dof}}$, and the vDOS tends to zero as $\omega\to 0$. We refer to this overconstrained variant as the overconstrained-network system, in contrast to the isostatic-network system defined above. As shown in Fig. 2, the excess constraints shift the low-frequency floppy modes of the isostatic-network system to higher frequencies in the overconstrained-network system (compare the cyan and green curves). By contrast, the high-frequency modes associated with stretching and compression of Si–O bonds above $\omega\approx 12.5$ THz are minimally affected, and the $N\,\delta e_{k,\text{Si--O}}^{2}$ data show little difference between the two network systems.

Finally, to transition from the overconstrained-network system back to the original atomistic system, we must reintroduce effects neglected so far: the short-range Si–Si interaction, the long-range Coulomb interaction, and the frustration introduced when unstressed springs are replaced by stressed (stretched or compressed) springs [52, 53, 54, 62, 58, 60, 67, 68]. While their consequences for low-frequency modes are discussed in the next section, here we focus on their impact at high frequencies. As seen in Fig. 2, the overconstrained-network system exhibits a clear separation: the band of modes dominated by stretching and compression of Si–O bonds (up to $4N_{\text{Si}}$ modes) is well separated from the remaining $5N_{\text{Si}}$ modes. In the original system, this separation is less pronounced. A key reason is that the Coulomb contribution, characterized by $\phi^{\prime\prime}_{L}$, reduces the effective stiffness of Si–O bonds: $\phi^{\prime\prime}_{S}=35.8$ eV Å*-2* is lowered to $\phi^{\prime\prime}_{S}+\phi^{\prime\prime}_{L}=35.8-11.1=24.7$ eV Å*-2* (Table 3). This weakening shifts the Si–O bond-stretching/compression band toward lower frequencies. Nevertheless, the trend persists in the original system: the $N\,\delta e_{k,\text{Si--O}}^{2}$ data show that bond stretching and compression continue to dominate the high-frequency region, whereas modes without such bond deformation dominate at low frequencies.

III.2 Marginal stability

We now examine marginal stability in silica glass. Since the concept of marginal stability is well established in packing-based glasses [49, 50, 51, 52, 53, 54, 55, 56], we discuss it alongside data for HSH and LJ glasses. Figure 3 presents $g(\omega)$, $g(\omega)/\omega^{2}$, and the participation ratio $\mathcal{P}_{k}$ (see Materials and Methods). The quantity $\mathcal{P}_{k}$ measures the fraction of particles participating in mode $k$ [89]. As limiting cases, $\mathcal{P}_{k}=1$ corresponds to a fully extended mode in which all particles vibrate equally, whereas $\mathcal{P}_{k}=1/N$ corresponds to a localized mode involving a single particle.

Let us first examine HSH glass in panels (d-f) and LJ glass in panels (g-i). For each system, we also analyze a overconstrained-network counterpart in which interacting pairs (contacts for HSH and neighbors within the cutoff for LJ) are connected by unstressed (neither pre-stretched nor pre-compressed) springs with stiffness set by the pairwise force constants $\phi^{\prime\prime}$. This construction removes pre-stress and the associated frustration present in the original systems and corresponds to the “unstressed” system in previous work [52, 53, 54, 62, 58, 60, 67, 68]. Examining panels (e,h), both overconstrained-network systems (green) exhibit a clear BP. Below the BP frequency $\omega_{\text{BP}}$, $g(\omega)/\omega^{2}$ converges to the Debye level $A_{D}$, i.e., $g(\omega)\to A_{D}\omega^{2}$, and the vibrational states are extended phonons with large $\mathcal{P}_{k}$ 111 Below $\omega_{\text{BP}}$, we also observe modes with small $\mathcal{P}_{k}$ indicative of localization, as seen in panels (f,i). These modes plausibly originate in regions of low local connectivity (low coordination number, i.e., few contacts or interacting neighbors), although a definitive identification requires further study. . Above $\omega_{\text{BP}}$, a band of soft modes of isostatic origin appears. In the isostatic limit, the low-frequency sector comprises both strictly floppy modes at zero frequency and additional soft, nonzero-frequency modes that likewise originate from isostaticity (though they are not classified as floppy). When excess constraints are introduced, both modes are lifted to finite frequencies and merge into a weakly dispersive, non-phononic band. This band remains spatially extended with large $\mathcal{P}_{k}$, and its accumulation produces the excess over the Debye law, i.e., the BP.

Then, when frustration is introduced by replacing unstressed springs with stressed (stretched or compressed) springs, the band of soft modes (of isostatic origin) shifts toward lower frequencies, and the BP accordingly moves downward with decreasing $\omega_{\text{BP}}$ (purple). In addition, QLVs with low $\mathcal{P}_{k}$ emerge at the low-frequency edge below $\omega_{\text{BP}}$. Following previous work [60, 63, 64], we partition modes into extended modes with $\mathcal{P}_{k}>\mathcal{P}_{\text{th}}$ and QLVs with $\mathcal{P}_{k}<\mathcal{P}_{\text{th}}$, and we compute their vDOSs, $g_{\text{EXT}}(\omega)$ and $g_{\text{QLV}}(\omega)$, separately. Here we take $\mathcal{P}_{\text{th}}=10^{-2}$ for HSH and LJ glasses. We observe that $g_{\text{EXT}}(\omega)$ converges to the Debye law $A_{D}\,\omega^{2}$ at a characteristic frequency $\omega_{0}$ (orange curves in panels (e,h)). The same behavior is found for SS glass in Fig. S1 of the SI. Furthermore, for HSH, LJ, and SS glasses, $g_{\text{QLV}}(\omega)$ follows $A_{0}(\omega/\omega_{0})^{4}\propto\omega^{4}$, as shown in Fig. 4, consistent with previous studies [60, 63, 64]. Importantly, the QLVs exhibit a gapless vDOS with a power-law dependence on $\omega$, indicating that the introduction of frustration drives the systems toward marginal stability.

The discussion of HSH and LJ glasses above extends to silica glass. In the isostatic-network system (cyan), the vDOS approaches a nonzero constant as $\omega\to 0$ (i.e., $g(\omega)\propto\omega^{0}$); accordingly, $g(\omega)/\omega^{2}\propto\omega^{-2}\to\infty$ and the BP collapses to zero frequency. The participation-ratio data show that these low-frequency states are spatially extended. Taken together, these results demonstrate the presence of extended floppy modes at zero frequency, together with additional soft, nonzero-frequency modes of isostatic origin [52, 53, 54], as discussed above. Turning to the overconstrained-network system (green), the excess constraints cause the vDOS to vanish as $\omega\to 0$. As seen in panel (b), $g(\omega)/\omega^{2}$ displays a clear BP, analogous to HSH and LJ glasses. Below $\omega_{\text{BP}}$, $g(\omega)/\omega^{2}\to A_{D}$ (hence $g(\omega)\to A_{D}\,\omega^{2}$), and the vibrational states are extended phonons with large $\mathcal{P}_{k}$. Above $\omega_{\text{BP}}$, an isostaticity-derived, weakly dispersive, non-phononic band appears at finite frequencies; it is spatially extended with large $\mathcal{P}_{k}$ and produces the excess over the Debye law, i.e., the BP.

Finally, to move from the overconstrained-network system to the original system, we reinstate all effects neglected in the overconstrained-network construction, including the short-range Si–Si interaction, the long-range Coulomb interaction, and the pre-stress that generates frustration 222 While all of the short-range Si–Si interaction, the long-range Coulomb interaction, and the pre-stress are needed for a faithful description, the essential ingredient that controls the approach to marginal stability is the pre-stress and the associated frustration; the short-range Si–Si and long-range Coulomb terms primarily set the magnitude and sign of the pre-stress and thereby tune the BP position. . With these ingredients restored, the BP shifts to lower frequencies with a reduced $\omega_{\text{BP}}$. Notably, QLVs also emerge at the low-frequency edge below $\omega_{\text{BP}}$. Following the procedure used above for HSH and LJ glasses [60, 63, 64], we compute $g_{\text{EXT}}(\omega)$ and $g_{\text{QLV}}(\omega)$ by setting $\mathcal{P}_{\text{th}}=5\times 10^{-2}$. As shown in Fig. 3(b) and Fig. 4, for $\omega<\omega_{0}$ we find $g_{\text{EXT}}(\omega)\to A_{D}\,\omega^{2}$ i.e., $g(\omega)$ follows the Debye law, while $g_{\text{QLV}}(\omega)$ follows the gapless form $A_{0}(\omega/\omega_{0})^{4}\propto\omega^{4}$. Previous simulations [70, 72] have also reported $g_{\text{QLV}}\propto\omega^{4}$ using different statistical analyses of QLVs 333 Although a quartic law $g_{\text{QLV}}(\omega)\propto\omega^{4}$ is often reported, the precise exponent remains under discussion [105, 106]. What matters here is that, irrespective of the exact exponent, $g_{\text{QLV}}(\omega)$ is gapless and follows a power law as $\omega\to 0$. Such a gapless power-law spectrum is a hallmark of marginal stability. . Taken together, these results show that, once the short-range Si–Si and long-range Coulomb interactions are restored and the associated pre-stress (frustration) is present, the system is driven toward marginal stability.

Silica glass, HS glass, LJ glass, and SS glass all share the characteristic that they become marginally stable states driven by frustration. The key difference lies in the origin of that frustration. In HS, LJ, and SS glasses, it arises from short-range interactions i.e., contact forces and van der Waals forces. By contrast, silica glass involves not only short-range interactions but also long-range Coulomb interactions, and frustration from both contributions collectively drives the approach to marginal stability. Because short-range frustration is mainly associated with repulsive forces, it tends to destabilize the system, whereas the attractive Coulomb interaction between Si and O atoms acts to stabilize and partly offset this effect. Consequently, in silica glass the difference in BP frequency between the overconstrained-network system and the original system is modest, about a factor of $1.5$ (compare the green and purple vertical lines in Fig. 3(b)). In HS, LJ, and SS glasses, the corresponding shift is larger, by a factor of $5$ to $8$.

III.3 Dynamical structure factor

Finally, we examine how isostaticity and marginal stability are reflected in the dynamical structure factor $S_{\alpha}(q,\omega)$ (with $\alpha=T,L$), which is accessible via inelastic scattering experiments. In the low-frequency regime, the vDOS can be estimated by integrating $S_{\alpha}(q,\omega)$ over wavenumber up to the Debye wavenumber $q_{D}$ [14, 27, 28, 30, 61, 55, 56],

[TABLE]

where $k_{B}$ is the Boltzmann constant and $M(\omega)$ is the effective mass (see Materials and Methods). Setting $\omega=\omega_{\text{BP}}$ allows us to resolve how excitations at each $q$ in $S_{T}(q,\omega_{\text{BP}})$ and $S_{L}(q,\omega_{\text{BP}})$ contribute to $g(\omega_{\text{BP}})/\omega_{\text{BP}}^{2}$, i.e., the BP. This viewpoint underlies several theoretical analyses [27, 28, 30, 55, 56, 61], which employ Green’s functions related to $S_{\alpha}(q,\omega)$ via the fluctuation–dissipation theorem.

Figure 5 presents $S_{T}(q,\omega)$ and $S_{L}(q,\omega)$, focusing on the low-frequency and low-wavenumber ($q\lesssim q_{D}$) regime. First, we examine the transverse $S_{T}(q,\omega)$, where data for the overconstrained-network system are shown alongside those for the original system. We begin with HSH glass in panels (d,e) and LJ glass in panels (g,h). For the overconstrained-network (unstressed) systems in panels (e,h), phonon excitations appear at low $q$ and low $\omega$ along the linear dispersion $\omega=c_{T}q$, where $c_{T}$ is the transverse sound speed. In addition, a broad, approximately wavenumber-independent band emerges around and above the BP and extends up to $q_{D}$, indicating non-phononic excitations. This band originates from isostaticity-derived modes that, in the presence of excess constraints, are lifted to finite frequencies (see Fig. 3). Therefore, the BP is built from two components: linearly dispersing phonons that follow $\omega=c_{T}q$, and a wavenumber-independent band of isostaticity-derived modes.

Turning to the original systems that include pre-stress (and hence frustration) in panels (d,g), the entire non-phononic band shifts to lower frequencies, with a corresponding decrease in $\omega_{\text{BP}}$. At the low-frequency edge below $\omega_{\text{BP}}$, a band of QLV excitations appears [60, 63, 64]. Phonon ridges along the linear dispersion $\omega=c_{T}q$ remain visible, but they are broader than in the overconstrained-network systems. This broadening is attributed to isostaticity-derived modes shifting into the low-frequency range and hybridizing with the phonons, which increases the phonon linewidths. Similar behavior is observed for HSL and SS glasses, as shown in Fig. S3 of the SI. In this way, marginal stability manifests as a wavenumber-independent band that is shifted downward by frustration, together with the emergence of QLV excitations at the low-frequency edge. Notably, a broad, wavenumber-independent band in the dynamical structure factor has also been reported recently in simulations of LJ glasses [104], consistent with the present results.

The discussion above for HSH and LJ glasses extends to silica glass in panels (a,b) of Fig. 5. For silica glass, additional data for the isostatic-network system are provided in Fig. S2 of the SI. First, in the isostatic-network system in Fig. S2, exact isostaticity generates a nearly wavenumber-independent band that reaches the zero-frequency limit; correspondingly, $S_{T}(q,\omega)$ accumulates substantial, nearly $q$-independent spectral weight as $\omega\to 0$, while $g(\omega)/\omega^{2}\to\infty$ as indicated by Eq. (21). In the overconstrained-network system in Fig. 5(b), excess constraints lift this isostaticity-derived band to finite frequencies around the BP, and a phonon branch appears along $\omega=c_{T}q$. In the original system in Fig. 5(a), restoring pre-stress from short-range and long-range Coulomb interactions shifts the entire non-phononic band to lower frequencies, and $\omega_{\text{BP}}$ decreases accordingly. At the low-frequency edge below $\omega_{\text{BP}}$, QLVs form a broad band. As discussed above, in silica glass the attractive Coulomb interaction between Si and O stabilizes the structure and partially offsets the destabilizing effect of short-range frustration. Consequently, the downward shift of the band is smaller in silica glass than in HSH, LJ, and SS glasses. Despite these quantitative differences, the underlying physics is the same across systems.

An intriguing aspect is the behavior of the longitudinal channel, $S_{L}(q,\omega)$. Across silica, HSH, LJ, and SS glasses, $S_{L}(q,\omega)$ is markedly weaker than $S_{T}(q,\omega)$ in the low-frequency regime, including the BP region, indicating that the BP predominantly reflects transverse motion [20, 29, 21, 61]. In silica glass, however, $S_{L}(q,\omega)$ shows a similar $\omega$–$q$ structure to $S_{T}(q,\omega)$: a broad, nearly wavenumber-independent band appears around the BP as shown in Fig. 5(c), whereas the corresponding feature is much less apparent in HSH and LJ glasses (Figs. 5(f) and (i)) and in SS glass (Fig. S3). This behavior likely reflects closer proximity to isostaticity in silica glass, which enhances the population of isostaticity-derived soft modes with both transverse and longitudinal character. Moreover, silica glass exhibits strong nonaffine elasticity not only under shear but also under volumetric deformation (see Materials and Methods), further reducing the contrast between longitudinal and transverse responses.

Since scattering experiments probe the longitudinal channel, we now compare our simulations with inelastic X-ray scattering data [12, 13, 14, 16]. Figure 6 shows the measured longitudinal dynamical structure factor for ambient-pressure silica glass at density $\rho=2.20$ g/cm3 and temperature $T=1620$ K. The plotted quantity is the inelastic part of the experimentally determined dynamical structure factor, $S_{L}(q,\omega)/\left[\hbar\,\omega\,(n(\omega,T)+1)\right]$ (which reduces to $S_{L}(q,\omega)/(k_{B}T)$ in the classical limit), in units of eV THz*-1*, where $n(\omega,T)$ is the Bose–Einstein population factor and $\hbar=h/2\pi$ with $h$ the Planck constant. Details on the procedure used to properly normalize the inelastic X-ray scattering data are provided in Ref. [16]. The data are limited to frequencies $\omega\gtrsim 1$ THz because subtracting the elastic line becomes delicate at lower frequencies. We also note that our simulation data in Fig.5(c) correspond to $T=0$ K; accordingly, the BP frequency is $\omega_{\text{BP}}\approx 1.21$ THz in simulation but $\omega_{\text{BP}}\approx 1.45$ THz in the experiment, consistent with the expected upward shift at elevated temperature $T=1620$ K. Despite this offset, the qualitative behavior is fully consistent between simulation and experiment: (i) at low wavenumber and frequency, longitudinal phonons follow the linear dispersion $\omega=c_{L}q$; and (ii) a broad, wavenumber-independent (dispersionless) non-phononic band is visible around the BP, extending from $q\approx 0.3$Å*-1* up to $\approx 1.2$ Å*-1*. By direct analogy with the simulations, we identify this dispersionless band in the inelastic X-ray scattering data as the band of isostaticity-derived modes.

IV Conclusion

We have explained the BP in silica glass in terms of isostaticity and marginal stability. When the tetrahedral network structure is extracted, it forms an isostatic-network system where the number of constraints equals the number of degrees of freedom. This isostatic network comprises both strictly zero-frequency floppy modes and many additional soft, low-frequency modes of isostatic origin, and the vDOS remains finite as $\omega\to 0$, i.e., $g(\omega)\propto\omega^{0}$. In practice, interunit couplings of van der Waals type add constraints beyond the degrees of freedom. The resulting spring-network is overconstrained and displays phonons at low frequency whose vDOS follows the Debye law, i.e., $g(\omega)\simeq A_{D}\omega^{2}$, while the isostaticity-derived soft modes are shifted to finite frequencies and form the non-phononic excess that constitutes the BP. This behavior of the overconstrained-network system is the same mechanism operative in HS, LJ, and SS glasses. Reinstating all effects neglected in the overconstrained-network construction, including short-range Si–Si interactions, long-range Coulomb interactions, and the associated pre-stress that generates frustration, shifts the non-phononic band downward, producing a low-frequency BP. Concurrently, QLVs appear at the low-frequency edge below the BP and exhibit a gapless power-law vDOS $\propto\omega^{4}$ as $\omega\to 0$. Together, these features indicate that frustration drives silica glass to a marginally stable state, as also observed in HS, LJ, and SS glasses.

Furthermore, we find that isostaticity and marginal stability are encoded in the dynamical structure factor as a broad, nearly wavenumber-independent band around the BP. In this frequency range, isostaticity-derived soft modes produce a dispersionless band. In silica glass, the band is sharply visible in both the transverse and longitudinal dynamical structure factor, whereas in HS, LJ, and SS glasses it is prominent in the transverse channel but much less apparent in the longitudinal one. This contrast is especially useful because scattering experiments generally cannot access the transverse component at low frequencies and instead probe the longitudinal channel. Indeed, a direct comparison between inelastic X-ray scattering measurements and our simulations of silica glass shows good agreement in the longitudinal dynamical structure factor, including both the low-wavenumber phonon branch and the broad dispersionless band around the BP. On this basis, we identify the dispersionless band observed in the inelastic X-ray scattering data as the band of isostaticity-derived modes. Thus, the BP can be analyzed directly via the longitudinal component of the dynamical structure factor, informing interpretations of existing measurements [4, 5, 6, 7, 8, 10, 11, 12, 13, 14, 15, 16] and guiding future investigations.

In conclusion, isostaticity and marginal stability are broadly applicable to both packing-type and network-forming glasses and capture fundamental aspects of the physics of amorphous solids. These principles manifest as a dispersionless, nearly wavenumber-independent band in the dynamical structure factor (including the longitudinal channel accessible to inelastic scattering), providing a direct, testable signature. At the end of the present work, we carried out an effective-medium mean-field analysis based on random spring networks, following Refs. [55, 56]. Within this framework we computed the vDOS and the dynamical structure factor; the results and explanations are reported in the SI. For the overconstrained-network system (often termed the “unstressed system” in the effective-medium literature), the analysis predicts an accumulation of isostaticity-derived soft modes at finite frequencies that produces the BP. When pre-stress is included to emulate the frustration present in the original system, these soft modes shift to lower frequencies and the BP moves downward. Crucially, the theory predicts that these soft modes imprint themselves on the dynamical structure factor as a broad, wavenumber-independent band around the BP. These effective-medium predictions are in excellent agreement with our simulations and with the inelastic X-ray scattering data on silica glass, establishing a consistent picture across theory, simulation, and experiment.

Acknowledgments

We thank Atsushi Ikeda, Kumpei Shiraishi, Osamu Yamamuro, and Maiko Kofu for useful discussions. This work was supported by JSPS KAKENHI Grant Numbers 22K03543, 23H04495, 25H01519.

Author contributions statement

H.M. designed the research with help from T.M. and E.M.; E.M. performed molecular dynamics simulations of silica glass; G.B. acquired and curated the experimental inelastic X-ray scattering data; H.M. conducted the theoretical analysis and wrote the paper; and all authors discussed the results and commented on the manuscript.

Author declarations

The authors declare no conflicts of interest.

Supplementary data

In the following, we report supplementary data, including vibrational states in HSL and SS glasses (Fig. S1), dynamical structure factors in the isostatic-network system of silica glass (Fig. S2), dynamical structure factors in HSL and SS glasses (Fig. S3), and physical quantities including elastic moduli and Debye values (Table S1).

Effective-medium mean-field analysis

Here, we carry out an effective-medium theory (EMT) analysis based on random spring networks, following Refs.[55, 56]. In this framework, particles (nodes) are connected by linear springs to form a random network. Two parameters control stability: the connectivity (contact number) $z$ and the level of pre-stress $e~(>0)$, which quantifies the internal forces carried by the springs. Mechanical stability requires $z$ to exceed the isostatic threshold $z_{c}=2d$ in $d$ dimensions, i.e., the excess contact number $\delta z=z-z_{c}>0$. In addition, if $e$ exceeds a critical value $e_{c}~(>0)$, the network becomes unstable; hence stability requires $e\leq e_{c}$. At $e=e_{c}$ the system lies exactly at the stability boundary, i.e., it is marginally stable. However, this “marginal stability” is in the sense of mean-field limit, and it does not match with the marginal stability of the quenched glasses studied in our simulations and experiments. Under an effective-medium (mean-field) approximation, one can compute the complex, frequency-dependent effective spring constant $k_{\rm eff}(\omega)$ that characterizes the elastic response of the network. For detailed derivations and the full formulation, see Refs. [55, 56].

Given $k_{\rm eff}(\omega)$, the vDOS $g(\omega)$ and the dynamical structure factor $S(q,\omega)$ are computed via

[TABLE]

where $q_{D}$ denotes the Debye wavenumber, and $\mathrm{Im}$ denotes the imaginary part. In this theory longitudinal and transverse polarizations are not distinguished, so $S(q,\omega)=S_{L}(q,\omega)=S_{T}(q,\omega)$. The two expressions obey

[TABLE]

which corresponds to Eq. (3) in the main text upon identifying $3S\!\to\!S_{T}+S_{L}$ (sum over two transverse and one longitudinal polarizations) and, in the EMT setting with identical nodes, $m$ playing the role of the effective mass. In what follows, we set the Debye wavenumber to $q_{D}=4$ (dimensionless units of the theory) and $m=1$.

First, to eliminate pre-stress we set $e=0$ and tune the connectivity to $z=z_{c}$. This corresponds to the isostatic-network system. Figure S4 shows $g(\omega)$ in (a) and $g(\omega)/\omega^{2}$ in (b) (cyan curves). As $\omega\to 0$, $g(\omega)$ remains finite and $g(\omega)/\omega^{2}\to\infty$. Thus, the low-frequency spectrum at isostaticity comprises strictly zero-frequency floppy modes together with many additional soft, low-frequency modes of isostatic origin.

Next, we introduce excess constraints by setting $\delta z=z-z_{c}=10^{-2}$ while keeping $e=0$. This corresponds to the overconstrained-network system (often termed the “unstressed system” in the effective-medium literature). In this case $g(\omega)\to 0$ as $\omega\to 0$, and $g(\omega)/\omega^{2}$ converges to the Debye level $A_{D}$, i.e., $g(\omega)\simeq A_{D}\,\omega^{2}$ (green curves). Introducing excess constraints lifts both the strictly floppy modes and the additional isostaticity-derived soft modes to finite frequencies, where they merge into a weakly dispersive, non-phononic band that forms the BP.

Finally, we include pre-stress to mimic frustration by setting $e=e_{c}(1-5\times 10^{-4})$ with $\delta z=10^{-2}$ held fixed. This choice indicates $1-e/e_{c}=5\times 10^{-4}$, i.e., the network is tuned arbitrarily close to the marginally stable point. As noted above, the EMT treats vibrational modes as spatially extended and therefore does not capture localization or the characteristic $\omega^{4}$ vDOS of quasi-localized vibrations; consequently, the EMT notion of “marginal stability” differs from that realized in quenched glasses in simulations and experiments. In particular, at exactly $e=e_{c}$ the mean-field BP collapses to $\omega_{\text{BP}}\to 0$, which disagrees with simulations and experiments. To enable a meaningful comparison with the original system in the main text (the quenched glasses studied in our simulations and experiments), we therefore adopt the slightly subcritical choice $e=e_{c}(1-5\times 10^{-4})$. Under this pre-stress, the isostaticity-derived soft-mode band shifts toward lower frequencies, and the BP correspondingly moves downward as $\omega_{\text{BP}}$ decreases (purple curves), in agreement with our simulations.

The same trends appear in the dynamical structure factor. Figure S5 shows $S(q,\omega)/k_{B}T$ for the original system ($\delta z=10^{-2}$, $e=e_{c}(1-5\times 10^{-4})$) in (a) and for the overconstrained-network system ($\delta z=10^{-2}$, $e=0$) in (b). At low wavenumbers and frequencies, phonons follow the linear dispersion $\omega=c\,q$, where $c$ is the sound speed. In addition, a broad, wavenumber-independent (dispersionless) non-phononic band appears around the BP in both cases. Compared with the overconstrained-network case, pre-stress in the original system shifts this band to lower frequencies; this downward shift enhances hybridization between isostaticity-derived soft modes and phonons, thereby broadening the phonon ridge along $\omega=c\,q$, consistent with our simulation results.

Taken together, these EMT predictions account for, and are consistent with, our simulations and the experimental observations, establishing a coherent picture across theory, simulation, and experiment.