# Fileset

[05_clean_manuscript_LP18921B_second_revision.pdf](https://mdr.nims.go.jp/filesets/908d3b42-08ec-4997-b490-d4780b24b364/download)

## Creator

[Tatsuya Mori](https://orcid.org/0000-0002-9443-594X), [Hideyuki Mizuno](https://orcid.org/0000-0002-2866-7400), [Yuzuki Motokawa](https://orcid.org/0009-0005-2439-3405), Dan Kyotani, [Soo Han Oh](https://orcid.org/0009-0002-9521-2606), [Yasuhiro Fujii](https://orcid.org/0000-0002-9522-6103), [Akitoshi Koreeda](https://orcid.org/0000-0002-3063-5563), [Shinji Kohara](https://orcid.org/0000-0001-9596-2680), [Seiji Kojima](https://orcid.org/0000-0003-1933-9269)

## Rights

©2026 American Physical Society[In Copyright](http://rightsstatements.org/vocab/InC/1.0/)

## Other metadata

[Continuum model for the terahertz dielectric response of glasses](https://mdr.nims.go.jp/datasets/278f2b5d-bc81-4f74-bd21-291b6c8cecc7)

## Fulltext

Continuum model for the terahertz dielectric response of glassesTatsuya Mori,1, ∗ Hideyuki Mizuno,2 Yuzuki Motokawa,1 Dan Kyotani,1 Soo HanOh,1 Yasuhiro Fujii,3, 4 Akitoshi Koreeda,5 Shinji Kohara,6 and Seiji Kojima11Department of Materials Science, University of Tsukuba,1-1-1 Tennodai, Tsukuba, Ibaraki 305-8573, Japan2Graduate School of Arts and Sciences, The University of Tokyo,3-8-1 Komaba, Meguro-ku, Tokyo 153-8902, Japan3Institute for Open and Transdisciplinary Research Initiatives,Osaka University, 2-1 Yamada-oka, Suita, Osaka 565-0871, Japan4Research Organization of Science and Technology, Ritsumeikan University,1-1-1 Noji-higashi, Kusatsu, Shiga 525-8577, Japan5Department of Physical Sciences, Ritsumeikan University,1-1-1 Noji-higashi, Kusatsu, Shiga 525-8577, Japan6Center for Basic Research on Materials, National Institute for MaterialsScience (NIMS), 1-2-1 Sengen, Tsukuba, Ibaraki 305-0047, JapanBoson peak dynamics in glasses produce a robust crossover in the terahertz (THz) dielectricresponse that standard Debye or Lorentz models do not capture. We develop a continuum descriptionof this THz response, coupling an infrared-effective charge fluctuation spectrum to a frequency-dependent shear modulus, and apply it to glycerol glass. The model reproduces the measuredcomplex dielectric function and the nearly linear infrared light-vibration coupling around the bosonpeak, and highlights the dominant role of transverse shear dynamics.I. INTRODUCTIONIn disordered solids, a universal excess of vibrationalmodes appears in the terahertz (THz) range—the bo-son peak (BP)—which is observed as a peak in the re-duced vibrational density of states (VDOS) g(ω)/ωD−1(where D is the spatial dimension) relative to the Debyelevel [1, 2]. Despite decades of work, the microscopicorigin of the BP remains an open problem that con-tinues to stimulate experiments [3–12], simulations [13–19], and theory [20–27]. The BP is tied to several hall-mark anomalies of glasses: unusually low thermal con-ductivity [3], nanoscale plasticity [14], and a markeddecrease in transmittance in the THz region above theBP frequency [11, 12]. Understanding the BP and itscoupling to light is essential for fundamental physicsand for predicting and controlling THz optical proper-ties. In particular, emerging THz communication tech-nologies demand window and substrate materials withlow permittivity and low loss [28]. Motivated by theseconsiderations, we focus on the dielectric response inthe THz regime and develop a framework that con-nects the observed absorption to the underlying vibra-tional dynamics. Empirically, the THz dielectric re-sponse ε(ω) in glasses shows a robust pattern that stan-dard models fail to capture: below the BP frequencyωBP the response is resonance-like (Lorentz-oscillator-like), whereas above ωBP it becomes relaxational (Debye-like) [11, 12]. Neither standard relaxational dielectricfunctions—Debye and common empirical generalizations(e.g., the Havriliak–Negami form [29])—nor an over-damped Lorentz oscillator reproduces this crossover; a∗ mori@ims.tsukuba.ac.jpunified analytic form for ε(ω) that captures this univer-sal crossover has been lacking.Several beyond-Debye/Lorentz approaches have beenpursued for the THz dielectric anomaly around the BP.Zaccone and co-workers formulated a phonon-polaritonicresponse by extending a Lorentz oscillator with awavenumber-dependent damping Γ(k), but the choice ofΓ(k) limits transferability across different glasses [30]. Adifferent route links dielectric relaxation to viscoelastic-ity: in the Gemant–DiMarzio–Bishop model the constantDebye viscosity is replaced by a frequency-dependentone, which can be expressed via a complex shear modu-lus [31, 32]. However, the model was devised for liquid-side relaxation and carries an explicit temperature depen-dence, making a direct application to low-temperatureglass spectra nontrivial. Along this line, Zeitler et al.discussed the scaling g(ω)/ω2 ∝ α(ω)/ω3 (with α(ω)the absorption coefficient) [33]; nevertheless, a unified,material-independent ε(ω) capturing the BP crossoverremains elusive.Light–vibration coupling coefficients in the BP re-gion are spectroscopic-channel dependent. Hyper-Ramanscattering (HRS) studies have shown an approximatelyconstant CHRS(ω) around the BP in vitreous silica [34,35] and have highlighted the tensor-rank dependence ofthe coupling in network glasses [35, 36], whereas ordi-nary Raman measurements have reported a nearly linearCRaman(ω) near the BP in a broad set of glasses [37];analysis of our Raman data for glycerol likewise yieldsan approximately linear CRaman(ω) in the BP region(see Appendix M). For the infrared (IR) channel con-sidered here, Taraskin et al., building on Maradudin’slinear-response theory for optical phonons in crystals [38]and decomposing the effective atomic-charge fluctuationsinto correlated (crystal-like) and uncorrelated (disorder-2induced) parts, derived the analytic form for the IR light–vibration coupling coefficient CIR(ω) = α(ω)/g(ω) =A + Bω2, which reproduces the IR absorption of sil-ica and several other oxide glasses [39]. However, forglycerol the experimentally extracted CIR(ω) displays apronounced linear dependence in the vicinity of the BP,clearly deviating from A + Bω2 (see Fig. 3(b)). Simi-lar linear trends have been reported for a broad rangeof amorphous materials, including hydrogen-bonded [11]and polymeric glasses [40], proteins [12], and even someinorganic glasses [41]. These IR observations call for aformulation that explains the linear CIR(ω) while remain-ing consistent with the measured complex dielectric func-tion ε(ω) (both its real and imaginary parts, ε′(ω) andε′′(ω)).In this work, we introduce a continuum model thatexplicitly incorporates the IR-effective charge fluctuationspectrum ∆q(k) and the effective complex shear modulusG(ω). With a low-k form ∆q(k) ≈ q0 + q2k2, the modelsimultaneously reproduces ε′(ω) and ε′′(ω) and accountsfor the linear frequency dependence of CIR(ω) in theBP region, in quantitative agreement with experimentaldata. Moreover, by increasing the constant-to-quadraticratio R ≡ q0/(q2k2D) (kD is the Debye wavenumber), theresponse continuously approaches a Taraskin-like form,so the present framework both explains the linear regimeand contains the Taraskin expression as a limiting case.II. CONTINUUM MODELA. Microscopic origin of the electric driving termBefore writing the continuum equation of motion, webriefly derive the electric source term from the micro-scopic rigid-ion dipole. In a rigid-ion picture, each atomor site i is characterized by a time-independent rigid-ioncharge qi, an equilibrium positionR0i , and a displacementui(t) from R0i . Its instantaneous position is thenRi(t) = R0i + ui(t), (1)and the microscopic dipole moment is the vectord(t) =∑iqiRi(t)=∑iqiR0i +∑iqiui(t) ≡ d0 + dvib(t).(2)Here d0 is the static dipole associated with the equilib-rium reference configuration, and dvib(t) is the vibra-tional contribution to the dipole. The first term is staticand does not generate a force on the vibrational coordi-nates; the field-induced vibrational response is governedbydvib(t) =∑iqiui(t). (3)The corresponding local IR-effective charge weight in theequilibrium reference configuration may be decomposedas q(r) = qav +∆q(r), where qav is its spatially uniformcomponent. For the charge-neutral rigid-ion system con-sidered here, qav = 0, leaving the fluctuation field ∆q(r).The continuum-field representation of this discrete sumis thereforedvib(t) ≃∫d3r∆q(r)u(r, t). (4)Here r denotes a position in the equilibrium referenceconfiguration, u(r, t) is the continuum displacement field,and ∆q(r) denotes the local IR-effective charge fluctua-tion field in that reference configuration. Equation (4)represents the discrete sum in Eq. (3) in continuum-fieldform: the spatial dependence of the IR-effective chargefluctuations is retained in ∆q(r), while the site displace-ments ui(t) are represented by u(r, t). The approximateequality in Eq. (4) reflects this discrete-to-continuum re-placement. The interaction energy between a spatiallyuniform external electric field E(t) and this vibrationaldipole isHvibint (t) = −E(t) · dvib(t)= −∫d3r∆q(r)E(t) · u(r, t).(5)Taking the functional derivative with respect to the dis-placement field u(r, t) gives the electric force densityfE(r, t) = − δHvibintδu(r, t)= ∆q(r)E(t). (6)After Fourier transforming in time, this becomesfE(r, ω) = ∆q(r)E(ω), where ω denotes angular fre-quency and E(ω) is the corresponding Fourier componentof the external electric field. This is the source term usedin Eq. (7). Appendix A gives the corresponding discrete-to-continuum derivation of this electric source term.B. Equation of motion and dielectric responseWe describe the THz dielectric response of a glassysolid using an isotropic viscoelastic continuum with bulkmodulus K, complex shear modulus G(ω), and the lo-cal IR-effective charge fluctuation field ∆q(r) introducedabove. In the frequency domain, the displacement fieldu(r, ω) driven by a spatially uniform THz electric fieldE(ω) = E(ω)ê obeys−ρω2u(r, ω) =(K +43G(ω))∇(∇·u(r, ω))−G(ω)∇×(∇× u(r, ω))+ ∆q(r)E(ω).(7)3Here ρ is the mass density. We use the Fourier represen-tationu(r, ω) =1√V∑k,αuα(k, ω) êα(k)eik·r,∆q(r) =1√V∑k∆q(k)eik·r,(8)where ∆q(r) is real, so that ∆q(−k) = ∆q(k)∗. Foreach wave vector, the three acoustic branches α =L, T1, T2 have polarization vectors êα(k) (êL(k) = k/|k|,êTi(k) · k = 0). The uniform field couples to branch αthrough the scalar projection Eα(ω) = êα(k) · E(ω) =(ê · êα(k))E(ω). Fourier transforming Eq. (7) and pro-jecting onto branch α give the mode equation(−ρω2 + Cα(ω)k2)uα(k, ω) = ∆q(k)Eα(ω), (9)with k = |k|, CL(ω) = M(ω) = K + 43G(ω), andCT1(ω) = CT2(ω) = G(ω).The field-parallel macroscopic polarization isP (ω) =1V∫d3r∆q(r) ê · u(r, ω). (10)Combining this definition with the mode response ob-tained from Eq. (9) gives the mode-sum formP (ω)E(ω)=1V∑k,α|∆q(k)|2(ê · êα(k))2−ρω2 + Cα(ω)k2. (11)We assume statistically homogeneous and isotropic IR-effective charge fluctuations. With the overbar denotinga disorder average, this gives∆q(k) = 0, ∆q(k)∆q(k′)∗ = δk,k′〈|∆q(k)|2〉dis,(12)where〈|∆q(k)|2〉disdepends only on k = |k|. We also de-fine the orientational factor γα =〈(ê · êα)2〉angles. Afterthe disorder and orientational averages,P (ω)E(ω)=1V∑k,α〈|∆q(k)|2〉disγα−ρω2 + Cα(ω)k2. (13)Finally, replacing the discrete sum by the isotropicDebye-sphere form1V∑kf(k) →∫ kD0f(k)3k2k3Ddk (14)and using ε(ω) = ε∞ + P (ω)/[ε0E(ω)] givesε(ω) = ε∞ +1ρε0∑α∫ kD0〈|∆q(k)|2〉disγα−ω2 + C̃α(ω)k23k2k3Ddk, (15)where C̃α(ω) = Cα(ω)/ρ. Appendices B–E give the cor-responding algebra and the branch-resolved form. Inan isotropic three-dimensional elastic medium, the threeacoustic polarizations form an orthonormal triad andare statistically equivalent, so that γL = γT1= γT2=1/3. Grouping the two degenerate transverse branchesand assuming that〈|∆q(k)|2〉disis branch-independent,Eq. (15) reduces toε(ω) = ε∞ +1ρε0∫ kD0[23〈|∆q(k)|2〉dis−ω2 + G̃(ω)k2+13〈|∆q(k)|2〉dis−ω2 + M̃(ω)k2]3k2k3Ddk.(16)where G̃(ω) = G(ω)/ρ and M̃(ω) = M(ω)/ρ, withM(ω) = K + (4/3)G(ω). Equation (16) is the expres-sion used to analyze the THz dielectric response and theBP-related excess contribution. It should be regardedas a factorized effective-medium closure: the IR-effectivecharge-coupling spectrum〈|∆q(k)|2〉disis introduced in-dependently of the effective mechanical response specifiedby G(ω) and the bulk modulus K. Possible connectedcorrelations between local charge-coupling environmentsand the microscopic mechanical disorder underlying theBP-related response are therefore not explicitly includedin this closure. Such correlations would generate non-factorizable corrections to Eq. (16) and are beyond thescope of the present minimal model.III. APPLICATION TO GLYCEROL GLASSFor quantitative comparison with experiment, themodel inputs are ρ, kD, G(ω), K, and ∆q(k). Forglycerol we use an effective shear modulus G(ω) ob-tained by applying the heterogeneous-elasticity-theorycoherent-potential-approximation (HET-CPA) proce-dure described in Ref. [50] to the experimental vibra-tional density of states g(ω) of glycerol glass at 80 K,and parameterize the IR-effective charge spectrum byretaining the leading terms of a Maclaurin expansion,∆q(k) ≈ q0 + q2k2. The constant term q0 is weakly con-strained by the glycerol dielectric spectra: when q0 andq2 are varied together, the fitted q0 remains consistentwith zero within the fitting uncertainty. In the mainanalysis we use the minimal form ∆q(k) = q2k2. Thecoefficient q2 is selected by a parameter scan to minimizethe misfit between calculated and measured ε′(ω) andε′′(ω). Details of the q0 uncertainty and the q2 scan aregiven in Appendix J. The parameters used in Fig. 1 aresummarized in Table I.The complex dielectric function ε(ω) of glycerol glasswas measured by terahertz time-domain spectroscopy(THz-TDS) [11, 12] using a commercial system (RT-10000, Tochigi Nikon Co.). The usable frequency range inthe present THz-TDS measurements was approximately0.2–2.5 THz. The THz beam path was enclosed in dryair, and the sample was mounted in a liquid-helium flowcryostat for temperature control between 14 and 295 K.41     0.0 0.5 1.0 1.5 2.0 2.52.83.03.23.43.63.84.04.2e'Frequency (THz)Glycerol (80 K)  THz-TDS model0.00.20.40.60.81.01.21.4e'' FIG. 1. Complex dielectric response of glycerol glass at 80 K.Symbols: THz-TDS measurements; solid curves: model.Plotted are ε′(ω) and ε′′(ω) versus frequency (THz). Themodel reproduces the crossover from resonance-like behav-ior below ωBP to a broad Debye-like response above ωBP.The weak convex-upward feature in ε′(ω) near ωBP is con-sistent with a shallow dip in G′(ω), i.e., a reduction ofVTA(ω) =√G′(ω)/ρ associated with the BP.TABLE I. Parameters used to calculate the dielectric responseof glycerol glass at 80 K in Fig. 1. ρ is the mass density,kD the Debye wavenumber, K the bulk modulus, and ε∞the high-frequency dielectric constant. q0 and q2 specify thequadratic parameterization of the IR-effective charge spec-trum ∆q(k) = q0+q2k2. The fixed mechanical inputs and theHET-CPA parameter set defining the frequency-dependentshear modulus G(ω) are described in Appendix G; the sound-velocity inputs are from Ref. [53], and G(ω) is shown in Fig. 2.ρ(g cm−3)kD(Å−1)K(GPa)ε∞q0(C cm−3)q2(C cm−3 Å2)Glycerol 1.26 1.90 10.58 2.765 0 1.6× 103Glycerol (≥ 99.5%, Sigma-Aldrich) was sealed in a liq-uid cell with a 0.2 mm optical thickness between twoz-cut α-quartz windows. The sample was vitrified in situby cooling through the glass transition (Tg ≃ 190 K [51]);the data analyzed correspond to the glassy state at 80 Kafter equilibration. Optical constants were extractedfrom a standard slab-transmission analysis to obtain thecomplex refractive index n(ω), and the complex dielectricfunction was then calculated via ε(ω) = n2(ω).Figure 1 compares the measured complex dielectricresponse of glycerol glass with the model. Over 0.3–2.5 THz, the model simultaneously reproduces both thereal and imaginary parts, ε′(ω) and ε′′(ω). Aroundthe BP frequency ωBP the response crosses over fromresonance-like behavior below ωBP to a broad, Debye-like lineshape above ωBP, in quantitative agreement with5     0.40.60.81.0Re[G(n)/G0] Re[G(n)/G0](a)0.00.20.40.6Im[G(n)/G0] Im[G(n)/G0]0.0 0.5 1.0 1.5 2.0 2.50.0000.0050.0100.0150.020g(n)/n2Frequency (THz)Glycerol (80 K)(b) INS CPA FIG. 2. (a) Real and imaginary parts of the effectiveshear modulus G(ω) from the HET-CPA analysis, plotted asRe[G(ω)/G0] (orange) and Im[G(ω)/G0] (blue), with G0 =5.43 GPa. (b) Reduced vibrational density of states g(ω)/ω2computed from G(ω) (magenta line) compared with INSdata [52] (filled circles). Here ν denotes the linear frequency,ν = ω/2π.experiment. This agreement is naturally explained if op-tical attenuation is mediated by the effective complexshear modulus G(ω): Rayleigh-like scattering of long-wavelength acoustic modes prevails below ωBP, whereasabove ωBP the dissipation associated with effective elas-tic heterogeneities is captured by a rapid increase of theviscoelastic loss modulus G′′(ω), which in turn producesthe strong IR absorption and broad Debye-like lineshape.The small convex-upward bump in ε′(ω) near ωBP isaccounted for by a shallow dip in the storage modu-lus G′(ω) (Fig. 2), which lowers the transverse soundspeed VTA(ω) =√G′(ω)/ρ, and slightly flattens the ω–kdispersion, yielding an excess VDOS in the BP region;accordingly, the ε′(ω) bump reflects the BP-related ex-cess. Within our isotropic three-dimensional formula-tion, Eq. (16), the transverse part of ε(ω) is governeddirectly by the shear modulus G(ω), whereas the longi-tudinal part involves the longitudinal modulus M(ω) =K+ 43G(ω). In the present model, K is a real, frequency-independent constant; thus all frequency dependence anddissipation of M(ω) originate from G(ω), while the con-stantK raises the longitudinal stiffness. A decompositionof the modeled dielectric function into transverse and lon-gitudinal contributions (Fig. 6) shows that, in the BP re-52   0.0000.0050.0100.0150.0200.025g(n)/n2 (THz-3)(a) INSGlycerol (80 K)0102030405060a(n)/n2 (cm-1THz-2) THz-TDS model0.0 0.5 1.0 1.5 2.00.00.51.01.52.02.5CIR(n) (103cm-1THz)Frequency (THz)(b)Glycerol (80 K)nBP-IR = 1.40 THznBP-INS = 1.00 THz Exp model FIG. 3. Reduced spectra and IR light-vibration coupling. (a)Reduced absorption α(ω)/ω2 (blue symbols) obtained fromFig. 1 via α(ω) = ωε′′(ω)/(c n′(ω)) (with c the speed of lightand n′(ω) the real part of the refractive index), together withthe reduced VDOS g(ω)/ω2 (orange symbols) from indepen-dent inelastic neutron scattering (INS) [52]; solid lines arethe model. (b) CIR(ω) = α(ω)/g(ω). The near-linear depen-dence around the BP is captured by the model, whereas thequadratic Taraskin form A+Bω2 does not capture this trend.Vertical dashed lines mark the BP positions from IR and INS,ωBP-IR/2π ≈ 1.40 THz and ωBP-INS/2π ≈ 1.00 THz, respec-tively.gion, both ε′(ω) and ε′′(ω) are almost entirely controlledby the transverse channel, with the longitudinal contribu-tion remaining comparatively small over the whole band.Thus the BP-related THz dielectric response is controlledprimarily by the frequency dependence of G(ω) and theassociated transverse-like vibrations. At the same time,a fully microscopic derivation of the dielectric responseremains a challenge.Using the modeled dielectric function, we convertε′′(ω) to the absorption coefficient α(ω) and, from thesame effective modulus G(ω), obtain g(ω); the reducedspectra α(ω)/ω2 and g(ω)/ω2 are plotted in Fig. 3(a).We then form CIR(ω) = α(ω)/g(ω) and compare withexperiment. As shown in Fig. 3(b), the resulting CIR(ω)agrees with the data in the vicinity of the BP and cap-tures the near-linear trend there, whereas the Taraskinform A+Bω2 [39] cannot account for this trend.IV. STATIC CHARGE CORRELATIONS ANDCOUPLING SPECTRUMIn Eq. (16) all microscopic information on the THz(IR) charge coupling is collected into the k-dependentoscillator strength ⟨|∆q(k)|2⟩dis. Since the k integralis restricted to the acoustic window 0 < k < kD, weonly require ∆q(k) over this range and adopt the mini-mal analytic form ∆q(k) ≈ q0 + q2k2. Microscopically,the mode oscillator strength is set by the mode effec-tive charge F(p) [54], i.e., the projection of the rigid-ion charge distribution onto each eigenmode p (equiv-alently by the charge–charge dynamic structure factorSZZ(k, ω) and its projection onto the eigenmodes [55,Chap. 10]). Appendix K relates this mode-resolvedquantity to static charge correlations in a limited sense.There, disordered eigenmodes are projected onto lon-gitudinal and transverse plane-wave components withmode-dependent k-space weights akα(ωp) [39, 56]. Be-fore further approximation, the squared projected modecharge contains the configurational four-point average⟨akα(ωp)a∗k′α′(ωp)ρZ(k)ρZ(k′)∗⟩conf . Under a diagonal,factorized approximation—neglecting cross/interferenceterms and representing the remaining charge correla-tions by the static structure factor SZZ(k)—this gives〈|ê · F(p)|2〉∝∑k,α |akα(ωp)|2SZZ(k) γα. This con-struction provides a microscopic route and diagnosticby which a static SZZ(k) can contribute to the effec-tive coupling spectrum represented by ⟨|∆q(k)|2⟩dis inEq. (16), rather than a full microscopic determination ofthat spectrum. This diagonal/factorized interpretationis expected to be most transparent in the low-frequencypropagon regime, where the k-space weights are relativelynarrow; above the BP, the weights broaden as the modesbecome diffuson-like, so the neglected off-diagonal termsand possible connected correlations may become moreimportant. A fully microscopic determination of this ef-fective charge-coupling spectrum would require dynamic,mode-resolved charge-correlation information, such asSZZ(k, ω) or an equivalent eigenmode-resolved effective-charge analysis, but the static SZZ(k) from molecular-dynamics (MD) suffices to diagnose whether medium-range correlations (e.g., the first sharp diffraction peak(FSDP)) should appear within the relevant acoustic win-dow.For glycerol, SZZ(k) exhibits no maximum at theFSDP (k ≃ 1.5 Å−1[57]) and varies smoothly through-out the acoustic window (Fig. 4). To examine whetherthe adopted minimal form ∆q(k) = q2k2 is compati-ble with the MD charge correlations, the inset of Fig. 4shows SZZ(k) on double-logarithmic axes. Because thesmallest-k values are affected by the finite simulationcell and the real-space truncation used in evaluating thepair-correlation integrals, we restrict the comparison to63     0 1 2 3 4 5 6 70.00.51.01.5Srr(k)SZZ(k)k (Å-1)kFSDPkD Glycerol (80K) SZZ(k)kD01234 Srr(k)1 2 310−410−310−210−1100 FIG. 4. Static charge and density correlations in glycerolglass (80 K). Shown are the charge–charge structure factorSZZ(k) (blue) and the mass-density static structure factorSρρ(k) (magenta) obtained from MD simulations. The FSDPand Debye wavenumbers (kFSDP, kD) are indicated. Sρρ(k)exhibits a pronounced FSDP at k ∼ 1.5 Å−1, whereas SZZ(k)has no FSDP maximum; instead its first (lowest-k) maximumappears only at higher k, on the interatomic-spacing scale.The inset shows SZZ(k) on double-logarithmic axes in thelow-k region. A power-law fit SZZ(k) = Akn over the lowestnumerically reliable interval 0.532 ≤ k ≤ 0.931 Å−1givesn ≃ 4.0, with coefficient of determination R2fit = 0.998. Thered line shows this fit, extended beyond the fitting interval asa guide to the eye. The fitted dependence is consistent withthe k4 oscillator-strength scaling implied by the minimal form∆q(k) = q2k2.the lowest numerically reliable interval, 0.532 ≤ k ≤0.931 Å−1. A free power-law fit SZZ(k) = Akn over thisinterval gives n ≃ 4.0, with coefficient of determinationR2fit = 0.998, showing that the accessible low-k behavioris compatible with a k4 dependence.Within the limited diagonal and factorized connectionestablished in Appendix K, this agreement provides astructural consistency check for the oscillator-strengthdependence ⟨|∆q(k)|2⟩dis ∝ k4 implied by ∆q(k) = q2k2.The absence of an FSDP maximum likewise providesno motivation for adding an explicit FSDP-like struc-ture to ∆q(k). Independently, the dielectric fit yieldsa value of q0 consistent with zero, as quantified in Ap-pendix J. Accordingly, for Fig. 1 we use the minimalform ∆q(k) = q2k2 and select q2 by the parameter scandescribed above. With this choice the model yields anear-linear CIR(ω) in the vicinity of the BP.We next examine how this near-linear behaviorchanges when a finite constant term q0 is introduced. Tothis end, we vary the dimensionless ratio R = q0/(q2k2D),which sets the balance between the k-independent (q0)and quadratic (q2k2) parts in ∆q(k). Increasing q0 en-hances the k-independent weight of ⟨|∆q(k)|2⟩dis—whitein k—which, upon Fourier synthesis with random phases4     0.0 0.2 0.4 0.6 0.8 1.00.00.20.40.60.81.0(a)Normalized Dq(k)k/kD 50.0 15.0 6.00 3.00 1.50 0.75 0.50 0.25 0.000.0 0.5 1.0 1.5 2.0 2.50.00.20.40.60.81.0Normalized CIR(w)Frequency (THz)(b) FIG. 5. Evolution with R = q0/(q2k2D) in ∆q(k) = q0 + q2k2:approaching the Taraskin-like regime. For each R, q0 andq2 are rescaled to keep∫ kD0dk k2|∆q(k)|2 constant, isolatingthe effect of the k-independent term q0 (see Appendix L).(a) Normalized ∆q(k) (scaled to unity at k = kD) plottedas a function of k/kD for several R. (b) Normalized CIR(ω)(scaled to unity at 2.5 THz) computed with the same G(ω).Increasing R strengthens the k-independent part of ∆q(k) anddrives CIR(ω) from a nearly linear form (small R) toward theTaraskin-like behavior CIR(ω) ≃ A+Bω2.over k < kD, corresponds to essentially white (spa-tially uncorrelated) charge fluctuations in real space. Asthe ratio R increases (Fig. 5(a)), the computed CIR(ω)evolves continuously toward a constant-term-dominatedform (Fig. 5(b)). Within our framework the mechanism isdirect: in the limit where the k-independent componentdominates the ε(ω) integral (∆q(k) ≈ q0), the absorp-tion α(ω) nearly follows the VDOS, α(ω) ∝ g(ω)/n′(ω);hence CIR(ω) becomes only weakly frequency dependent.V. RELATION TO THE TARASKIN COUPLINGFORMFinally, we turn to the question of why Taraskin’smodel does not reproduce the nearly linear CIR(ω) ob-served near the BP. Taraskin’s approach starts from afactorization ansatz α(ω) ≈ CIR(ω)g(ω) [39], motivatedby the mode-resolved linear-response expression for crys-tals derived by Maradudin [38]. Microscopically, how-ever, linear response gives α(ω) =∑p Cpδ(ω − ωp),i.e., a sum over modes with mode-dependent couplingsCp, not a product of a smooth coupling and the VDOSg(ω) =∑p δ(ω − ωp). Identifying α(ω) with CIR(ω)g(ω)amounts to a modeling assumption: it replaces themode-resolved couplings Cp by a single smooth func-tion of frequency CIR(ω). This limitation was empha-sized early on for disordered solids, where absorption in-volves mode- (or subband-) resolved couplings that canvary strongly with frequency (see Galeener and Sen [58]).7A related non-factorizable convolution form appears inlow-frequency Raman scattering (Schmid and Schirma-cher [59]). Under Taraskin’s decomposition into uncor-related and correlated charge components, this factoriza-tion yields the quadratic form CIR(ω) = A + Bω2 usedin their analysis. Within our framework, Eq. (16) canbe interpreted as an effective continuum k-space repre-sentation of Maradudin’s mode-resolved linear-responseformula [38], in which the discrete mode sum is re-placed by an integral over k and the mode oscillatorstrengths are represented by the k-dependent spectrum⟨|∆q(k)|2⟩dis (see Appendices B–E and Appendix K). Inthis effective representation, the spectral factor is theviscoelastic Green’s function of the acoustic branches,[−ω2 + C̃α(ω)k2]−1, determined by the complex mod-uli G(ω) and M(ω). A fully microscopic determinationof these effective inputs is left for future work.VI. CONCLUSIONSA continuum model based on IR-effective charge fluc-tuations ∆q(k) and a frequency-dependent shear mod-ulus G(ω) reproduces the measured ε′(ω), ε′′(ω), andthe near-linear CIR(ω) around the BP in glycerol glass.The single ratio R = q0/(q2k2D) tunes the IR couplingfrom this linear regime toward the constant/quadratic-coupling behavior widely used to describe THz absorp-tion in disordered materials, CIR(ω) = A + Bω2. Withmaterial-specific inputs, the framework is readily extend-able to other glasses and provides a practical route toassess THz dielectric losses.ACKNOWLEDGMENTSThis work was supported by JSPS KAKENHI GrantNos. 23H01139 and 23K25836 (to T.M.), 25H01519 and22K03543 (to H.M.), and 24K08045 (to Y.F.); by theAsahi Glass Foundation (to T.M.); and by support fromGIC & NGF (to T.M.).Appendix A: Microscopic origin of the electricdriving termHere we derive the electric source term in Eq. (7) fromthe microscopic rigid-ion dipole. The mechanical part ofEq. (7) is represented in the continuum description bythe elastic bulk modulus K and the effective, generallydissipative shear modulus G(ω).We first work in the time domain. The external electricfield is assumed to be spatially uniform over the lengthscale of the vibrational displacement field. When a po-larization direction is specified, we writeE(t) = E(t)ê, (A1)where ê is a time-independent unit vector. The corre-sponding frequency-domain convention is E(ω) = E(ω)ê,as used in the main text.For a microscopic configuration, the displaced posi-tion and the decomposition of the microscopic dipole intostatic and vibrational parts are given in Eqs. (1)–(3) ofthe main text. In the rigid-ion approximation the chargesqi are fixed to the moving sites. The static term d0is independent of the vibrational coordinates and there-fore does not contribute to the vibrational driving force,whereas dvib(t) is the dynamical part relevant for infraredabsorption.To pass to the continuum model, introduce the mi-croscopic rigid-ion charge distribution in the referenceconfiguration,qri(r) =∑iqi δ(r−R0i ), (A2)and identify the discrete displacement with the contin-uum displacement at the corresponding material point,ui(t) ≃ u(R0i , t). (A3)Then the vector vibrational dipole can be written asdvib(t) ≃∫d3r qri(r)u(r, t). (A4)In the continuum-field description used in the maintext, the reference-configuration site distribution is rep-resented by the local IR-effective charge fluctuation field∆q(r),qri(r) → ∆q(r). (A5)Applying this replacement to Eq. (A4) gives thecontinuum-field vibrational dipole in Eq. (4) of the maintext. The position vector r labels a material point inthe equilibrium reference configuration, and the spatialdependence of the IR-effective charge fluctuations is re-tained in ∆q(r). The time dependence of the vibrationaldipole is carried by u(r, t), which represents the site dis-placements in Eq. (A3). If one writes the instantaneousmicroscopic rigid-ion charge density in real space, with xdenoting a position vector in the instantaneous configu-ration,ρri(x, t) =∑iqi δ[x−R0i − ui(t)], (A6)that density describes the positions of the movingcharged sites. This should be distinguished from thereference-configuration field ∆q(r) in Eq. (4), which mul-tiplies the displacement field in the linear vibrationaldipole. For a neutral microscopic cell, the zero-wave-number component of the bare rigid-ion charge distri-bution qri(r) is constrained by charge neutrality. Thecharge-coupling quantity entering the dielectric responsebelow is the effective disorder variance ⟨|∆q(k)|2⟩dis.The interaction energy and the corresponding electricforce density are given in Eqs. (5) and (6) of the maintext. After Fourier transformation in time, the lattersupplies the electric driving term ∆q(r)E(ω) in Eq. (7).8Appendix B: Continuum equation and modedecompositionWe start from the elastodynamic equation in Eq. (7)and use the Fourier representation and the associated re-ality condition, acoustic polarization basis, and branchprojection given in Eq. (8) and the accompanying text.Here uα(k, ω) is a complex response amplitude at a fixedfrequency; the amplitude at −k is obtained by applyingthe same mode equation to the wave vector −k. Withthese definitions, Fourier transforming Eq. (7) and pro-jecting onto branch α gives the mode equation, Eq. (9)of the main text. Solving Eq. (9) for uα(k, ω) givesuα(k, ω) =∆q(k)−ρω2 + Cα(ω)k2Eα(ω). (B1)Appendix C: Macroscopic polarization and theappearance of |∆q(k)|2Starting from the field-parallel macroscopic polariza-tion defined in Eq. (10) and inserting the Fourier repre-sentation in Eq. (8) yieldsP (ω) =1V∫d3r[1√V∑k∆q(k)eik·r]× 1√V∑k′,αuα(k′, ω)(ê · êα(k′))eik′·r=1V∑k,k′,α∆q(k)uα(k′, ω)(ê · êα(k′))×∫d3r ei(k+k′)·r.(C1)The spatial integral produces a Kronecker delta (for adiscrete k-mesh),∫d3r ei(k+k′)·r = V δk′,−k, (C2)so thatP (ω) =1V∑k,α∆q(k)uα(−k, ω)(ê · êα(−k)). (C3)The amplitude at −k follows by applying the mode equa-tion, Eq. (9), to the wave vector −k:uα(−k, ω) =∆q(−k)−ρω2 + Cα(ω)k2[êα(−k) ·E(ω)]. (C4)Substituting Eq. (C4) into Eq. (C3) givesP (ω) =E(ω)V∑k,α∆q(k)∆q(−k)−ρω2 + Cα(ω)k2(ê · êα(−k))2=E(ω)V∑k,α|∆q(k)|2−ρω2 + Cα(ω)k2(ê · êα(k))2.(C5)In the second line we used the reality of the charge-coupling field in the reference configuration, ∆q(−k) =∆q(k)∗. The squared projection is unchanged by thepossible sign change of êα between k and −k. Thus thefactor |∆q(k)|2 arises from the convolution of the charge-coupling field at k with the displacement response at −k,together with the reality condition for ∆q(r). At thisstage we have not yet made any statistical assumptionsabout ∆q(r) or about the orientation of k and êα.Appendix D: Disorder average and orientationalcoupling factorsWe use the statistically homogeneous and isotropic dis-order correlation defined in Eq. (12) of the main text.Averaging Eq. (C5) over the disorder using that relationgivesP (ω) =E(ω)V∑k,α⟨|∆q(k)|2⟩dis−ρω2 + Cα(ω)k2(ê · êα)2. (D1)We now perform an orientational average over the rela-tive angle between the external field direction ê and thepolarization basis {êα(k)}. We assume that the statisticsof ∆q(k) are independent of these orientations; hence thedisorder and orientational averages factorize. We definethe orientational coupling factorγα ≡ ⟨(ê · êα)2⟩angles, (D2)and replace (ê · êα)2 in Eq. (D1) by its angular aver-age, giving Eq. (13) of the main text. Using ε(ω) =ε∞ + P (ω)/[ε0E(ω)] with that result yields the generaldiscrete-k expression for the dielectric function:ε(ω) = ε∞ +1ρε0V∑k,α⟨|∆q(k)|2⟩dis γα−ω2 + C̃α(ω)k2, (D3)where C̃α(ω) = Cα(ω)/ρ. Next we apply the isotropicDebye-sphere replacement in Eq. (14) of the main textto Eq. (D3). This gives the general branch-resolved di-electric function in Eq. (15) of the main text. The re-placement follows the standard isotropic Debye prescrip-tion, whereby the discrete k sum is approximated bya continuum integral over the isotropic Debye sphere0 < k < kD, with kD providing an interatomic-scale(short-wavelength) cutoff for the present continuum for-mulation.Appendix E: Isotropic limit: γ = 1/3 and the finalexpressionIn an isotropic three-dimensional elastic medium thethree polarization vectors {êL, êT1, êT2} form an or-thonormal basis for any given k. For a fixed ê we alwayshave(ê · êL)2 + (ê · êT1)2 + (ê · êT2)2 = |ê|2 = 1. (E1)9Because the three directions are statistically equivalent,their orientational averages are identical:γL = γT1= γT2≡ γ. (E2)Taking the orientational average of Eq. (E1) immediatelygives3γ = 1 ⇒ γ =13. (E3)Substituting Eq. (E3) into Eq. (15) and grouping the twodegenerate transverse branches (CT1= CT2= G) and thesingle longitudinal branch (CL = M) gives Eq. (16) ofthe main text, with G̃(ω) and M̃(ω) defined immediatelybelow that equation.Appendix F: Branch-resolved polarizations PT (ω)and PL(ω) and the dielectric contributionsWe define the branch-resolved contributions of thetransverse and longitudinal acoustic branches to the field-parallel polarization PT (ω) and PL(ω), as follows:PT (ω) =E(ω)V∑k|∆q(k)|2−ρω2 +G(ω)k2(ê · êT )2, (F1)PL(ω) =E(ω)V∑k|∆q(k)|2−ρω2 +M(ω)k2(ê · êL)2, (F2)where in Eq. (F1) êT represents one of the two trans-verse polarizations. In the isotropic limit each acousticbranch carries weight 1/3; combining the two degeneratetransverse branches then givesP (ω) =23PT (ω) +13PL(ω), (F3)and Eq. (16) follows from ε(ω) = ε∞ + P (ω)/[ε0E(ω)]after applying the Debye-sphere replacement in Eq. (14).Equation (16) can therefore be viewed as a sum of trans-verse and longitudinal propagators. Using this branch-resolved decomposition, we define the corresponding con-tributions to the modeled dielectric function, εT (ω) andεL(ω), which are shown in Fig. 6. Accordingly, we writeε(ω) = ε∞ + εT (ω) + εL(ω), (F4)withεT (ω) =1ρε0∫ kD0[23⟨|∆q(k)|2⟩dis−ω2 + G̃(ω)k2]3k2k3Ddk, (F5)εL(ω) =1ρε0∫ kD0[13⟨|∆q(k)|2⟩dis−ω2 + M̃(ω)k2]3k2k3Ddk. (F6)Here G̃(ω) and M̃(ω) are the mass-density-normalizedmoduli defined immediately below Eq. (16). Figure 66    0.00.20.40.60.81.01.2e'(a)0.0 0.5 1.0 1.5 2.0 2.50.00.10.20.30.40.50.6e"Frequency (THz)(b) Total Transverse Longitudinal FIG. 6. Transverse and longitudinal contributions to the mod-eled complex dielectric function. (a) Real part and (b) imag-inary part. Red: transverse contribution εT (ω); blue: lon-gitudinal contribution εL(ω); black: total εT (ω) + εL(ω) =ε(ω) − ε∞. The functions εT (ω) and εL(ω) are defined inEqs. (F5) and (F6).summarizes the transverse and longitudinal dielectriccontributions defined in Eqs. (F5) and (F6) by plot-ting their real and imaginary parts together with theirsum, εT (ω) + εL(ω) = ε(ω) − ε∞. Near the boson peak(BP) frequency ωBP/(2π) ≃ 1.0 THz, both ε′(ω) andε′′(ω) are dominated by the transverse channel. Thissuppression of the longitudinal contribution follows di-rectly from the structure of the longitudinal modulus,M(ω) = K + (4/3)G(ω). In the present model, K is areal, frequency-independent constant; thus all frequencydependence and dissipation ofM(ω) originate fromG(ω),while the constant K raises the longitudinal stiffness andsuppresses the longitudinal dielectric response.Appendix G: Effective shear modulus G(ω) used forglycerolOur dielectric model requires the complex shear mod-ulus G(ω) as an input function. Importantly, the present10formulation does not rely on a specific microscopic originof the BP: any physically reasonable G(ω) that repro-duces the BP part of the reduced vibrational density ofstates (VDOS) is sufficient for the calculations discussedin the main text.For glycerol, the effective shear modulus G(ω) usedthroughout the dielectric calculations was obtained byapplying the heterogeneous-elasticity-theory coherent-potential-approximation (HET-CPA) procedure de-scribed in Ref. [50] to the experimental INS VDOS ofglycerol glass at 80 K [52]. The fixed mechanical in-puts used to set the elastic and wavenumber scales wereρ = 1.26 g cm−3, vT = 1871 m s−1, vL = 3614 m s−1,and kD = 1.897 Å−1. This numerical input set is listedfor glycerol in Ref. [50], and the sound velocities orig-inate from Brillouin light-scattering measurements [53].The corresponding bulk modulus was K = 10.58 GPa,calculated as K = ρ(v2L − 4v2T /3). The HET-CPA fitto the BP part of the 80 K INS VDOS yielded ke/kD =0.399 (ke = 0.757 Å−1), σ = 0.766 (σ2 = 0.587), andG0 = 5.43 GPa. These parameters define the VDOS-constrained mechanical input G(ω) entering Eq. (16).Figure 2(a) shows the resulting real and imagi-nary parts of G(ω), plotted in the normalized formRe[G(ω)/G0] and Im[G(ω)/G0]. Using the same Green-function-based relation between G(ω) and the VDOS asin Ref. [50], we compute the reduced VDOS g(ω)/ω2shown in Fig. 2(b). The agreement with the INS spec-trum demonstrates that the effective modulus used inthe main-text dielectric calculations is consistent withthe experimentally observed BP spectrum. Once the log-normal HET-CPAG(ω) used in the main analysis was ob-tained from the VDOS analysis, it and K were kept fixedduring the subsequent (q0, q2) charge-parameter scansand the calculation of the dielectric response reportedin the main text; no HET-CPA parameter was varied inthose steps.Appendix H: Molecular dynamics glassconfiguration (glycerol)The atomistic glycerol glass structure shown in Fig. 7was generated by classical molecular–dynamics (MD)simulation using LAMMPS. The system contains N =12096 atoms, corresponding to Nmol = 864 glycerolmolecules in an all-atom representation (14 atoms permolecule), in a cubic periodic simulation cell of sidelength L = 47.16 Å (mass density ρ = 1.26 g cm−3).The interatomic interactions were described by an all-atom AMBER-based glycerol force field following thereparameterized model of Blieck et al. [42], which is basedon the AMBER glycerol model introduced by Chelli etal. [43, 44]. The specific parameter set employed hereis summarized in Ref. [45]. Long-range electrostaticswere treated using the particle–particle particle–mesh(PPPM) method in LAMMPS. Starting from an equi-librated liquid configuration at T = 350 K, the system10   FIG. S3. Atomic configuration of glycerol glass at 80 K (MD). Representative MD snapshot of the 80 K glycerol glass used for the structure-factor analysis (basis for Fig. 3 in the main text).    FIG. 7. Atomic configuration of glycerol glass at 80 K (MD).Representative MD snapshot of the 80 K glycerol glass usedfor the structure-factor analysis (basis for Fig. 4 in the maintext).was cooled to T = 80 K in the canonical (NVT) ensem-ble at fixed volume with a cooling rate of 0.4 K/ps toobtain a glassy state, subsequently equilibrated at 80 K.Figure 7 shows a representative snapshot from the equi-librated 80 K glass used for the structural analysis. Thesame 80 K configurations were used as input for the cal-culation of Sρρ(k) and SZZ(k) described in the followingsubsection.Appendix I: Static structure factors Sρρ(k) andSZZ(k)From the MD-generated glass configuration, we eval-uate the mass-density static structure factor Sρρ(k) andthe (normalized) charge–charge structure factor SZZ(k).We classify atomic sites into species (types) α, β, withnumber fractions cα ≡ Nα/N and total number densityn ≡ N/V . From the partial pair correlation functionsgαβ(r) we define the total correlation functions:hαβ(r) ≡ gαβ(r)− 1. (I1)Assuming isotropy, we use the spherically averagedFourier transform of hαβ(r). In practice the real-spaceintegral is truncated at a maximum distance rmax (typi-cally rmax = L/2 for a cubic box of side length L), andwe apply the standard Lorch modification function fL(r)11to reduce termination artifacts [46]:hαβ(k) = 4π∫ rmax0dr r2hαβ(r)sin(kr)krfL(r),fL(r) =sin(πr/rmax)πr/rmax.(I2)The mass-density static structure factor is then evaluatedin the Faber–Ziman form [47] asSρρ(k) = 1 +nm̄2∑α,βcαcβmαmβhαβ(k),m̄ ≡∑γcγmγ .(I3)where mα is the atomic mass of species α (note thatnm̄ = ρ).For the charge–charge correlations, we use the samepartial functions but weight them by the (dimensionless)partial charges zα (i.e., the atomic partial charge is ezα).We follow the conventional multicomponent definition ofthe charge–charge structure factor [55, Chap. 10] and usethe normalized form:SZZ(k) = 1 +nz̄2∑α,βcαcβzαzβhαβ(k),z̄2 ≡∑γcγz2γ .(I4)Charge neutrality implies∑α cαzα = 0; the abovenormalization corresponds to dividing the conventionalSZZ(k) by z̄2 and ensures SZZ(k) → 1 at large k.Because the pair-correlation integrals are evaluated ina finite simulation cell and truncated at rmax = L/2, thesmallest-k values are sensitive to finite-size and termina-tion effects. We therefore exclude the points below k =0.532 Å−1from the power-law analysis and use the low-est numerically reliable interval, 0.532 ≤ k ≤ 0.931 Å−1.A free log–log fit SZZ(k) = Akn over this interval givesn ≃ 4.0, with coefficient of determination R2fit = 0.998.The red line in the inset of Fig. 4 shows this fit, extendedbeyond the fitting interval as a guide to the eye.Appendix J: Parameter scan for q2In the main text, the IR-effective charge fluctuationspectrum is represented by the minimal analytic param-eterization ∆q(k) = q0 + q2k2 over the k range enteringthe dielectric integral, 0 < k < kD. For the glycerolcalculations reported here, we adopt the minimal form∆q(k) = q2k2 and determine q2 by a one-parameter scan.The role of the constant term q0 is assessed separatelyby the two-parameter check described below. For eachtrial value of q2, we compute the complex permittivityε(ω) = ε′(ω) + iε′′(ω) with all other model inputs fixed,including the HET-CPA-derived mechanical input G(ω)and the values listed in Table I. We then evaluate a scalarmisfit (sum of squared residuals)S(q2) =N∑i=1[(ε′calc(ωi; q2)− ε′exp(ωi))2+(ε′′calc(ωi; q2)− ε′′exp(ωi))2],(J1)where {ωi} are the experimental frequencies within thefitting window. We focus on the BP region (around∼ 1 THz for glycerol) where the vibrational contri-bution is prominent. The scan was performed overq2 = (1.25–1.90)× 103 Ccm−3 Å2with a step of 0.005×103 Ccm−3 Å2, and Fig. 8 shows the normalized mis-fit S(q2)/Smin. A clear minimum is found at q2 =1.595 × 103 Ccm−3 Å2which is reported in Table I asq2 = 1.6× 103 Ccm−3 Å2(rounded). Repeating the scanfor slightly shifted fitting windows within the experimen-tal overlap changes the optimum by no more than onescan step (Fig. 8(a)).As a check on the minimal choice q0 = 0, we also min-imized the same scalar misfit with respect to both q0and q2, while keeping all other model inputs fixed, in-cluding the HET-CPA-derived mechanical input G(ω).For the main fitting window 0.5–1.5 THz, this two-parameter minimization gives q0 = 31 ± 39C cm−3 andq2 = (1.578 ± 0.019) × 103 Ccm−3 Å2. The uncertain-ties are formal one-standard-error estimates for the least-squares misfit defined above. Thus q0 = 0 lies within thefitting uncertainty. When q0 is fixed to zero and q2 is re-optimized, the misfit increases by only 0.58% relative tothe two-parameter minimum. The value used in Table I,q2 = 1.6 × 103 Ccm−3 Å2, is therefore consistent withthe two-parameter minimization. We consequently re-tain q0 = 0 as the minimal parameterization for glycerolglass at 80 K in the present frequency window.Appendix K: Plane-wave expansion and amicroscopic route to the static SZZ(k)This section provides a microscopic interpretation ofthe k-dependent oscillator-strength spectrum used in thecontinuum expression for ε(ω) in the main text. Westart from the standard linear-response mode-sum rep-resentation for the vibrational dielectric loss in the har-monic approximation, expressed in terms of mode effec-tive charges. We then project disordered eigenmodesonto longitudinal and transverse plane-wave componentswith mode-dependent k-space weights. After displayingthe four-point configurational average that appears be-fore approximation, we introduce a diagonal, factorizedapproximation in which cross/interference terms are ne-glected and the remaining charge correlations are rep-resented by the static charge–charge structure factorSZZ(k). This construction gives a microscopic route and1211   FIG. S4. Panel (a) shows log10[𝑆(𝑞2) 𝑆min⁄ ] for three fitting windows within the experimental overlap, demonstrating that the minimum is robust. Panel (b) shows a zoom near the minimum for the main fitting window (0.5–1.5 THz). The dashed line marks the adopted value 𝑞2 =1.595 × 103 C cm−3 Å2.   1.3 1.4 1.5 1.6 1.7 1.8 1.9100101102103 0.5-1.5 THz 0.5-2.5 THz 0.32-2.53 THzlog10[S(q2)/Smin]q2 [103 C cm-3 Å2](a)(b)1.585 1.590 1.595 1.600 1.60512345 0.5-1.5 THzS(q2)/Sminq2 [103 C cm-3 Å2]  FIG. 8. Panel (a) shows log10[S(q2)/Smin] for three fittingwindows within the experimental overlap, demonstrating thatthe minimum is robust. Panel (b) shows a zoom near the min-imum for the main fitting window (0.5–1.5 THz). The dashedline marks the adopted value q2 = 1.595× 103 Ccm−3 Å2.diagnostic by which static charge correlations can con-tribute to the effective coupling spectrum ⟨|∆q(k)|2⟩disappearing in Eq. (16); it is not a full microscopic deter-mination of that spectrum over the entire experimentalfrequency window.(i) Dielectric loss expressed in vibrational eigenmodes.We recall a standard mode-sum expression for the di-electric loss. For an isotropic system, it can be expressedas [48, 49, 54]:ε′′(ω) =4π23V∑p∑j=x,y,z|F (p)j |22ωpδ(ω − ωp), (K1)where p labels eigenmodes of frequency ωp and F(p)j is theCartesian component of the mode charge-coupling vector(whose squared magnitude gives the oscillator strength).In the most general form (Born effective charges) one has[48, 49]F(p)j ≡∑i,kZi,jk√mie(p)i,k , (K2)with Zi,jk the Born effective charge tensor and e(p)i thenon-mass-weighted polarization vector. For the presentanalysis, we adopt the rigid-ion approximation Zi,jk ≃qiδjk, which reduces Eq. (K2) to the mode charge-coupling vectorF(p) ≡∑iqi√mie(p)i ,∑j|F (p)j |2 = |F(p)|2. (K3)Equation (K1) then becomesε′′(ω) =4π23V∑p|F(p)|22ωpδ(ω − ωp). (K4)(ii) Plane-wave expansion with explicit longitudi-nal/transverse polarizations. We project the displace-ment pattern of each eigenmode onto plane-wave compo-nents. The coefficients akα(ωp) specify the correspond-ing k-space weight distribution. Taraskin et al. explicitlystate that the low-frequency disordered modes can be ap-proximately expanded in plane waves and that the distri-bution of the plane-wave coefficients shows two relativelynarrow peaks associated with transverse/longitudinalacoustic (TA/LA) hybridization (while broadening athigher frequencies as disorder-induced attenuation be-comes stronger) [39, 56]. To make the polarization con-tent explicit, we writee(p)i ≃∑k∑α∈{L,T1,T2}√mim̄akα(ωp) êα(k) eik·ri , (K5)where m̄ ≡ N−1∑i mi. For each k, the three polariza-tion vectors form an orthonormal triad: êL(k) ∥ k andêTa(k) · k = 0 (a = 1, 2).(iii) From plane-wave expansion to a static SZZ(k).Substituting Eq. (K5) into Eq. (K3) givesF(p) ≃ 1√m̄∑k,αakα(ωp) êα(k) ρZ(k),ρZ(k) ≡∑iqieik·ri .(K6)Consider the squared projection onto the unit vector êspecifying the direction of the applied electric field. Be-fore the diagonal approximation is made, substitutingEq. (K6) and performing the orientational average gives〈|ê · F(p)|2〉conf,angles≃ 1m̄∑k,α∑k′,α′Γαα′(k,k′)× ⟨akα(ωp)a∗k′α′(ωp)×ρZ(k)ρZ(k′)∗⟩conf , (K7)13whereΓαα′(k,k′) ≡ ⟨(ê · êα(k))×(ê · êα′(k′))⟩angles=13êα(k) · êα′(k′). (K8)Equation (K7) shows explicitly that the general averagedoscillator strength contains a four-point configurationalaverage involving both the plane-wave weights and thecharge-density Fourier components. To obtain the static-SZZ(k) form used below, we make a diagonal, factorizedapproximation: cross/interference terms with (k, α) ̸=(k′, α′) are neglected, and the remaining diagonal part isfactorized as⟨akα(ωp)a∗k′α′(ωp)ρZ(k)ρZ(k′)∗⟩conf≃ δk,k′δα,α′ |akα(ωp)|2〈|ρZ(k)|2〉conf. (K9)Here |akα(ωp)|2 denotes the configuration-averaged orrepresentative k-space weight; connected correlations be-tween the mode weights and the charge correlations arenot retained in this approximation. Using Eq. (K9) anddefining the orientational factor γα = Γαα(k,k), we ob-tain〈|ê · F(p)|2〉≃ 1m̄∑k,α|akα(ωp)|2〈|ρZ(k)|2〉confγα=Nq̄ 2m̄∑k,α|akα(ωp)|2 SZZ(k) γα. (K10)where q̄ 2 ≡ N−1∑i q2i , andSZZ(k) ≡1Nq̄ 2⟨ρZ(k)ρZ(−k)⟩is the normalized static charge–charge structure factor,consistent with the definition adopted in Appendix I(note that qi = ezi implies q̄ 2 = e2z̄ 2 in the notationused there), andγα ≡〈(ê · êα)2〉angles. (K11)In an isotropic three-dimensional medium the three po-larizations are statistically equivalent, giving γL = γT1=γT2= 1/3 (hence a transverse-to-longitudinal weight ra-tio 2 : 1). If one further assumes that the two trans-verse polarizations are equivalent on average, |akT1|2 ≃|akT2|2 ≡ |akT |2, then Eq. (K10) can be grouped as〈|ê · F(p)|2〉≃ Nq̄ 2m̄∑k[23|akT (ωp)|2+13|akL(ωp)|2]SZZ(k).(K12)(iv) Back to ε′′(ω) and comparison with Eq. (16) inthe main text. Using 13∑j |F(p)j |2 = ⟨|ê · F(p)|2⟩angles,Eq. (K4) can be viewed as a mode sum over the orien-tationally averaged projection. Substituting Eq. (K12)yields the approximate formε′′(ω) ∝∑pδ(ω − ωp)ωp∑k[23|akT (ωp)|2+13|akL(ωp)|2]SZZ(k),(K13)up to an overall prefactor. Equation (K13) makes ex-plicit that (i) within the diagonal, factorized approxi-mation of Eq. (K9), a static SZZ(k) can contribute tothe effective oscillator-strength spectrum via the plane-wave expansion, and (ii) the transverse and longitudi-nal contributions appear in the 2/3 vs 1/3 ratio charac-teristic of an isotropic medium, matching the polariza-tion weighting used in the continuum expression of themain text [Eq. (16)] and in Appendices D–F. This static-SZZ(k) construction should therefore be understood asa diagonal/factorized microscopic route and diagnostic,not as a quantitative microscopic determination of theeffective coupling spectrum over the full 0.3–2.5 THzexperimental window. This approximation is expectedto be most transparent in the low-frequency propagonregime, where the k-space weights have relatively narrowacoustic peaks. Above the BP, the weights broaden asthe modes become diffuson-like, and the neglected off-diagonal terms and possible connected correlations be-tween mode weights and charge correlations may becomemore important. A full microscopic determination ofthe effective charge-coupling spectrum would require adynamic, mode-resolved charge-correlation description,such as SZZ(k, ω) or an equivalent eigenmode-resolvedeffective-charge analysis.Finally, the physical meaning of the plane-waveweights is analogous to a k-space spectral function:Taraskin et al. report that at low frequencies thedistribution of akα(ω) is dominated by relatively nar-row TA/LA peaks (hybridization at fixed ω), while athigher frequencies the distribution broadens in k spaceas disorder-induced attenuation becomes strong [39, 56].In our continuum formalism, the corresponding broad-ening is captured in the imaginary part of the effectivepropagators (−ω2 + C̃α(ω)k2)−1 appearing in Eq. (16).Appendix L: Rescaling procedure for the R-scan inFig. 5We parametrize the k-dependent effective charge fluc-tuation as∆q(k) = q0 + q2k2, (L1)and define the dimensionless ratioR ≡ q0q2k2D, (L2)14which measures the relative weight of the k-independentcomponent at the Debye wavenumber kD. When scan-ning R, it is useful to isolate the effect of the shape of∆q(k) without allowing the overall fluctuation strengthto grow artificially at large R. To this end, we rescale(q0, q2) for each R so that the Debye-sphere integratedstrength is kept constant:∫ kD0dk k2 |∆q(k)|2 = const. (L3)Using ∆q(k) = q0 + q2k2, the integral can be evaluatedanalytically as∫ kD0dk k2 (q0+ q2k2)2 = q20k3D3+2q0q2k5D5+ q22k7D7. (L4)Substituting q0 = Rq2k2D, we obtain∫ kD0dk k2 |∆q(k)|2 = q22k7D J(R),J(R) =R23+2R5+17.(L5)We choose the constant to be equal to the R = 0 value,for which q0 = 0 and q2 = q(0)2 . In the calculationsfor Fig. 5 we used kD = 1.90 Å−1and q(0)2 = 1.595 ×103 Ccm−3 Å2(rounded to 1.6 × 103 in Table I). Thisyields the rescaling formulasq2(R) = q(0)2√(1/7)J(R), q0(R) = Rq2(R) k2D. (L6)This procedure keeps the overall fluctuation strengthfixed while varying the relative weight of the k-independent component, thereby avoiding an artificial in-crease in the dielectric strength at large R and enablinga cleaner assessment of the R-dependence.Appendix M: Raman coupling coefficient of glycerolLow-frequency Raman spectra of glycerol glass weremeasured in backscattering geometry with the con-focal micro-Raman setup described in Refs. [11, 12].A single-longitudinal-mode, frequency-doubled diode-pumped solid-state Nd:YAG laser at 532 nm was usedfor excitation; the scattered light was collected through ahome-built microscope with ultranarrow-band notch fil-ters and analyzed by an HR320 single monochromator(1200 grooves/mm) equipped with a CCD camera (An-dor DU420). Both the parallel HH and crossed HV po-larization geometries were recorded, and the depolarizedHV spectrum was used in the present analysis.The Raman spectra were measured at 80 K. Follow-ing Refs. [11, 12], the measured Raman intensity in theHV geometry, IHV(ω), and the corresponding imaginary7  0.0000.0050.0100.0150.0200.025g(n)/n2 (THz-3)(a) INSGlycerol (80 K)01020c'' HV(n)/n (arb. units) Raman0.0 0.5 1.0 1.5 2.00.00.51.0CRaman(n) (arb. units)Frequency (THz)(b) Glycerol (80 K) FIG. 9. Raman coupling coefficient of glycerol. (a) ReducedVDOS g(ω)/ω2 from the 80 K INS data [52] (black symbols)and the reduced depolarized Raman response χ′′HV(ω)/ω ob-tained from the HV spectrum measured at 80 K (red sym-bols). (b) Raman coupling coefficient CRaman(ω) obtained bydividing the reduced Raman response by the reduced VDOS.The Raman coupling coefficient increases approximately lin-early in the BP region.Raman susceptibility χ′′HV(ω) are related by IHV(ω) ∝[nB(ω, T )+1]χ′′HV(ω), where nB(ω, T ) = [exp(ℏω/kBT )−1]−1 is the Bose–Einstein occupation number, ℏ is thereduced Planck constant, kB is the Boltzmann constant,and T denotes temperature. The reduced Raman re-sponse shown in Fig. 9(a) was therefore evaluated asχ′′HV(ω)/ω ∝ IHV(ω)/{ω[nB(ω, T ) + 1]}.The Raman light–vibration coupling coefficientCRaman(ω) was obtained by combining this reducedRaman response with the 80 K INS VDOS [52], usingχ′′HV(ω)/ω = CRaman(ω)g(ω)/ω2. Thus, over the com-mon frequency range of the two data sets, CRaman(ω)is proportional to [χ′′HV(ω)/ω]/[g(ω)/ω2]. As shownin Fig. 9(b), the Raman coupling coefficient increasesapproximately linearly in the BP region.15[1] T. Nakayama, Boson peak and terahertz frequency dy-namics of vitreous silica, Rep. Prog. Phys. 65, 1195(2002).[2] M. A. Ramos, Low-Temperature Thermal and VibrationalProperties of Disordered Solids (World Scientific, Lon-don, 2022).[3] R. C. Zeller and R. O. Pohl, Thermal conductivity andspecific heat of noncrystalline solids, Phys. Rev. B 4, 2029(1971).[4] U. Buchenau, M. Prager, N. Nücker, A. J. Dianoux,N. Ahmad, and W. A. Phillips, Low-frequency modes invitreous silica, Phys. Rev. B 34, 5665 (1986).[5] B. Rufflé, D. A. Parshin, E. Courtens, and R. Vacher,Boson peak and its relation to acoustic attenuation inglasses, Phys. Rev. Lett. 100, 015501 (2008).[6] G. Baldi, V. M. Giordano, G. Monaco, and B. Ruta,Sound attenuation at terahertz frequencies and the bo-son peak of vitreous silica, Phys. Rev. Lett. 104, 195501(2010).[7] A. I. Chumakov, G. Monaco, A. Monaco, W. A. Crich-ton, A. Bosak, R. Rüffer, A. Meyer, F. Kargl, L. Comez,D. Fioretto, H. Giefers, S. Roitsch, G. Wortmann,M. H. Manghnani, A. Hushur, Q. Williams, J. Balogh,K. Parliński, P. Jochym, and P. Piekarz, Equivalence ofthe boson peak in glasses to the transverse acoustic vanHove singularity in crystals, Phys. Rev. Lett. 106, 225501(2011).[8] W. A. Phillips, Amorphous Solids: Low-TemperatureProperties (Springer, Berlin, 1981).[9] V. K. Malinovsky and A. P. Sokolov, The nature of bosonpeak in Raman scattering in glasses, Solid State Com-mun. 57, 757 (1986).[10] S. Kojima, Low-frequency Raman investigation of theliquid-glass transition in glycerol, Phys. Rev. B 47, 2924(1993).[11] M. Kabeya, T. Mori, Y. Fujii, A. Koreeda, B. W.Lee, J.-H. Ko, and S. Kojima, Boson peak dynamics ofglassy glucose studied by integrated terahertz-band spec-troscopy, Phys. Rev. B 94, 224204 (2016).[12] T. Mori, Y. Jiang, Y. Fujii, S. Kitani, H. Mizuno, A. Ko-reeda, L. Motoji, H. Tokoro, K. Shiraki, Y. Yamamoto,and S. Kojima, Detection of boson peak and fractaldynamics of disordered systems using terahertz spec-troscopy, Phys. Rev. E 102, 022502 (2020).[13] H. Shintani and H. Tanaka, Universal link between theboson peak and transverse phonons in glass, Nat. Mater.7, 870 (2008).[14] A. Tanguy, B. Mantisi, and M. Tsamados, Vibrationalmodes as a predictor for plasticity in a model glass, Eu-rophys. Lett. 90, 16004 (2010).[15] S. Franz, G. Parisi, P. Urbani, and F. Zamponi, Universalspectrum of normal modes in low-temperature glasses,Proc. Natl. Acad. Sci. U.S.A. 112, 14539 (2015).[16] E. Lerner, G. Düring, and E. Bouchbinder, Statistics andproperties of low-frequency vibrational modes in struc-tural glasses, Phys. Rev. Lett. 117, 035501 (2016).[17] H. Mizuno, H. Shiba, and A. Ikeda, Continuum limitof the vibrational properties of amorphous solids, Proc.Natl. Acad. Sci. U.S.A. 114, E9767 (2017).[18] Y.-C. Hu and H. Tanaka, Origin of the boson peak inamorphous solids, Nat. Phys. 18, 669 (2022).[19] E. Lerner and E. Bouchbinder, Boson-peak vibrationalmodes in glasses feature hybridized phononic and quasilo-calized excitations, J. Chem. Phys. 158, 194503 (2023).[20] W. Schirmacher, Thermal conductivity of glassy materi-als and the “boson peak”, Europhys. Lett. 73, 892 (2006).[21] W. Schirmacher, M. Paoluzzi, F. C. Mocanu,D. Khomenko, G. Szamel, F. Zamponi, and G. Ruocco,The nature of non-phononic excitations in disorderedsystems, Nat. Commun. 15, 3107 (2024).[22] M. Wyart, Scaling of phononic transport with connec-tivity in amorphous solids, Europhys. Lett. 89, 64001(2010).[23] E. DeGiuli, A. Laversanne-Finot, G. Düring, E. Lerner,and M. Wyart, Effects of coordination and pressure onsound attenuation, boson peak and elasticity in amor-phous solids, Soft Matter 10, 5628 (2014).[24] U. Buchenau, Y. M. Galperin, V. L. Gurevich, and H. R.Schober, Anharmonic potentials and vibrational localiza-tion in glasses, Phys. Rev. B 43, 5039 (1991).[25] V. L. Gurevich, D. A. Parshin, and H. R. Schober, An-harmonicity, vibrational instability, and the boson peakin glasses, Phys. Rev. B 67, 094203 (2003).[26] C. Rainone, P. Urbani, F. Zamponi, E. Lerner, andE. Bouchbinder, Mean-field model of interacting quasilo-calized excitations in glasses, SciPost Phys. Core 4, 008(2021).[27] A. Moriel, E. Lerner, and E. Bouchbinder, Boson peakin the vibrational spectra of glasses, Phys. Rev. Research6, 023053 (2024).[28] R. T. Ako, A. Upadhyay, W. Withayachumnankul,M. Bhaskaran, and S. Sriram, Dielectrics for terahertzmetasurfaces: Material selection and fabrication tech-niques, Adv. Opt. Mater. 8, 1900750 (2020).[29] S. Havriliak and S. Negami, A complex plane representa-tion of dielectric and mechanical relaxation processes insome polymers, Polymer 8, 161 (1967).[30] L. Casella, M. Baggioli, T. Mori, and A. Zaccone, Physicsof phonon-polaritons in amorphous materials, J. Chem.Phys. 154, 014501 (2021).[31] E. A. DiMarzio and M. Bishop, Connection between themacroscopic electric and mechanical susceptibilities, J.Chem. Phys. 60, 3802 (1974).[32] K. Niss, B. Jakobsen, and N. B. Olsen, Dielectric andshear mechanical relaxations in glass-forming liquids: atest of the Gemant-DiMarzio-Bishop model, J. Chem.Phys. 123, 234510 (2005).[33] J. Kölbel, W. Schirmacher, E. Shalaev, and J. A. Zeitler,Terahertz dynamics in the glycerol-water system, Phys.Rev. B 107, 104203 (2023).[34] B. Hehlen, E. Courtens, R. Vacher, A. Yamanaka,M. Kataoka, and K. Inoue, Hyper-Raman scattering ob-servation of the boson peak in vitreous silica, Phys. Rev.Lett. 84, 5355 (2000).[35] B. Hehlen and G. Simon, The vibrations of vitreous sil-ica observed in hyper-Raman scattering, J. Raman Spec-trosc. 43, 1941 (2012).[36] G. Simon, B. Hehlen, E. Courtens, E. Longueteau, andR. Vacher, Hyper-Raman scattering from vitreous boronoxide: Coherent enhancement of the boson peak, Phys.Rev. Lett. 96, 105502 (2006).[37] N. V. Surovtsev and A. P. Sokolov, Frequency behavior16of Raman coupling coefficient in glasses, Phys. Rev. B66, 054205 (2002).[38] A. A. Maradudin and R. F. Wallis, Lattice anharmonic-ity and optical absorption in polar crystals. II. classicaltreatment in the linear approximation, Phys. Rev. 123,777 (1961).[39] S. N. Taraskin, S. I. Simdyankin, S. R. Elliott, J. R.Neilson, and T. Lo, Universal features of terahertz ab-sorption in disordered materials, Phys. Rev. Lett. 97,055504 (2006).[40] Y. Hashimoto, T. Mori, and S. Kojima, Broadbandterahertz time-domain spectroscopy and low-frequencyRaman scattering of glassy polymers: Boson peak ofPMMA, Mol. Cryst. Liq. Cryst. 629, 258 (2016).[41] S. N. Taraskin, S. I. Simdyankin, and S. R. Elliott,The atomic charge distribution in glasses obtained byterahertz spectroscopy, J. Phys.: Condens. Matter 19,455216 (2007).[42] J. Blieck, F. Affouard, P. Bordat, A. Lerbret, andM. Descamps, Molecular dynamics simulations of glyc-erol glass-forming liquid, Chem. Phys. 317, 253 (2005).[43] R. Chelli, P. Procacci, G. Cardini, R. G. Della Valle, andS. Califano, Glycerol condensed phases part I. A molec-ular dynamics study, Phys. Chem. Chem. Phys. 1, 871(1999).[44] R. Chelli, P. Procacci, G. Cardini, and S. Califano,Glycerol condensed phases part II. A molecular dynam-ics study of the conformational structure and hydrogenbonding, Phys. Chem. Chem. Phys. 1, 879 (1999).[45] A. V. Egorov, A. P. Lyubartsev, and A. Laaksonen,Molecular dynamics simulation study of glycerol-waterliquid mixtures, J. Phys. Chem. B 115, 14572 (2011).[46] E. Lorch, Neutron diffraction by germania, silica andradiation-damaged silica glasses, J. Phys. C: Solid StatePhys. 2, 229 (1969).[47] T. E. Faber and J. M. Ziman, A theory of the electricalproperties of liquid metals. III. The resistivity of binaryalloys, Philos. Mag. 11, 153 (1965).[48] A. Pasquarello and R. Car, Dynamical charge tensors andinfrared spectrum of amorphous SiO2, Phys. Rev. Lett.79, 1766 (1997).[49] L. Pedesseau, S. Ispas, and W. Kob, First-principlesstudy of a sodium borosilicate glass-former. II. the glassstate, Phys. Rev. B 91, 134202 (2015).[50] D. Kyotani, S. H. Oh, S. Kitani, Y. Fujii, H. Hi-jiya, H. Mizuno, S. Kohara, A. Koreeda, A. Masuno,H. Kawaji, S. Kojima, Y. Yamamoto, and T. Mori, Rela-tionship between the boson peak and first sharp diffrac-tion peak in glasses, Sci. Rep. 15, 9617 (2025).[51] Y. E. Ryabov, Y. Hayashi, A. Gutina, and Y. Feldman,Features of supercooled glycerol dynamics, Phys. Rev. B67, 132202 (2003).[52] J. Wuttke, W. Petry, G. Coddens, and F. Fujara, Fastdynamics of glass-forming glycerol, Phys. Rev. E 52, 4026(1995).[53] F. Scarponi, L. Comez, D. Fioretto, and L. Palmieri,Brillouin light scattering from transverse and longitudi-nal acoustic waves in glycerol, Phys. Rev. B 70, 054203(2004).[54] M. F. Thorpe and S. W. de Leeuw, Coulomb effects indisordered solids, Phys. Rev. B 33, 8490 (1986).[55] J.-P. Hansen and I. R. McDonald, Theory of Simple Liq-uids, 3rd ed. (Academic Press, Amsterdam, 2006).[56] S. N. Taraskin and S. R. Elliott, Propagation of plane-wave vibrational excitations in disordered systems, Phys.Rev. B 61, 12017 (2000).[57] D. C. Champeney, R. N. Joarder, and J. C. Dore, Struc-tural studies of liquid D-glycerol by neutron diffraction,Mol. Phys. 58, 337 (1986).[58] F. L. Galeener and P. N. Sen, Theory for the first-ordervibrational spectra of disordered solids, Phys. Rev. B 17,1928 (1978).[59] B. Schmid and W. Schirmacher, Raman scattering andthe low-frequency vibrational spectrum of glasses, Phys.Rev. Lett. 100, 137402 (2008).