###### Abstract

Glasses show vibrational properties that are markedly different to those of
crystals which are known as phonons. For example, excess low-frequency modes
(the so-called boson peak), vibrational localization, and strong scattering of
phonons have been the most discussed topics, and a theoretical understanding
of these phenomena is challenging. To address this problem, computational
simulations are a powerful tool, which have been employed by many previous
works. In this chapter, we describe simulation methods for studying the
vibrational properties of glasses (and any solid-state materials). We first
present a method for studying vibrational eigenmodes. Since vibrational
motions of particles are excited along eigenmodes, the eigenmodes are
fundamental to descriptions of vibrational properties. The eigenmodes in
glasses are non-phonon modes in general, and some of them are even localized
in space. We next present a method of analysing phonon transport, which is
also crucial for understanding vibrational properties. Since phonons are not
eigenmodes in glasses, they are decomposed into several different, non-phonon
eigenmodes. As a result, phonons in glasses are strongly scattered. In
addition, we describe how to analyse the elastic response. The elastic
response of glasses is also anomalous with respect to that of crystals.
Finally, we briefly introduce recent advances that have been achieved by means
of large-scale computational simulations.

##  Chapter 0 Computational simulations of
the vibrational properties of glasses

###  1 Introduction

Currently, there is no doubt that computational simulations play an important
role in the development of fundamental science as well as engineering
applications. Molecular Dynamics (MD) simulations and Monte Carlo (MC)
simulations have been established to simulate the behaviours of materials in
dense gas, liquid, and solid states at the microscopic, molecular level [1,
2]. The importance of these molecular simulations comes from the fact that
they provide exact, quasi-experimental data on well-defined simulation models
of materials. The usefulness of simulations is also based on the fact that
they can access data that cannot be obtained through experiments. From the
theoretical point of view, exact data on prototypical models are valuable for
understanding the fundamental mechanisms of phenomena as well as for testing
the validity of proposed theories. In the past, molecular simulations have
been employed to solve many important problems. For example, in statistical
physics, phase transitions, such as the gas-liquid transition and the
paramagnetic-ferromagnetic transition, are among the topics that are most
widely studied by means of computer simulations [3]. Additionally, to address
the problem of glass transition, many previous works have relied on
simulations of, e.g., the drastically slowed dynamics near the glass
transition [4, 5].

Computational simulations are also a powerful tool for studying the
vibrational properties of glasses. For crystals, thanks to their periodicity
and symmetry, analytical formulations can be obtained for the vibrational
motions of the molecules, which give the concept of phonons (lattice waves)
[6, 7] 111 Phonons are quantized lattice waves, but here we use the term
“phonons” for general lattice waves. We can also consider phonons in glasses,
which are sinusoidal waves propagating in disordered structures. Such phonons
in glasses will be discussed in Section 4. . Particularly, the Debye theory
has been established to explain the behaviour of the vibrational density of
states (vDOS) in a crystal. In contrast, for glasses, due to the lack of
periodicity and symmetry, analytical calculations are much difficult to
perform. Although some mean-field theories, such as replica theory [8] and
effective medium theory [9, 10, 11, 12], have been proposed, it is crucial to
test their validity and extend them to the 3-dimensional case. Also, in many
past works, e.g., Refs. 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,
experiments have been performed to study the thermal and vibrational
properties of glasses. Although the vDOS and dynamic structure factor can be
measured in experiments by monitoring the responses of various probes, such as
light, X-rays, and neutrons, it is generally difficult to directly observe the
vibrational motions of molecules. Considering this situation of separate
purely theoretical and purely experimental works, computational simulations
take advantage of the ability to directly observe and understand the
vibrational motions of particles in well-defined simulation models. There are
many relevant problems regarding the vibrational properties of glasses, e.g.,
excess low-frequency modes (the so-called boson peak), vibrational
localization, and the strong scattering of phonons. A theoretical
understanding of these phenomena remains to be developed, and doing so will be
challenging. To address these problems, computer simulations have been
employed in many previous works, e.g., Refs. 26, 27, 28, 29, 30, 31, 32, 33,
34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52,
53, 54, 55, 56, 57, 58, 59, 60, 61, 62.

In this chapter, we review simulation methods for analysing the vibrational
properties of glasses; these methods can also be applied for any solid-state
materials. In Section 2, we describe simulation models of glasses. We first
present a brief description of MD simulations, and we next explain how to
obtain glass systems through simulations. Section 3 introduces a method of
analysing vibrational eigenmodes. Since the vibrational motions of particles
are excited along eigenmodes, the eigenmodes are fundamental to descriptions
of vibrational properties. Section 4 provides a method of analysing phonon
transport, which is also a fundamental property. In addition, we also present
a method of analysing elastic moduli in Section 5. Finally, Section 6 briefly
introduces recent advances that have been achieved by means of large-scale
computational simulations.

###  2 Simulation models

Several excellent books (e.g., Refs. 1, 2) have already been published that
describe molecular simulations, i.e., MD and MC simulations. Here, we provide
only a brief description of MD simulations and explain how we can obtain glass
systems through MD simulations (we can also use MC simulations to obtain
glasses in the same way). In any solid-state material, including glasses and
crystals, the constituent particles vibrate around the so-called inherent
structure 222 The inherent structure of a crystal is a periodic lattice
structure. . We also explain how to obtain this inherent structure through
computational simulations.

####  1 Molecular Dynamics (MD) simulations

In an MD simulation, we define a system composed of an enormous number of
particles and numerically solve the equations of motion to evolve the dynamics
of those constituent particles. Here, let us consider a 3-dimensional system
composed of $N$ particles in a cubic box of length $L$ and volume $V\equiv
L^{3}$. The mass of particle $i$ ($i=1,2,\cdots,N$) is denoted by $m_{i}$, and
the position of particle $i$ is $\mathbf{r}_{i}=[r_{ix},r_{iy},r_{iz}]^{T}$ (a
3-dimensional vector), where $T$ denotes transposition. To represent the
configuration of all the constituent particles, we introduce a
$3N$-dimensional vector
$\mathbf{r}\equiv\left[\mathbf{r}_{1}^{T},\mathbf{r}_{2}^{T},\cdots,\mathbf{r}_{N}^{T}\right]^{T}$.
In this chapter, we treat vectorial quantities, such as $\mathbf{r}_{i}$ and
$\mathbf{r}$, as vertical vectors. We suppose that the total potential energy
of the system, $\Phi$, is a function of the positions of the particles:
$\Phi(\mathbf{r})\equiv\Phi(\mathbf{r}_{1},\mathbf{r}_{2},\cdots,\mathbf{r}_{N})$.
For the case of a pair-wise potential, where particles $i$ and $j$ interact
through the potential $\phi_{ij}$, $\Phi(\mathbf{r})=\sum_{i<j}\phi_{ij}$. The
equation of motion is then

$m_{i}\frac{d^{2}\mathbf{r}_{i}}{dt^{2}}=-\frac{\partial\Phi}{\partial\mathbf{r}_{i}}\quad(i=1,2,\cdots,N),$
(1)

or, in terms of the $3N$-dimensional vector $\mathbf{r}$,

$\mathcal{M}\frac{d^{2}\mathbf{r}}{dt^{2}}=-\frac{\partial\Phi}{\partial\mathbf{r}},$
(2)

where $t$ is time and we introduce the mass matrix $\mathcal{M}$, which is a
$3N\times 3N$ diagonal matrix:

$\mathcal{M}\equiv\text{diag}(m_{1},m_{1},m_{1},m_{2},m_{2},m_{2},\cdots,m_{N},m_{N},m_{N}).$
(3)

We start with some initial configuration of particles, $\mathbf{r}(t=0)$, and
velocity, $d\mathbf{r}(t=0)/dt$, and solve Eq. (1) (or Eq. (2)). Here, we need
to implement appropriate boundary conditions in the space. In this chapter, we
consider periodic boundary conditions in all directions. Thermodynamic
quantities such as temperature $T$ and pressure $p$ can be calculated from the
trajectory $\mathbf{r}(t)$: for example, the temperature is calculated from
the total kinetic energy of the particles. We may use a heat bath and/or a
pressure bath to control the temperature and/or the pressure, respectively. As
the simulation runs, the system approaches a steady, equilibrium state.

Several potentials $\Phi(\mathbf{r})$ have been proposed for modelling
different types of glass systems. Here, we briefly introduce the potentials
for modelling atomic glasses (packed glasses) and covalently bonded glasses
(network glasses).

  1. (I)

Atomic glasses (packed glasses):
To model an atomic glass composed of rare gas atoms such as argons, we can
employ the Lennard-Jones (LJ) potential [63]:

$\Phi_{\text{LJ}}(\mathbf{r})=\sum_{i<j}4\epsilon_{ij}\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right],$
(4)

where $r_{ij}\equiv|\mathbf{r}_{i}-\mathbf{r}_{j}|$ is the distance between
particles $i$ and $j$ and $\sigma_{ij}$ and $\epsilon_{ij}$ represent the
length and energy scales, respectively. We also often employ the soft-core
(SC) potential, which is the repulsive part of the LJ potential:

$\Phi_{\text{SC}}(\mathbf{r})=\sum_{i<j}\epsilon_{ij}\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}.$
(5)

In addition, as the simplest model of a glass, a finite-range harmonic
potential has been employed:

$\Phi_{\text{HA}}(\mathbf{r})=\sum_{i<j}\frac{\epsilon_{ij}}{2}\left(1-\frac{r_{ij}}{\sigma_{ij}}\right)^{2}H(\sigma_{ij}-r_{ij}),$
(6)

where $H(x)$ is the Heaviside step function: $H(x)=1$ for $x\geq 0$ and
$H(x)=0$ for $x<0$. This harmonic potential was originally proposed for
modelling granular materials, emulsions, foams, etc. [64].

  2. (II)

Covalently bonded glasses (network glasses):
To model silica (SiO2) glass, which is one example of a covalently bonded
glass, the Beest-Kramer-Santen potential [65] has been proposed:

$\Phi_{\text{BKS}}(\mathbf{r})=\sum_{i<j}\left[\frac{q_{i}q_{j}e^{2}}{r_{ij}}+A_{ij}\exp(-B_{ij}r_{ij})-\frac{C_{ij}}{r_{ij}^{6}}\right],$
(7)

where particles $i$ and $j$ are Si or O atoms, $q_{i}e$ is the electric charge
($e$ is elementary charge), and $A_{ij}$, $B_{ij}$, and $C_{ij}$ are
constants. The first term in Eq. (7) represents the Coulomb interaction. In
addition, the Stillinger-Weber potential [66], $\Phi_{\text{SW}}(\mathbf{r})$,
has been proposed for modelling amorphous silicon.
$\Phi_{\text{SW}}(\mathbf{r})$ includes three-body interactions to represent
the bending rigidity. For the explicit equation for $\Phi_{\text{SW}}$, please
see Ref. 66.

In addition, for metallic glasses, Daw and Baskes have proposed the embedded-
atom method (EAM) based on density functional theory [67]. Additionally, for
polymer glasses, we can perform coarse-grained simulations by using the
Kremer–Grest model [68], which treats polymer chains as linear series of
monomer beads (particles). In the following section, we take a system with the
SC potential [Eq. (5)] as an example. We present numerical results not only
for the glass but also for the crystal, which are obtained from Refs. 40, 41,
42. By comparing the glass and crystal, we discuss the characteristic features
of glasses. The values of quantities are presented in units of the mass ($m$),
length ($\sigma$), and energy ($\epsilon$) scales that are typical of the
constituent particles.

(a) Glass.

(b) Crystal.

Figure 1:  Radial distribution function $g(r)$ of a soft-core (SC) system at a
low temperature of $T=10^{-2}$, well below $T_{g}$ and $T_{m}$. (a) Glass. (b)
Face-centred cubic (FCC) crystal.

####  2 Glass and inherent structure

By means of an MD simulation, we obtain the glass as follows [4, 5]. We first
equilibrate the system in the liquid phase above the melting temperature,
$T_{m}$. We next rapidly quench the system below $T_{m}$. If we avoid
crystallization 333 To avoid crystallization, we may employ a poly-disperse
system composed of a number of different kinds of particles with different
mass, length, and energy scales. , we can maintain the liquid phase even below
$T_{m}$; such a phase is called a supercooled liquid phase. When we further
quench the system below the glass transition temperature, $T_{g}$, the system
is finally frozen in a disordered configuration, without crystallization. To
study the structural properties of this system, we can measure the radial
distribution function $g(r)$ or the static structure factor $S(q)$ [63].
Figure 1 presents the $g(r)$ functions for two SC systems, (a) glass and (b)
crystal, for comparison. We can see long-range order in the crystal, while the
glass shows only short-range order. Although glasses are solids like crystals
(i.e., they exhibit rigidity), they also possess disordered structures as
liquids do [4, 5].

In the glass phase below $T_{g}$, the particles vibrate around the inherent
structure, which is denoted by
$\mathbf{R}\equiv\left[\mathbf{R}_{1}^{T},\mathbf{R}_{2}^{T},\cdots,\mathbf{R}_{N}^{T}\right]^{T}$
444 Glasses may show ageing phenomena at a time scale that is much longer than
that of vibrations. Ageing induces rearrangements of some particles and alters
the inherent structure. In this chapter, we will not consider this effect. .
In the inherent structure, the system is in a state of mechanical equilibrium,
and the potential energy $\Phi(\mathbf{r})$ takes its minimum value in the
$3N$-dimensional space:

$\left.\frac{\partial\Phi}{\partial\mathbf{r}}\right|_{\mathbf{r}=\mathbf{R}}=0.$
(8)

Thus, we can numerically obtain $\mathbf{R}$ by minimizing $\Phi(\mathbf{r})$.
For this purpose, several numerical techniques have been established, e.g.,
the steepest descent method [69], the conjugate gradient method [69], and the
fire algorithm [70]. We note that minimizing the potential corresponds to
quenching the system to zero temperature, $T=0$, and the inherent structure is
the configuration of the particles at $T=0$.

###  3 Vibrational eigenmodes

In a glass (or any solid-state material), the constituent particles vibrate
around the inherent structure. These vibrational motions of particles are
excited along vibrational eigenmodes. The eigenmodes are therefore fundamental
to understanding the vibrational properties. In this section, we explain how
to obtain the eigenmodes by means of computational simulations. We then
describe the vDOS, which represents the statistics of the eigenmodes. We also
explain the order parameters used to characterize each vibrational eigenmode.
In particular, we introduce the participation ratio, which measures the extent
of localization, and the phonon order parameter, which measures the extent of
phonon-like vibrations.

####  1 General description

Let us introduce the displacement vector ($3N$-dimensional vector) of the
particles relative to the inherent structure $\mathbf{R}$ (multiplied by the
mass factor $\sqrt{\mathcal{M}}$) as follows:

$\mathbf{u}=\sqrt{\mathcal{M}}(\mathbf{r}-\mathbf{R}).$ (9)

We then expand the potential $\Phi(\mathbf{r})$ around $\mathbf{R}$ in a power
series of $\mathbf{u}$ as follows ($T$ denotes transposition):

$\Phi(\mathbf{r})=\Phi(\mathbf{R})+\frac{1}{2}\mathbf{u}^{T}\mathcal{D}\mathbf{u}+\mathcal{O}(|\mathbf{u}|^{3}).$
(10)

The first derivative of $\Phi(\mathbf{r})$ is zero (mechanical equilibrium),
as expressed in Eq. (8). $\mathcal{D}$ is the so-called dynamical matrix (a
$3N\times 3N$ matrix) which is the second derivative of $\Phi(\mathbf{r})$ [6,
7]:

$\mathcal{D}=\left.\frac{1}{\sqrt{\mathcal{M}}}\left[\frac{\partial^{2}\Phi}{\partial\mathbf{r}\partial\mathbf{r}^{T}}\right|_{\mathbf{r}=\mathbf{R}}\right]\frac{1}{\sqrt{\mathcal{M}}}.$
(11)

$\mathcal{D}$ contains $3\times 3$ matrix elements, $(\mathcal{D})_{ij}$,
corresponding to particles $i$ and $j$:

$(\mathcal{D})_{ij}=\left.\frac{1}{\sqrt{m_{i}m_{j}}}\left[\frac{\partial^{2}\Phi}{\partial\mathbf{r}_{i}\partial\mathbf{r}_{j}^{T}}\right|_{\mathbf{r}=\mathbf{R}}\right].$
(12)

In the case of a pair-wise potential, $\phi_{ij}(r_{ij})$, which depends only
on the distance $r_{ij}$ between particles $i$ and $j$, $(\mathcal{D})_{ij}$
is formulated as

$(\mathcal{D})_{ij}=\left\\{\begin{aligned} &\frac{1}{m_{i}}\sum_{k=1,k\neq
i}^{N}\left[\frac{d^{2}\phi_{ik}}{dr_{ik}^{2}}\frac{\mathbf{r}_{ik}\mathbf{r}^{T}_{ik}}{r_{ik}^{2}}+\frac{1}{r_{ik}}\frac{d\phi_{ik}}{dr_{ik}}\left(I_{3}-\frac{\mathbf{r}_{ik}\mathbf{r}^{T}_{ik}}{r_{ik}^{2}}\right)\right]&(i=j),\\\
&-\frac{1}{\sqrt{m_{i}m_{j}}}\left[\frac{d^{2}\phi_{ij}}{dr_{ij}^{2}}\frac{\mathbf{r}_{ij}\mathbf{r}^{T}_{ij}}{r_{ij}^{2}}+\frac{1}{r_{ij}}\frac{d\phi_{ij}}{dr_{ij}}\left(I_{3}-\frac{\mathbf{r}_{ij}\mathbf{r}^{T}_{ij}}{r_{ij}^{2}}\right)\right]&(i\neq
j),\end{aligned}\right.$ (13)

where $I_{3}$ is the $3\times 3$ unit matrix,
$\mathbf{r}_{ij}\equiv\mathbf{r}_{i}-\mathbf{r}_{j}$, and
$r_{ij}\equiv\left|\mathbf{r}_{ij}\right|$. Since $|\mathbf{u}|\ll 1$ at low
temperatures, we neglect higher-order terms of $\mathcal{O}(|\mathbf{u}|^{3})$
in $\Phi(\mathbf{r})$ as expressed in Eq. (10) and substitute this expression
into Eq. (2) to obtain the linearized equation of motion (harmonic
approximation):

$\frac{d^{2}\mathbf{u}(t)}{dt^{2}}=-\mathcal{D}\mathbf{u}(t).$ (14)

To solve Eq. (14), we perform a Fourier transform and obtain

$\omega^{2}\tilde{\mathbf{u}}(\omega)=\mathcal{D}\tilde{\mathbf{u}}(\omega),$
(15)

where $\omega$ is the frequency and $\tilde{\mathbf{u}}(\omega)$ is the
Fourier transform of $\mathbf{u}(t)$:

$\tilde{\mathbf{u}}(\omega)=\int\mathbf{u}(t)\exp\left(-\text{i}\omega
t\right)dt.$ (16)

$u(t)$ is then obtained by taking the inverse Fourier transform of
$\tilde{\mathbf{u}}(\omega)$:

${\mathbf{u}}(t)=\frac{1}{2\pi}\int\tilde{\mathbf{u}}(\omega)\exp\left(\text{i}\omega
t\right)d\omega.$ (17)

We thus encounter the eigenvalue problem of the matrix $\mathcal{D}$ [Eq.
(15)].

Since $\mathcal{D}$ is a symmetric matrix, it can be diagonalized by an
orthonormal matrix ${X}$, and its eigenvalues $\lambda^{k}$
($k=1,2,\cdots,3N$) are all real numbers:

${X}^{T}\mathcal{D}{X}=\text{diag}(\lambda^{1},\lambda^{2},\cdots,\lambda^{3N}),$
(18) $X^{T}X=I_{3N},$ (19)

where $I_{3N}$ is the $3N\times 3N$ unit matrix. We write the matrix in the
form of
$X\equiv\left[\mathbf{e}^{1},\mathbf{e}^{2},\cdots,\mathbf{e}^{3N}\right]$,
thereby introducing $3N$-dimensional eigenvectors,
$\mathbf{e}^{k}\equiv\left[\mathbf{e}_{1}^{kT},\mathbf{e}_{2}^{kT},\cdots,\mathbf{e}_{N}^{kT}\right]^{T}$
($k=1,2,\cdots,3N$). The orthonormality condition, Eq. (19), then gives

$\mathbf{e}^{k}\cdot\mathbf{e}^{l}=\sum_{i=1}^{N}\mathbf{e}^{k}_{i}\cdot\mathbf{e}^{l}_{i}=\delta_{kl},$
(20)

where $\delta_{kl}$ is the Kronecker delta 555 The set of $3N$ eigenvectors
$\left\\{\mathbf{e}^{1},\mathbf{e}^{2},\cdots,\mathbf{e}^{3N}\right\\}$ can be
treated as an orthonormal basis in $3N$-dimensional space for $\mathbf{u}$. .
We thus obtain $3N$ sets of eigenvalues $\lambda^{k}$ and eigenvectors
$\mathbf{e}^{k}$ for the matrix $\mathcal{D}$:

$\lambda^{k}\mathbf{e}^{k}=\mathcal{D}\mathbf{e}^{k}\quad(k=1,2,\cdots,3N),$
(21)

and $\mathcal{D}$ can be described as

$\mathcal{D}=\sum_{k=1}^{3N}\lambda^{k}\mathbf{e}^{k}\mathbf{e}^{kT}.$ (22)

If the system is stable, the eigenvalues are positive, $\lambda^{k}>0$, and
the eigenfrequencies, $\omega^{k}\equiv\sqrt{\lambda^{k}}$, are real numbers.
A general solution for $\tilde{\mathbf{u}}(\omega)$ in Eq. (15) is thus
obtained as a superposition of the $\mathbf{e}^{k}$:

$\tilde{\mathbf{u}}(\omega)=\sum_{k=1}^{3N}A^{k}\delta(\omega-\omega^{k})\mathbf{e}^{k},$
(23)

where the $A^{k}$ are complex constants, and $\mathbf{u}(t)$ is

$\mathbf{u}(t)=\frac{1}{2\pi}\sum_{k=1}^{3N}A^{k}\exp\left(\text{i}\omega^{k}t\right)\mathbf{e}^{k}.$
(24)

The values of the $A^{k}$ are determined from the initial conditions
$\mathbf{u}=\mathbf{u}_{0}$ and $d\mathbf{u}/dt=\dot{\mathbf{u}}_{0}$ at
$t=0$, as
$\text{Re}(A^{k})=2\pi\left(\mathbf{e}^{k}\cdot\mathbf{u}_{0}\right)$ and
$\text{Im}(A^{k})=-2\pi\left(\mathbf{e}^{k}\cdot\dot{\mathbf{u}}_{0}\right)/\omega^{k}$
(where Re and Im denote the real and imaginary parts, respectively), and we
finally obtain $\mathbf{u}(t)$ as follows:

$\displaystyle\mathbf{u}(t)$
$\displaystyle=\sum_{k=1}^{3N}u^{k}(t)\mathbf{e}^{k},$ (25) $\displaystyle
u^{k}(t)$
$\displaystyle=\left(\mathbf{e}^{k}\cdot\mathbf{u}_{0}\right)\cos\left(\omega^{k}t\right)+\left(\mathbf{e}^{k}\cdot\dot{\mathbf{u}}_{0}\right)\frac{\sin\left(\omega^{k}t\right)}{\omega^{k}},$

which is a superposition of terms of the form $u^{k}(t)\mathbf{e}^{k}$, which
represents the vibrations of particles along the eigenvector $\mathbf{e}^{k}$
(vibration of particle $i$ along $\mathbf{e}_{i}^{k}$) with frequency
$\omega^{k}$ and which we call the vibrational eigenmodes $k$.

The kinetic energy, $K$, and the potential energy,
$\Delta\Phi\equiv\Phi-\Phi(\mathbf{R})$, are described as sums of the energies
of the eigenmodes $k$:

$\displaystyle K(t)$
$\displaystyle=\frac{1}{2}\frac{d\mathbf{u}(t)}{dt}\cdot\frac{d\mathbf{u}(t)}{dt}=\frac{1}{2}\sum_{k=1}^{3N}\left(\frac{d{u}^{k}(t)}{dt}\right)^{2},$
(26) $\displaystyle\Delta\Phi(t)$
$\displaystyle=\frac{1}{2}\mathbf{u}(t)^{T}\mathcal{D}\mathbf{u}(t)=\frac{1}{2}\sum_{k=1}^{3N}\left(\omega^{k}{u}^{k}(t)\right)^{2}.$

The total energy,
$K(t)+\Delta\Phi(t)=\left(1/2\right)\sum_{k=1}^{3N}\left[{\omega^{k}}^{2}\left(\mathbf{e}^{k}\cdot\mathbf{u}_{0}\right)^{2}+\left(\mathbf{e}^{k}\cdot\dot{\mathbf{u}}_{0}\right)^{2}\right]$,
is constant, independent of time (energy conservation).

####  2 Vibrational eigenmodes in elastic media (elastic waves)

As the simplest case, let us consider the eigenmodes in an elastic medium,
i.e., elastic waves. Elastic media are continuum systems (i.e., not
particulate systems), and their eigenfrequencies and eigenvectors are
determined by continuum mechanics [6, 7].

An elastic wave is specified by the wavevector $\mathbf{q}$ and the
polarization $\alpha$. In the 3-dimensional case, there are three
polarizations: two transverse ($\alpha=T_{1},\ T_{2}$) and one longitudinal
($\alpha=L$). An eigenvector of elastic waves is a continuous function of the
position $\mathbf{r}\in V$ (where $V$ is the volume of the system) 666 The
continuous variable $\mathbf{r}$ for an elastic medium corresponds to the
particle index “$i$” in a particulate system. :

$\mathbf{e}^{\mathbf{q},\alpha}_{\text{el}}\equiv\left[\mathbf{e}^{\mathbf{q},\alpha}_{\text{el}}(\mathbf{r});\
\mathbf{r}\in
V\right]\equiv\left[\mathbf{s}_{\alpha}(\hat{\mathbf{q}})\frac{\exp(\text{i}\mathbf{q}\cdot\mathbf{r})}{\sqrt{V}};\
\mathbf{r}\in V\right].$ (27)

$\mathbf{s}_{\alpha}(\hat{\mathbf{q}})$ is the polarization vector, which
depends on the direction of $\mathbf{q}$, i.e.,
$\hat{\mathbf{q}}\equiv\mathbf{q}/q$ (where $q\equiv|\mathbf{q}|$ is the
wavenumber), and can be determined from the elastic equation of motion that
corresponds to the eigen equation (21):

$\omega^{2}\mathbf{s}_{\alpha}(\hat{\mathbf{q}})=\left(\frac{\mathcal{C}(\hat{\mathbf{q}})}{\rho}\right)q^{2}\mathbf{s}_{\alpha}(\hat{\mathbf{q}}),$
(28)

where $\rho$ is the mass density and $\mathcal{C}(\hat{\mathbf{q}})$ (a
$3\times 3$ tensor) is the elastic modulus tensor in the $\hat{\mathbf{q}}$
direction 777 Given the elastic modulus tensor $C_{\alpha\beta\gamma\delta}$
(a $3\times 3\times 3\times 3$ tensor), where
$\alpha,\beta,\gamma,\delta=x,y,z$, $\mathcal{C}(\hat{\mathbf{q}})$ is defined
as
$\left(\mathcal{C}(\hat{\mathbf{q}})\right)_{\alpha\delta}=\sum_{\beta,\gamma}C_{\alpha\beta\gamma\delta}\hat{q}_{\beta}\hat{q}_{\gamma}$,
where $\hat{\mathbf{q}}\equiv(\hat{q}_{x},\hat{q}_{y},\hat{q}_{z})$. We will
describe the measurement of $C_{\alpha\beta\gamma\delta}$ in Section 5. .

Since $\left(\mathcal{C}(\hat{\mathbf{q}})/\rho\right)q^{2}$ (a $3\times 3$
matrix) is positive and symmetric, its eigenvalues $\omega^{2}$ are real,
positive numbers, and the eigenvectors $\mathbf{s}_{\alpha}(\hat{\mathbf{q}})$
are orthonormalized as
$\mathbf{s}_{\alpha}(\hat{\mathbf{q}})\cdot\mathbf{s}_{\alpha^{\prime}}(\hat{\mathbf{q}})=\delta_{\alpha\alpha^{\prime}}$.
The eigenvectors $\mathbf{e}^{\mathbf{q},\alpha}_{\text{el}}(\mathbf{r})$ are
then orthonormalized as 888 For an elastic medium, we need to replace
$\sum_{i=1}^{N}$ in Eq. (20) with $\int_{V}d^{3}\mathbf{r}$.

$\mathbf{e}^{\mathbf{q},\alpha}_{\text{el}}\cdot\mathbf{e}^{\mathbf{q}^{\prime},\alpha^{\prime}\ast}_{\text{el}}=\int_{V}d^{3}\mathbf{r}\
\mathbf{e}^{\mathbf{q},\alpha}_{\text{el}}(\mathbf{r})\cdot\mathbf{e}^{\mathbf{q}^{\prime},\alpha^{\prime}\ast}_{\text{el}}(\mathbf{r})=\delta_{\mathbf{q}\mathbf{q}^{\prime}}\delta_{\alpha\alpha^{\prime}},$
(29)

where $\ast$ denotes the complex conjugate. Eq. (28) determines the
eigenfrequency $\omega_{\alpha}(\mathbf{q})$ as a linear function of $q$:

$\omega_{\alpha}(\mathbf{q})=c_{\alpha}(\hat{\mathbf{q}})q,$ (30)

where $c_{\alpha}(\hat{\mathbf{q}})$ is the sound speed, which is the square
root of the eigenvalue of $\mathcal{C}(\hat{\mathbf{q}})/\rho$. We note that
Eq. (28) is the so-called Christoffel equation, which can be solved
analytically following the Every algorithm [71].

In the case of an isotropic medium, the elastic modulus tensor contains two
independent moduli: the shear modulus $G$ and the bulk modulus $K$.
$\mathbf{s}_{T_{1}}(\hat{\mathbf{q}})$ and
$\mathbf{s}_{T_{2}}(\hat{\mathbf{q}})$ are perpendicular to $\mathbf{q}$,
while $\mathbf{s}_{L}(\hat{\mathbf{q}})$ is parallel to $\mathbf{q}$:
$\hat{\mathbf{q}}\cdot\mathbf{s}_{T_{1}}(\hat{\mathbf{q}})=\hat{\mathbf{q}}\cdot\mathbf{s}_{T_{2}}(\hat{\mathbf{q}})=0$
and $\mathbf{s}_{L}(\hat{\mathbf{q}})=\hat{\mathbf{q}}$.
$c_{\alpha}(\hat{\mathbf{q}})$ does not depend on $\hat{\mathbf{q}}$, and it
is calculated as

$c_{T_{1}}=c_{T_{2}}=\sqrt{\frac{G}{\rho}},\qquad
c_{L}=\sqrt{\frac{K+4G/3}{\rho}}.$ (31)

Figure 2:  Schematic illustration of the polarization vector
$\mathbf{s}_{\alpha}(\hat{\mathbf{q}})$ in a medium with cubic symmetry. We
show $\mathbf{s}_{T_{1}}(\hat{\mathbf{q}})$,
$\mathbf{s}_{T_{2}}(\hat{\mathbf{q}})$, and $\mathbf{s}_{L}(\hat{\mathbf{q}})$
in three directions of $\hat{\mathbf{q}}$: $[100]$, $[110]$, and $[111]$
(Miller indices).

On the other hand, in the case of an anisotropic medium,
$\mathbf{s}_{\alpha}(\hat{\mathbf{q}})$ and $c_{\alpha}(\hat{\mathbf{q}})$
depend on $\hat{\mathbf{q}}$. The simplest example is a medium with cubic
symmetry, which has three independent moduli: the pure shear modulus $G_{p}$,
the simple shear modulus $G_{s}$, and the bulk modulus $K$ 999 $G_{p}$,
$G_{s}$, and $K$ correspond to the elastic moduli for pure shear, simple
shear, and bulk deformations, respectively. These three deformations are
illustrated in Fig. 11 of Section 5. . Figure 2 shows schematic illustrations
of $\mathbf{s}_{\alpha}(\hat{\mathbf{q}})$ in three directions, namely,
$[100]$, $[110]$, and $[111]$, where we use Miller indices 101010
$\hat{\mathbf{q}}=(1,0,0)$ and all equivalent vectors for $[100]$,
$\hat{\mathbf{q}}=(1/\sqrt{2},1/\sqrt{2},0)$ and all equivalent vectors for
$[110]$, and $\hat{\mathbf{q}}=(1/\sqrt{3},1/\sqrt{3},1/\sqrt{3})$ and all
equivalent vectors for $[111]$. . Additionally, $c_{\alpha}(\hat{\mathbf{q}})$
is calculated as

$c_{T_{1}}=c_{T_{2}}=\sqrt{\frac{G_{s}}{\rho}},\quad
c_{L}=\sqrt{\frac{K+4G_{p}/3}{\rho}}\quad\text{for $[100]$},$ (32)
$c_{T_{1}}=\sqrt{\frac{G_{s}}{\rho}},\quad
c_{T_{2}}=\sqrt{\frac{G_{p}}{\rho}},\quad
c_{L}=\sqrt{\frac{K+G_{p}/3+G_{s}}{\rho}}\quad\text{for $[110]$},$ (33)
$c_{T_{1}}=c_{T_{2}}=\sqrt{\frac{2G_{p}+G_{s}}{3\rho}},\quad
c_{L}=\sqrt{\frac{K+4G_{s}/3}{\rho}}\quad\text{for $[111]$}.$ (34)

A general solution for $c_{\alpha}(\hat{\mathbf{q}})$ as a function of
$\hat{\mathbf{q}}$ is given in Ref. 72.

(a) Glass.

(b) Crystal.

(c) Glass.

(d) Crystal.

Figure 3: Visualization of the eigenmodes in SC systems. (a),(c) Glass.
(b),(d) FCC crystal. We show the eigenvector field,
$\mathbf{e}^{k}\equiv\left[\mathbf{e}_{1}^{kT},\mathbf{e}_{2}^{kT},\cdots,\mathbf{e}_{N}^{kT}\right]^{T}$,
as arrows. In (c),(d), the eigenvectors on some plane are presented in a two-
dimensional plot.

####  3 Vibrational eigenmodes in crystals (phonons)

We next consider the eigenmodes in crystals. The inherent structure
$\mathbf{R}$ of a crystal is a periodic lattice structure (we show the radial
distribution function $g(r)$ in Fig. 1(b)). In this case, the eigenmodes have
been established to be phonons (lattice waves) [6, 7]. Here, we consider
single-component crystals, for which the eigenmodes are phonons of the
acoustic type. For multi-component crystals, phonons of the optical type also
appear, which are not considered below.

A phonon is specified by the wavevector $\mathbf{q}$ and the polarization
$\alpha$, as is an elastic wave. There are $N$ wavevectors in the first
Brillouin zone,
$\mathbf{q}=\mathbf{q}_{1},\mathbf{q}_{2},\cdots,\mathbf{q}_{N}$, and three
polarizations, $\alpha=T_{1},\ T_{2},\ L$; therefore, there are $3N$ phonons
in total. The eigenvector $\mathbf{e}_{\text{ph}}^{\mathbf{q},\alpha}$ is
described as

$\mathbf{e}^{\mathbf{q},\alpha}_{\text{ph}}\equiv\left[\mathbf{s}_{\alpha}^{T}(\mathbf{q})\frac{\exp(\text{i}\mathbf{q}\cdot\mathbf{R}_{1})}{\sqrt{N}},\mathbf{s}_{\alpha}^{T}(\mathbf{q})\frac{\exp(\text{i}\mathbf{q}\cdot\mathbf{R}_{2})}{\sqrt{N}},\cdots,\mathbf{s}_{\alpha}^{T}(\mathbf{q})\frac{\exp(\text{i}\mathbf{q}\cdot\mathbf{R}_{N})}{\sqrt{N}}\right]^{T},$
(35)

where the factor $1/\sqrt{N}$ is necessary for the orthonormality condition,
$\mathbf{e}^{\mathbf{q},\alpha}_{\text{ph}}\cdot\mathbf{e}^{\mathbf{q}^{\prime},\alpha^{\prime}\ast}_{\text{ph}}=\delta_{\mathbf{q}\mathbf{q}^{\prime}}\delta_{\alpha\alpha^{\prime}}$.
The eigen equation (21) then reduces to

$\omega^{2}\mathbf{s}_{\alpha}({\mathbf{q}})={D}({\mathbf{q}})\mathbf{s}_{\alpha}({\mathbf{q}}),$
(36)

where ${D}({\mathbf{q}})$ (a $3\times 3$ matrix) is the dynamical matrix of
the unit cell of the lattice structure:

${D}({\mathbf{q}})=\frac{1}{N}\sum_{i,j=1}^{N}\left(\mathcal{D}\right)_{ij}\exp(-\text{i}\mathbf{q}\cdot(\mathbf{R}_{i}-\mathbf{R}_{j})).$
(37)

Since ${D}({\mathbf{q}})$ is positive and symmetric, Eq. (36) gives the real
number of eigenfrequencies $\omega_{\alpha}(\mathbf{q})$ and polarization
vectors $\mathbf{s}_{\alpha}(\mathbf{q})$ with the orthonormality condition,
$\mathbf{s}_{\alpha}({\mathbf{q}})\cdot\mathbf{s}_{\alpha^{\prime}}({\mathbf{q}})=\delta_{\alpha\alpha^{\prime}}$.
The eigenvectors of the phonons in Eq. (35) are of the same form as those of
the elastic waves in Eq. (27); however, we remark that they are discretized
waves in a particulate system. As an example, we visualize a transverse phonon
in Figs. 3(b) and 3(d). We can see the vortex structure that is characteristic
of transverse phonons.

An important feature is that the dynamical matrix, ${D}({\mathbf{q}})$,
converges to that of elastic waves,
$(\mathcal{C}(\hat{\mathbf{q}})/\rho)q^{2}$, at low wavenumbers $q\ll 1$ with
111111 We can show that
$\left(\mathcal{C}(\hat{\mathbf{q}})\right)_{\alpha\delta}=\sum_{\beta,\gamma}\left(C^{B}_{\alpha\beta\gamma\delta}+C^{C}_{\alpha\beta\gamma\delta}\right)\hat{q}_{\beta}\hat{q}_{\gamma}$
with $\hat{\mathbf{q}}\equiv(\hat{q}_{x},\hat{q}_{y},\hat{q}_{z})$, where
$C^{B}_{\alpha\beta\gamma\delta}$ and $C^{C}_{\alpha\beta\gamma\delta}$ are
respectively the Born term and the correction term of elastic modulus tensor,
given in Eq. (85) of Section 5.

$\frac{\mathcal{C}(\hat{\mathbf{q}})}{\rho}=-\frac{1}{2N}\sum_{i,j=1}^{N}\left(\mathcal{D}\right)_{ij}\left[\hat{\mathbf{q}}\cdot\left(\mathbf{R}_{i}-\mathbf{R}_{j}\right)\right]^{2},$
(38)

and thus, the phonons smoothly converge to elastic waves [6, 7] 121212
Acoustic phonons (the type considered here) converge to elastic waves at $q\ll
1$; however, optical phonons (not considered here) do not. Optical phonons
represent vibrations of particles within a unit cell (vibrations at
microscopic scales). . In particular, the phonons in cubic crystals, such as
FCC crystals, converge to the elastic waves in media with cubic symmetry. The
eigenfrequency $\omega_{\alpha}(\mathbf{q})$ of the phonons converges to a
linear function of $q$ as shown in Eq. (30).

####  4 Vibrational eigenmodes in glasses

Unlike for the phonons in crystals, we do not have any established formulation
for the eigenmodes in glasses. In this case, computational simulations are a
powerful tool. We can numerically diagonalize the dynamical matrix
$\mathcal{D}$ (solve the eigen equation (21)) to obtain the eigenfrequencies,
$\omega^{k}\equiv\sqrt{\lambda^{k}}$, and the eigenvectors,
$\mathbf{e}^{k}\equiv\left[\mathbf{e}_{1}^{kT},\mathbf{e}_{2}^{kT},\cdots,\mathbf{e}_{N}^{kT}\right]^{T}$
($k=1,2,\cdots,3N$). Note that when we diagonalize the $\mathcal{D}$ of a
crystal, we obtain phonons, as presented in Figs. 3(b) and 3(d). We also note
that for a 3-dimensional system under periodic boundary conditions, there
appear three zero-frequency modes that represent uniform translations of
particles. These modes are known as Goldstone modes and are usually discarded.

Figures 3(a) and 3(c) present an example of eigenmodes in an SC glass. We can
see localized vibration, in which some particles vibrate considerably while
other particles vibrate much less. This localization is a characteristic
feature of glasses (disordered systems) [26, 27, 28]. Although phonon-like
modes exist in glasses, the eigenmodes are non-phonon modes in general. We may
expect that at low frequencies, glasses will behave as uniform elastic media
and their eigenmodes will converge to elastic waves, as in the case of
crystals and phonons. However, recent simulations [50, 51, 52, 53, 54, 55, 56,
57, 58, 59] have demonstrated that this is not the case, as will be discussed
in Section 6.

(a) Glass.

(b) Crystal.

Figure 4:  The vDOS $g(\omega)$ of an SC system. (a) Glass. (b) FCC crystal.
The dashed line shows the Debye vDOS, $g_{D}(\omega)=A_{D}\omega^{2}$. The
Debye frequency $\omega_{D}$, the boson peak (BP) frequency
$\omega_{\text{BP}}$, and the position of the lowest-frequency, van Hove
singularity $\omega_{T_{2}S}$ are indicated by arrows.

(a) Glass.

(b) Crystal.

Figure 5:  The reduced vDOS $g(\omega)/\omega^{2}$ of an SC system. (a) Glass.
(b) FCC crystal. The dashed line shows the Debye level $A_{D}$. See also the
caption of Fig. 4.

####  5 Vibrational density of states

Now that we have obtained information on the eigenfrequencies $\omega^{k}$ and
eigenvectors $\mathbf{e}^{k}$, we can next analyse them to understand the
vibrational properties of glasses (or any solids). One of the most important
quantities is the vDOS $g(\omega)$. $g(\omega)$ measures how many eigenmodes
exist around the frequency $\omega$:

$g(\omega)=\frac{1}{3N}\sum_{k=1}^{3N}\delta(\omega-\omega^{k}),$ (39)

where $\delta(x)$ is the Dirac delta function. Figure 4(a) presents
$g(\omega)$ for an SC glass. For comparison, we also plot $g(\omega)$ for an
SC crystal in Fig. 4(b). Before looking at the results in these figures, let
us define the Debye vDOS $g_{D}(\omega)$ [6, 7], which is the vDOS calculated
by assuming glasses or crystals to be elastic media and is a useful reference
for capturing the features of the vDOS.

The eigenfrequencies of elastic waves are given by
$\omega_{\alpha}(\hat{\mathbf{q}})=c_{\alpha}(\hat{\mathbf{q}})q$ as in Eq.
(30). We can then calculate the vDOS by using Eq. (39) as follows (where we
replace the summation $\sum_{k=1}^{3N}$ with the integral
$(L/2\pi)^{3}\sum_{\alpha=T_{1},T_{2},L}\int d\mathbf{q}$):

$g_{D}(\omega)=A_{D}\omega^{2}=\frac{3}{\omega_{D}^{3}}\omega^{2},$ (40)

where $A_{D}\equiv{3}/{\omega_{D}^{3}}$ is the Debye level and $\omega_{D}$ is
the Debye frequency. $\omega_{D}\equiv\bar{c}k_{D}$, with
$k_{D}\equiv\left(6\pi^{2}{N}/{V}\right)^{1/3}$ being the Debye wavenumber and
$\bar{c}$ being the average sound speed:

$\bar{c}=\left(\frac{1}{3}\sum_{\alpha=T_{1},T_{2},L}\int\frac{d\Omega(\hat{\mathbf{q}})}{4\pi}\frac{1}{c_{\alpha}(\hat{\mathbf{q}})^{3}}\right)^{-1/3},$
(41)

where $d\Omega$ denotes the solid angle. In the case of an isotropic system
such as a glass, $\omega_{D}$ is given in the following simple form:

$\omega_{D}=\left[\frac{18\pi^{2}(N/V)}{2{c_{T}}^{-3}+{c_{L}}^{-3}}\right]^{1/3}.$
(42)

For an anisotropic system, we need to solve Eq. (28) to obtain
$c_{\alpha}(\hat{\mathbf{q}})$ as a function of $\hat{\mathbf{q}}$ [71]. An
analytical solution $c_{\alpha}(\hat{\mathbf{q}})$ for cubic crystals is given
in Ref. 72. Note that to calculate the Debye vDOS $g_{D}(\omega)$ (with Debye
level $A_{D}$ and Debye frequency $\omega_{D}$), we need to measure the
elastic moduli of the system, as described later in Section 5.

In Fig. 4, the Debye vDOS $g_{D}(\omega)$ is plotted as a dashed line. Since
the phonons in a crystal converge to elastic waves at low $\omega$, the
$g(\omega)$ function of the crystal converges to $g_{D}(\omega)$, as
demonstrated in Fig. 4(b). We can also clearly observe van Hove singularities
[6, 7]. In Fig. 4(b), we use $\omega_{T_{2}S}$ to indicate the position of the
lowest-frequency, van Hove singularity which corresponds to that of $T_{2}$
transverse phonons.

Compared to the crystal, the glass shows smoother variations and broader
distributions of the eigenmodes. Remarkably, the low-$\omega$ portion of
$g(\omega)$ shows an enhancement in eigenmodes over the Debye vDOS. To
demonstrate this point more clearly, Figure 5 presents the reduced vDOS, which
is $g(\omega)$ divided by $\omega^{2}$, i.e., $g(\omega)/\omega^{2}$. Note
that $g_{D}(\omega)/\omega^{2}$ is equal to the Debye level $A_{D}$, i.e., a
constant value, as shown by the dashed line in Fig. 5. For the crystal,
$g(\omega)/\omega^{2}$ coincides with $A_{D}$ at low $\omega$, whereas for the
glass, $g(\omega)/\omega^{2}$ clearly shows an excess peak over $A_{D}$, which
we call the boson peak (BP) [16, 17, 25] 131313 The excess eigenmodes result
in enhancement of the specific heat of a glass. This excess specific heat is
also called the boson peak [13, 14, 15]. . The frequency at which
$g(\omega)/\omega^{2}$ reaches a maximum is called the BP frequency and is
denoted by $\omega_{\text{BP}}$ ($\omega_{\text{BP}}\approx 0.9$ in Fig.
5(a)). We remark that $g(\omega)/\omega^{2}$ also shows a peak for the
crystal, as shown in Fig. 5(b); however, this peak corresponds to a van Hove
singularity and appears at a frequency much higher than $\omega_{\text{BP}}$
($\omega_{T_{2}S}\approx 7\gg\omega_{\text{BP}}\approx 0.9$). The van Hove
singularity picks up vibrations at a microscopic length scale on the order of
the lattice constant [6, 7], whereas the boson peak exhibits vibrations at a
much longer, mesoscopic length scale [32, 33].

As $\omega$ decreases below the BP frequency $\omega_{\text{BP}}$, the
$g(\omega)/\omega^{2}$ of the glass approaches the Debye level $A_{0}$, as
shown in Fig. 5(a). We might expect that $g(\omega)/\omega^{2}$ will converge
to $A_{0}$ at some frequency, where the eigenmodes will converge to elastic
waves. However, this is not the case, as will be discussed in detail in
Section 6.

####  6 Characterization of vibrational eigenmodes

We employ order parameters to characterize the vibrations of particles in each
eigenmode $k$. Here, we introduce two order parameters. One is the
participation ratio, which measures the extent of localization [26, 27, 28].
The other is the phonon order parameter, which measures the extent of phonon-
like vibrations [51, 52].

(a) Glass.

(b) Crystal.

Figure 6:  The participation ratio $\mathcal{P}^{k}$ of an SC system. (a)
Glass. (b) FCC crystal. The dashed line corresponds to $\mathcal{P}^{k}=2/3$,
the value for elastic waves. See also the caption of Fig. 4.

##### Participation ratio

As already seen in Fig. 3, vibrational localization occurs in glasses. The
participation ratio quantitatively measures the extent of localization [26,
27, 28]. Given an eigenvector,
$\mathbf{e}^{k}\equiv\left[\mathbf{e}_{1}^{kT},\mathbf{e}_{2}^{kT},\cdots,\mathbf{e}_{N}^{kT}\right]^{T}$,
for eigenmode $k$, its participation ratio $\mathcal{P}^{k}$ is calculated as

$\mathcal{P}^{k}\equiv\frac{1}{N}\left[\sum_{i=1}^{N}\mathbf{e}^{k}_{i}\cdot\mathbf{e}^{k}_{i}\right]^{2}\left[\sum_{i=1}^{N}\left(\mathbf{e}^{k}_{i}\cdot\mathbf{e}^{k}_{i}\right)^{2}\right]^{-1}=\frac{1}{N}\left[\sum_{i=1}^{N}\left(\mathbf{e}^{k}_{i}\cdot\mathbf{e}^{k}_{i}\right)^{2}\right]^{-1},$
(43)

where the last equality comes from the orthonormality condition in Eq. (20).
$\mathcal{P}^{k}$ measures the fraction of particles that participate in the
vibrations. As extreme cases,
$\mathcal{P}^{k}=\left(\sum_{i=1}^{N}m_{i}\right)^{2}/\left(N\sum_{i=1}^{N}m_{i}^{2}\right)$
for an ideal mode in which all constituent particles vibrate equally 141414
For the case of identical mass $m_{i}\equiv m$,
$\mathcal{P}^{k}=\left(\sum_{i=1}^{N}m_{i}\right)^{2}/\left(N\sum_{i=1}^{N}m_{i}^{2}\right)=1$.
, $\mathcal{P}^{k}=1/N\ll 1$ for an ideal mode involving only one particle,
and $\mathcal{P}^{k}=2/3$ for elastic waves,
$\mathbf{e}^{\mathbf{q},\alpha}_{\text{el}}$, as given in Eq. (27) 151515 For
elastic waves, we can calculate their participation ratio as
$\mathcal{P}^{k}\equiv(1/V)\left[\int_{V}d\mathbf{r}^{3}\
\text{Re}\left(\mathbf{e}^{\mathbf{q},\alpha}_{\text{el}}(\mathbf{r})\right)\cdot\text{Re}\left(\mathbf{e}^{\mathbf{q},\alpha}_{\text{el}}(\mathbf{r})\right)\right]^{2}\left[\int_{V}d\mathbf{r}^{3}\
\left\\{\text{Re}\left(\mathbf{e}^{\mathbf{q},\alpha}_{\text{el}}(\mathbf{r})\right)\cdot\text{Re}\left(\mathbf{e}^{\mathbf{q},\alpha}_{\text{el}}(\mathbf{r})\right)\right\\}^{2}\right]^{-1}=2/3$.
.

Figure 6 plots $\mathcal{P}^{k}$ as a function of the frequency $\omega^{k}$
for (a) an SC glass and (b) an SC crystal. We observe that the values of
$\mathcal{P}^{k}$ in the crystal fluctuate around $\mathcal{P}^{k}=2/3$,
indicating that the vibrations are similar to elastic waves. Indeed, the
eigenvectors of phonons, $\mathbf{e}^{\mathbf{q},\alpha}_{\text{ph}}$, have
the same form as those of elastic waves,
$\mathbf{e}^{\mathbf{q},\alpha}_{\text{el}}$ (compare Eqs. (35) and (27)).
However, phonons are not continuum waves but discrete waves, and their
discrete nature causes fluctuations in $\mathcal{P}^{k}$ around
$\mathcal{P}^{k}=2/3$.

On the other hand, the values of $\mathcal{P}^{k}$ in the glass are much
different to those in the crystal. Remarkably, the glass shows localization at
the low-$\omega$ and high-$\omega$ edges, which is a characteristic feature of
disordered matter. For the low-frequency case, localization starts to occur
near the BP frequency $\omega_{\text{BP}}$. In addition to the participation
ratio, we can also employ, for example, an effective mass [26] to measure the
extent of localization.

##### Phonon order parameter

As the frequency decreases, the eigenmodes in glasses tend to show phonon-like
vibrations. The phonon order parameter evaluates the extent to which an
eigenmode $\mathbf{e}^{k}$ is similar to phonons
$\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}$ [51, 52]:

$\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}\equiv\left[\mathbf{s}_{\alpha}^{T}(\mathbf{q})\frac{\exp(\text{i}\mathbf{q}\cdot\mathbf{R}_{1})}{\sqrt{N}},\mathbf{s}_{\alpha}^{T}(\mathbf{q})\frac{\exp(\text{i}\mathbf{q}\cdot\mathbf{R}_{2})}{\sqrt{N}},\cdots,\mathbf{s}_{\alpha}^{T}(\mathbf{q})\frac{\exp(\text{i}\mathbf{q}\cdot\mathbf{R}_{N})}{\sqrt{N}}\right]^{T}.$
(44)

Here we define $\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}$ to extend the
concept of phonons to any solid-states including amorphous states. The
eigenvector $\mathbf{e}^{k}$ can always be expanded in terms of a series of
$3N$ phonons $\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}$:
$\mathbf{e}^{k}=\sum_{\mathbf{q},\alpha}A^{k}_{\mathbf{q},\alpha}\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}$.
In the case of a crystal where the inherent structure $\mathbf{R}$ is a
lattice structure, $\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}$ are exactly
the same as the eigenvectors $\mathbf{e}^{\mathbf{q},\alpha}_{\text{ph}}$ in
Eq. (35), which are orthonormalized as
$\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}\cdot\mathbf{u}^{\mathbf{q}^{\prime},\alpha^{\prime}\ast}_{\text{ph}}=\delta_{\mathbf{q}\mathbf{q}^{\prime}}\delta_{\alpha\alpha^{\prime}}$.
For the case of a glass where $\mathbf{R}$ is a disordered structure, the
phonons $\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}$ (which are not the
eigenvectors in general) are not exactly identical to but can be approximately
orthonormalized as
$\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}\cdot\mathbf{u}^{\mathbf{q}^{\prime},\alpha^{\prime}\ast}_{\text{ph}}\approx\delta_{\mathbf{q}\mathbf{q}^{\prime}}\delta_{\alpha\alpha^{\prime}}$.
Thus, we can calculate the projection of $\mathbf{e}^{k}$ onto one particular
phonon, $\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}$, as follows:

$\displaystyle
O^{k}_{\mathbf{q},\alpha}=\left|A^{k}_{\mathbf{q},\alpha}\right|^{2}\approx\left|\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}\cdot\mathbf{e}^{k}\right|^{2}.$
(45)

Here, note that $\sum_{\mathbf{q},\alpha}O^{k}_{\mathbf{q},\alpha}\approx 1$
since $\mathbf{e}^{k}\cdot\mathbf{e}^{k}=1$.

If eigenmode $k$ shows the phonon vibration, then $\mathbf{e}^{k}$ can be
described as the sum of a finite number of phonons with large overlap:
$\mathbf{e}^{k}=\sum_{\mathbf{q},\alpha;\ O^{k}_{\mathbf{q},\alpha}\geq
N_{m}/3N}A^{k}_{\mathbf{q},\alpha}\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}$.
Here, we define “large overlap” to mean that $O^{k}_{\mathbf{q},\alpha}\geq
N_{m}/3N$, i.e., the phonons overlap by more than $N_{m}$ modes. Considering
this, we define the phonon order parameter $O^{k}$ as

$\displaystyle O^{k}$ $\displaystyle=\sum_{\mathbf{q},\alpha;\
O^{k}_{\mathbf{q},\alpha}\geq N_{m}/3N}O^{k}_{\mathbf{q},\alpha}.$ (46)

$O^{k}=1$ for a phonon, whereas $O^{k}=0$ for a mode that is considerably
different to a phonon. To calculate the phonon order parameter in Eq. (46), we
need to set an appropriate value of $N_{m}$; the obtained results and
conclusions should not depend on the choice of $N_{m}$. $N_{m}=100$ might be
appropriate, but this depends on the considered system. Simulation data for
the phonon order parameter will be presented in Section 6.

###  4 Phonon transport

Phonon transport is a fundamental property of solid-state materials [6, 7].
Phonons play a role in carrying heat energy. In crystals, phonons are the
vibrational eigenmodes, and they propagate without any attenuation in the
zero-temperature, harmonic limit. The scattering of phonons is induced through
anharmonic effects at finite temperatures. This situation is drastically
changed in glasses, where the eigenmodes are non-phonon modes, as seen in
Section 3. In this case, a phonon will be decomposed into several different
eigenmodes in some frequency range [29], and as a result, it will be strongly
scattered. The disordered structures therefore scatter phonons propagating in
glasses.

A phonon is specified by the wavevector $\mathbf{q}$ (wavenumber
$q\equiv|\mathbf{q}|$ and $\hat{\mathbf{q}}\equiv\mathbf{q}/q$) and the
polarization $\alpha$, where $\alpha=T_{1},\ T_{2}$ (transverse) or $L$
(longitudinal). Its transport is characterized by two quantities: the
propagation frequency $\Omega(\mathbf{q})$ and the attenuation rate
$\Gamma_{\alpha}(\mathbf{q})$. (For the frequency $\Omega(\mathbf{q})$, we do
not explicitly write the polarization index $\alpha$.) The propagation speed
is given by $c_{\alpha}(\mathbf{q})=\Omega(\mathbf{q})/q$ 161616
$c_{\alpha}(\mathbf{q})=\Omega(\mathbf{q})/q$ is the phase speed, while the
group speed is given by
$v_{\alpha}(\mathbf{q})=\left.d\Omega/dq\right|_{\hat{\mathbf{q}}}$. In the
low-wavenumber region, these two speeds coincide with each other. . Below, we
will describe how to measure $\Omega(\mathbf{q})$ or $c_{\alpha}(\mathbf{q})$,
and $\Gamma_{\alpha}(\mathbf{q})$ through simulations.

####  1 Dynamic structure factor

Let us introduce the dynamic structure factor to study phonon transport. The
dynamic structure factor can be measured through scattering experiments
involving light, X-rays, or neutrons. In the computational simulation
approach, we first perform an MD simulation at a finite temperature $T>0$ and
obtain the corresponding trajectory data,
$\mathbf{r}(t)\equiv\left[\mathbf{r}_{1}^{T}(t),\mathbf{r}_{2}^{T}(t),\cdots,\mathbf{r}_{N}^{T}(t)\right]^{T}$.
From this trajectory, the dynamic structure factor,
$S_{\alpha}(\mathbf{q},\omega)$ (where $\alpha=T$ or $L$), is calculated as
the Fourier transform of the current-current correlation function [37, 38, 39,
41]:

$S_{\alpha}(\mathbf{q},\omega)=\left(\frac{q}{\omega}\right)^{2}\frac{1}{2\pi}\int\frac{1}{N}\left<\mathbf{j}_{\alpha}(\mathbf{q},t)\cdot\mathbf{j}_{\alpha}^{\ast}(\mathbf{q},0)\right>_{0}\exp\left({\text{i}\omega
t}\right)dt,$ (47)

where $\left<\right>_{0}$ denotes the ensemble average over configurations at
$t=0$. $\mathbf{j}_{\alpha}(\mathbf{q},t)$ is the transverse ($\alpha=T$) or
longitudinal ($\alpha=L$) current:

$\displaystyle\mathbf{j}_{T}(\mathbf{q},t)$
$\displaystyle=\sum_{i=1}^{N}\left[\mathbf{v}_{i}(t)-\left(\mathbf{v}_{i}(t)\cdot\hat{\mathbf{q}}\right)\hat{\mathbf{q}}\right]\exp\left(\text{i}\mathbf{q}\cdot\mathbf{r}_{i}(t)\right),$
(48) $\displaystyle\mathbf{j}_{L}(\mathbf{q},t)$
$\displaystyle=\sum_{i=1}^{N}\left(\mathbf{v}_{i}(t)\cdot\hat{\mathbf{q}}\right)\hat{\mathbf{q}}\exp\left(\text{i}\mathbf{q}\cdot\mathbf{r}_{i}(t)\right),$

where $\mathbf{v}_{i}=d\mathbf{r}_{i}/dt$ is the velocity of particle $i$. For
the longitudinal case, the conservation law relates the current
$\mathbf{j}_{L}(\mathbf{q},t)$ and the number density
$\hat{\rho}(\mathbf{q},t)\equiv\sum_{i=1}^{N}\exp\left(\text{i}\mathbf{q}\cdot\mathbf{r}_{i}(t)\right)$
in the form

$\frac{\partial\hat{\rho}(\mathbf{q},t)}{\partial
t}=\text{i}\mathbf{q}\cdot\mathbf{j}_{L}(\mathbf{q},t),$ (49)

and $S_{L}(\mathbf{q},\omega)$ corresponds to the density-density correlation
function 171717 The factor $(q/\omega)^{2}$ in $S_{\alpha}(\mathbf{q},\omega)$
as expressed in Eq. (47) is necessary for $S_{L}(\mathbf{q},\omega)$ to
correspond to the density-density correlation function. .

We can understand the dynamic structure factor $S_{\alpha}(\mathbf{q},\omega)$
as follows. The phonon displacement vector,
$\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}$, is given in Eq. (44). When we
take the polarization vectors to be
$\mathbf{s}_{T}(\mathbf{q})=\hat{\mathbf{q}}^{\perp}$ (where $\perp$ denotes a
perpendicular vector) and $\mathbf{s}_{L}(\mathbf{q})=\hat{\mathbf{q}}$, the
current-current correlation function
$\left<\mathbf{j}_{\alpha}(\mathbf{q},t)\cdot\mathbf{j}^{\ast}_{\alpha}(\mathbf{q},0)\right>_{0}$
is formulated as

$\frac{1}{N}\left<\mathbf{j}_{\alpha}(\mathbf{q},t)\cdot\mathbf{j}^{\ast}_{\alpha}(\mathbf{q},0)\right>_{0}=\left<\left(\mathbf{v}(t)\cdot\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}\right)\left(\mathbf{v}(0)\cdot\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}\right)^{\ast}\right>_{0}\times\left[1+\mathcal{O}\left(\sqrt{T}\right)\right],$
(50)

where
$\mathbf{v}(t)\equiv\left[\mathbf{v}_{1}^{T}(t),\mathbf{v}_{2}^{T}(t),\cdots,\mathbf{v}_{N}^{T}(t)\right]^{T}$
is the velocity vector of the system. Therefore,
$S_{\alpha}(\mathbf{q},\omega)$ at the low-$T$ harmonic limit converges to the
correlation function of the velocity projected onto
$\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}$. Thermal energy activates many
different phonons with different wavevectors and different polarizations. From
among these activated phonons, $S_{\alpha}(\mathbf{q},\omega)$ picks up a
specific phonon $\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}$ of wavevector
$\mathbf{q}$ and polarization $\alpha$. We remark that the actual polarization
vectors of the phonons can deviate from
$\mathbf{s}_{T}(\mathbf{q})=\hat{\mathbf{q}}^{\perp}$ and
$\mathbf{s}_{L}(\mathbf{q})=\hat{\mathbf{q}}$. However, they still
considerably overlap with
$\mathbf{s}_{T}(\mathbf{q})=\hat{\mathbf{q}}^{\perp}$ and
$\mathbf{s}_{L}(\mathbf{q})=\hat{\mathbf{q}}$, and we can properly analyse
phonon transport by means of $S_{\alpha}(\mathbf{q},\omega)$.

(a) Glass, $\alpha=T$.

(b) Crystal, $\alpha=T$.

(c) Glass, $\alpha=L$.

(d) Crystal, $\alpha=L$.

Figure 7:  The dynamic structure factor $S_{\alpha}(\mathbf{q},\omega)$ in the
$[110]$ direction for an SC system at a low temperature of $T=10^{-2}$, well
below $T_{g}$ and $T_{m}$. (a) Glass, $\alpha=T$. (b) FCC crystal, $\alpha=T$.
(c) Glass, $\alpha=L$. (d) FCC crystal, $\alpha=L$.

Figure 7 shows $S_{\alpha}(\mathbf{q},\omega)$ in the $[110]$ direction
(Miller indices) for an SC system. The temperature is very low, $T=10^{-2}$,
well below the glass transition temperature $T_{g}$ and the melting
temperature $T_{m}$. We plot $S_{\alpha}(\mathbf{q},\omega)$ for a glass and
for an FCC crystal. We can clearly observe the Brillouin peaks for both the
glass and the crystal. For the case of the FCC crystal,
$S_{T}(\mathbf{q},\omega)$ for transverse phonons shows two peaks at a fixed
$q$, as shown in Fig. 7(b), which correspond to a $T_{2}$ phonon at the lower
$\omega$ and a $T_{1}$ phonon at the higher $\omega$ (see Fig. 2 for the
polarization vectors $\mathbf{s}_{T_{2}}$ and $\mathbf{s}_{T_{1}}$). On the
other hand, the longitudinal $S_{L}(\mathbf{q},\omega)$ shows a single
Brillouin peak, as shown in Fig. 7(d).

Unlike the FCC crystal, since the glass is an isotropic system, it shows a
single Brillouin peak in both the transverse and longitudinal dynamic
structure factors, as shown in Figs. 7(a) and 7(c). Notably, the Brillouin
peaks are rather broad compared to those in the crystal. This finding
indicates that phonon attenuation is considerably enhanced in glasses, as we
will see below.

Based on the data for $S_{\alpha}(\mathbf{q},\omega)$, the propagation
frequency $\Omega(\mathbf{q})$ and the attenuation rate
$\Gamma_{\alpha}(\mathbf{q})$ can be extracted by fitting the spectral region
around the Brillouin peak to the damped harmonic oscillator model [37, 38, 39,
41]:

$S_{\alpha}(\mathbf{q},\omega)\propto\frac{\Gamma_{\alpha}(\mathbf{q})\Omega^{2}(\mathbf{q})}{\left(\omega^{2}-\Omega^{2}(\mathbf{q})\right)^{2}+\omega^{2}\Gamma_{\alpha}^{2}(\mathbf{q})}.$
(51)

This fitting procedure in the frequency domain is equivalent to fitting the
correlation function
$\left<\mathbf{j}_{\alpha}(\mathbf{q},t)\cdot\mathbf{j}^{\ast}_{\alpha}(\mathbf{q},0)\right>_{0}$
to the following function in the time domain:

$\frac{1}{N}\left<\mathbf{j}_{\alpha}(\mathbf{q},t)\cdot\mathbf{j}^{\ast}_{\alpha}(\mathbf{q},0)\right>_{0}\propto\cos\left(\Omega(\mathbf{q})t\right)\exp\left(-\frac{\Gamma_{\alpha}(\mathbf{q})}{2}t\right).$
(52)

In the following, we will show data for $\Omega(\mathbf{q})$,
$c_{\alpha}(\mathbf{q})\equiv\Omega(\mathbf{q})/q$, and
$\Gamma_{\alpha}(\mathbf{q})$ that are extracted from the
$S_{\alpha}(\mathbf{q},\omega)$ presented in Fig. 7.

(a) Glass, $\alpha=T$.

(b) Crystal, $\alpha=T_{2}$ in $[110]$.

(c) Glass, $\alpha=L$.

(d) Crystal, $\alpha=L$ in $[110]$.

Figure 8:  The propagation speed $c_{\alpha}(\Omega)$ for an SC system. (a)
Glass, $\alpha=T$. (b) FCC crystal, $\alpha=T_{2}$ in $[110]$. (c) Glass,
$\alpha=L$. (d) FCC crystal, $\alpha=L$ in $[110]$. For the glass [(a),(c)]
which is an isotropic system, we plot the data for the different directions
($[100]$, $[110]$, $[111]$) and for $\alpha=T_{1},\ T_{2}$ all together, which
coincide with each other. The dashed line indicates the macroscopic sound
speed calculated from the elastic moduli, as $c_{T0}=\sqrt{G/\rho}$ and
$c_{L0}=\sqrt{(K+4G/3)/\rho}$ as in Eq. (31). For the crystal [(b),(d)], the
macroscopic sound speed for $[110]$ is plotted by the dashed line, as
$c_{T_{2}0}=\sqrt{G_{p}/\rho}$ or $c_{L0}=\sqrt{(K+G_{p}/3+G_{s})/\rho}$ as in
Eq. (33). The BP frequency $\omega_{\text{BP}}$ and the position of the
lowest-frequency, van Hove singularity $\omega_{T_{2}S}$ are indicated by
arrows.

(a) Glass, $\alpha=T,\ L$.

(b) Crystal, $\alpha=T_{2},\ L$ in $[110]$.

Figure 9:  The attenuation rate $\Gamma_{\alpha}(\Omega)$ for an SC system.
(a) Glass, $\alpha=T,\ L$. (b) FCC crystal, $\alpha=T_{2},\ L$ in $[110]$. The
dashed line indicates the power-law scaling with the frequency $\Omega$. See
also the caption of Fig. 8.

####  2 Propagation frequency and attenuation rate

Figures 8 and 9 plot the sound speed $c_{\alpha}(\Omega)$ and the attenuation
rate $\Gamma_{\alpha}(\Omega)$, respectively, as functions of the frequency
$\Omega$. Here, we treat $c_{\alpha}$ and $\Gamma_{\alpha}$ as functions of
$\Omega$ by transforming $\mathbf{q}$ into $\Omega$ via the relation
$\Omega=\Omega(\mathbf{q})$. For the FCC crystal, we focus on the lower-
frequency transverse ($\alpha=T_{2}$) and longitudinal ($\alpha=L$) phonons in
the $[110]$ direction. On the other hand, for the glass, which is an isotropic
system, there are no distinctions among the different directions or between
$\alpha=T_{1}$ and $T_{2}$, so we show the results for all three different
directions ($[100]$, $[110]$, $[111]$) and for $\alpha=T_{1},\ T_{2}$ all
together; these results indeed coincide with each other, as seen in Figs. 8
and 9.

Let us first consider the data for the crystal. Since the phonons in crystals
converge to elastic waves at low $\Omega$, $c_{\alpha}(\Omega)$ converges to
the macroscopic value $c_{\alpha 0}$, i.e., $c_{T_{2}0}=\sqrt{G_{p}/\rho}$ or
$c_{L0}=\sqrt{(K+G_{p}/3+G_{s})/\rho}$ (see Eq. (33)), as is indeed shown in
Figs. 8(b) and 8(d). In this frequency region, the dispersion curve is a
straight line, $\Omega(\mathbf{q})=c_{\alpha 0}q$, as in Eq. (30). It starts
to deviate from linearity as $\Omega$ (or $q$) increases towards the edge of
the first Brillouin zone, where the van Hove singularity occurs [6, 7]. In
addition, both the transverse and longitudinal $\Gamma_{\alpha}(\Omega)$ show
an $\Omega^{2}$ dependence, $\Gamma_{\alpha}\propto\Omega^{2}$, as seen in
Fig. 9(b). This damping originates from purely anharmonic couplings between
phonons through the Umklapp process at finite $T>0$. The $\omega^{2}$
dependence can be explained by thermo-elasticity theory or the Boltzmann
equation analysis [73, 74]. At $T=0$, there are no anharmonic effects, and the
attenuation becomes exactly zero.

The sound speed $c_{\alpha}(\Omega)$ of the glass also converges to its
macroscopic value, $c_{T0}=\sqrt{G/\rho}$ or $c_{L0}=\sqrt{(K+4G/3)/\rho}$
(see Eq. (31)), at low $\Omega$, as shown in Figs. 8(a) and 8(c) 181818 The
longitudinal sound speed does not converge to $c_{L0}$ in the frequency range
presented in Fig. 8(c). However, recent work has shown that it converges as
$\Omega$ decreases [57]. . As $\Omega$ increases, $c_{\alpha}(\Omega)$
decreases and reaches a minimum. This behaviour has been observed in
experiments [19, 20] and simulations [38, 39, 41], and it is called sound
softening. The softening of transverse phonons occurs near the BP frequency
$\omega_{\text{BP}}$. Above the softening region, $c_{\alpha}(\Omega)$
increases with increasing $\Omega$; this behaviour is called sound hardening
[44].

Notably, the attenuation rate $\Gamma_{\alpha}(\Omega)$ of the glass behaves
much differently from that of the crystal, as shown in Fig. 9(a). First,
$\Gamma_{\alpha}(\Omega)$ takes much larger values than in the crystal.
Second, $\Gamma_{\alpha}(\Omega)$ shows Rayleigh scattering,
$\Gamma_{\alpha}\propto\Omega^{4}$, at low $\Omega$ below
$\omega_{\text{BP}}$, and it depends on $\Omega^{2}$, i.e.,
$\Gamma_{\alpha}\propto\Omega^{2}$, at high $\Omega$ above
$\omega_{\text{BP}}$ 191919 The longitudinal phonons do not show Rayleigh
scattering behaviour in the frequency range presented in Fig. 9(b). However,
recent work has shown that a Rayleigh scattering regime emerges as $\Omega$
decreases [57]. . These anomalous attenuation behaviours have also been
observed in experiments [18, 19, 20] and simulations [38, 39, 41]. There are
two sources of scattering in glasses: one is anharmonic effects, as in
crystals, and the other is the disordered structure. The latter effects are
dominant over the former and induce much stronger sound damping compared to
that in crystals; in particular, Rayleigh scattering at low $\Omega$ and an
$\Omega^{2}$ dependence at high $\Omega$ are induced. Although we cannot
observe any evidence of them in Fig. 9(a), the former, anharmonic effects
should become visible as some dependence on $\Omega$, possibly a dependence on
$\Omega^{2}$ as in the crystal [73, 74], in the lower $\Omega$ regime. Indeed,
both experimental works [22, 21] and theoretical works [75, 76, 77] have
studied these anharmonic effects, which add to the effects of the disordered
structure. Notably, close to the glass transition, the fractal frequency
dependence of the damping, $\Gamma_{\alpha}\propto\Omega^{3/2}$, has been
reported by experiment [21] and predicted by theory [77].

At low $\Omega$ below the boson peak, we observe a macroscopic sound speed and
Rayleigh scattering in the attenuation behaviour. These phonon transport
characteristics can be understood by considering the glass as an elastic
medium with point defects. This analysis indicates that the disordered
structure is not uniformly coarse-grained even at macroscopic length scales,
but rather, it still plays a defect-type role. On the other hand, at shorter
length scales, phonons are more strongly influenced by the disordered
structure, which induces sound softening and an $\Omega^{2}$ dependence of
damping. $\Gamma_{\alpha}\propto\Omega^{2}$ exhibits dynamics characteristic
of viscous damping [37]. In this regime, the scattering is so strong that a
phonon does not propagate as a plane wave but rather immediately attenuates to
become diffusive. Refs. 78, 79 refer to this vibrational behaviour as
diffusion (nonpropagating, delocalized vibration). We therefore expect the BP
frequency $\omega_{\text{BP}}$ to represent the upper bound on the frequency
at which phonons can propagate as plane waves, as will be discussed next.

(a) Glass, $\alpha=T,\ L$.

(b) Crystal, $\alpha=T_{2},\ L$ in $[110]$.

Figure 10:  Ioffe-Regel (IR) frequency for an SC system. (a) Glass,
$\alpha=T,\ L$. (b) FCC crystal, $\alpha=T_{2},\ L$ in $[110]$. We plot
$\pi\Gamma_{\alpha}/\Omega$ versus $\Omega$. The data for $\Omega$ and
$\Gamma_{\alpha}$ are taken from Figs. 8 and 9. The point of intersection with
$\pi\Gamma_{\alpha}/\Omega=1$ as indicated by the dashed line gives the value
of $\Omega_{\alpha\text{IR}}$. See also the caption of Fig. 8.

####  3 Ioffe-Regel (IR) frequency

We introduce the Ioffe-Regel (IR) frequency, $\Omega_{\alpha\text{IR}}$, which
corresponds to an upper bound on the frequency at which a phonon can propagate
as a plane wave [37, 38, 41]. For $\alpha$ phonons, $\Omega_{\alpha\text{IR}}$
is defined as
$\pi\Gamma_{\alpha}(\Omega=\Omega_{\alpha\text{IR}})/\Omega_{\alpha\text{IR}}=1$.
Above $\Omega_{\alpha\text{IR}}$, the phonon decay time
($=\Gamma_{\alpha}^{-1}$) becomes shorter than half of the vibrational period
($=\pi/\Omega$), i.e., the phonon decays within half of the duration of one
period. Figure 10 plots $\pi\Gamma_{\alpha}(\Omega)/\Omega$ as a function of
$\Omega$. The point of intersection with $\pi\Gamma_{\alpha}(\Omega)/\Omega=1$
gives the value of $\Omega_{\alpha\text{IR}}$. Note that the fitting functions
in Eqs. (51) and (52) (the damped harmonic oscillator model) may not be
appropriate for measuring phonon transport properties above
$\Omega_{\alpha\text{IR}}$ [60].

For the crystal, $\pi\Gamma_{\alpha}(\Omega)/\Omega\ll 1$ in the whole
$\Omega$ region, which means that phonons propagate as plane waves. In
contrast, the glass shows a finite value of $\Omega_{\alpha\text{IR}}$. In
particular, the $\Omega_{T\text{IR}}$ value for transverse phonons nearly
coincides with the BP frequency,
$\Omega_{T\text{IR}}\approx\omega_{\text{BP}}$. This coincidence has been
observed in many glasses [37, 38, 41, 45]. Moreover, this is consistent with
the expectation that $\omega_{\text{BP}}$ should serve as the upper bound on
the frequency at which phonons can propagate as plane waves.

However, the $\Omega_{L\text{IR}}$ for longitudinal phonons in the glass is
located at a much higher frequency,
$\Omega_{L\text{IR}}\gg\Omega_{T\text{IR}}\approx\omega_{\text{BP}}$ [37, 38,
41, 45], meaning that longitudinal phonons can propagate even above
$\omega_{\text{BP}}$. This result originates from the fact that the
attenuation of longitudinal phonons ($\Gamma_{L}(\Omega)$) is much lower than
that of transverse phonons ($\Gamma_{T}(\Omega)$) at a fixed $\Omega$, as
shown in Fig. 9(a), which can, in turn, be understood in terms of
heterogeneous elasticity theory [9, 10]. This theory assumes that shear
modulus heterogeneity dominates compared with bulk modulus heterogeneity,
which is true in the present SC system [40, 41, 42], and predicts that the
shear modulus heterogeneity induces anomalous behaviours for both transverse
and longitudinal phonons. In this theoretical framework, we can understand
that $\Gamma_{T}(\Omega)$ becomes larger than $\Gamma_{L}(\Omega)$ because
transverse phonons are more sensitive to the shear modulus heterogeneity.

####  4 Zero-temperature measurement

As we have discussed so far, the disordered structure has a significant impact
on the phonon transport in a glass. To study these structure-induced effects,
it is most straightforward to measure phonon transport at zero temperature
$T=0$, without thermal fluctuations. We will present two methods of doing so
below.

##### Dynamic structure factor in the zero-temperature limit

We first take the zero-temperature limit (as $T\rightarrow 0$) of the dynamic
structure factor $S_{\alpha}(\mathbf{q},\omega)$. As $T\rightarrow 0$ in the
harmonic limit, the velocity vector $\mathbf{v}(t)$ converges as follows:

$\mathbf{v}(t)=\frac{1}{\sqrt{\mathcal{M}}}\frac{d\mathbf{u}(t)}{dt}\
\longrightarrow\
\sum_{k=1}^{3N}\frac{du^{k}(t)}{dt}\left(\frac{\mathbf{e}^{k}}{\sqrt{\mathcal{M}}}\right),$
(53)

where we apply Eq. (25) for $\mathbf{u}(t)$ and

$\frac{du^{k}(t)}{dt}=-\left(\mathbf{e}^{k}\cdot\mathbf{u}_{0}\right)\omega^{k}\sin\left(\omega^{k}t\right)+\left(\mathbf{e}^{k}\cdot\dot{\mathbf{u}}_{0}\right)\cos\left(\omega^{k}t\right).$
(54)

Eqs. (47) and (50) then lead to

$S_{\alpha}(\mathbf{q},\omega)\ \longrightarrow\
\frac{T}{2}\left(\frac{q}{\omega}\right)^{2}\sum_{k=1}^{3N}\left|\left(\frac{\mathbf{e}^{k}}{\sqrt{\mathcal{M}}}\right)\cdot\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}\right|^{2}\left[\delta\left(\omega-\omega^{k}\right)+\delta\left(\omega+\omega^{k}\right)\right],$
(55)

where we use
$\left<\left(\mathbf{e}^{k}\cdot\dot{\mathbf{u}}_{0}\right)^{2}\right>_{0}=T$,
which is the equipartition law for energy. By taking the polarization vector
to be $\mathbf{s}_{T}(\mathbf{q})=\hat{\mathbf{q}}^{\perp}$ or
$\mathbf{s}_{L}(\mathbf{q})=\hat{\mathbf{q}}$, we can also write Eq. (55) in a
more familiar form:

$S_{\alpha}(\mathbf{q},\omega)\ \longrightarrow\
\frac{T}{2N}\left(\frac{q}{\omega}\right)^{2}\sum_{k=1}^{3N}F_{\alpha}^{k}(\mathbf{q})\left[\delta\left(\omega-\omega^{k}\right)+\delta\left(\omega+\omega^{k}\right)\right],$
(56) $\displaystyle F_{T}^{k}(\mathbf{q})$
$\displaystyle=\left|\sum_{i=1}^{N}\left(\frac{\mathbf{e}^{k}_{i}}{\sqrt{m_{i}}}\times\hat{\mathbf{q}}\right)\exp(\text{i}\mathbf{q}\cdot\mathbf{R}_{i})\right|^{2},$
(57) $\displaystyle F_{L}^{k}(\mathbf{q})$
$\displaystyle=\left|\sum_{i=1}^{N}\left(\frac{\mathbf{e}^{k}_{i}}{\sqrt{m_{i}}}\cdot\hat{\mathbf{q}}\right)\exp(\text{i}\mathbf{q}\cdot\mathbf{R}_{i})\right|^{2}.$

For the case of single-component crystals with mass $m_{i}\equiv m$,
$\mathbf{e}^{k}=\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}=\mathbf{e}^{\mathbf{q},\alpha}_{\text{ph}}$,
and we obtain

$S_{\alpha}(\mathbf{q},\omega)\ \longrightarrow\
\frac{T}{2m}\left(\frac{q}{\omega}\right)^{2}\left[\delta\left(\omega-\omega_{\alpha}(\mathbf{q})\right)+\delta\left(\omega+\omega_{\alpha}(\mathbf{q})\right)\right].$
(58)

We thus confirm that $S_{\alpha}(\mathbf{q},\omega)$ is the delta function as
$T\rightarrow 0$; that the propagation frequency is the eigenfrequency,
$\Omega(\mathbf{q})=\omega_{\alpha}(\mathbf{q})$; and that the attenuation
rate is exactly zero, $\Gamma_{\alpha}(\mathbf{q})=0$. For the case of a
glass, it has been numerically verified that $S_{\alpha}(\mathbf{q},\omega)$
at low $T$ coincides with the value at $T\rightarrow 0$ [30].

##### Direct measurement at zero temperature

The second method is to directly measure phonon transport at $T=0$ [46, 57].
We start with the inherent structure
$\mathbf{R}\equiv\left[\mathbf{R}_{1}^{T},\mathbf{R}_{2}^{T},\cdots,\mathbf{R}_{N}^{T}\right]^{T}$
and excite a phonon of wavevector $\mathbf{q}$ and polarization $\alpha$ by
setting the initial conditions at $t=0$ to
${\mathbf{u}}\equiv\mathbf{u}_{0}=\mathbf{0}$ and
$d{\mathbf{u}}/dt\equiv\dot{\mathbf{u}}_{0}=\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}$
in Eq. (44). We next solve the linearized equation of motion, Eq. (14):

$\frac{d^{2}\mathbf{u}(t)}{dt^{2}}=-\mathcal{D}\mathbf{u}(t)+\dot{\mathbf{u}}_{0}\delta(t).$
(59)

From $\mathbf{u}(t)$, we calculate the (normalized) velocity-velocity time
correlation function:

$C(t)\equiv\left(\frac{\dot{\mathbf{u}}(t)}{\sqrt{\mathcal{M}}}\cdot\frac{\dot{\mathbf{u}}_{0}}{\sqrt{\mathcal{M}}}\right)\left(\frac{\dot{\mathbf{u}}_{0}}{\sqrt{\mathcal{M}}}\cdot\frac{\dot{\mathbf{u}}_{0}}{\sqrt{\mathcal{M}}}\right)^{-1}.$
(60)

The function $C(t)$ represents the propagation and attenuation behaviours of
the initially excited phonon $\dot{\mathbf{u}}_{0}$ with $\mathbf{q}$ and
$\alpha$. $\Omega(\mathbf{q})$ and $\Gamma_{\alpha}(\mathbf{q})$ are then
extracted by fitting the simulated data for $C(t)$ to the damped harmonic
oscillator model in Eq. (52) [46, 57]:
$C(t)\equiv\cos\left(\Omega(\mathbf{q})t\right)\exp\left(-{\Gamma_{\alpha}(\mathbf{q})}t/{2}\right)$.

###  5 Elastic deformation

In this section, we describe how to measure the elastic response, particularly
the elastic moduli, of glasses and any solid-state materials. The elastic
moduli are necessary for calculating the Debye vDOS and the macroscopic sound
speed, which are useful references for understanding the vibrational
properties of a material. There are two methods of measuring the elastic
moduli. One is the so-called fluctuation formulation. In this method, we do
not apply any external deformation but rather perform an equilibrium,
molecular simulation and use the formulation developed based on linear
response theory. The other method is to measure the elastic moduli directly.
We apply an external strain to the system and measure the stress as a function
of the strain, i.e., the stress-strain curve. The slope of this curve gives
the value of the elastic modulus.

####  1 General description

We first provide a general explanation of the elastic response and elastic
moduli.

##### Elastic modulus tensor

The free energy is given by [63]

$F=-T\ln
Z=-T\ln\left\\{\text{Tr}\exp\left[-\frac{1}{T}\left(\sum_{i=1}^{N}\frac{\mathbf{p}_{i}^{2}}{2m_{i}}+\Phi(\mathbf{r})\right)\right]\right\\},$
(61)

where we set the Boltzmann constant to $k_{B}=1$ and
$Z\equiv\text{Tr}\exp\left(-{H}/{T}\right)$ is the partition function
($H\equiv\sum_{i=1}^{N}{\mathbf{p}_{i}^{2}}/{m_{i}}+\Phi$ is the Hamiltonian,
and Tr denotes the trace operator). Here, we consider the Helmholtz free
energy and the canonical ensemble; however, the values of the elastic moduli
are insensitive to the ensemble at low $T$ 202020 The isothermal moduli (in
canonical ensemble) and the adiabatic moduli (in micro-canonical ensemble)
coincide in the low-$T$ harmonic limit [6, 7]. . The stress tensor
$\sigma_{\alpha\beta}$ is defined as the first derivative of the free energy
$F$ with respect to the strain tensor $\epsilon_{\alpha\beta}$ as [80, 81, 82,
83, 84]

$\sigma_{\alpha\beta}\equiv\frac{1}{V}\left.\frac{\partial
F}{\partial\epsilon_{\alpha\beta}}\right|_{\epsilon_{\alpha\beta}\to 0},$ (62)

where $\alpha,\beta~{}(\gamma,\delta)=x,y,z$. There are two different
definitions of $\epsilon_{\alpha\beta}$ 212121 The stress tensors defined by
$\epsilon_{\alpha\beta}=e_{\alpha\beta}$ and $\eta_{\alpha\beta}$ in Eq. (62)
coincide with each other. : the linear strain (infinitesimal strain) tensor,
$e_{\alpha\beta}$, and the Green-Lagrange strain (finite strain) tensor,
$\eta_{\alpha\beta}$:

$\displaystyle e_{\alpha\beta}$
$\displaystyle\equiv\frac{1}{2}\left(\frac{\partial u_{\alpha}}{\partial
r_{\beta}}+\frac{\partial u_{\beta}}{\partial
r_{\alpha}}\right)\equiv\frac{1}{2}\left(u_{\alpha\beta}+u_{\beta\alpha}\right),$
(63) $\displaystyle\eta_{\alpha\beta}$
$\displaystyle\equiv\frac{1}{2}\left(\frac{\partial u_{\alpha}}{\partial
r_{\beta}}+\frac{\partial u_{\beta}}{\partial
r_{\alpha}}+\sum_{\gamma=x,y,z}\frac{\partial u_{\gamma}}{\partial
r_{\alpha}}\frac{\partial u_{\gamma}}{\partial
r_{\beta}}\right)\equiv\frac{1}{2}\left(u_{\alpha\beta}+u_{\beta\alpha}+\sum_{\gamma=x,y,z}u_{\gamma\alpha}u_{\gamma\beta}\right),$

where $[r_{x},r_{y},r_{z}]^{T}$ and $[u_{x},u_{y},u_{z}]^{T}$ represent the
spatial coordinates and the displacement field, respectively, and
$u_{\alpha\beta}\equiv{\partial u_{\alpha}}/{\partial r_{\beta}}$ is the
displacement gradient tensor.

We then define the elastic modulus tensor in two different ways [80, 81, 82,
83, 84]. Firstly the modulus tensor $\tilde{C}_{\alpha\beta\gamma\delta}$ is
defined as the second derivative of $F$ with respect to the Green-Lagrange
strain tensor $\eta_{\alpha\beta}$:

$\tilde{C}_{\alpha\beta\gamma\delta}\equiv\frac{1}{V}\left.\frac{\partial^{2}F}{\partial\eta_{\alpha\beta}\partial\eta_{\gamma\delta}}\right|_{\eta_{\alpha\beta}\to
0,\ \eta_{\gamma\delta}\to 0}.$ (64)

Secondly the modulus tensor ${C}_{\alpha\beta\gamma\delta}$ is defined as the
first derivative of $\sigma_{\alpha\beta}$ with respect to the linear strain
tensor $e_{\gamma\delta}$:

${C}_{\alpha\beta\gamma\delta}\equiv\left.\frac{\partial\sigma_{\alpha\beta}}{\partial
e_{\gamma\delta}}\right|_{e_{\gamma\delta}\to 0}.$ (65)

These two modulus tensors do not coincide for the system under initial stress
$\sigma^{0}_{\alpha\beta}\equiv\sigma_{\alpha\beta}(e_{\gamma\delta}=0)$
222222 $\tilde{C}_{\alpha\beta\gamma\delta}$ is defined in a Lagrangian
framework, whereas ${C}_{\alpha\beta\gamma\delta}$ is based on an Eulerian
framework [84, 81]. . $\tilde{C}_{\alpha\beta\gamma\delta}$ and
${C}_{\alpha\beta\gamma\delta}$ are related as follows (please see Refs. 84,
81 for details):

$\displaystyle C_{\alpha\beta\gamma\delta}$
$\displaystyle=\tilde{C}_{\alpha\beta\gamma\delta}+C^{C}_{\alpha\beta\gamma\delta},$
(66) $\displaystyle C^{C}_{\alpha\beta\gamma\delta}$
$\displaystyle=-\frac{1}{2}\left(2\sigma^{0}_{\alpha\beta}\delta_{\gamma\delta}-\sigma^{0}_{\alpha\gamma}\delta_{\beta\delta}-\sigma^{0}_{\alpha\delta}\delta_{\beta\gamma}-\sigma^{0}_{\beta\gamma}\delta_{\alpha\delta}-\sigma^{0}_{\beta\delta}\delta_{\alpha\gamma}\right).$

We refer to $C^{C}_{\alpha\beta\gamma\delta}$ as the correction term in the
following. If the initial stress tensor $\sigma^{0}_{\alpha\beta}$ is zero,
then $C^{C}_{\alpha\beta\gamma\delta}$ is zero, and
$C_{\alpha\beta\gamma\delta}$ and $\tilde{C}_{\alpha\beta\gamma\delta}$
exactly coincide. Since $C_{\alpha\beta\gamma\delta}$ is necessary for
calculations of the Debye vDOS and the macroscopic sound speed, we focus on
$C_{\alpha\beta\gamma\delta}$ below.

Here we note that strictly speaking, for calculations of the Debye vDOS and
the macroscopic sound speed, we should apply the modulus tensor
${C}^{\prime}_{\alpha\beta\gamma\delta}$ which is defined as the first
derivative of $\sigma_{\alpha\beta}$ with respect to the displacement gradient
tensor $u_{\gamma\delta}$ (instead of $e_{\gamma\delta}$):

${C}^{\prime}_{\alpha\beta\gamma\delta}\equiv\left.\frac{\partial\sigma_{\alpha\beta}}{\partial
u_{\gamma\delta}}\right|_{u_{\gamma\delta}\to 0}.$ (67)

However, as long as we consider the situation under hydrostatic pressure where
the initial stress is $\sigma^{0}_{\alpha\beta}=-p\delta_{\alpha\beta}$,
${C}_{\alpha\beta\gamma\delta}$ and ${C}^{\prime}_{\alpha\beta\gamma\delta}$
are identical [84, 81].

Figure 11:  Schematic illustrations of (a) bulk, (b) pure shear, and (c)
simple shear deformations. The bulk ($K$), pure shear ($G_{p}$), and simple
shear ($G_{s}$) moduli correspond to these deformations.

##### Bulk, pure shear, and simple shear elastic moduli

We consider three types of deformations in particular, namely, a volume-
changing bulk deformation and pure shear (plane or triaxial strains) and
simple shear deformations, which are illustrated in Fig. 11. Since there are
two independent pure shear deformations and three independent simple shear
deformations, we consider six independent deformations in total, including one
bulk deformation [85]. The strain tensors for these six deformations are
232323 We can also consider other pure shear deformations, which are described
as superpositions of the present two independent pure shear deformations.

$\epsilon_{xx}=\epsilon_{yy}=\epsilon_{zz}=\epsilon_{b}\quad\text{(bulk)},$
(68)
$\displaystyle\epsilon_{xx}=\epsilon_{p1},\quad\epsilon_{yy}=-\epsilon_{p1}$
$\displaystyle\text{(pure shear1)},$ (69)
$\displaystyle\epsilon_{xx}=\epsilon_{p2},\quad\epsilon_{yy}=\epsilon_{p2},\quad\epsilon_{zz}=-2\epsilon_{p2}$
$\displaystyle\text{(pure shear2)},$
$\displaystyle\epsilon_{xy}=\epsilon_{yx}=\epsilon_{s1}$
$\displaystyle\text{(simple shear1)},$ (70)
$\displaystyle\epsilon_{xz}=\epsilon_{zx}=\epsilon_{s2}$
$\displaystyle\text{(simple shear2)},$
$\displaystyle\epsilon_{yz}=\epsilon_{zy}=\epsilon_{s3}$
$\displaystyle\text{(simple shear3)},$

where no written components in $\epsilon_{\alpha\beta}$ are zero, and
$\epsilon_{b}$, $\epsilon_{p1}$, $\epsilon_{p2}$, $\epsilon_{s1}$,
$\epsilon_{s2}$, $\epsilon_{s3}$ represent the applied strains corresponding
to the deformations. $\delta V/V\equiv 3\epsilon_{b}$ represents the volume
change, while $\gamma_{p1}\equiv 2\epsilon_{p1}$, $\gamma_{p2}\equiv
3\epsilon_{p2}$, $\gamma_{s1}\equiv 2\epsilon_{s1}$, $\gamma_{s2}\equiv
2\epsilon_{s2}$, $\gamma_{s3}\equiv 2\epsilon_{s3}$ present the shear strains
for the corresponding shear deformations.

When we apply a deformation to the system, the corresponding stress varies
with the applied strain. For the bulk deformation, the two pure shear
deformations, and the three simple shear deformations, the applied stresses
are the pressure, $p\equiv-(\sigma_{xx}+\sigma_{yy}+\sigma_{zz})/3$; two pure
shear stresses, $\sigma_{p1}\equiv{(\sigma_{xx}-\sigma_{yy})}/{2}$ and
$\sigma_{p2}\equiv{(\sigma_{xx}+\sigma_{yy}-2\sigma_{zz})}/{4}$; and three
simple shear stresses, $\sigma_{s1}\equiv\sigma_{xy}$,
$\sigma_{s2}\equiv\sigma_{xz}$, and $\sigma_{s3}\equiv\sigma_{yz}$,
respectively. Accordingly, the bulk modulus $K$, the two pure shear moduli
$G_{p1}$ and $G_{p2}$, and the three simple shear moduli $G_{s1}$, $G_{s2}$,
and $G_{s3}$ are defined as follows:

$K\equiv\left.-\frac{\partial p}{\partial(\delta V/V)}\right|_{(\delta
V/V)\rightarrow
0}=\left.\frac{\partial(\sigma_{xx}+\sigma_{yy}+\sigma_{zz})/3}{\partial(3\epsilon_{b})}\right|_{\epsilon_{b}\rightarrow
0}\quad\text{(bulk)},$ (71) $\displaystyle
G_{p1}\equiv\left.\frac{\partial\sigma_{p1}}{\partial\gamma_{p1}}\right|_{\gamma_{p1}\rightarrow
0}=\left.\frac{\partial(\sigma_{xx}-\sigma_{yy})/2}{\partial(2\epsilon_{p1})}\right|_{\epsilon_{p1}\rightarrow
0}$ $\displaystyle\text{(pure shear 1)},$ (72) $\displaystyle
G_{p2}\equiv\left.\frac{\partial\sigma_{p2}}{\partial\gamma_{p2}}\right|_{\gamma_{p2}\rightarrow
0}=\left.\frac{\partial(\sigma_{xx}+\sigma_{yy}-2\sigma_{zz})/4}{\partial(3\epsilon_{p2})}\right|_{\epsilon_{p2}\rightarrow
0}$ $\displaystyle\text{(pure shear 2)},$ $\displaystyle
G_{s1}\equiv\left.\frac{\partial\sigma_{s1}}{\partial\gamma_{s1}}\right|_{\gamma_{s1}\rightarrow
0}=\left.\frac{\partial\sigma_{xy}}{\partial(2\epsilon_{s1})}\right|_{\epsilon_{s1}\rightarrow
0}$ $\displaystyle\text{(simple shear 1)},$ (73) $\displaystyle
G_{s2}\equiv\left.\frac{\partial\sigma_{s2}}{\partial\gamma_{s2}}\right|_{\gamma_{s2}\rightarrow
0}=\left.\frac{\partial\sigma_{xz}}{\partial(2\epsilon_{s2})}\right|_{\epsilon_{s2}\rightarrow
0}$ $\displaystyle\text{(simple shear 2)},$ $\displaystyle
G_{s3}\equiv\left.\frac{\partial\sigma_{s3}}{\partial\gamma_{s3}}\right|_{\gamma_{s3}\rightarrow
0}=\left.\frac{\partial\sigma_{yz}}{\partial(2\epsilon_{s3})}\right|_{\epsilon_{s3}\rightarrow
0}$ $\displaystyle\text{(simple shear 3)}.$

$K$, $G_{p1}$, $G_{p2}$, $G_{s1}$, $G_{s2}$, and $G_{s3}$ can be formulated in
terms of the elastic modulus tensor $C_{\alpha\beta\gamma\delta}$ [85]:

$K=\frac{(C_{xxxx}+C_{yyyy}+C_{zzzz}+C_{xxyy}+C_{yyxx}+C_{xxzz}+C_{zzxx}+C_{yyzz}+C_{zzyy})}{9},$
(74) $\displaystyle G_{p1}$
$\displaystyle=\frac{(C_{xxxx}+C_{yyyy}-C_{xxyy}-C_{yyxx})}{4},$ (75)
$\displaystyle G_{p2}$
$\displaystyle=\frac{(C_{xxxx}+C_{yyyy}+4C_{zzzz}+C_{xxyy}+C_{yyxx}-2C_{xxzz}-2C_{zzxx}-2C_{yyzz}-2C_{zzyy})}{12},$
$G_{s1}={C_{xyxy}},\qquad G_{s2}={C_{xzxz}},\qquad G_{s3}={C_{yzyz}}.$ (76)

Note that $C_{xyxy}=C_{xyyx}$, $C_{xzxz}=C_{xzzx}$, $C_{yzyz}=C_{yzzy}$. For
cubic crystals such as FCC crystals, the two pure shear moduli coincide, i.e.,
$G_{p}\equiv G_{p1}=G_{p2}$, and the three simple shear moduli also coincide,
i.e., $G_{s}\equiv G_{s1}=G_{s2}=G_{s3}$; however, the pure shear and simple
shear moduli generally take different values, i.e., $G_{p}\neq G_{s}$. On the
other hand, for glasses, which are isotropic systems, all five shear moduli
coincide: $G\equiv G_{p1}=G_{p2}=G_{s1}=G_{s2}=G_{s3}$.

##### Affine and non-affine elastic moduli

The elastic response can be decomposed into affine and non-affine components
[31, 84, 86, 87]. Affine deformation refers to particles that follow the
applied affine strain field and are displaced affinely at all scales. More
specifically, for a linear strain $\epsilon_{\alpha\beta}=e_{\alpha\beta}$,
the affine deformation causes the displacement of the position of particle
$i$, $\mathbf{r}_{i}=[r_{ix},r_{iy},r_{iz}]^{T}$, as follows:

$r_{i\kappa}\rightarrow
r_{i\kappa}+\frac{1}{2}\left(u_{\alpha\beta}r_{i\beta}\delta_{\alpha\kappa}+u_{\beta\alpha}r_{i\alpha}\delta_{\beta\kappa}\right),$
(77)

where $\kappa=x,y,z$. Similarly, for the Green-Lagrange strain
$\epsilon_{\alpha\beta}=\eta_{\alpha\beta}$, particles are displaced as
follows:

$r_{i\kappa}\rightarrow
r_{i\kappa}+\frac{1}{2}\left(u_{\alpha\beta}r_{i\beta}\delta_{\alpha\kappa}+u_{\beta\alpha}r_{i\alpha}\delta_{\beta\kappa}+\sum_{\gamma=x,y,z}u_{\gamma\alpha}u_{\gamma\beta}r_{i\alpha}r_{i\beta}\delta_{\gamma\kappa}\right).$
(78)

The affine strain tensors $\epsilon^{A}_{\alpha\beta}=e^{A}_{\alpha\beta}$ and
$\eta^{A}_{\alpha\beta}$ are defined as displacing particles as expressed in
Eqs. (77) and (78).

By applying the affine strain tensors to Eqs. (64) and (65), we can define the
affine elastic modulus tensors as

$\tilde{C}^{A}_{\alpha\beta\gamma\delta}\equiv\frac{1}{V}\left.\frac{\partial^{2}F}{\partial\eta^{A}_{\alpha\beta}\partial\eta^{A}_{\gamma\delta}}\right|_{\eta^{A}_{\alpha\beta}\to
0,\ \eta^{A}_{\gamma\delta}\to 0},$ (79)
${C}^{A}_{\alpha\beta\gamma\delta}\equiv\left.\frac{\partial\sigma_{\alpha\beta}}{\partial
e^{A}_{\gamma\delta}}\right|_{e^{A}_{\gamma\delta}\to 0}.$ (80)

The elastic response of a crystal is characterized mostly by affine
deformation and the corresponding affine modulus [80, 82]. However,
importantly, the elastic response of a glass does not follow solely affine
deformation. We additionally need to consider non-affine deformation, which
results in additional particle displacements at the microscopic scale and
causes the displacements to deviate from the applied affine field [31, 84, 86,
87]. Non-affine deformation contributes negatively to the overall modulus;
thus, the corresponding modulus is defined as the non-affine modulus
$-C_{\alpha\beta\gamma\delta}^{N}$. We therefore describe the elastic modulus
tensor as

${C}_{\alpha\beta\gamma\delta}={C}^{A}_{\alpha\beta\gamma\delta}-C_{\alpha\beta\gamma\delta}^{N}.$
(81)

As we will see below, the non-affine component is important in glasses, being
comparable in magnitude to the affine modulus [31, 84, 86, 87].

####  2 Fluctuation formulation

We now present how to measure the elastic moduli by means of computational
simulations. We first introduce the fluctuation formulation. In the following,
we consider a system with a pair-wise potential, where particles $i$ and $j$
interact through the potential $\phi_{ij}(r_{ij})$ ($r_{ij}$ is the distance
between them) and the total potential is $\Phi=\sum_{i<j}\phi_{ij}(r_{ij})$.
We also employ the linear strain tensor
$\epsilon_{\alpha\beta}=e_{\alpha\beta}$. Please see, e.g., Refs. 83, 84, 86,
for detailed derivations of the following formulations.

##### Finite temperature $T>0$

By using Eqs. (61) and (62), we can obtain the stress tensor in the form of

$\displaystyle\sigma_{\alpha\beta}$
$\displaystyle=\left<\hat{\sigma}_{\alpha\beta}\right>,$ (82)
$\displaystyle\hat{\sigma}_{\alpha\beta}$
$\displaystyle=\frac{1}{V}\left[-\sum_{i=1}^{N}m_{i}v_{i\alpha}v_{i\beta}+\sum_{i<j}\left(r_{ij}{\displaystyle\frac{\partial\phi_{ij}}{\partial
r_{ij}}}\right){n_{ij\alpha}n_{ij\beta}}\right],$

where $\left<\right>$ denotes the canonical ensemble average,
$\hat{\rho}={N}/{V}$ is the number density, and
$\mathbf{n}_{ij}=\left[n_{ijx},n_{ijy},n_{ijz}\right]^{T}\equiv(\mathbf{r}_{i}-\mathbf{r}_{j})/r_{ij}$
is the unit vector connecting particle $j$ to particle $i$. The elastic
modulus tensor $C_{\alpha\beta\gamma\delta}$ is formulated using Eqs. (64) and
(66), which includes the affine modulus and the non-affine modulus as in Eq.
(81). The affine modulus tensor $C^{A}_{\alpha\beta\gamma\delta}$ [Eq. (80)]
is formulated as

$\displaystyle C^{A}_{\alpha\beta\gamma\delta}$
$\displaystyle=C^{B}_{\alpha\beta\gamma\delta}+C^{K}_{\alpha\beta\gamma\delta}+C^{C}_{\alpha\beta\gamma\delta},$
(83) $\displaystyle C^{B}_{\alpha\beta\gamma\delta}$
$\displaystyle=\frac{1}{V}\left<\sum_{i<j}\left(r_{ij}^{2}{\displaystyle\frac{\partial^{2}\phi_{ij}}{\partial{r_{ij}}^{2}}}-{r_{ij}}{\displaystyle\frac{\partial\phi_{ij}}{\partial
r_{ij}}}\right){n_{ij\alpha}n_{ij\beta}n_{ij\gamma}n_{ij\delta}}\right>,$
$\displaystyle C^{K}_{\alpha\beta\gamma\delta}$
$\displaystyle=2\hat{\rho}T(\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma}),$
$\displaystyle C^{C}_{\alpha\beta\gamma\delta}$
$\displaystyle=-\frac{1}{2}\left(2\left<\hat{\sigma}_{\alpha\beta}\right>\delta_{\gamma\delta}-\left<\hat{\sigma}_{\alpha\gamma}\right>\delta_{\beta\delta}-\left<\hat{\sigma}_{\alpha\delta}\right>\delta_{\beta\gamma}-\left<\hat{\sigma}_{\beta\gamma}\right>\delta_{\alpha\delta}-\left<\hat{\sigma}_{\beta\delta}\right>\delta_{\alpha\gamma}\right),$

where $C^{B}_{\alpha\beta\gamma\delta}$ is the so-called Born term and
$C^{K}_{\alpha\beta\gamma\delta}$ is the kinetic contribution. Since these two
terms are formulated by Eq. (79) based on the Green-Lagrange strain, we need
the correction term $C^{C}_{\alpha\beta\gamma\delta}$ to obtain the modulus
tensor based on the linear strain, as shown in Eq. (66) [84, 81]. In addition,
the non-affine modulus is formulated as

$C^{N}_{\alpha\beta\gamma\delta}=\frac{V}{T}\left[\left<\hat{\sigma}_{\alpha\beta}\hat{\sigma}_{\gamma\delta}\right>-\left<\hat{\sigma}_{\alpha\beta}\right>\left<\hat{\sigma}_{\gamma\delta}\right>\right]=\frac{V}{T}\left<\delta\hat{\sigma}_{\alpha\beta}\delta\hat{\sigma}_{\gamma\delta}\right>,$
(84)

where
$\delta\hat{\sigma}_{\alpha\beta}\equiv\hat{\sigma}_{\alpha\beta}-\left<\hat{\sigma}_{\alpha\beta}\right>$
is the stress fluctuation. $C^{N}_{\alpha\beta\gamma\delta}$ is therefore
formulated as the correlation function for stress fluctuations.

The elastic moduli, $K$, $G_{p}$, and $G_{s}$, of SC systems: a glass and an
FCC crystal. We present values obtained using two methods: the fluctuation
formulation (FF) at a finite temperature of $T=10^{-2}$ (almost at the zero-
temperature limit) and direct measurement (DM) at zero temperature, $T=0$. We
present the average values of pure shear modulus, $G_{p}=(G_{p1}+G_{p2})/2$,
and simple shear modulus, $G_{s}=(G_{s1}+G_{s2}+G_{s3})/3$. The affine and
non-affine moduli are also presented.

Method $K$ $K^{A}$ $K^{N}$ $G_{p}$ $G_{p}^{A}$ $G_{p}^{N}$ $G_{s}$ $G_{s}^{A}$
$G_{s}^{N}$ Glass FF $40.9$ $41.0$ $0.1$ $6.2$ $14.7$ $8.5$ $6.3$ $14.8$ $8.5$
DM $40.5$ $40.5$ $0.0$ $6.3$ $14.4$ $8.1$ $6.4$ $14.7$ $8.3$ Crystal FF $32.9$
$33.0$ $0.1$ $5.1$ $5.2$ $0.1$ $16.1$ $16.3$ $0.2$ DM $32.6$ $32.6$ $0.0$
$5.4$ $5.4$ $0.0$ $16.2$ $16.2$ $0.0$

In our computer simulation, we first perform an equilibrium ($NVT$) MD or MC
simulation at a finite $T>0$ and generate configurations
$\mathbf{r}\equiv\left[\mathbf{r}_{1}^{T},\mathbf{r}_{2}^{T},\cdots,\mathbf{r}_{N}^{T}\right]^{T}$.
From these configurations, we evaluate $C_{\alpha\beta\gamma\delta}$ by using
Eqs. (81) to (84). From $C_{\alpha\beta\gamma\delta}$, we then calculate the
elastic moduli, $K$, $G_{p}$, and $G_{s}$, using Eqs. (74) to (76). Table 2
presents the values of $K$, $G_{p}$, and $G_{s}$ for a glass and an FCC
crystal with the SC potential. Here, we also present the values of the affine
moduli $K^{A}$, $G_{p}^{A}$, and $G_{s}^{A}$ and the non-affine moduli
$K^{N}$, $G_{p}^{N}$, and $G_{s}^{N}$. The temperature is very low,
$T=10^{-2}$, almost at the zero-temperature limit. We confirm that the crystal
shows negligible values of the non-affine components in both the bulk and
shear moduli. The elastic response of the crystal is therefore characterized
by affine deformation, which is due to the symmetry of the lattice structure
[80, 82].

In contrast to the crystal, the glass shows large values of the non-affine
components in the shear moduli, comparable to those of the affine components.
Therefore, non-affine deformation plays an important role in the elastic
response of the glass, making it distinct from the elastic response of the
crystal [31, 84, 86, 87]. Note that $G_{p}$ and $G_{s}$ coincide in the glass
due to its isotropic structure. On the other hand, for bulk deformation, the
glass also shows a very small non-affine component, which is due to the
isotropic nature of bulk deformation. In particular, we can demonstrate that a
mono-disperse system with an inverse-power-law potential, such as the present
SC potential [Eq. (5)], shows a non-affine modulus of zero under bulk
deformation at zero temperature.

##### Zero-temperature limit $T\rightarrow 0$

It is useful to take the limit as $T\rightarrow 0$ in the fluctuation
formulation given in Eqs. (81) to (84) [83, 84]. As $T\rightarrow 0$, the
configuration converges to the inherent structure,
$\mathbf{r}\rightarrow\mathbf{R}=\left[\mathbf{R}_{1}^{T},\mathbf{R}_{2}^{T},\cdots,\mathbf{R}_{N}^{T}\right]^{T}$,
and we obtain $C^{B}_{\alpha\beta\gamma\delta}$,
$C^{K}_{\alpha\beta\gamma\delta}$, and $C^{C}_{\alpha\beta\gamma\delta}$ as
follows:

$\displaystyle C^{B}_{\alpha\beta\gamma\delta}$ $\displaystyle\
\longrightarrow\
\frac{1}{V}\sum_{i<j}\left.\left(r_{ij}^{2}{\displaystyle\frac{\partial^{2}\phi_{ij}}{\partial{r_{ij}}^{2}}}-{r_{ij}}{\displaystyle\frac{\partial\phi_{ij}}{\partial
r_{ij}}}\right){n_{ij\alpha}n_{ij\beta}n_{ij\gamma}n_{ij\delta}}\right|_{\mathbf{r}=\mathbf{R}},$
(85) $\displaystyle C^{K}_{\alpha\beta\gamma\delta}$ $\displaystyle\
\longrightarrow\ 0,$ $\displaystyle C^{C}_{\alpha\beta\gamma\delta}$
$\displaystyle\ \longrightarrow\ -\frac{1}{2}\left(2\sigma_{\alpha\beta
0}\delta_{\gamma\delta}-\sigma_{\alpha\gamma
0}\delta_{\beta\delta}-\sigma_{\alpha\delta
0}\delta_{\beta\gamma}-\sigma_{\beta\gamma
0}\delta_{\alpha\delta}-\sigma_{\beta\delta 0}\delta_{\alpha\gamma}\right),$

where $\sigma_{\alpha\beta 0}$ is the stress tensor $\sigma_{\alpha\beta}$ at
$T\rightarrow 0$,

$\sigma_{\alpha\beta}\ \longrightarrow\ \sigma_{\alpha\beta
0}\equiv\frac{1}{V}\left.\sum_{i<j}\left(r_{ij}{\displaystyle\frac{\partial\phi_{ij}}{\partial
r_{ij}}}\right){n_{ij\alpha}n_{ij\beta}}\right|_{\mathbf{r}=\mathbf{R}}.$ (86)

The non-affine term $C_{\alpha\beta\gamma\delta}^{N}$ is formulated in terms
of the eigenfrequencies $\omega^{k}$ and the eigenvectors $\mathbf{e}^{k}$ as
follows. The stress fluctuation is

$\displaystyle\delta\hat{\sigma}_{\alpha\beta}\ \longrightarrow\
{\displaystyle\frac{\partial\sigma_{\alpha\beta
0}}{\partial\mathbf{R}}}\cdot\frac{\mathbf{u}(t)}{\sqrt{\mathcal{M}}}=\sum_{k=1}^{3N}u^{k}(t)\left({\displaystyle\frac{\partial\sigma_{\alpha\beta
0}}{\partial\mathbf{R}}}\cdot\frac{\mathbf{e}^{k}}{\sqrt{\mathcal{M}}}\right),$
(87)

where we use $\mathbf{u}(t)$ in Eq. (25). We then obtain

$C_{\alpha\beta\gamma\delta}^{N}\ \longrightarrow\
V\sum_{k=1}^{3N}\frac{1}{{\omega^{k}}^{2}}\left({\displaystyle\frac{\partial\sigma_{\alpha\beta
0}}{\partial\mathbf{R}}}\cdot\frac{\mathbf{e}^{k}}{\sqrt{\mathcal{M}}}\right)\left({\displaystyle\frac{\partial\sigma_{\gamma\delta
0}}{\partial\mathbf{R}}}\cdot\frac{\mathbf{e}^{k}}{\sqrt{\mathcal{M}}}\right),$
(88)

where we use
$\left<\left(\mathbf{e}^{k}\cdot\dot{\mathbf{u}}_{0}\right)^{2}\right>={\omega^{k}}^{2}\left<\left(\mathbf{e}^{k}\cdot{\mathbf{u}}_{0}\right)^{2}\right>=T$,
which is the equipartition law for energy. It is worth noting that the affine
component, $C^{A}_{\alpha\beta\gamma\delta}\equiv
C^{B}_{\alpha\beta\gamma\delta}+C^{K}_{\alpha\beta\gamma\delta}+C^{C}_{\alpha\beta\gamma\delta}$,
is determined by static structural properties, whereas the non-affine term,
$C^{N}_{\alpha\beta\gamma\delta}$, reflects the vibrational properties. We
also note that the non-affine term $C_{\alpha\beta\gamma\delta}^{N}$ in Eq.
(88) can be written as

$C_{\alpha\beta\gamma\delta}^{N}\ \longrightarrow\
V{\displaystyle\frac{\partial\sigma_{\alpha\beta
0}}{\partial\mathbf{R}^{T}}}\left[\left.\frac{\partial^{2}\Phi}{\partial\mathbf{r}\partial\mathbf{r}^{T}}\right|_{\mathbf{r}=\mathbf{R}}\right]^{-1}{\displaystyle\frac{\partial\sigma_{\alpha\beta
0}}{\partial\mathbf{R}}}.$ (89)

From Eqs. (85), (86), and (89), we can confirm that elastic modulus at zero
temperature does not depend on mass of particles but rather it depends on only
the potential.

(a) Glass (bulk).

(b) Crystal (bulk).

(c) Glass (pure shear$1$).

(d) Crystal (pure shear$1$).

(e) Glass (simple shear$1$).

(f) Crystal (simple shear$1$).

Figure 12:  The elastic response of an SC system. (a),(c),(e) Glass.
(b),(d),(f) FCC crystal. We plot $p$ versus $\delta V/V$ for bulk deformation
in (a),(b), $\sigma_{p1}$ versus $\gamma_{p1}$ for pure shear $1$ deformation
in (c),(d), and $\sigma_{s1}$ versus $\gamma_{s1}$ for simple shear $1$
deformation in (e),(f). The dashed line represents affine deformation, while
the solid line represents the overall response, including non-affine
deformation. The slopes of the stress-strain curves yield the values of the
elastic moduli, $K$, $G_{p1}$, and $G_{s1}$, as in Eqs. (71) to (73). The
obtained values are presented in Table 2.

####  3 Direct measurement

We can apply an external strain to the system to directly measure the elastic
moduli. By increasing the strain, we obtain the stress-strain curve, and the
slope of this curve yields the corresponding elastic modulus. Specifically,
the bulk, pure shear $1$, and simple shear $1$ deformations, which are
illustrated in Fig. 11, yield the $p-\delta V/V$, $\sigma_{p1}-\gamma_{p1}$,
and $\sigma_{s1}-\gamma_{s1}$ curves, respectively. Figure 12 (solid lines)
plots these stress-strain curves, which we obtained by applying the linear
strain $\epsilon_{\alpha\beta}=e_{\alpha\beta}$. The slopes of these curves
give the bulk modulus $K$, the pure shear modulus $G_{p1}$, and the simple
shear modulus $G_{s1}$, as shown in Eqs. (71) to (73). We note that other pure
shear, $G_{p2}$, and simple shear, $G_{s2}$ and $G_{s3}$, moduli can be also
obtained in the same manner.

The obtained values of the elastic moduli are presented in Table 2. From Table
2, we can confirm that these values coincide well with those obtained from the
fluctuation formulation. We note that the fluctuation formulation and direct
measurement method were implemented at a finite temperature of $T=10^{-2}$ and
at zero temperature ($T=0$), respectively. However, almost negligible
differences are observed at these different $T$s, since $T=10^{-2}$ is almost
at the zero-temperature limit and the elastic moduli are insensitive to $T$ in
the low-$T$ regime.

To obtain the response under affine deformation, we can apply an affine
strain, $\epsilon^{A}_{\alpha\beta}=e^{A}_{\alpha\beta}$, to displace the
particles according to Eq. (77). After the application of an affine strain,
the particles are generally not in mechanical equilibrium; however, we do not
permit the particles to relax to the equilibrium state. By doing so, we obtain
the stress-strain curve under affine deformation, which is shown by the dashed
lines in Fig. 12. On the other hand, when we allow the particles to relax, we
obtain the overall stress-strain curve, as shown by the solid lines in Fig.
12. This relaxation process, i.e., non-affine deformation, generally reduces
the stress and causes a reduction in the elastic modulus.

As shown in Figs. 12(b), 12(d), and 12(f), for the crystal, the stress-strain
curves of the overall response (solid line) and the affine response (dashed
line) coincide. The elastic response of the crystal is therefore determined by
the affine deformation without relaxation, as discussed above. In contrast to
the crystal, the glass shows a large non-affine contribution in its shear
deformations, as shown in Figs. 12(c) and 12(e). For the bulk deformation of
the glass, the non-affine relaxation is negligible, as shown in Fig. 12(a);
this is due to the isotropic nature of bulk deformation, as also discussed
above.

###  6 Recent advances concerning the vibrational properties of glasses

At the end of this chapter, we will introduce recent advances in the
understanding of the vibrational properties of glasses. Currently, the
capability of computers is rapidly growing. This enables us to perform large-
scale computational simulations and to understand the vibrations of the
particles in glasses in greater detail. In particular, recent simulations [50,
51, 52, 53, 54, 55, 56, 57, 58, 59] have revealed the vibrational properties
in the low-frequency continuum limit. In this section, we present simulation
data for vibrational eigenmodes and phonon transport, giving particular
attention to the low-frequency regime. The data presented below were obtained
in Refs. 51, 57 by simulating an atomic glass (packed glass) with a harmonic
potential, $\Phi_{\text{HA}}$, as expressed in Eq. (6).

As described in Section 3, the eigenmodes in a crystal, i.e., phonons,
smoothly converge to elastic waves at low frequencies, and their vDOS
converges to the Debye vDOS [6, 7]. This means that crystals behave as elastic
media at long length scales, where the microscopic lattice structures are
uniformly coarse-grained. Similarly, we might expect that the disordered
structures of glasses are also uniformly coarse-grained at long length scales
and that they would therefore also behave as elastic media. The eigenmodes and
vDOS of a glass would then be expected to converge to elastic waves and the
Debye vDOS at low frequencies. However, in contrast to this expectation, we
will show below that the disordered structures of glasses are not completely
uniformly coarse-grained but rather act as defects even at macroscopic scales.
Glasses therefore behave not as uniform elastic media but rather as elastic
media with defects.

Figure 13:  Vibrational eigenmodes in a glass with a harmonic potential. (a)
Plots of the reduced vDOS $g(\omega)/\omega^{2}$ and of the phonon order
parameter $O^{k}$ and participation ratio $\mathcal{P}^{k}$ for each eigenmode
$k$ as functions of the frequency $\omega$. (b) The vDOS $g(\omega)$. (c)
$O^{k}$ versus $\mathcal{P}^{k}$ for the low-frequency eigenmodes with
$\omega^{k}<\omega_{\text{ex0}}$. Regarding the vDOS, data are also presented
for $g_{\text{ex}}(\omega)$, the vDOS of the extended modes
($\mathcal{P}^{k}>10^{-2}$), and for $g_{\text{loc}}(\omega)$, the vDOS of the
localized modes ($\mathcal{P}^{k}<10^{-2}$). The horizontal lines at the top
and bottom of (a) represent the Debye level $A_{D}$ and $\mathcal{P}=10^{-2}$,
respectively. In the top panel of (a), we use arrows to indicate the
characteristic frequencies $\omega_{\text{ex0}}$, $\omega_{\text{BP}}$, and
$\omega_{\ast}$, as discussed in the main text.

####  1 Vibrational eigenmodes in glasses

Figure 13 presents the results for the vibrational eigenmodes in a glass with
a harmonic potential. The top panel of Fig. 13(a) (circles) shows the reduced
vDOS $g(\omega)/\omega^{2}$ together with the Debye level $A_{0}$. The reduced
vDOS clearly exhibits a maximum, i.e., the boson peak (the arrow indicates the
BP frequency, $\omega_{\text{BP}}$). As the frequency decreases below
$\omega_{\text{BP}}$, $g(\omega)/\omega^{2}$ decreases toward but does not
reach $A_{0}$ in the present frequency region. We will carefully discuss this
point below.

To enable characterization of the eigenmodes, Figure 13(a) also presents data
on the phonon order parameter $O^{k}$ defined in Eq. (46) (middle panel) and
the participation ratio $\mathcal{P}^{k}$ defined in Eq. (43) (bottom panel).
Let us first consider $O^{k}$, which measures the extent to which the
eigenmode $k$ exhibits phonon-like vibrations [51, 52]: $O^{k}$ takes values
from $1$ (phonon) to $0$ (non-phonon). At high $\omega$, $O^{k}$ is nearly
zero, which confirms that these eigenmodes are considerably different to
phonons. We may define the frequency, $\omega_{\ast}$, at which $O^{k}$
converges to zero 242424 The present glass system shows a characteristic
plateau in the vDOS [35, 34, 36]. The value of $\omega_{\ast}$ is usually
defined as the onset frequency of this plateau, which is plotted in Fig.
13(a). However, we can confirm in Fig. 13(a) that the frequency at which
$O^{k}$ converges to zero is consistent with the onset frequency
$\omega_{\ast}$. . It has been reported that at the high
$\omega>\omega_{\ast}$, the eigenmodes are disordered and extended [35, 34,
36]; this report is consistent with the present result of $O^{k}\approx 0$. As
$\omega$ decreases from $\omega_{\ast}$ to $\omega_{\text{BP}}$, $O^{k}$
smoothly increases to $\approx 0.3$. This result indicates that the eigenmodes
near $\omega_{\text{BP}}$ show phonon-like vibrations to some extent. Notably,
as $\omega$ further decreases below $\omega_{\text{BP}}$, the modes can be
divided into two groups: $O^{k}$ increases with decreasing $\omega$ in one
group, whereas $O^{k}$ decreases in the other group. In the former group,
$O^{k}$ converges to almost $1$ at $\omega_{\text{ex0}}$ (we will provide the
precise definition of $\omega_{\text{ex0}}$ later), indicating that these
modes are phonon modes 252525 The values of $O^{k}$ for these phonon modes are
close to but not exactly $1$, which indicates that they are weakly perturbed.
An exact value of $O^{k}=1$ may be realized only in the limit of $\omega\to
0$. .

We next consider $\mathcal{P}^{k}$, which evaluates the extent of spatial
localization of eigenmode $k$ [26, 27, 28]. $\mathcal{P}^{k}$ takes values
from 1 (extended over all particles equally) to $1/N\ll 1$ (localized to one
particle). $\mathcal{P}^{k}$ more clearly exhibits the division of the modes
into two groups: one group approaches $\mathcal{P}^{k}=\mathcal{O}(1)$ with
decreasing $\omega$, and the other approaches
$\mathcal{P}^{k}=\mathcal{O}(1/N)$. Figure 13(c) shows that the non-phonon
modes (small $O^{k}$) are localized (small $\mathcal{P}^{k}$) at low
$\omega<\omega_{\text{ex0}}$, whereas the phonon modes (large $O^{k}$) are
extended (large $\mathcal{P}^{k}$). Therefore, the data for $O^{k}$ and
$\mathcal{P}^{k}$ unambiguously demonstrate that phonon modes and non-phonon
localized modes coexist at the low $\omega<\omega_{\text{ex0}}$.

The distinction between phonon modes and localized modes enables us to
separately consider the vDOSs for these two types of modes. We define
$g_{\text{ex}}(\omega)$ as the vDOS for modes with
$\mathcal{P}^{k}>\mathcal{P}_{c}$, and we define $g_{\text{loc}}(\omega)$ as
the vDOS for modes with $\mathcal{P}^{k}<\mathcal{P}_{c}$. Here, we set a
reasonable threshold value of $\mathcal{P}_{c}=10^{-2}$, as in the data for
$\mathcal{P}^{k}$ in Fig. 13(a) (bottom panel) 262626 However, the results are
insensitive to the choice of $\mathcal{P}_{c}$ for $5\times
10^{-3}<\mathcal{P}_{c}<2\times 10^{-2}$ [51]. . We plot the reduced versions
of $g_{\text{ex}}(\omega)$ and $g_{\text{loc}}(\omega)$ in Fig. 13(a) (top
panel). We also plot the vDOS itself in Fig. 13(b). $g_{\text{ex}}(\omega)$
converges exactly to the Debye vDOS $A_{D}\omega^{2}$ at a finite value of
$\omega$, which we define as $\omega_{\text{ex0}}$. On the other hand,
$g_{\text{loc}}(\omega)$ follows a different scaling law,
$g_{\text{loc}}(\omega)\propto\omega^{4}$. Thus, we conclude that the phonon
modes that follow the Debye law ($g_{\text{ex}}(\omega)=A_{D}\omega^{2}$) and
the localized modes that follow the other, non-Debye law
($g_{\text{loc}}(\omega)\propto\omega^{4}$) coexist at
$\omega<\omega_{\text{ex0}}$. Note that the $\omega^{4}$ scaling is the same
law proposed in the soft-potential model [88, 89, 90, 91, 92, 93]. The total
vDOS $g(\omega)$ can therefore be described as

$g(\omega)=g_{\text{ex}}(\omega)+g_{\text{loc}}(\omega)=A_{D}\omega^{2}+A_{\text{loc}}\omega^{4}.$
(90)

Since $g_{\text{loc}}(\omega)$ always takes finite values at
$\omega<\omega_{\text{ex0}}$, $g(\omega)$ coincides with the Debye vDOS only
at zero frequency. This means that the disordered structures of glasses are
not uniformly coarse-grained even at macroscopic scales but rather continue to
produce localized vibrations.

Based on the vDOS in Eq. (90), the heat capacity $C(T)$ can be predicted
within the harmonic approximation as follows [7, 6]:

$C(T)=3\int\hbar\omega\left(\frac{\partial f(\omega,T)}{\partial
T}\right)g(\omega)d\omega,$ (91)

where $f(\omega,T)=\left[\exp(\hbar\omega/k_{B}T)-1\right]^{-1}$ is the Bose-
Einstein distribution, $k_{B}$ is the Boltzmann constant, and $\hbar=h/2\pi$
with $h$ being the Plank constant. At the low temperature,
$k_{B}T\lesssim\hbar\omega_{\text{ex0}}$, $C(T)$ consists of two terms: the
Debye term $C_{\text{ex}}(T)\propto T^{3}$ and the non-Debye term
$C_{\text{loc}}(T)\propto T^{5}$. This result for $C(T)$ in the harmonic
approximation cannot correctly capture experimental observations [13, 14, 15],
in which a linear $T$ dependence appears at low temperatures. The present
results therefore demonstrate that anharmonicities should play an important
role in the low-$T$ properties of glasses, e.g., a two-level system might give
rise to the linear $T$ dependence [94, 23, 24]. This is completely different
to the case of crystals, which are well described within the harmonic
approximation [7, 6].

(a) Transverse phonon transport.

(b) Longitudinal phonon transport.

Figure 14:  Phonon transport in a glass with a harmonic potential. (a)
Transverse phonon transport. (b) Longitudinal phonon transport. We plot the
sound speed $c_{\alpha}(\Omega)$ (top panel) and the attenuation rate
$\Gamma_{\alpha}(\Omega)$ (bottom) as functions of the frequency $\Omega$. The
dashed line indicates the macroscopic sound speed (top panel) or the power-law
scaling with $\Omega$ (bottom). The inset in the bottom panel plots
$\Gamma_{\alpha}/\Omega^{4}$ versus $\Omega$ to clearly demonstrate the
Rayleigh scattering law. We use arrows to indicate the characteristic
frequencies $\omega_{\text{ex0}}$, $\omega_{\text{BP}}$, and $\omega_{\ast}$
that are obtained from the vibrational eigenmodes as shown in Fig. 13.

####  2 Phonon transport in glasses

We next present the results for phonon transport in a glass with a harmonic
potential. Figure 14 shows the sound speeds, $c_{\alpha}(\Omega)$, and the
attenuation rates, $\Gamma_{\alpha}(\Omega)$, for transverse ($\alpha=T$, Fig.
14(a)) and longitudinal ($\alpha=L$, Fig. 14(b)) phonon transport. These
results were obtained through direct measurement at zero temperature (see
Section 4). In this figure, arrows indicate the characteristic frequencies,
$\omega_{\ast}$, $\omega_{\text{BP}}$, and $\omega_{\text{ex0}}$, all of which
are obtained from the eigenmode results in Fig. 13. In the boson peak regime,
$\Omega\sim\omega_{\text{BP}}$, $c_{\alpha}$ takes its minimum value,
corresponding to sound softening. $\Gamma_{\alpha}$ shows an $\Omega^{2}$
dependence, $\Gamma_{\alpha}\propto\Omega^{2}$. We also find that the IR limit
for transverse phonons is $\Omega_{T\text{IR}}\approx\omega_{\text{BP}}$.
These results indicate that such a phonon does not propagate as a plane wave
but rather exhibits dynamics characteristic of viscous damping [37]. As
$\Omega$ decreases to $\Omega\lesssim\omega_{\text{ex0}}$, we observe a clear
crossover to Rayleigh scattering behaviour. $c_{\alpha}$ converges to its
macroscopic value: $c_{T0}=\sqrt{G/\rho}$ or $c_{L0}=\sqrt{(K+4G/3)/\rho}$.
Additionally, $\Gamma_{\alpha}$ shows an $\Omega^{4}$ dependence,
$\Gamma_{\alpha}\propto\Omega^{4}$.

All of the above observations were previously discussed in Section 4. The
crossover in the nature of the phonon transport has been observed in both
experiments [18, 19, 20] and simulations [38, 39, 41, 60, 45]. This crossover
can be predicted by mean-field theories [9, 10, 11, 12]. In particular, both
transverse and longitudinal phonons show similar transport properties: we
indeed observe the same crossover frequency $\omega_{\text{ex0}}$ for both of
them. This result is consistent with the prediction of heterogeneous
elasticity theory [9, 10]: the shear modulus heterogeneity induces anomalous
behaviours for both transverse and longitudinal phonon transport.

We can understand the crossover in phonon transport in terms of the underlying
vibrational eigenmodes as follows. Because the eigenmodes at
$\Omega\sim\omega_{\text{BP}}$ correspond to disordered and extended
vibrations and an initially excited phonon is decomposed into these
eigenmodes, it immediately attenuates to become diffusive. On the other hand,
the eigenmodes at $\Omega\lesssim\omega_{\text{ex0}}$ consist of phonon modes
and localized modes. Here, we measure the overlap $O_{\text{loc}}$ between the
initially excited phonon,
$\dot{\mathbf{u}}_{0}=\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}$ [Eq. (44)],
and the localized modes, $\mathbf{e}^{k}$, with participation ratio
$\mathcal{P}^{k}<10^{-2}$, by calculating

$O_{\text{loc}}=\sum_{k;\
\mathcal{P}^{k}<10^{-2}}\left|\mathbf{u}^{\mathbf{q},\alpha}_{\text{ph}}\cdot\mathbf{e}^{k}\right|^{2},$
(92)

and we find that only a few $\%$ (at most) of the initially excited phonon is
made up of localized modes. From this observation, we can understand that the
initially excited phonon is decomposed mainly into phonon modes, so it
attenuates only slowly. Therefore, the present results unambiguously link
these two types of phonon transport to the eigenmodes in the corresponding
frequency regimes. We are able to identify the crossover frequency as
$\omega_{\text{ex0}}$, at which the nature of the underlying eigenmodes
changes.

As was already discussed in Section 3, we note that the IR limit for
longitudinal phonons, $\Omega_{L\text{IR}}$, is much higher than that for
transverse phonons:
$\Omega_{L\text{IR}}\gg\Omega_{T\text{IR}}\approx\omega_{\text{BP}}$.
Therefore, a longitudinal phonon can propagate even at
$\Omega>\omega_{\text{BP}}$. This result indicates that although the
disordered eigenmodes are dominant at $\Omega>\omega_{\text{BP}}$,
longitudinal phonon modes also exist in this $\Omega$ regime, which support
the propagation of longitudinal phonons.

In summary, we have shown that even in the low-frequency limit, the glass
exhibits localized eigenmodes and Rayleigh scattering in phonon transport. The
same conclusion has also been obtained for an LJ glass [38, 52]. These results
indicate that the disordered structures of glasses are not uniformly coarse-
grained even at macroscopic scales but rather play the role of defects,
influencing the vibrational properties of the material. Below the frequency
$\omega_{\text{ex0}}$, which corresponds to the continuum limit frequency,
glasses behave as elastic media with defects.

As a final remark, the above results are $T=0$ properties, in the harmonic
approximation limit. Here we mention anharmonic properties of vibrational
excitations. At finite temperatures $T>0$, systems undergo anharmonic
processes. In crystals phonons excited by thermal fluctuations couple through
phonon-phonon interactions [6, 7]. For example, in three-phonon processes, one
phonon splits into two different ones, or conversely, two phonons combine into
a single excitation. Similarly, in glasses, vibrational eigenmodes also show
anharmonic processes due to mode-mode interactions [95, 96].

In contrast to crystals, however, it has been demonstrated that an additional
anharmonic channel of different origin emerges in glasses, which induces
unconventional intermittent rearrangements of particles [97, 98, 99, 100]. The
intermittent rearrangements can be considered as relics of the liquid state
which survive the complete dynamic arrest taking place at the glass transition
temperature [101, 102, 103]. It was also demonstrated that glasses exhibit
anomalous temperature dependence of phonon attenuation at finite temperatures
[21, 22, 104, 105, 106], which can be attributed to this additional anharmonic
channel.

## References

  * 1. M. P. Allen and D. J. Tildesley, _Computer Simulation of Liquids_. Oxford University Press, Oxford (1986).
  * 2. D. Frenkel and B. Smit, _Understanding Molecular Simulation_ , second edn. Academic Press, New York (2002).
  * 3. M. Newman and G. Barkema, _Monte Carlo Methods in Statistical Physics_. Oxford University Press, U.S.A. (1999).
  * 4. A. Cavagna, Supercooled liquids for pedestrians, _Physics Reports_. 476(4), 51 – 124 (2009).  ISSN 0370-1573.  https://doi.org/10.1016/j.physrep.2009.03.003.  URL http://www.sciencedirect.com/science/article/pii/S0370157309001112.
  * 5. L. Berthier and G. Biroli, Theoretical perspective on the glass transition and amorphous materials, _Rev. Mod. Phys._ 83, 587–645 (Jun, 2011).  10.1103/RevModPhys.83.587.  URL https://link.aps.org/doi/10.1103/RevModPhys.83.587.
  * 6. N. W. Ashcroft and N. D. Mermin, _Solid State Physics_. Harcourt College Publishers, New York (1976).
  * 7. C. Kittel, _Introduction to Solid State Physics_ , 7th edn. John Wiley and Sons, New York (1996).
  * 8. S. Franz, G. Parisi, P. Urbani, and F. Zamponi, Universal spectrum of normal modes in low-temperature glasses, _Proceedings of the National Academy of Sciences_. 112, 14539 (2015).
  * 9. W. Schirmacher, Thermal conductivity of glassy materials and the “boson peak”, _Europhys. Lett._ 73, 892 (2006).
  * 10. W. Schirmacher, G. Ruocco, and T. Scopigno, Acoustic attenuation in glasses and its relation with the boson peak, _Phys. Rev. Lett._ 98, 025501 (2007).
  * 11. M. Wyart, Scaling of phononic transport with connectivity in amorphous solids, _EPL (Europhysics Letters)_. 89(6), 64001 (2010).  URL http://stacks.iop.org/0295-5075/89/i=6/a=64001.
  * 12. E. DeGiuli, A. Laversanne-Finot, G. During, E. Lerner, and M. Wyart, Effects of coordination and pressure on sound attenuation, boson peak and elasticity in amorphous solids, _Soft Matter_. 10, 5628–5644 (2014).  10.1039/C4SM00561A.  URL http://dx.doi.org/10.1039/C4SM00561A.
  * 13. R. C. Zeller and R. O. Pohl, Thermal conductivity and specific heat of noncrystalline solids, _Phys. Rev. B_. 4, 2029–2041 (Sep, 1971).  10.1103/PhysRevB.4.2029.  URL http://link.aps.org/doi/10.1103/PhysRevB.4.2029.
  * 14. W. A. Phillips, _Amorphous Solids: Low Temperature Properties_ , third edn. Springer, Berlin (1981).
  * 15. S. R. Elliott, _Physics of Amorphous Materials_ , 2nd edn. Longman, New York (1990).
  * 16. U. Buchenau, N. Nücker, and A. J. Dianoux, Neutron scattering study of the low-frequency vibrations in vitreous silica, _Phys. Rev. Lett._ 53, 2316 (1984).
  * 17. O. Yamamuro, T. Matsuo, K. Takeda, T. Kanaya, T. Kawaguchi, and K. Kaji, Inelastic neutron scattering study of low energy excitations in glassy 1‐butene, _The Journal of Chemical Physics_. 105(2), 732–737 (1996).  10.1063/1.471928.  URL https://doi.org/10.1063/1.471928.
  * 18. C. Masciovecchio, G. Baldi, S. Caponi, L. Comez, S. Di Fonzo, D. Fioretto, A. Fontana, A. Gessini, S. C. Santucci, F. Sette, G. Viliani, P. Vilmercati, and G. Ruocco, Evidence for a crossover in the frequency dependence of the acoustic attenuation in vitreous silica, _Phys. Rev. Lett._ 97, 035501 (Jul, 2006).  10.1103/PhysRevLett.97.035501.  URL http://link.aps.org/doi/10.1103/PhysRevLett.97.035501.
  * 19. G. Monaco and V. M. Giordano, Breakdown of the debye approximation for the acoustic modes with nanometric wavelengths in glasses, _Proc. Natl. Acad. Sci. USA_. 106, 3659 (2009).
  * 20. G. Baldi, V. M. Giordano, G. Monaco, and B. Ruta, Sound attenuation at terahertz frequencies and the boson peak of vitreous silica, _Phys. Rev. Lett._ 104, 195501 (May, 2010).  10.1103/PhysRevLett.104.195501.  URL https://link.aps.org/doi/10.1103/PhysRevLett.104.195501.
  * 21. C. Ferrante, E. Pontecorvo, G. Cerullo, A. Chiasera, G. Ruocco, W. Schirmacher, and T. Scopigno, Acoustic dynamics of network-forming glasses at mesoscopic wavelengths, _Nature communications_. 4, 1793 (2013).
  * 22. G. Baldi, V. M. Giordano, B. Ruta, R. Dal Maschio, A. Fontana, and G. Monaco, Anharmonic damping of terahertz acoustic waves in a network glass and its effect on the density of vibrational states, _Phys. Rev. Lett._ 112, 125502 (Mar, 2014).  10.1103/PhysRevLett.112.125502.  URL http://link.aps.org/doi/10.1103/PhysRevLett.112.125502.
  * 23. T. Pérez-Castañeda, C. Rodríguez-Tinoco, J. Rodríguez-Viejo, and M. A. Ramos, Suppression of tunneling two-level systems in ultrastable glasses of indomethacin, _Proceedings of the National Academy of Sciences_. 111(31), 11275–11280 (2014).  10.1073/pnas.1405545111.
  * 24. T. Pérez-Castañeda, R. J. Jiménez-Riobóo, and M. A. Ramos, Two-level systems and boson peak remain stable in 110-million-year-old amber glass, _Phys. Rev. Lett._ 112, 165901 (Apr, 2014).  10.1103/PhysRevLett.112.165901.  URL https://link.aps.org/doi/10.1103/PhysRevLett.112.165901.
  * 25. M. Kabeya, T. Mori, Y. Fujii, A. Koreeda, B. W. Lee, J.-H. Ko, and S. Kojima, Boson peak dynamics of glassy glucose studied by integrated terahertz-band spectroscopy, _Phys. Rev. B_. 94, 224204 (Dec, 2016).  10.1103/PhysRevB.94.224204.  URL https://link.aps.org/doi/10.1103/PhysRevB.94.224204.
  * 26. H. R. Schober and B. B. Laird, Localized low-frequency vibrational modes in glasses, _Phys. Rev. B_. 44, 6746–6754 (Oct, 1991).  10.1103/PhysRevB.44.6746.  URL http://link.aps.org/doi/10.1103/PhysRevB.44.6746.
  * 27. V. Mazzacurati, G. Ruocco, and M. Sampoli, Low-frequency atomic motion in a model glass, _Europhys. Lett._ 34, 681 (1996).
  * 28. S. N. Taraskin and S. R. Elliott, Anharmonicity and localization of atomic vibrations in vitreous silica, _Phys. Rev. B_. 59, 8572–8585 (Apr, 1999).  10.1103/PhysRevB.59.8572.  URL http://link.aps.org/doi/10.1103/PhysRevB.59.8572.
  * 29. S. N. Taraskin and S. R. Elliott, Propagation of plane-wave vibrational excitations in disordered systems, _Phys. Rev. B_. 61, 12017 (2000).
  * 30. G. Ruocco, F. Sette, R. Di Leonardo, G. Monaco, M. Sampoli, T. Scopigno, and G. Viliani, Relaxation processes in harmonic glasses?, _Phys. Rev. Lett._ 84, 5788–5791 (Jun, 2000).  10.1103/PhysRevLett.84.5788.  URL https://link.aps.org/doi/10.1103/PhysRevLett.84.5788.
  * 31. A. Tanguy, J. P. Wittmer, F. Leonforte, and J.-L. Barrat, Continuum limit of amorphous elastic bodies: A finite-size study of low-frequency harmonic vibrations, _Phys. Rev. B_. 66, 174205 (2002).
  * 32. F. Leonforte, R. Boissière, A. Tanguy, J. P. Wittmer, and J.-L. Barrat, Continuum limit of amorphous elastic bodies. iii. three-dimensional systems, _Phys. Rev. B_. 72, 224206 (Dec, 2005).  10.1103/PhysRevB.72.224206.  URL https://link.aps.org/doi/10.1103/PhysRevB.72.224206.
  * 33. F. Léonforte, A. Tanguy, J. P. Wittmer, and J.-L. Barrat, Inhomogeneous elastic response of silica glass, _Phys. Rev. Lett._ 97, 055501 (Jul, 2006).  10.1103/PhysRevLett.97.055501.  URL https://link.aps.org/doi/10.1103/PhysRevLett.97.055501.
  * 34. M. Wyart, L. E. Silbert, S. R. Nagel, and T. A. Witten, Effects of compression on the vibrational modes of marginally jammed solids, _Phys. Rev. E_. 72, 051306 (Nov, 2005).  10.1103/PhysRevE.72.051306.  URL http://link.aps.org/doi/10.1103/PhysRevE.72.051306.
  * 35. L. E. Silbert, A. J. Liu, and S. R. Nagel, Vibrations and diverging length scales near the unjamming transition, _Phys. Rev. Lett._ 95, 098301 (Aug, 2005).  10.1103/PhysRevLett.95.098301.  URL http://link.aps.org/doi/10.1103/PhysRevLett.95.098301.
  * 36. L. E. Silbert, A. J. Liu, and S. R. Nagel, Normal modes in model jammed systems in three dimensions, _Phys. Rev. E_. 79, 021308 (Feb, 2009).  10.1103/PhysRevE.79.021308.  URL http://link.aps.org/doi/10.1103/PhysRevE.79.021308.
  * 37. H. Shintani and H. Tanaka, Universal link between the boson peak and transverse phonons in glass, _Nature Mater._ 7, 870 (2008).
  * 38. G. Monaco and S. Mossa, Anomalous properties of the acoustic excitations in glasses on the mesoscopic length scale, _Proceedings of the National Academy of Sciences_. 106, 16907 (2009).
  * 39. A. Marruzzo, W. Schirmacher, A. Fratalocchi, and G. Ruocco, Heterogeneous shear elasticity of glasses: the origin of the boson peak, _Scientific Reports_. 3, 1407 (2013).
  * 40. H. Mizuno, S. Mossa, and J.-L. Barrat, Elastic heterogeneity, vibrational states, and thermal conductivity across an amorphisation transition, _EPL (Europhysics Letters)_. 104(5), 56001 (2013).  URL http://stacks.iop.org/0295-5075/104/i=5/a=56001.
  * 41. H. Mizuno, S. Mossa, and J.-L. Barrat, Acoustic excitations and elastic heterogeneities in disordered solids, _Proceedings of the National Academy of Sciences_. 111(33), 11949–11954 (2014).  10.1073/pnas.1409490111.  URL http://www.pnas.org/content/111/33/11949.abstract.
  * 42. H. Mizuno, S. Mossa, and J.-L. Barrat, Relation of vibrational excitations and thermal conductivity to elastic heterogeneities in disordered solids, _Phys. Rev. B_. 94, 144303 (Oct, 2016).  10.1103/PhysRevB.94.144303.  URL https://link.aps.org/doi/10.1103/PhysRevB.94.144303.
  * 43. A. Ikeda, L. Berthier, and G. Biroli, Dynamic criticality at the jamming transition, _J. Chem. Phys._ 138(12), 12A507 (2013).  http://dx.doi.org/10.1063/1.4769251.  URL http://scitation.aip.org/content/aip/journal/jcp/138/12/10.1063/1.4769251.
  * 44. D. Crespo, P. Bruna, A. Valles, and E. Pineda, Phonon dispersion relation of metallic glasses, _Phys. Rev. B_. 94, 144205 (Oct, 2016).  10.1103/PhysRevB.94.144205.  URL https://link.aps.org/doi/10.1103/PhysRevB.94.144205.
  * 45. Y. M. Beltukov, C. Fusco, D. A. Parshin, and A. Tanguy, Boson peak and ioffe-regel criterion in amorphous siliconlike materials: The effect of bond directionality, _Phys. Rev. E_. 93, 023006 (Feb, 2016).  10.1103/PhysRevE.93.023006.  URL https://link.aps.org/doi/10.1103/PhysRevE.93.023006.
  * 46. S. Gelin, H. Tanaka, and A. Lemaitre, Anomalous phonon scattering and elastic correlations in amorphous solids, _Nature Mater._ 15(11), 1177–1181 (Nov., 2016).  ISSN 1476-1122.  URL http://dx.doi.org/10.1038/nmat4736.
  * 47. R. Milkus and A. Zaccone, Local inversion-symmetry breaking controls the boson peak in glasses and crystals, _Phys. Rev. B_. 93, 094204 (Mar, 2016).  10.1103/PhysRevB.93.094204.  URL https://link.aps.org/doi/10.1103/PhysRevB.93.094204.
  * 48. J. Krausser, R. Milkus, and A. Zaccone, Non-affine lattice dynamics of defective fcc crystals, _Soft Matter_. 13, 6079–6089 (2017).  10.1039/C7SM00843K.  URL http://dx.doi.org/10.1039/C7SM00843K.
  * 49. K. Baumgarten, D. Vågberg, and B. P. Tighe, Nonlocal elasticity near jamming in frictionless soft spheres, _Phys. Rev. Lett._ 118, 098001 (Feb, 2017).  10.1103/PhysRevLett.118.098001.  URL https://link.aps.org/doi/10.1103/PhysRevLett.118.098001.
  * 50. E. Lerner, G. Düring, and E. Bouchbinder, Statistics and properties of low-frequency vibrational modes in structural glasses, _Phys. Rev. Lett._ 117, 035501 (Jul, 2016).  10.1103/PhysRevLett.117.035501.  URL http://link.aps.org/doi/10.1103/PhysRevLett.117.035501.
  * 51. H. Mizuno, H. Shiba, and A. Ikeda, Continuum limit of the vibrational properties of amorphous solids, _Proceedings of the National Academy of Sciences_. 114(46), E9767–E9774 (2017).  10.1073/pnas.1709015114.
  * 52. M. Shimada, H. Mizuno, and A. Ikeda, Anomalous vibrational properties in the continuum limit of glasses, _Phys. Rev. E_. 97, 022609 (Feb, 2018).  10.1103/PhysRevE.97.022609.  URL https://link.aps.org/doi/10.1103/PhysRevE.97.022609.
  * 53. E. Lerner and E. Bouchbinder, Effect of instantaneous and continuous quenches on the density of vibrational modes in model glasses, _Phys. Rev. E_. 96, 020104 (Aug, 2017).  10.1103/PhysRevE.96.020104.  URL https://link.aps.org/doi/10.1103/PhysRevE.96.020104.
  * 54. E. Lerner and E. Bouchbinder, Frustration-induced internal stresses are responsible for quasilocalized modes in structural glasses, _Phys. Rev. E_. 97, 032140 (Mar, 2018).  10.1103/PhysRevE.97.032140.  URL https://link.aps.org/doi/10.1103/PhysRevE.97.032140.
  * 55. L. Angelani, M. Paoluzzi, G. Parisi, and G. Ruocco, Probing the non-debye low-frequency excitations in glasses through random pinning, _Proceedings of the National Academy of Sciences_. 115(35), 8700–8704 (2018).
  * 56. E. Bouchbinder and E. Lerner, Universal disorder-induced broadening of phonon bands: from disordered lattices to glasses, _New Journal of Physics_. 20(7), 073022 (2018).  URL http://stacks.iop.org/1367-2630/20/i=7/a=073022.
  * 57. H. Mizuno and A. Ikeda, Phonon transport and vibrational excitations in amorphous solids, _Phys. Rev. E_. 98, 062612 (Dec, 2018).  10.1103/PhysRevE.98.062612.  URL https://link.aps.org/doi/10.1103/PhysRevE.98.062612.
  * 58. M. Shimada, H. Mizuno, M. Wyart, and A. Ikeda, Spatial structure of quasilocalized vibrations in nearly jammed amorphous solids, _Phys. Rev. E_. 98, 060901 (Dec, 2018).  10.1103/PhysRevE.98.060901.  URL https://link.aps.org/doi/10.1103/PhysRevE.98.060901.
  * 59. L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel, and E. Flenner, Low-frequency vibrational modes of stable glasses, _Nature Communications_. 10, 26 (2019).