# Quantifying and Using System Uncertainty
in UAV Navigation

Fabio Arnez  Université Paris-Saclay, CEA, List
F-91120, Palaiseau, France
<EMAIL_ADDRESS>Ansgar Radermacher  Université Paris-Saclay, CEA, List
F-91120, Palaiseau, France
<EMAIL_ADDRESS>Huascar Espinoza  KDT JU
Brussels, Belgium
<EMAIL_ADDRESS>
###### Abstract

As autonomous systems increasingly rely on Deep Neural Networks (DNN) to
implement the navigation pipeline functions, uncertainty estimation methods
have become paramount for estimating confidence in DNN predictions. Bayesian
Deep Learning (BDL) offers a principled approach to model uncertainties in
DNNs. However, DNN components from autonomous systems partially capture
uncertainty, or more importantly, the uncertainty effect in downstream tasks
is ignored. This paper provides a method to capture the overall system
uncertainty in a UAV navigation task. In particular, we study the effect of
the uncertainty from perception representations in downstream control
predictions. Moreover, we leverage the uncertainty in the system’s output to
improve control decisions that positively impact the UAV’s performance on its
task.

††Preprint version. Accepted at the ICRA 2022 Workshop on Releasing Robots
into the Wild: Simulations, Benchmarks, and Deployment. Copyright 2022 by the
author(s).

##  I Introduction

Navigation in complex environments still represents a big challenge for
automated systems (AS). Particular instances of this problem are autonomous
driving and autonomous aerial navigation in the context of Unmanned Aerial
Vehicles (UAVs). In both cases, the navigation task is addressed by first
acquiring rich and complex raw sensory information (e.g., from camera, radar,
LiDAR, etc.), which is then processed to drive the agent towards its goal.
Usually, this process is done in sequence, where tasks and specific software
components are linked together in the so-called perception-planning-control
software pipeline [1, 2]. In the last decade, Deep Neural Networks (DNNs) have
become a popular choice to implement navigation components. For this purpose,
three paradigms exist to develop and train DNN-based components: Modular
learning (isolated) [3], End-to-End (E2E) learning [4, 5, 6], and a mixed
approach [7, 8].

Despite the remarkable progress in representation learning, DNNs should also
represent the confidence in their predictions to deploy them in safety-
critical systems. McAllister et al. [2] proposed using Bayesian Deep Learning
(BDL) to implement the components from navigation pipelines. Bayesian methods
offer a principled framework to model and capture uncertainty. Nevertheless,
if the Bayesian approach is followed, all the components in the system should
use BDL methods to enable uncertainty propagation in the navigation pipeline.
Hence, BDL components should admit uncertainty information as an input to
account for the uncertainty from the outputs of preceding BDL components (See
Fig. 1 bottom).

In recent years, a large body of literature has employed uncertainty
estimation methods in robotic tasks thanks to its potential to improve the
safety of automated functions [9], and the capacity to increase the task
performance [10, 11]. However, uncertainty is captured partially in navigation
pipelines that utilize DNNs. BDL methods are used mainly in perception task,
and downstream components (e.g., planning and control) usually ignore the
uncertainty from the preceding components or do not capture uncertainty in
their predictions. Although some works propagate downstream perceptual
uncertainty from intermediate representations [12, 13, 14], the overall system
output does not take into account all the uncertainty sources from DNN
components in the pipeline.

Figure 1: UAV Uncertainty-Aware Aerial Navigation Pipeline

Quantifying uncertainty in a BDL-based system (i.e., a pipeline of BDL
components) still remains a challenging task. Uncertainties from BDL
components must be assembled in a principled way to provide a reliable measure
of overall system uncertainty, based on which safe decisions can be made [2,
15, 16] (See Fig. 1 top). In this paper, we propose to capture and use the the
overall system uncertainty in the navigation pipeline to improve the UAV
performance on its task.

##  II System Task Formulation

In this paper, we address the problem of autonomous aerial navigation. The
goal of the UAV is to navigate through a set of gates with unknown locations
disposed in a circular track. Following prior work from [7, 12], the
navigation architecture consists of two neural networks: one for perception
and the other for control, as shown in Fig. 2. To achieve goal, the navigation
task is formulated as a sequential-decision making problem, where a control
action is produced given an environment observation. In this regard, the
simulation environment provides at time $t$ an observation comprised of an RGB
image $\mathbf{x}$ acquired from a front-facing camera on the UAV. The
perception component defines an encoder function
$q_{\phi}:\mathcal{X}\rightarrow\mathcal{Z}$ that maps the input image
$\mathbf{x}$ to a rich low dimensional representation
$\mathbf{z}\in\mathbb{R}^{10}$. Next, a control policy
$\pi_{w}:\mathcal{Z}\rightarrow\mathcal{Y}$ maps the compact representation
$\mathbf{z}$, to control commands
$\mathbf{y}=[\dot{x},\dot{y},\dot{z},\dot{\psi}]\in\mathbb{R}^{4}$,
corresponding to linear and yaw velocities in the UAV body frame.

In the perception component, a cross-modal variational autoencoder (CMVAE)
[17, 7] is used to learn a rich and robust compact representation. A CMVAE is
a variant of traditional variational autoencoders (VAE) [18] that learns a
single latent representation for multiple data modalities. In this case, two
data modalities are used: the RGB images and the pose of the gate relative to
the UAV body-frame. The CMVAE encoder $q_{\phi}$ maps then an input image
$\mathbf{x}$ to a noisy representation with mean $\mu_{\phi}(\mathbf{x})$ and
variance $\sigma^{2}_{\phi}(\mathbf{x})$ in the latent space, from where
latent vectors $z$ are sampled,
$\mathbf{z}\sim\mathcal{N}(\mu_{\phi},\sigma^{2}_{\phi})$. The encoder
$q_{\phi}$ is based on the Dronet architecture [4], and additional constraints
on the latent space are imposed to promote the learning of robust disentangled
representations. For the downstream control task (control policy $\pi$), a
feed-forward network is used to operate on the latent vectors $\mathbf{z}$.
For more information about the general architecture for aerial navigation, we
refer the reader to [7, 12].

##  III Methodology

###  III-A Uncertainty from Perception Representations

Although the CMVAE encoder $q_{\phi}$ employs Bayesian inference to obtain
latent vectors $\mathbf{z}$, CM-VAE does not capture epistemic uncertainty
since the encoder lacks a distribution over parameters $\phi$. To capture
uncertainty in the perception encoder we follow prior work from [19, 20] that
attempts to capture epistemic uncertainty in VAEs. We adapt the CM-VAE to
capture the posterior $q(\mathbf{z}\mid\mathbf{x},\mathcal{D}_{p})$ as shown
in eq. 1.

$q_{\Phi}(\mathbf{z}\mid\mathbf{x},\mathcal{D}_{p})=\int{q(\mathbf{z}\mid\mathbf{x},\phi)p(\phi\mid\mathcal{D}_{p})d\phi}$
(1)

To approximate eq. 1 we take a set ${\Phi=\\{\phi_{m}\\}^{M}_{m}}$ of encoder
parameters samples $\phi_{m}\sim p(\phi\mid\mathcal{D}_{p})$, to obtain a set
of latent samples $\\{\mathbf{z}_{m}\\}^{M}_{m=1}\sim
q_{\Phi}(\mathbf{z}\mid\mathbf{x},\mathcal{D}_{p})$ at the output of the
encoder. In practice, we modify CMVAE by adding a dropout layer in the
encoder. Then, we use Monte Carlo Dropout (MCD) [21] to approximate the
posterior on the encoder weights $p(\phi\mid\mathcal{D}_{p})$. Finally, for a
given input image $\mathbf{x}$ we perform $M$ stochastic forward passes (with
dropout turned on) to compute a set of $M$ latent vector samples $\mathbf{z}$
at runtime.

Figure 2: System architecture for aerial navigation

###  III-B Input Uncertainty for Control

In BDL, downstream uncertainty propagation assumes that a neural network
component is able to handle or admit uncertainty at the input. In our case,
this implies that the neural network for control is able to handle the
uncertainty coming from the perception component. When BNNs are used with LV
(BNN+LV), both types of uncertainty can be captured [22, 23, 24, 12]. To
capture the overall system uncertainty at the output of the controller, we
compute the posterior predictive distribution for target variable
$\mathbf{y^{*}}$ associated with a new input image $\mathbf{x^{*}}$, as shown
in eq. 2.

$p(\mathbf{y^{*}}\mid\mathbf{x^{*}},\mathcal{D}_{c},\mathcal{D}_{p})=\\\
\iint{\pi(\mathbf{y}\mid\mathbf{z},\mathbf{w})p(\mathbf{w}\mid\mathcal{D}_{c})q_{\Phi}(\mathbf{z}\mid\mathbf{x^{*}},\mathcal{D}_{p})d\mathbf{w}d\mathbf{z}}$
(2)

The above integrals are intractable, and we rely on approximations to obtain
an estimation of the predictive distribution. The posterior
$p(\mathbf{w}\mid\mathcal{D}_{c})$ is difficult to evaluate, thus we can
approximate the inner integral using an ensemble of neural networks [25]. In
practice, we train an ensemble of $N$ probabilistic control policies
$\\{\mathbf{\pi}_{n}(\mathbf{y}\mid\mathbf{z},\mathbf{w}_{n})\\}_{n=1}^{N}$,
with weights $\\{\mathbf{w}_{n}\\}^{N}_{n=1}\sim p(\mathbf{w|\mathcal{D}})$,
and where each control policy $\pi_{n}$ in the ensemble predicts the mean
$\mu$ and variance $\sigma^{2}$ for each velocity command, i.e.,
$\mathbf{y}_{\mu}=[\mu_{\dot{x}},\mu_{\dot{y}},\mu_{\dot{z}},\mu_{\dot{\psi}}]$
and
$\mathbf{y}_{\sigma^{2}}=[\sigma^{2}_{\dot{x}},\sigma^{2}_{\dot{y}},\sigma^{2}_{\dot{z}},\sigma^{2}_{\dot{\psi}}]$.
For training each probabilistic control policy we use imitation learning and
the heteroscedastic loss function, as suggested by [26, 27].

The outer integral is approximated by taking a set of samples from the
perception component latent space. In [12] latent samples are drawn using the
encoder mean and variance
$\mathbf{z}\sim\mathcal{N}(\mu_{\phi},\sigma^{2}_{\phi})$. For the sake of
simplicity, we directly use the samples obtained in the perception component
$\\{\mathbf{z}_{m}\\}^{M}_{m}\sim
q_{\Phi}(\mathbf{z}\mid\mathbf{x},\mathcal{D}_{p})$ to take into account the
epistemic uncertainty from the previous stage. Finally, the predictions that
we get from passing each latent vector $\mathbf{z}$ through each ensemble
member are used to estimate the posterior predictive distribution in eq. 2.
From the control policy perspective, using multiple latent samples
$\mathbf{z}$ can be seen as taking a better picture of the latent space
(perception representation) to gather more information about the environment.
Interestingly, we can also make a connection between our sampling approach and
MEMO [28], a method that performs augmentations on test dataset points to
improve prediction robustness. Our method is illustrated in Fig. 2.

###  III-C Decision-Making Under Uncertainty

For control, assigning the same weight to each ensemble member prediction to
compute the ensemble mixture predictive mean and variance [27, 25] can result
in sub-optimal solutions when facing multimodal predictions that could arise
from ambiguous inputs. For example, the predictions from one control policy
network could attempt to move the UAV to the left, while the predictions from
another control policy could try to move the UAV to the right. This can result
in the UAV moving straight and hence can lead to a tragic result.

To overcome this problem, we take inspiration from techniques on active
learning for Bayesian Deep Learning [29, 30], that use mutual information.
However, in our approach, given the samples from the latent space, we propose
to choose the control predictions (the density) from the ensemble member that
minimizes the mutual information, as presented in eq. 3.

${\pi}^{*}=\underset{\pi}{\arg\min}\big{(}\mathbb{I}(y;z,w)\big{)},\forall
y\in\mathbf{y}=[\dot{x},\dot{y},\dot{z},\dot{\psi}]$ (3)

In our navigation architecture, the mutual information $\mathbb{I}$ is
formulated as follows,

$\mathbb{I}(y;z,w)=\int_{z}{\mathbf{KL}\big{(}\pi(y\mid
z,w)\parallel\pi(y)\big{)}p(z)dz}$ (4)

To estimate the mutual information from the previous integral, we use the
variational lower bound approximation from [31] taking the latent samples
$\mathbf{z}$ and each ensemble member predictions ($\hat{\mu}$ and
$\hat{\sigma}$), as presented in eq. 5:

$\mathbb{I}(y;z,w)\geqslant\\\
\mathbb{E}_{z_{1:M}}\Big{[}\sum_{i=1}^{M}{\mathbf{KL}\big{(}\pi(y\mid
z,w)\parallel m(y;z_{1:M},w_{1:N})\big{)}}\Big{]}$ (5)

Once an ensemble member is chosen, we can use the control policy predicted
densities by taking the mean or the modes, as presented in Fig.2.

##  IV Experiments

For our experiments, we study the impact of the uncertainty from perception
representations in a downstream control component. We seek to answer the
following research questions. RQ1. How does uncertainty from perception
representations affect the downstream control predictions? RQ2. Can we improve
the UAV performance using perception and control uncertainties?

(a) Circular track view

(b) Noisy circular track view

Figure 3: Aerial navigation tracks: birds-eye view (left), and view from UAV
perspective (right)

###  IV-A Experimental setup

####  IV-A1 Navigation Model Baselines

All the navigation architectures are based on [7] and are implemented using
PyTorch. Table I shows the uncertainty-aware navigation architectures used in
our experiments, and the type of perception component, the number latent
variable samples (LVS), the type of control policy, and the number of control
prediction samples (CPS) at the output of the system. For instance,
$\mathcal{M}_{0}$ represents our Bayesian navigation pipeline.
$\mathcal{M}_{0}$ perception component captures epistemic uncertainty using
MCD with 32 forward passes for each input to get 32 latent variable
predictions. For the sake of simplicity, perception predictions are directly
used as latent variable samples in downstream control. The control component
uses a ensemble of 5 probabilistic control policies obtaining 160 control
prediction samples.

TABLE I: Uncertainty-aware Navigation models in the experiments Model | Perception $(q_{\phi})$ | LVS | Control Policy $(\pi)$ | CPS
---|---|---|---|---
$\mathcal{M}_{0}$ | MCD-CMVAE | 32 | Ensemble (5) Prob. | 160
$\mathcal{M}_{1}$ | CMVAE | 32 | Ensemble (5) Prob. | 160
$\mathcal{M}_{2}$ | CMVAE | 1 | Ensemble (5) Prob. | 5
$\mathcal{M}_{3}$ | CMVAE | 32 | Deterministic | 32
$\mathcal{M}_{4}$ | CMVAE | 1 | Prob. | 1

(a) Navigation model input, Perception prediction densities, and Control
policy ensemble prediction densities ($\hat{\mu},\hat{\sigma}$)

(b) Control prediction densities
($\hat{\mu}_{\dot{y}},\hat{\sigma}_{\dot{y}}$) and
($\hat{\mu}_{\dot{\psi}},\hat{\sigma}_{\dot{\psi}}$) for each ensemble member

Figure 4: Bayesian Navigation model: Component prediction densities

####  IV-A2 Datasets

We use two independent datasets for each component in the navigation pipeline.
The CMVAE uses a dataset of 300k images with gate pose labels. The control
component uses a dataset of 17k images with drone velocity labels. The
perception dataset split divided into 80% for training, and the remaining 20%
for validation and testing. The control dataset uses a split of 90% for
training and the remaining for validation and testing. In both cases the image
size is 64x64 pixels.

####  IV-A3 Evaluation Procedure

First, we observe $\mathcal{M}_{0}$ predictions with an out-of-distribution
input sample that introduces ambiguity to the model (double-gate). Then, we
evaluate the navigation model using AirSim. We use a circular track with eight
equally spaced gates positioned initially in a radius of 8 m and constant
height. To assess the system performance to perturbations in the environment,
we generate new tracks adding random noise to each gate radius and height. In
the context of the AirSim [32] simulation environment, a track is entirely
defined by a set of gates, their poses in three-dimensional space, and the
agent navigation direction. For perception-based navigation, the complexity of
a track resides in the gate-visibility difficulty [33, 34], i.e., how well the
camera Field-of-View (FoV) captures the gate. A natural way to increase track
complexity is by adding a random displacement to the position of each gate. A
track without random displacement in the gates has a circular fashion. Gate
position randomness alters the shape of the track, affecting the gate
visibility, i.e., gates are: not visible, partially visible or multiple gates
can be captured in the UAV FoV as presented in Fig. 3. We measure the system
performance by taking the average number of gates passed in six different
tracks. The UAV mission considers a maximum of 32 gates which is equivalent to
4 laps (8 gates / lap), and two levels of noise for gate offset: Gate Radius
Noise (GRN) and Gate Height Noise (GHN). Each navigation model has two trials
on each track. In addition, we consider two control decision-making strategies
(CDMS), to use control predictions. The first strategy uses the deep ensemble
mean. The second control strategy uses mutual information lower bound to
select the ensemble member predicted density from where we choose the lowest
velocity modes for $\dot{y},\dot{z},\dot{\psi}$, while for $\dot{x}$ we choose
a conservative strategy by selecting the lowest predicted velocity.

TABLE II: Uncertainty-Aware Navigation Models Performance:
Avg. number of gates passed Model | CDMS | Passed Gates with Noise Level
---|---|---
| $GRN=1.0$
---
$GHN=2.0$
| $GRN=1.5$
---
$GHN=3.0$
$\mathcal{M}_{0}$ | MI-mode | 28.88 | 23.05
$\mathcal{M}_{0}$ | DE-mean | 19.77 | 9.22
$\mathcal{M}_{1}$ | DE-mean | 17.67 | 6.0
$\mathcal{M}_{2}$ | DE-mean | 17.33 | 4.0
$\mathcal{M}_{3}$ | DE-mean | 8.33 | 5.0
$\mathcal{M}_{4}$ | DE-mean | 15.16 | 4.38

###  IV-B Results

Fig. 4(a) shows $\mathcal{M}_{0}$ predictions at the output of the perception
($\mathbf{z}$) and control ($\hat{\mu}$, $\hat{\sigma}$) components.
Predictions are made using an input sample image with two gates to see if the
model is able to capture the ambiguity from the sample. Interestingly, the
control densities, using all the CPS from the ensemble, show that
$\mathcal{M}_{0}$ is able to represent the ambiguity in the input. The
predicted velocity densities shows complex multimodal distributions (two
peaks) for $\hat{\mu}_{\dot{y}}$, $\hat{\mu}_{\dot{\psi}}$ and
$\hat{\sigma}_{\dot{\psi}}$. the multimodal predictions are observed more in
detail in Fig. 4(b), in the prediciton densities for each ensemble member for
$\dot{y}$ and $\dot{\psi}$ velocity commands.

Finally, Table II presents the navigation performance results. In general,
learning to predict uncertainty in the control components can boost the
performance significantly. In the case of $\mathcal{M}_{0}$, the strategy to
select control predictions clearly impacts the navigation performance. In
general, using the control ensemble mean can lead to similar results. However,
exploiting the ensemble members predicted distributions (as seen 4) can boost
the model performance.

##  V Conclusion

We presented a method to capture and propagate uncertainty along a navigation
pipeline implemented with Bayesian deep learning components for UAV
navigation. We analyzed the effect of uncertainty propagation regarding system
component predictions and performance. Our experiments show that our approach
for capturing and propagating uncertainty along the system can provide
valuable predictions and uncertainty estimates to build dependable systems.
However, a proper use of component predictions and uncertainty estimates are
needed to positively impact the system performance. In future work, we aim to
explore sampling-free methods for uncertainty estimation [35] to reduce the
computational budget and memory footprint in our approach.

## Acknowledgment

This work has received funding from the COMP4DRONES project, under Joint
Undertaking (JU) grant agreement N°826610. The JU receives support from the
European Union’s Horizon 2020 research and innovation programme and from
Spain, Austria, Belgium, Czech Republic, France, Italy, Latvia, Netherlands.

## References

  * [1] R. Siegwart, I. R. Nourbakhsh, and D. Scaramuzza, _Introduction to autonomous mobile robots_. MIT press, 2011\.
  * [2] R. McAllister, Y. Gal, A. Kendall, M. Van Der Wilk, A. Shah, R. Cipolla, and A. Weller, “Concrete problems for autonomous vehicle safety: Advantages of bayesian deep learning,” in _Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence_. International Joint Conferences on Artificial Intelligence, Inc., 2017.
  * [3] S. Grigorescu, B. Trasnea, T. Cocias, and G. Macesanu, “A survey of deep learning techniques for autonomous driving,” _Journal of Field Robotics_ , vol. 37, no. 3, pp. 362–386, 2020.
  * [4] A. Loquercio, A. I. Maqueda, C. R. Del-Blanco, and D. Scaramuzza, “Dronet: Learning to fly by driving,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 2, pp. 1088–1095, 2018.
  * [5] F. Codevilla, E. Santana, A. M. López, and A. Gaidon, “Exploring the limitations of behavior cloning for autonomous driving,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2019, pp. 9329–9338.
  * [6] W. Zeng, W. Luo, S. Suo, A. Sadat, B. Yang, S. Casas, and R. Urtasun, “End-to-end interpretable neural motion planner,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 8660–8669.
  * [7] R. Bonatti, R. Madaan, V. Vineet, S. Scherer, and A. Kapoor, “Learning visuomotor policies for aerial navigation using cross-modal representations,” _arXiv preprint arXiv:1909.06993_ , 2019.
  * [8] M. Mueller, A. Dosovitskiy, B. Ghanem, and V. Koltun, “Driving policy transfer via modularity and abstraction,” in _Conference on Robot Learning_. PMLR, 2018, pp. 1–15.
  * [9] R. Michelmore, M. Kwiatkowska, and Y. Gal, “Evaluating uncertainty quantification in end-to-end autonomous driving control,” _arXiv preprint arXiv:1811.06817_ , 2018.
  * [10] F. Nozarian, C. Müller, and P. Slusallek, “Uncertainty quantification and calibration of imitation learning policy in autonomous driving.” in _TAILOR_ , 2020, pp. 146–162.
  * [11] E. Ohn-Bar, A. Prakash, A. Behl, K. Chitta, and A. Geiger, “Learning situational driving,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 11 296–11 305.
  * [12] F. Arnez, H. Espinoza, A. Radermacher, and F. Terrier, “Improving robustness of deep neural networks for aerial navigation by incorporating input uncertainty,” in _International Conference on Computer Safety, Reliability, and Security_. Springer, 2021, pp. 219–225.
  * [13] B. Ivanovic, K.-H. Lee, P. Tokmakov, B. Wulfe, R. McAllister, A. Gaidon, and M. Pavone, “Heterogeneous-agent trajectory forecasting incorporating class uncertainty,” _arXiv preprint arXiv:2104.12446_ , 2021.
  * [14] S. Casas, C. Gulino, S. Suo, K. Luo, R. Liao, and R. Urtasun, “Implicit latent variable model for scene-consistent motion forecasting,” in _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XXIII 16_. Springer, 2020, pp. 624–641.
  * [15] A. Lavin, C. M. Gilligan-Lee, A. Visnjic, S. Ganju, D. Newman, S. Ganguly, D. Lange, A. G. Baydin, A. Sharma, A. Gibson _et al._ , “Technology readiness levels for machine learning systems,” _arXiv preprint arXiv:2101.03989_ , 2021.
  * [16] H. Rueß and S. Burton, “Safe ai–how is this possible?” _arXiv preprint arXiv:2201.10436_ , 2022.
  * [17] A. Spurr, J. Song, S. Park, and O. Hilliges, “Cross-modal deep variational hand pose estimation,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 89–98.
  * [18] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” _arXiv preprint arXiv:1312.6114_ , 2013.
  * [19] E. Daxberger and J. M. Hernández-Lobato, “Bayesian variational autoencoders for unsupervised out-of-distribution detection,” _arXiv preprint arXiv:1912.05651_ , 2019.
  * [20] A. Jesson, S. Mindermann, U. Shalit, and Y. Gal, “Identifying causal-effect inference failure with uncertainty-aware models,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 11 637–11 649, 2020.
  * [21] Y. Gal and Z. Ghahramani, “Dropout as a bayesian approximation: Representing model uncertainty in deep learning,” in _international conference on machine learning_ , 2016, pp. 1050–1059.
  * [22] S. Depeweg, J. Hernández-Lobato, F. Doshi-Velez, and S. Udluft, “Learning and policy search in stochastic dynamical systems with bayesian neural networks,” in _5th International Conference on Learning Representations, ICLR 2017-Conference Track Proceedings_ , 2017.
  * [23] S. Depeweg, J.-M. Hernandez-Lobato, F. Doshi-Velez, and S. Udluft, “Decomposition of uncertainty in bayesian deep learning for efficient and risk-sensitive learning,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 1184–1193.
  * [24] M. Henaff, A. Canziani, and Y. LeCun, “Model-predictive policy learning with uncertainty regularization for driving in dense traffic,” in _International Conference on Learning Representations_ , 2018.
  * [25] F. K. Gustafsson, M. Danelljan, and T. B. Schön, “Evaluating scalable bayesian deep learning methods for robust computer vision,” _arXiv preprint arXiv:1906.01620_ , 2019.
  * [26] A. Kendall and Y. Gal, “What uncertainties do we need in bayesian deep learning for computer vision?” in _Advances in neural information processing systems_ , 2017, pp. 5574–5584.
  * [27] B. Lakshminarayanan, A. Pritzel, and C. Blundell, “Simple and scalable predictive uncertainty estimation using deep ensembles,” in _Advances in neural information processing systems_ , 2017, pp. 6402–6413.
  * [28] M. Zhang, S. Levine, and C. Finn, “Memo: Test time robustness via adaptation and augmentation,” _arXiv preprint arXiv:2110.09506_ , 2021.
  * [29] Y. Gal, R. Islam, and Z. Ghahramani, “Deep bayesian active learning with image data,” in _International Conference on Machine Learning_. PMLR, 2017, pp. 1183–1192.
  * [30] A. Kirsch, J. Van Amersfoort, and Y. Gal, “Batchbald: Efficient and diverse batch acquisition for deep bayesian active learning,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [31] B. Poole, S. Ozair, A. Van Den Oord, A. Alemi, and G. Tucker, “On variational bounds of mutual information,” in _International Conference on Machine Learning_. PMLR, 2019, pp. 5171–5180.
  * [32] S. Shah, D. Dey, C. Lovett, and A. Kapoor, “Airsim: High-fidelity visual and physical simulation for autonomous vehicles,” in _Field and service robotics_. Springer, 2018, pp. 621–635.
  * [33] R. Madaan, N. Gyde, S. Vemprala, M. Brown, K. Nagami, T. Taubner, E. Cristofalo, D. Scaramuzza, M. Schwager, and A. Kapoor, “Airsim drone racing lab,” _arXiv preprint arXiv:2003.05654_ , 2020.
  * [34] Y. Song, M. Steinweg, E. Kaufmann, and D. Scaramuzza, “Autonomous drone racing with deep reinforcement learning,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2021, pp. 1205–1212.
  * [35] B. Charpentier, O. Borchert, D. Zügner, S. Geisler, and S. Günnemann, “Natural posterior network: Deep bayesian predictive uncertainty for exponential family distributions,” _arXiv preprint arXiv:2105.04471_ , 2021\.


short-format= TTA short = TTA, long = transport triggered architecture
VLIWshort = VLIW, long = very long instruction word FU short = FU, long =
functional unit RF short = RF, long = register file LSU short = LSU, long =
load-store unit GCU short = GCU, long = general control unit TCE short = TCE,
long = TTA-based Co-Design Environment TUT short = TUT, long = Tampere
University of Technology AG short = AG, long = address generator TFG short =
TFG, long = twiddle factor generator CADDshort = CADD, long = complex adder
FFT short = FFT, long = fast Fourier transform LSB short = LSB, long = least
significant bit MSB short = MSB, long = most significant bit LUT short = LUT,
long = lookup table ROM short = ROM, long = read-only memory OFDMshort = OFDM,
long = orthogonal frequency division multiplexing ASICshort = ASIC, long =
application specific integrated circuit FPGAshort = FPGA, long = field-
programmable gate array CGRAshort = CGRA, long = coarse-grained reconfigurable
array GPP short = GPP, long = general purpose processor GPU short = GPU, long
= graphics processung unit PE short = PE, long = processing element DIT short
= DIT, long = decimation-in-time DIF short = DIF, long = decimation-in-
frequency SNR short = SNR, long = signal-to-noise ratio

# Low-Power Programmable Processor for Fast Fourier Transform Based on
Transport Triggered Architecture

###### Abstract

This paper describes a low-power processor tailored for fast Fourier transform
computations where transport triggering template is exploited. The processor
is software-programmable while retaining an energy-efficiency comparable to
existing fixed-function implementations. The power savings are achieved by
compressing the computation kernel into one instruction word. The word is
stored in an instruction loop buffer, which is more power-efficient than
regular instruction memory storage. The processor supports all power-of-two
FFT sizes from 64 to 16384 and given 1 mJ of energy, it can compute 20916
transforms of size 1024.

Index Terms—  Fast Fourier Transform, Transport Triggered Architecture,
Application-Specific Instruction-Set Processor

##  1 Introduction

FFT is one of the most widely used signal processing algorithms thanks to its
ability to represent a time-domain signal in a frequency domain. For example,
FFT is used in OFDM systems, which are employed in wireless communication
devices. Due to the popularity of embedded and battery-powered systems,
minimizing power consumption is a major objective.

Regarding power consumption, ASIC implementations are considered as more
efficient compared to reconfigurable hardware (such as FPGA, CGRA) or GPP.
However, ASIC-based FFT processors are mostly fixed-function and lack
programmability. While the reconfigurable fabric offers more silicon
reusability than ASIC, their functionality can be only modified by a hardware
design process similar to ASIC. The goal of this work is to propose a
software-programmable mixed radix-4/2 FFT processor with an energy-efficiency
comparable to fixed-function ASIC implementations. Other software
implementations of FFT are implemented on either GPP [1] or a GPU [2].
However, both of these approaches aim for the best performance and do not
provide sufficiently low-power solutions.

Fixed-function ASIC FFT processors can be divided into two categories \-
pipelined and memory based. Pipelined architectures ([3], [4], [5]) rely on a
cascade of PE processing the input data stream. The intermediate results are
stored in a distributed memory system. Due to a higher number of PE, pipelined
architectures consume more power and occupy larger silicon area than memory
based architectures. However, they usually have higher throughput, which can
lead to a high energy-efficiency.

Memory based architectures ([6], [7]) typically have one PE and data is
processed in a sequential fashion. They typically use only one or two global
memory elements, thus a conflict-free memory access has to be maintained. The
proposed architecture is memory based, using a single-port data memories as it
allows for a convenient software-programmable implementation.

The processor was designed using a TTA template [8]. It improves a previous
FFT processor [9] by further increasing its energy-efficiency. Several
optimizations were applied to allow compressing the computation kernel into
only one - repeatedly executed - instruction word that can be executed in a
more energy-efficient way.

##  2 FFT Algorithm

The proposed processor supports all the power-of-two FFT sizes from $2^{6}$ to
$2^{14}$. A mixed radix-4/2 algorithm was used, following a DIT approach [10]
as it provides better SNR compared to DIF approach[11]. Otherwise, they share
the same arithmetic complexity. Radix-4 is used in a majority of the stages
because is requires less operations per FFT than the radix-2 algorithm [12].
At the same time, radix-4 butterfly operation requires only trivial
operations. Higher radices require more complicated operations. However, using
only radix-4 would restrict the processor to only power-of-four FFT sizes.
Therefore, in the last stage of the computation, radix-2 butterflies are used
for FFT sizes which can not be computed using radix-4 algorithm (i.e. for FFT
sizes $2^{k}$ where $k$ is odd). The computation follows an in-place approach
where output samples are written back to the same memory locations from which
the operands were read. This allows to utilize only one memory module of the
size equal to the computed FFT size.

##  3 Transport Triggered Architecture

TTA [8] is a processor template, which exposes its internal datapaths to a
programmer. Similarly to VLIW [13], it utilizes long instruction words and
instruction-level parallelism. The difference is that TTA gives a programmer
the control over the data flow. It is possible to bypass accesses to RF by
feeding results from one FU directly to the input of another. Register
bypassing reduces the required RF size and hardware complexity leading to
significant power savings [14].

The data transports are defined by a _move_ instruction - the only instruction
of the TTA’s instruction set. _Moving_ data into a _trigger_ port of a FU
triggers the desired operation. FU can also have _operand_ ports for
additional data that can be loaded anytime without triggering the operation.
Memory access is performed by LSU in a similar way as any other instruction. A
control unit responsible for instruction fetching, decoding, and executing is
also implemented as one of the FU. Data moves are distributed over an
interconnection network consisting of several parallel buses. The number of
the parallel buses determines the maximum number of instructions that can be
executed in parallel, i.e., the maximum number of simultaneous data moves.

##  4 Proposed Processor Architecture

The proposed architecture is shown in Fig. 1. The architecture consists of ten
32-bit wide buses (B0–B9) and one 1-bit bus (b), represented by horizontal
lines. FU and one RF are connected to the buses. Vertical lines represent
sockets, which connect input/output ports of FU to the interconnection
network. The connections are marked as dots.

Fig. 1:  Architecture of the proposed TTA processor. _TFG_ : twiddle factor
generator: _AG_ : address generator. _DLY_ : rotating register as a delay
unit. _CMUL_ : complex multiplier. _CADD_ : complex adder. _ADD_ : adder.
_LSU_ : load-store unit. _SH_ : shifter. _RF_ : 8x32 register file. _GCU_ :
general control unit. Gray color denotes unused connections. Trigger ports are
marked with ‘*’.

Two LSU are connected to a data memory system that behaves like a dual-port
memory. In fact, two single-port memories are used and connected to the LSU
via an added logic, which provides a conflict-free memory access. The parallel
memory system was chosen due to a lower power consumption of single-port
memories compared to multi-port memories [15].

The parity of the address determines which one of the two single-port memories
is accessed. In the case when both LSU are trying to access an address with
the same parity (i.e. the same memory module), the processor is temporarily
locked and the accesses are resolved sequentially. However, the conflict-free
memory access is guaranteed for the FFT addressing scheme.

The streamlined instruction schedule (see Section 6) implies generation of two
parallel streams of addresses - read and write. In order to guarantee a
different parity for any two parallel addresses (thus conflict-free memory
access), a special scheduler module was put between the LSU and the parallel
memory logic described in a previous paragraph. The scheduler internally
buffers and reschedules the LSU data in a way that always two parallel read
addresses or two parallel write addresses are loaded into the parallel memory
logic. Because the address generator preserves parity (see Section 5.1), the
scheduler guarantees a conflict-free memory access. The internal buffering is
not recognized by a high-level compiler and, therefore, the programming is
only possible by low-level assembly. However, it is possible to provide a
software-exposed switch in a form of another port of LSU or a special FU that
toggles the scheduler on and off, thus preserving a full compiler support for
generic applications.

Loop buffer [16] \- a critical component of the design - is implemented as a
part of the GCU. It is a small instruction memory cache used for storing
frequently repeated instruction words, e.g., loops. Reading from a loop buffer
consumes significantly less power than reading an instruction directly from
the instruction memory.

Each single-port data memory is composed of increasingly sized memory blocks
(32, 32, 64, 128, …, 4096 - summing up to total 8192). Based on the access
address, only one block is selected at a time while the other do not receive
any control signals. This significantly decreases dynamic power consumption
when computing smaller FFT sizes.

The processor was designed using TCE toolset developed at TUT [17]. TCE
provides a comprehensive set of tools for designing TTA processors including a
retargetable compiler and a hardware description of the most common FU. For
data and instruction memories, low power Cacti-P models were used [18].

##  5 Special Function Units

This section describes special FU developed specifically for this work. All
the other FU were taken from TCE component libraries.

Complex numbers are represented by two 16-bit fixed point numbers sharing one
32-bit data word. The real part occupies LSB of the data word while the
imaginary part takes the MSB.

In order to prevent overflow, each addition is divided by two. When summed up,
the complex adder divides the result by four in case of radix-4 and by two in
case of radix-2 butterfly. The complex multiplier divides the result by two.

###  5.1 Address Generator

The AG is responsible for computing the memory addresses for butterfly
operands. It is generated from a linear counter by a bit pair permutation
following the same pattern as the reference implementation [19]. An example of
an address generation for a 128-point and 256-point FFT is illustrated in Fig.
2. Each $b_{i}$ represents an $i$-th bit of a linear counter. The ‘index’ bits
are sufficient to represent the index within one stage while ‘stage’
determines the current stage of the computation. The position of the LSB bit
pair is determined by the ‘stage’ part of the linear counter.

The address generator preserves the parity of the linear counter. Thus, any
two consecutive addresses have a different parity and if fed in parallel into
the parallel memory logic (described in 4), a conflict-free memory access is
guaranteed.

Fig. 2:  Address generation from a linear counter for 128-point (above) and
256-point (below) FFT.

###  5.2 Twiddle Factor Generator

The generation of twiddle factors is based on a LUT implemented as a single-
port synchronous ROM of pre-computed values. It follows the same approach as
the one described in [20]. The address for the LUT ROM is computed from the
linear index by a bit permutation and scaling based on the current FFT size.
Only $N/8+1$ complex coefficients need to be stored in the LUT [20]. All the
remaining coefficients can be reconstructed by a trivial manipulation
(negating and swapping the real and imaginary parts) of the stored
coefficients. Therefore, in order to support the maximum 16384-point FFT, the
LUT has to contain 2049 coefficients. A side function of the TFG FU is
determining whether the current stage is radix-4 or radix-2. This information
is then used by the CADD.

rx2 | cnt | result
---|---|---
0 | 00 | a + b + c + d
0 | 01 | a - i*b + c + i*d
0 | 10 | a - b + c - d
0 | 11 | a + i*b - c - i*d
1 | 00 | a + b
1 | 01 | a - b
1 | 10 | c + d
1 | 11 | c - d
Fig. 3:  An operation performed by a complex adder based on the value of its
‘rx2’ input and an internal counter (‘cnt’). Four operands (a, b, c, d) and
rx2 are constant until the next reset of the counter. i denotes an imaginary
unit.

###  5.3 Complex Adder

The CADD performs a butterfly operation on four inputs. Based on its ‘rx2’
input, it performs either one radix-4 or two radix-2 butterflies.

Traditionally, the CADD would be implemented as a four-input FU with the four
inputs buffered in register files before feeding them in parallel into the
CADD’s ports. However, due to the single-instruction kernel requirement, the
register file buffering is not possible since the data can be moved only to a
single location. Therefore, the proposed CADD FU has one serial data input
port and performs the buffering internally. This makes the FU unusable for
high-level programming since this mode of operation can not be recognized by a
high-level compiler.

Figure 3 shows the CADD’s results based on its ‘rx2’ input. The ‘cnt’ column
is an internal counter that increments each time a data sample is loaded into
the FU’s _trigger_ port. Both signals form an opcode selecting the operation
of the complex adder.

###  5.4 Complex Multiplier

The complex multiplier performs generic complex multiplication of two
operands. The proposed implementation requires four multipliers and two
adders.

###  5.5 Rotating Register

Rotating register is used to delay the address of a butterfly’s input sample
for the in-place computation. After the butterfly operation is complete, the
output of the rotating register is used as an address for the results to store
them back to the memory.

##  6 Instruction Schedule

| type | tech. | volt. | freq | WL | t | power | FFT/mJ | FFT/mJ norm. | programmable
---|---|---|---|---|---|---|---|---|---|---
|  | (nm) | (V) | (MHz) | (bits) |  ($\mu$s) | (mW) |  |  |
[4] | pipelined | 65 | 1.10 | 50 | 16 | 21.5 | 17.60 | 2641 | 3196 | no
proposed | memory based | 28 | 0.60 | 450 | 16 | 11.4 | 4.19 | 20916 | 3243 | yes
[21] | memory based | 65 | 1.20 | 500 | 16 | 2.6 | 170.0 | 2287 | 3292 | no
[9] | memory based | 130 | 1.50 | 250 | 16 | 20.6 | 60.40 | 802 | 3609 | yes
[6] | memory based | 600 | 3.30 | 173 | 20 | 30.0 | 845.0 | 39 | 6058 | no
proposed | memory based | 65 | 1.00 | 450 | 16 | 11.4 | 12.21 | 7171 | 7171 | yes
[7] | memory based | 90 | 1.00 | 160 | 16 | 2.6 | 29.00 | 13360 | 18498 | no
[5] | pipelined | 55 | 0.90 | 18 | 16 | 1.5 | 8.88 | 77131 | 52865 | no
Table 1: Comparison of various FFT processor architectures (1024-point FFT).

The computation of one radix-4 butterfly can be visualized with the aid of a
reservation table in Fig. 4. Each column represents one clock cycle. Buses are
represented by rows and their names (on the right) correspond to the ones
shown in Fig. 1. Gray square denotes that an instruction, i.e., data transfer,
is executed on the bus during the clock cycle. The instruction (data _move_)
transferred on each bus is shown on the left. The syntax respects the
following pattern: source.port $\rightarrow$ destination.port. Source and
destination are FU. Port can be t (trigger), o (operand), r (result) and rx2
(output port of TFG signalizing whether the butterfly is radix-4 or radix-2).

Fig. 4: Bus reservation table of computing one radix-4 butterfly

Full FFT is computed by repeating the above pattern multiple times every four
clock cycles. At 13 clock cycle, the bus utilization reaches 100% and the
instruction word becomes constant until no new samples need to be computed.
Thus, the execution can be separated into three stages: prologue (first 13
cycles), kernel (length depends on FFT size) and epilogue (last 13 cycles).
The size of the prologue and epilogue is constant for all FFT sizes. Because
the kernel consists of only one repeated instruction word, it can be loaded
into the loop buffer from where it can be fetched consuming minimal power.

Apart from the prologue, kernel and epilogue, a setup code consisting of 6
instructions is present to distribute static parameters between FU. Thus, the
size of the complete code is 33 ($6+13+1+13$) instructions. The architecture
uses 51-bit wide instruction words.

##  7 Evaluation

The processor was synthesized using Synopsys Design Compiler and two IC
technologies were used - a 28 nm FDSOI low-power technology and another 65 nm
technology.

In order to be able to compare different technologies, the energy was
normalized according to the following formula [6], [22]:

$E_{n}=E\frac{L_{ref}U^{2}_{ref}\left(\frac{1}{3}W^{2}_{ref}+\frac{2}{3}W_{ref}\right)}{LU^{2}\left(\frac{1}{3}W^{2}+\frac{2}{3}W\right)},$
(1)

where $E_{n}$ is the normalized energy; $E$, $L$, $U$ and $W$ are parameters
of the proposed architecture (energy, technology size, voltage and word
length, respectively); the ref suffix marks the reference technology (65 nm,
1.0 V, 16 bits).

Table 1 compares the proposed architecture with selected state-of-the-art
solutions and traditional architectures. The chosen focus point is a
1024-point FFT as a mid-point between the smallest and largest supported FFT
sizes. The frequency 450 MHz is close to the maximum achievable frequency (500
MHz for 28nm/0.60V and 550 for 65nm/1.00V). The maximum achievable frequency
is 1150 MHz with 28nm/1.10V technology.

##  8 Conclusion

In this paper, a low-power software-programmable FFT processor was proposed,
which is is based on a TTA template. The key contribution is reducing the
computation kernel into only one repeated instruction word and executing it
from a loop buffer instead of fetching from an instruction memory every clock
cycle. This reduces the power consumption of an instruction memory to a
negligible value. In order to achieve the instruction word compression, an
internal buffering was introduced in case of a complex adder and a memory
access, which renders them unusable by a high level language compiler.
However, it is possible to provide a software-accessible switch to disable the
memory buffering for more generic applications. Additional functionality can
be introduced by adding other functional units. Synthesis power evaluation
performed at two different ASIC technologies (28 nm and 65 nm) shows that the
processor can provide an energy efficiency comparable with fixed-function ASIC
processors.

##  9 Acknowledgments

The authors thank the following sources of financial support: Tampere
University of Technology Graduate School, Business Finland (FiDiPro Program
funding decision 40142/14), and ECSEL JU project FitOptiVis (project number
783162).

## References

  * [1] M. Khelifi, D. Massicotte, and Y. Savaria,  “Parallel independent FFT implementation on intel processors and Xeon phi for LTE and OFDM systems,”  in IEEE Nordic Circ. Syst. Conf. & Intl. Symp. SoC, 2015, pp. 1–4.
  * [2] X. Lyu, J. Zuo, and H. Xie,  “Non-equispaced FFT computation with CUDA and GPU,”  in Intl. Conf. on Virtual Reality and Visualization (ICVRV), 2016, pp. 227–234.
  * [3] T. H. Tran, S. Kanagawa, D. P. Nguyen, and Y. Nakashima,  “ASIC design of MUL-RED radix-2 pipeline FFT circuit for 802.11ah system,”  in Proc. IEEE Low-Power and High-Speed Chips Symp., 2016, pp. 1–3.
  * [4] M. Garrido, R. Andersson, F. Qureshi, and O. Gustafsson,  “Multiplierless unity-gain SDF FFTs,”  IEEE T. Very Large Scale Integration Syst., vol. 24, no. 9, pp. 3003–3007, 2016.
  * [5] M. Garrido, S. Huang, and S. Chen,  “Feedforward FFT hardware architectures based on rotator allocation,”  IEEE T. Circ. Syst. I: Regular Papers, vol. 65, no. 2, pp. 581–592, 2018.
  * [6] B. M. Baas,  “A low-power, high-performance, 1024-point FFT processor,”  IEEE J. Solid-State Circuits, vol. 34, no. 3, pp. 380–387, 1999\.
  * [7] S. Huang and S. Chen,  “A high-parallelism memory-based FFT processor with high SQNR and novel addressing scheme,”  in Proc. IEEE ISCAS, 2016, pp. 2671–2674.
  * [8] H. Corporaal,  Microprocessor Architectures: From VLIW to TTA,  John Wiley & Sons, Inc., 1997.
  * [9] T. Pitkänen and J. Takala,  “Low-power application-specific processor for FFT computations,”  J. Signal Process. Syst., vol. 63, no. 1, pp. 165–176, 2011.
  * [10] A. Saidi,  “Decimation-in-time-frequency FFT algorithm,”  in Proc. IEEE ICASSP, 1994, vol. 3, pp. III–453.
  * [11] W. H. Chang and T. Q. Nguyen,  “On the fixed-point accuracy analysis of FFT algorithms,”  IEEE T. Signal Processing, vol. 56, no. 10, pp. 4673–4682, 2008\.
  * [12] T. Pitkänen,  Fast Fourier Transforms on Energy-Efficient Application-Specific Processors,  Ph.D. thesis, Tampere University of Technology, Finland, 2014.
  * [13] J. A. Fisher,  “Very long instruction word architectures and the ELI-512,”  in Proc. of the 10th Annual Intl. Symp. on Computer Arch., Stockholm, Sweden, 1983, pp. 140–150.
  * [14] P. Jääskeläinen, H. Kultala, T. Viitanen, and J. Takala,  “Code density and energy efficiency of exposed datapath architectures,”  J. Signal Process. Syst., vol. 80, no. 1, pp. 49–64, 2015.
  * [15] T. Pitkänen, J. K. Tanskanen, R. Mäkinen, and J. Takala,  “Parallel memory architecture for application-specific instruction-set processors,”  J. Signal Process. Syst., vol. 57, no. 1, pp. 21–32, 2009.
  * [16] V. Guzma, T. Pitkänen, and J. Takala,  “Reducing instruction memory energy consumption by using instruction buffer and after scheduling analysis,”  in Intl. Symp. on SoC Proc., 2010, pp. 99–102.
  * [17] P. Jääskeläinen, T. Viitanen, J. Takala, and H. Berg,  “HW/SW co-design toolset for customization of exposed datapath processors,”  in Computing Platforms for Software-Defined Radio, Waqar Hussain, Jari Nurmi, Jouni Isoaho, and Fabio Garzia, Eds., pp. 147–164. Springer, 2017.
  * [18] S. Li, K. Chen, J. H. Ahn, J. B. Brockman, and N. P. Jouppi,  “CACTI-P: Architecture-level modeling for SRAM-based structures with advanced leakage reduction techniques,”  in Proc. Intl. Conf. on Computer-Aided Design, 2011, pp. 694–701.
  * [19] T. Pitkänen, R. Mäkinen, J. Heikkinen, T. Partanen, and J. Takala,  “Low-power, high-performance TTA processor for 1024-point fast Fourier transform,”  in Embedded Computer Systems: Architectures, Modeling, and Simulation, Berlin, Heidelberg, 2006, pp. 227–236, Springer.
  * [20] T. Pitkänen, T. Partanen, and J. Takala,  “Low-power twiddle factor unit for FFT computation,”  in Embedded Computer Systems: Architectures, Modeling, and Simulation, Berlin, Heidelberg, 2007, pp. 65–74, Springer.
  * [21] M. A. Shami, M. A. Tajammul, and A. Hemani,  “Configurable FFT processor using dynamically reconfigurable resource arrays,”  J. Signal Process. Syst., 2018.
  * [22] Y. Chen, Y. W. Lin, Y. C. Tsao, and C. Y. Lee,  “A 2.4-gsample/s DVFS FFT processor for MIMO OFDM communication systems,”  IEEE Journal of Solid-State Circuits, vol. 43, no. 5, pp. 1260–1273, 2008.


# Properties of stable ensembles of Euclidean random matrices

Philipp Baumgärtel  Florian Vogel  Matthias Fuchs  Fachbereich Physik,
Universität Konstanz, 78457 Konstanz, Germany

###### Abstract

We study the spectrum of a system of coupled disordered harmonic oscillators
in the thermodynamic limit. This Euclidean random matrix ensemble has been
suggested as model for the low-temperature vibrational properties of glass.
Exact numerical diagonalization is performed in three and two spatial
dimensions. It reveals a low-frequency regime of sound waves that are damped
by Rayleigh scattering. At large frequencies localized modes exist. In
between, the central peak in the vibrational density of states is well
described by Wigner’s semicircle law expected for random matrix systems. We
compare our results with predictions from two recent self-consistent field
theories.

##  I Introduction

The nature of the vibrational excitations in athermal amorphous solids remains
an open and important topic affecting inter alia thermal properties of glasses
at low temperatures [1]. The vibrational properties of glassy materials differ
strongly from the ones of crystals [2] as is well established by scattering
experiments [3, 4, 5, 6, 7]. Computer simulations of various particle models
have shown that the preparation of the glass state matters and a number of
diverse phenomena have been discovered [8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19]. As a consistent comprehensive theoretical picture is still lacking
[20], a simple idealized model appears desirable, where a part of the
phenomena can be studied in detail.

In 1999, Mezard, Parisi, and Zee introduced ensembles of Euclidean random
matrices (ERM) [21], which were studied as simple models for low temperature
glasses by Grigera and colleagues [22, 23, 24, 25], and by Schirmacher and
colleagues [26, 27, 28, 29]. Random matrix ensembles have successfully been
employed in a wide variety of physical systems with disorder [30, 31, 32], and
the special ERM ensemble may arguably be considered the most idealized model
for the vibrations in glass. Particles perform harmonic motion around random
positions. The restoring forces depend on the distances between the positions
via a positive spring function, and translational invariance is postulated.
Only properties averaged over the random positions are studied. Based on
diagrammatic perturbation expansions in field theoretic approaches, self-
consistent theories for the central Greens functions have been developed [23,
24, 26, 27, 33]. They allow the (approximate) calculation of the vibrational
density of states (vDOS) and of the dynamical structure factor.

While ERM had been studied intensively up to a decade ago, open questions
remained on the sound damping and on the spatial characteristics of the
eigenmodes. In the present contribution, we investigate the most simple ERM
model, already considered in [23, 27], an homogeneous and isotropic system
with a positive Gaussian spring function. The only state parameter turns out
to be the rescaled density $n$, which encodes the amount of disorder. We
perform large scale numerical investigations including studies of finite size
corrections, in order to reveal the complete characteristics of this specific
ERM ensemble. We also compare with predictions from two self-consistent
theories [23, 33], where different series of diagrams in the perturbation
expansion [25] were re-summed.

##  II Model

In this work we study an Euclidean random matrix (ERM) model [21, 31, 29] in
which we consider $N$ particles which are randomly placed in a box of Volume
$V=L^{d}$. Here, $d$ denotes the dimension of the system, and $L$ is its
length. We apply periodic boundary conditions to the system and study an
uniform distribution of particles. The set of random positions
$\\{\pmb{r}_{i}\\}$ will be called inherent positions. We consider a harmonic
motion of the particles around their inherent positions which leads us to the
definition of a random matrix $\pmb{M}$ via the interaction potential $U$

$\displaystyle U(\\{\phi_{i}\\})$
$\displaystyle=\frac{1}{4}\sum_{i,j}f(\pmb{r}_{i}-\pmb{r}_{j})(\phi_{i}-\phi_{j})^{2}$
$\displaystyle=\frac{1}{2}\sum_{i,j}M_{ij}\phi_{i}\phi_{j},$ (1)

with

$\displaystyle
M_{ij}=-f(\pmb{r}_{i}-\pmb{r}_{j})+\delta_{ij}\sum_{k}f(\pmb{r}_{i}-\pmb{r_{k}}).$
(2)

Here, $f(\pmb{r})$ is called spring function and $\phi_{i}$ is a small scalar
displacement of particle $i$ from its inherent site. In this work we consider
the simple case, where the spring function is isotropic and given by the
Gaussian

$\displaystyle f(r)=\exp(-r^{2}/2),$ (3)

with $r$ the dimensionless distance. For positive spring functions, $f(r)>0$,
the potential $U$ is positive and thus the matrix $\bf M$ is positive semi-
definite.
In the limit of large systems, a single state parameter $n=N/V$ determines the
properties of quantities averaged over the disorder.

In the harmonic approximation, the equations of motion of the system are given
by

$\displaystyle\ddot{\phi}_{i}=-\sum_{j=1}^{N}M_{ij}\phi_{j}\quad,\;\mbox{for
}1\leq i\leq N$ (4)

Translational invariance and hence momentum conservation follow immediately
from the potential $U(\\{\phi_{i}\\})$. Consequently, $\bf M$ has the
eigenvalue zero. The associated eigenvector $\pmb{e}^{0}$ corresponds to the
uniform shift $\pmb{e}^{0}=(1,1,....,1)/\sqrt{N}$.
Note that we already set the length scale of our system to one by the
definition of the spring function in Eq. (3), and frequency and time are also
chosen dimensionless quantities in Eq. (4).

##  III Methods

We use two methods to study the characteristics of the system. The first one,
in which we diagonalize the random matrix $M_{ij}$, will be called normal mode
analysis. The second one, where we solve the equations of motion, will be
called excited wave analysis. In both cases, averages over the disorder are
finally performed by sampling different inherent positions.

###  III.1 Normal mode analysis

In the normal mode analysis [12, 8, 3, 10, 11, 13, 14], we calculate the
eigenvalues $\lambda^{k}$, corresponding to the eigenfrequencies
$\omega^{k}=\sqrt{\lambda^{k}}$, and the eigenvectors $\pmb{e}^{k}$ of the
random matrix $\pmb{M}$. For this, we use the standard diagonalization routine
of matlab and a routine called lobpcg [34] which can handle sparse matrices
efficiently. Note, that the symmetry and semi-positivity of $\pmb{M}$ assures
$\lambda^{k}\geq 0$ and that the $\pmb{e}^{k}$ form an orthonormal basis.

The density of states per particle in the energy domain is calculated by [35]

$\displaystyle
g_{\lambda}(\lambda)=\frac{1}{N}\overline{\sum_{k}\delta(\lambda-\lambda^{k})}$
(5)

and can be transformed into the frequency domain with $\lambda=\omega^{2}$
leading to

$\displaystyle g(\omega)=2\omega\;g_{\lambda}(\lambda(\omega)).$ (6)

Here, the overbar denotes the average over disorder. As we expect discrete
eigenfrequencies in finite systems, the density of states $g(\omega)$ can be
subject to inaccuracies occurring due to the binning process if the bin size
is chosen incorrectly. A quantity which resolves this issue is called
integrated density of states [36, 37, 38] and can be calculated by

$\displaystyle
I(\omega)=\int_{0}^{\omega}g(\omega^{\prime})\;\mathrm{d}\omega^{\prime}.$ (7)

The integrated density of states counts the number of eigenfrequencies up to a
frequency $\omega$.

We calculate the dynamical structure factor by
$S(q,\omega)=2\omega\,S_{\lambda}(p,\lambda)$, where

$\displaystyle S_{\lambda}(q,\lambda)$
$\displaystyle=\overline{\sum_{k}Q^{k}(q)\delta(\lambda-\lambda^{k})},$ (8)
$\displaystyle Q^{k}(q)$
$\displaystyle=\frac{1}{N}\left|\sum_{j}e^{k}_{j}\;\exp(i\,qx_{j})\right|^{2}.$
(9)

Here, we have assumed an excitation along the $x$-axis with
$\pmb{q}=q\,\pmb{\hat{e}}_{x}$. Note that throughout all analyses we only
allow discrete wavevectors $q_{l}=l\,2\pi/L$ with $l=\pm 1,\pm 2,\ldots$,
satisfying the periodic boundary conditions. The dynamic structure factor can
be used to extract the dispersion relation $\Omega(q)$ and the damping
$\Gamma(q)$ by fitting $S(q,\omega)$ to a damped harmonic oscillator model
[39, 40]

$\displaystyle
S(q,\omega)\propto\frac{\Omega^{2}(q)\Gamma(q)}{(\omega^{2}-\Omega^{2}(q))^{2}+\omega^{2}\Gamma^{2}(q)}$
(10)

Another quantity which is used to characterise the eigenmodes of the systems
is the participation ratio [10, 11, 13, 15]

$\displaystyle P^{k}=\frac{1}{N}\frac{1}{\sum_{i}(e_{i}^{k}\,e_{i}^{k})^{2}}.$
(11)

The participation ratio $P^{k}=1/N\ll 1$ is small for an ideal localized mode
involving only one particle, while $P^{k}=2/3$ for an ideal plane wave.

###  III.2 Excited wave analysis

If the eigenvectors $\pmb{e}^{k}$ of the system are known, we can also easily
solve the equations of motion by

$\displaystyle\phi_{i}(t)=\sum_{k}u^{k}(t)e_{i}^{k},$ (12)

where

$\displaystyle
u^{k}(t)=\pmb{e}^{k}\cdot\pmb{\phi}(0)\cos(\omega^{k}t)+\pmb{e}^{k}\cdot\dot{\pmb{\phi}}(0)\frac{\sin(\omega^{k}t)}{\omega^{k}}.$
(13)

The scalar product abbreviates the sum over particles,
$\pmb{e}^{k}\cdot\pmb{a}=\sum_{i}e^{k}_{i}\,a_{i}$. We can excite a standing
wave with wavenumber $q$ in our inherent structure by choosing for example
$\phi_{i}(0)=0$ and $\dot{\phi}_{i}(0)=~{}\sin(qx_{i}+\Phi)$ as initial
conditions [41, 12, 14, 9], with $\Phi=0,\pi/2$. Equivalent to fitting the
damped harmonic oscillator model to the dynamic structure factor in the
frequency domain, we can calculate the correlation function

$\displaystyle C(q,t)$ $\displaystyle=\overline{R(q,t)},$ (14) $\displaystyle
R(q,t)$
$\displaystyle=\frac{\sum_{i}\dot{\phi}_{i}(0)\,\dot{\phi}_{i}(t)}{\sum_{i}\dot{\phi}_{i}(0)\,\dot{\phi}_{i}(0)}\,$
(15)

and fit it with

$\displaystyle C(q,t)=\exp(-\Gamma(q)\,t/2)\cos(\Omega(q)t).$ (16)

Note that the average now also includes both phases $\Phi=0,\pi/2$ for each
set of inherent positions. This again allows us to extract the dispersion
relation $\Omega(q)$ and the attenuation $\Gamma(q)$. While the dynamic
structure factor includes a binning process in its calculation, the
correlation function does not, yielding better results for small wavevectors
$q$.

We can rewrite the correlation function in terms of ”hybridization”
coefficients

$\displaystyle\xi^{k}(q)=\frac{\pmb{e}^{k}\cdot\dot{\pmb{\phi}}(0)}{\sqrt{\dot{\pmb{\phi}}(0)\cdot\dot{\pmb{\phi}}(0)}}$
(17)

by

$\displaystyle R(q,t)$
$\displaystyle=\frac{\sum_{k}\pmb{e}^{k}\cdot\dot{\pmb{\phi}}(0)\cos(\omega^{k}t)\pmb{e}^{k}\cdot\dot{\pmb{\phi}}(0)}{\dot{\pmb{\phi}}(0)\cdot\dot{\pmb{\phi}}(0)}$
(18) $\displaystyle=\sum_{k}(\xi^{k}(q))^{2}\cos(\omega^{k}t),$ (19)

where $\sum_{k}(\xi^{k}(q))^{2}=1$ must hold.

##  IV Numerical details

We study systems with periodic boundary conditions. In our case this means
that the periodic boundaries influence the calculation of the dynamical matrix
$\vec{M}$. A particle at the boundary of the simulation box also interacts
with the particles of the periodic copies of the simulation box. This enforces
the selection of wavevectors $q_{l}$ introduced above.
In order to use sparse matrices (storing less then $N^{2}$ entries) for the
large systems, we introduce a cut-off radius $\sigma$ at which the spring
function $f(r)$ is truncated. We choose $\sigma=4$. See below for a discussion
of the accruing errors.
In order to obtain the randomly generated inherent structures we use the
standard random number generator of matlab which is the Mersenne Twister
algorithm [42].
The main simulation parameter we will vary is the dimensionless density $n$.
At a given number of particles $N$ the density determines the size of the
periodic simulation box via $L=(N/n)^{1/d}$. Hence, if we increase the number
of particles $N$ we also increase the size of the simulation box $L$. At the
same time we also get access to smaller wave vectors $q$.

We calculate the full set of eigenvalues for system sizes up to $N=4\times
10^{4}$ particles and the smallest $2000$ eigenvalues for larger systems. This
drastically reduces the computation time and especially the storage
consumption. See below for a discussion of the accruing errors.

For each density $n$ and system size $N$ we perform calculations in $250$
realizations of the random inherent positions and present the averages. When
calculating the 3D participation ratios, we tested an ensemble size of
$5\times 10^{5}$ simulations and found qualitatively the same result.

##  V Results and discussions

First we show the results obtained via the normal mode analysis and afterwards
we show the excited wave analysis.
Fig. 1a) shows the density of states for systems of size $N=4\times 10^{4}$
for varying densities $n$. As is often done, it is divided by the Debye
behaviour
$g_{\mathrm{D}}(\omega)=A_{\mathrm{D}}\,\omega^{2}=\omega^{2}/\omega_{\mathrm{D}}^{3}$.
Here $A_{\mathrm{D}}$ is the Debye level and $\omega_{\mathrm{D}}$ the Debye
frequency.

Figure 1: (a) Density of states $g(\omega)$ and (b) integrated density of
states $I(\omega)$ divided by the Debye behaviour for different densities $n$
(see legend) and as function of the rescaled frequency. The dashed lines
indicate the Debye levels $nA_{\mathrm{D}}/3$.

We can observe a Debye spectrum for $\omega\to 0$. It is expected because of
the breaking of translational invariance in the solid inherent structures
which leads to the existence of sound waves for small wavevector. For lowering
$n$, viz. increasing average separation of the inherent positions, the elastic
restoring forces become weaker and the Debye-frequency decreases. We observe
an excess in the vDOS above the Debye-law for larger frequencies $\omega$,
which we call the boson peak of the ERM. This interpretation will be discussed
in the Sect. VI. The location of the boson peak $\omega_{\mathrm{BP}}$ scales
with $\sqrt{n}$. Fig. 1b) shows the integrated density of states $I(\omega)$
divided by the Debye behaviour for the same densities $n$. The dashed lines
indicate the Debye levels $A_{\mathrm{D}}/3$, which can most directly be
extracted from these rescaled $I(\omega)$ data.

Figure 2: Density of states $g(\omega)$ for different densities $n$ (see
legend) plotted as function of $\omega-\omega_{\mathrm{BP}}$. The dashed line
shows the predicted Wigner semicircle law [33].

Fig. 2 shows the vDOS presented as function of $\omega-\omega_{\mathrm{BP}}$.
In this plot the curves coincide quite well [43], highlighting that the
majority of eigenfrequencies is accumulated in the boson peak. The overall
shape of the vDOS is well described by Wigner’s semicircle law, which would
hold if the entries of $\bf M$ were independent and identically distributed
zero-mean Gaussian entries [31]. The semicircle law is shown by the dashed
lines in Fig. 2; note that for simplicity its normalization to unity is not
adjusted to fit the data best. The semicircle is (in energy space) located
around $\omega_{\mathrm{BP}}^{2}$ and has a radius of $R=n\sqrt{2a}$, where
$a=\hat{f}(0)/(\sqrt{8}n)$. In frequency space the explicit form is given by
[33]

$\displaystyle g_{\mathrm{wigner}}(\omega)=\frac{4\,\omega}{\pi
R^{2}}\sqrt{R^{2}-(\omega^{2}-\omega_{\mathrm{BP}}^{2})^{2}}.$ (20)

Note that
$g_{\mathrm{wigner}}(\omega_{\mathrm{BP}})=4\,\omega_{\mathrm{BP}}/(\pi
R)=4/\pi\,2^{1/4}\propto n^{0}$. A shift of the boson peak frequency is the
dominant effect when changing the disorder via changing $n$. The plot hides
the $n$-dependence of the vDOS for small frequencies, as the Debye-level is
very low compared to the boson peak amplitude for the considered $n$.

From $g(\omega)$ and $I(\omega)$ we extract the density dependence of the
relevant frequencies, $\omega_{\mathrm{D}}$ and $\omega_{\mathrm{BP}}$ shown
in Fig. 3. The uncertainty of $\omega_{\mathrm{D}}$ is calculated from the
confidence interval of the fit of the Debye-level $A_{\mathrm{D}}$ and the
uncertainty of $\omega_{\mathrm{BP}}$ is estimated by evaluating the
frequencies of the bins left and right of the maximum of the vDOS. Since the
uncertainties in $\omega_{\mathrm{BP}}$ are very small they are omitted
throughout this work. The energy scale of the boson peak arises from the
pairwise interaction among all particles. This explains the scaling
$\omega_{\mathrm{BP}}\propto n^{1/2}$. For high enough $n$, the boson peak
position is the square-root of the mean value of the diagonal entries
$\omega_{\mathrm{BP}}^{2}=\overline{M_{ii}}=n\hat{f}(0)$ [33].

Figure 3:  Rescaled Debye frequency
$\hat{\omega}_{\mathrm{D}}=\omega_{\mathrm{D}}\,n^{-5/6}$ (black circles) and
boson peak frequency
$\hat{\omega}_{\mathrm{BP}}=\omega_{\mathrm{BP}}\,n^{-1/2}$ (black crosses) as
functions of the density $n$. The $\hat{\omega}_{\mathrm{D}}$ and
$\hat{\omega}_{\mathrm{BP}}$ are $n$-independent in stable states.
Additionally, the ratio $\hat{\omega}_{\mathrm{BP}}/\hat{\omega}_{\mathrm{D}}$
extracted from the simulation data (red diamonds) is shown and compared to
theory (line) from Ref. [33].

Considering the amplitude of the Debye-law, $A_{D}=1/\omega_{D}^{3}$, we
observe that $A_{D}\propto n^{-5/2}$ and thus $\omega_{\mathrm{D}}$ scales
with $n^{5/6}$. The increase of the Debye frequency with decreasing $n$ is a
non-trivial effect of the increasing disorder, which will be explained based
on the dispersion relations shown in Fig. 5 below. While the obtained boson
peak positions $\omega_{\mathrm{BP}}$ perfectly match with the values from the
theory of [33], we find small deviations in the Debye frequency
$\omega_{\mathrm{D}}$ and therefore also in the ratios
$\omega_{\mathrm{BP}}/\omega_{\mathrm{D}}$. Still, the scaling of
$\omega_{\mathrm{D}}$ and $\omega_{\mathrm{BP}}/\omega_{\mathrm{D}}$ with the
density $n$ matches well with the theory.

Fig. 4 shows the dynamic structure factors $S(q,\omega)$ at $n=1.0$ for 4
different $q$ values. The structure factors are characterised by a pronounced
peak which shifts to higher frequencies $\omega$ with increasing wave vector
$q$. At the same time the peak broadens and its height decreases. In the limit
of $q\to\infty$, $S(q,\omega)$ approaches the vDOS $g(\omega)$ shown by the
dashed line.

Figure 4: Dynamic structure factor $S(q,\omega)$ for $q_{4},q_{6},q_{8}$ and
$q_{25}$ (marked in Fig. 5; $q$ increases from left to right) at $n=1.0$. The
density of states is shown by the black dashed line, which agrees with
$S(q_{25},\omega)$.

We use Eq. (10) to extract the dispersion relation $\Omega(q)$ from the
structure factor. The results are shown in Fig. 5a). We observe a linear
dispersion relation $\Omega(q)=c_{\mathrm{T}}q$ for small wave vectors and
that $\Omega(q)\to\omega_{BP}$ saturates for $q\to\infty$. Here,
$c_{\mathrm{T}}$ denotes the speed of sound, it inherits the scaling
$c_{T}\propto\sqrt{n}$ from the dispersion relation. For increasing density
$n$ the slope of the dispersion relation, and therefore, the speed of sound
$c_{\mathrm{T}}$ gets larger; the system gets stiffer with the average
separation of particles getting smaller. The dispersion relation obtained via
the excited wave analysis is shown exemplary for $n=1.0$ by the purple circles
and agrees perfectly with the dispersion relation from the structure factor.
We compare the dispersion relation with the bare dispersion,
$\Omega_{0}(q)=\sqrt{\epsilon_{0}(q)n}=\sqrt{(\hat{f}(0)-\hat{f}(q))n}$, at
the smallest density $n=0.5$. The largest deviations from the bare dispersion
should be observed at the biggest disorder, i.e. the smallest density. Yet, at
the largest simulated disorder we still observe a good agreement between the
bare dispersion and $\Omega(q)$. Comparing with the self-consistent theories
[24, 33], one notices that both overestimate the effects of disorder; see Fig.
1 in Ref. 33, where dressed and bare dispersion relations are shown. Both
predict a stronger change of $\Omega(q)$ than observed, yet the more elaborate
resummation including non-planar diagrams lies closer to the data.
In Fig. 5b) we show the rescaled dispersion relation
$\Omega(q)\,\omega_{\mathrm{BP}}^{-1}$ as function of the rescaled wave vector
$q/q_{BP}$ with $q_{BP}=\omega_{\mathrm{BP}}/\,c_{\mathrm{T}}$. The rescaled
data collapse very well indicating that the boson peak frequency
$\omega_{\mathrm{BP}}$ is the relevant frequency for the rescaling. Using the
Debye-frequency does not lead to a comparable rescaling (not shown).

Figure 5:  (a) Dispersion relation $\Omega(q)$ for varying densities $n$ (see
legend in Fig. 2). The black arrows at $n=1.0$ mark the $q$ values at which
the dynamic structure factors in Fig. 4 are shown. The black circles at
$n=0.5$ indicate the bare dispersion
$\sqrt{\epsilon_{0}(q)n}=\sqrt{(\hat{f}(0)-\hat{f}(q))n}$ and the purple
diamonds show the dispersion relation obtained from the excited wave analysis.
(b) Rescaled dispersion relation $\Omega(q)\,\omega_{\mathrm{BP}}^{-1}$ in
dependency of the rescaled wave vector
$q\,c_{\mathrm{T}}/\omega_{\mathrm{BP}}$. Additionally, the prediction from
Ref. [33] (black dashed line) and from Ref. [23] (green dashed line) both at
$n=1.0$ are shown.

The scaling of the sound velocity with density, $c_{T}\propto\sqrt{n}$, also
explains the scaling of the Debye frequency, as $A_{D}\approx
1/(2n\pi^{2}c_{T}^{3})$ approximately holds for high enough densities [22,
33]. The Debye amplitude $A_{D}$ grows for increasing disorder because the
sound velocity softens. Additionally, $A_{D}$ is inversely proportional to the
number of degrees of freedom, which become fewer with lowering $n$ at fixed
volume. Both effects together cause the dependence $\omega_{D}\propto n^{5/6}$
in $d=3$ which is appreciably stronger than the dependence of the boson peak
frequency on disorder.
Fig. 6 shows the participation ratios $P^{k}$ for a single realisation of the
inherent positions at different system sizes $N$ at $n=1.0$. We can clearly
see frequencies with large participation ratios of magnitude close to the
ideal plane wave at small frequencies $\omega^{k}$. These frequencies have
distinct gaps between them and belong to the phonon bands which have discrete
frequencies due to the finite size of our system. The first participation
ratios are therefore located at $\omega=c_{\mathrm{T}}2\pi/L$.

Figure 6: Participation ratio $P^{k}$ in dependency of the eigenfrequencies
$\omega^{k}$ for system sizes $N=10^{4}$ (purple), $N=3\times 10^{4}$ (green)
and $N=3\times 10^{5}$ (cyan). Note that in the last case, only the 2000
lowest participation ratios can be calculated. The black dashed line is
located at the boson peak frequency. The red lines show the average, maximum
and minimum participation ratio of an large ensemble of $N=10^{4}$ systems.

With increasing $\omega$ the participation ratio gets smaller and at the boson
peak it distinctively drops. Above the boson peak we observe localized modes.
The behaviour is qualitatively the same for different system sizes. We have
run additional simulations for an increased ensemble size of $5\times 10^{5}$
systems each with $N=10^{4}$. We calculate the average, minimum and maximum
participation ratio in distinct frequency bins. The maximum and minimum
participation ratio per bin is shown in Fig. 6 by the red lines dashed dotted
lines and the average participation ratio by the full line. The large ensemble
confirms the overall shape of the participation ratio distribution. Contrary
to the interpretation in [33], we do not find (quasi-) localised modes. This
may be a consequence of studying scalar excitations [20, 44].

We can now turn to the results obtained by the excited wave analysis. Fig. 7
shows the veloctiy correlation function $C(q,t)$ for $q_{4}$ and $q_{7}$. One
can estimate the uncertainties of $C(q,t)$ by the standard deviation. For
better visibility we do not shown them in Fig. 7 because they are smaller than
the symbols.
Clearly, a damped oscillation can be observed where the damping and the
frequency becomes larger for large $q$. For some $q$ there are beats with
large magnitude for larger times $t$ which can be somewhat eliminated by
averaging over different inherent structures. However, we get a more reliable
result for the damping if we fit $\exp(-\Gamma(q)t)$ to the envelope
$C_{\mathrm{env}}(q,t)$ of the correlation function.

Figure 7: Velocity correlation function $C(q,t)$ of $N=4\times 10^{5}$ systems
at $n=1.0$ for $q_{4}$ (a) and $q_{7}$ (b). The black dashed lines correspond
to a fitted damped oscillation.

Fig. 8a) and Fig.8b) show the envelope $C_{\mathrm{env}}(q,t)$ for the
wavevectors $q_{4}$ and $q_{7}$ at $N=4\times 10^{5}$. Here, the uncertainties
are those of $C(q,t)$ evaluated at the respective maxima and minima. The
initial decay of the envelope is exponential (a fit is shown by the dashed
black lines) but large deviations from the exponential decay are visible. This
effect was already observed and discussed in [12]. They argued that this is a
finite size effect and that the deviations start at a system size dependent
time and are stronger for smaller wavevector. Our results confirm this
observation as can be seen in Fig. 8c) where the envelope of the correlation
function for three different systems sizes at a similar wavevector $q$ are
shown. The uncertainty of the damping $\Gamma(q)$ therefore is due to the
uncertainty of the range in which the envelope can by fitted by an
exponential.

Figure 8:  Envelope of the velocity correlation functions shown in Fig.7 for
$q_{4}$ (a) and $q_{7}$ at $n=1.0$ (b). An exponential decay is fitted and is
drawn by the black dashed line. (c) Envelope $C_{\mathrm{env}}(q,t)$ of three
different system sizes at a similar wavevector.

So far, we neglected the fact that only the $2000$ smallest eigenvectors are
calculated for systems with $N>4\times 10^{4}$. We will now argue that the
$2000$ smallest eigenvectors are sufficient to capture the small $q$ behaviour
of the velocity correlation function. For this we look at the hybridization
coefficients shown in Fig. 9.
As can be seen, for small values $q$ we get very large hybridization
coefficients only for a narrow band of frequencies. As $q$ increases, a wider
range of frequencies is involved in the response of the system to a standing
wave. At a certain value of $q$ the frequencies with large hybridization
coefficients start to lie outside of the smallest $2000$ eigenvalues. At this
point the approximation of using only the $2000$ smallest eigenvectors
necessarily fails. In Fig. 9b) we show the sum over the hybridization
coefficients for different values of $q$. As long as this sum is $1$ our
approximation for $C(q,t)$ is valid. If the sum distinctively differs from $1$
we have not included enough frequencies and the approximation fails. For the
shown system this is the case at $q>q_{7}$.

Figure 9: (a) Hybdridization coefficients for the $7$ smallest values of $q$
in a system at size $N=4\times 10^{5}$. The larger $q$ becomes, the wider the
range becomes over which non-vanishing hybridization coefficients are spread.
The corresponding values of $\Omega(q)$ are indicated by the dashed grey
lines. (b) The sum over the smallest $2000$ hybridization coefficients in this
system drastically drops as soon as $q$ exceeds $q_{7}$.

Finally, we show the obtained damping $\Gamma(q)$ in Fig. 10. We observe that
the attenuation becomes larger for increasing wave vectors $q$. The damping
saturates for large $q$. For very small $q$ we observe a weak quartic
(Rayleigh) damping which can be fitted by $\Gamma(q)\sim q^{4}$. The
theoretically predicted prefactor $B_{R}$ [33] lies within a factor of two
relative to the exact prefactor. The theory prediction $B_{r}q^{4}$ is
included in Fig. 10. In the high density regime holds
$B_{R}\approx\frac{7}{48\pi}\frac{\omega_{BP}^{4}}{nc_{T}^{3}}$ around the
sound pole.

Figure 10: Sound attenuation $\Gamma(q)$ in dependency of the wave vector $q$
at $n=1.0$. The data obtained from different systems sizes is combined. Note
that the small systems are used to obtain the large $q$ behaviour, while the
large systems are used to obtain the small $q$ behaviour. The magenta colored
line corresponds to Rayleigh damping $\Gamma(q)=B_{\mathrm{R}}q^{4}$ where the
strength $B_{\mathrm{R}}$ is taken form the theory [33]. The blue line shows
the attenuation as calculated from the dynamic structure factor $S(q,\omega)$.

We have also studied 2D systems. Compared to their $3D$ counterparts the box
size $L$ of the 2D systems at our maximum particle number of $N=10^{6}$ is
larger. Hence, smaller wavevectors $q$ are accessible in the 2D systems. In
general, the 2D systems behave quite similar to their 3D counterparts. Note
that the Debye behaviour becomes $g(\omega)=A_{\mathrm{D}}\,\omega$ in 2D with
$A_{D}=1/\omega_{D}^{2}$.
In Fig. 11 the density of states $g(\omega)$, the dependences of the
frequencies on density, $\hat{\omega}(n)$, the damping $\Gamma(q)$, and the
participation ratios $P^{k}$ of the 2D systems are shown. We observe a Debye
spectrum for $\omega\to 0$. Lowering the density $n$ again increases the Debye
level i.e. decreases the Debye frequency. The Boson peak frequency
$\omega_{\mathrm{BP}}$ again scales with $\sqrt{n}$, while the Debye frequency
scales with $n^{1}$. We observe Rayleigh damping for small wavevectors $q$.
Note that in 2D the Rayleigh damping becomes $\Gamma(q)\propto q^{3}$ [41].
In Fig. 11d) the corresponding participation ratios are shown. Again as in 3D,
a crossover from extended to localized modes can be observed, which happens
already at lower frequencies than $\omega_{BP}$ in 2D. An ensemble of $500$
systems allows to characterize the distribution of $P^{k}$, which lies lower
at the boson peak frequency in 2D than in 3D; compare with Fig. 11d).

Fig. 12 shows the contact number $\overline{z}=M/N-1$ where $M$ is the number
of nonzero entries of $\vec{M}$. Recall from Sect. IV, that a finite cut-off
$\sigma$ was chosen in order to speed-up the matrix diagonalizations. Thus the
contact number, which would be infinite for the Gaussian spring function in
Eq. (3), becomes finitie.

Figure 11:  2D Systems. (a) Density of states $g(\omega)$ for different
densitites $n$ devided by the Debye behaviour. (b) Obtained and rescaled Debye
frequency $\hat{\omega}_{\mathrm{D}}=\omega_{\mathrm{D}}\,n^{-1}$ and boson
peak frequency $\hat{\omega}_{\mathrm{BP}}=\omega_{\mathrm{BP}}\,n^{-1/2}$ in
dependency of the density $n$. Additionally, the ratio
$\hat{\omega}_{\mathrm{BP}}/\hat{\omega}_{\mathrm{D}}$. (c) Sound attenuation
for the systems sizes $N=10^{4}$ and $N=4\times 10^{5}$ at $n=1$. A $q^{3}$
fit is shown for the small $q$ regime. (d) Participation ratio $P^{k}$ in
dependency of the eigenfrequencies $\omega^{k}$ at $n=1$ for $N=3\times
10^{4}$ (green) $N=4\times 10^{5}$ (purple). Only the 2000 lowest
participation ratios can be calculated for $N=4\times 10^{5}$. The black
dashed line is located at $\omega_{BP}$. The red lines show the average,
maximum and minimum participation ratio of an ensemble of systems with
$N=3\times 10^{4}$.

The contact number depends linear on the density $n$. This dependency becomes
evident if one considers a $d$ dimensional sphere with the cut-off radius
$\sigma$ and volume $V_{d}$ around a test particle. The contact number $z$ of
the test particle is given by the number of other particles in the sphere and
thus simply reads $V_{d}\,n$. In Fig. 12 the expected contact number is shown
by the dashed black lines. At all densities $n$ the contact number is above
Maxwell’s isostatic stability criterion $\overline{z}_{c}=2\,d$ which is
indicated by the dashed dotted lines [45]. Thus we conclude that the cut-off
leads to negligible errors only.

Figure 12: Contact number $\overline{z}$ in dependency of the density $n$ for
the 3D (red circles) and 2D (purple diamonds) systems. The dashed lines show
the expected contact numbers and the dashed dotted lines are located at the
isostatic stability criterion $\overline{z}_{c}=2\,d$

##  VI Conclusions and outlook

Aim of the present contribution is to argue that the ERM system at large
contact numbers captures the pertinent phenomena of transverse vibrations in
stable glass at vanishing temperature [12, 13, 41, 3]. Two spatial dimensions,
$d=3$ and $d=2$, were studied. We used dimensionless frequencies and
wavevectors, which can be mapped onto experimental systems in the following
way: The frequency or time scale of the ERM model is given by the boson peak
frequency which was denoted $\omega_{BP}$. The corresponding length scale is
obtained using the (transverse) sound velocity, $q_{BP}=\omega_{BP}/c_{T}$.
Both scales allow to map the vDOS and the dynamical structure factor onto
measured data. The strength of the disorder then is the only free parameter.
Within the ERM model with a Gaussian spring function it is quantified by the
dimensionless density $n$. When comparing to real systems, the amplitude of
the Debye-level at the position of the boson peak, viz.
$g_{D}(\omega_{BP})=A_{D}\,\omega_{BP}^{2}$ relative to the amplitude of the
boson peak, viz. $g(\omega_{BP})$, can be used to match $n$. This gives the
dominant variation with disorder in the ERM system, when rescaled variables
$\omega/\omega_{BP}$ and $q/q_{BP}$ are employed. Importantly, all parameters
of the ERM model are thus accessible by experiment easily and by well defined
procedures.

Our numerical results should be compared to earlier numerical work on the same
ERM system using diagonalizations and approximate calculations with the method
of moments [22, 24]. We extend the range of accessible wavevectors so that
clear statements on the damping of sounds modes becomes possible. Considering
the slow crossover to $\Gamma(q)\propto q^{d+1}$ exhibited in Fig. 10 and
panel c) of Fig 11 for $q\to 0$, we conclude that previously no statements on
Rayleigh damping had been possible in the numerical ERM solutions.

We also compared our numerical results to predictions obtained from two self-
consistent theories, one where all diagrams of first order in a diagrammatic
perturbative expansion in $1/n$ were re-summed [23, 24], and a more recent
one, where all diagrams of second order in $1/n$ were re-summed [33]. Both
theories qualitatively agree in the predictions of a Wigner-semicircle law for
the boson peak and of a Debye-law at low frequencies, with the crucial
difference of the sound damping. Non-planar diagrams in the diagrammatic
perturbation expansion are required to correctly predict Rayleigh damping of
sound. The non-planar diagrams capture non-local correlations of elastic
fluctuations and arise first in second order in $1/n$. Thus they were missed
in the older approach which consequently predicts hydrodynamic damping
$\Gamma(q)\propto q^{2}$. See Ref. [33] for more details of the comparison of
both theories.

Finally, let us address the interpretation of the main peak in the vDOS of the
ERM model. When discussing Fig. 1 and Fig. 2, we called this peak boson peak
and suggested to consider it as simple model of the (transverse contribution
to the) boson peak in the vDOS of real glasses of simple particle systems.
This interpretation, which differs from the older literature on the ERM system
[23, 24, 26], rests on the following arguments:
$(i)$ The boson peak in the vDOS of real glasses survives in the zero
temperature limit and thus arguably should be contained in an harmonic
approach such as the ERM model.
$(ii)$ The reported universality of the boson peak [46, 44] is mirrored in its
origin in the disorder, which is the single conceptual extension of the
harmonic approach to disordered solids beyond the classical Born-Debye theory
of crystals: The nature of the dominant normal modes changes in the frequency
region of the boson peak. While the vibrational modes are mainly extended
below $\omega_{BP}$, they are localised above of the peak [47, 48, 49], which
can be clearly seen in Fig. 6. This transition in the nature of modes has been
related to the anomalous thermal properties of glasses [50, 51, 47].
$(iii)$ It was argued in [48, 49] that the normal modes in the region of the
Boson peak follow the statistics of a Gaussian orthogonal random matrix
ensemle (GOE). This is why the shape of the boson peak resembles Wigner’s
semicircle law which is one of the hallmarks of a GOE matrix [52]. Clearly,
our peak shown and analysed in the figures 1,2,6, and 11 is in accordance with
this characteristic of the boson peak. This also hints at why the planar ERM
model [24, 23] describes this part of the spectrum rather well, both
quantitatively and qualitatively. It is known, that only the planar diagrams
survive in thermodynamic limit for GOE matrices [52]. Again, see [33] for a
detailed comparison of the planar and non-planar self-consistent model.
$(iv)$ The density of states of a stable ERM system exhibits only a single
peak as seen in simulations of stable glasses [13, 41]. Arguably, the
vibrational phenomena we presented agree qualitatively and even semi-
quantitatively with these simulations.
$(v)$ An alternative explanation for the origin of the prominent peak in the
ERM model is, that it is a smeared out van Hove singularity [53]. But this can
not rationalise the observed transition from extended to localised eigenmodes
at $\omega_{BP}$. Almost all eigenmodes above $\omega_{BP}$ are localized (for
the considered simple ERM spring function) and while the dispersion relation
approaches $\omega_{BP}$ in the limit of $q\to\infty$, the number of extended
modes in the boson peak region is far too small in order to justify its
interpretation as smeared van Hove singularity. This is further supported by
the analytic derivation of Wigner’s semicircle law. It rests on the solution
of the self-consistency equation for $g(\omega)$ neglecting the coupling to
other modes [23, 33]. Again, this is in accordance with observed transition in
the nature of normal modes.

The present study considered the simplest ERM system of a Gaussian spring
function. Work is under way to extend it to richer spring functions.
Additionally, recent works have shown that the vectorial character of the
displacement field in glass is important as vortex like eigenmodes become
possible [11, 20, 54, 44]. These stuctures can be studied in appropriately
generalized ERM systems.

##  VII Acknowledgments

We thank Giancarlo Ruocco, Walter Schirmacher, and Grzegorz Szamel for
fruitful discussions. The work was supported by the Deutsche
Forschungsgemeinschaft (DFG) via SFB 1432 project CO7.

## References

  * Binder and Kob [2005] K. Binder and W. Kob,  _Glassy Materials And Disordered Solids: An Introduction To Their Statistical Mechanics_ (World Scientific Publishing Company, 2005).
  * Schirmacher and Ruocco [2020] W. Schirmacher and G. Ruocco, Heterogeneous elasticity: The tale of the boson peak (2020).
  * Monaco and Mossa [2009] 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).
  * Baldi  _et al._ [2010] 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 (2010).
  * Baldi  _et al._ [2013] G. Baldi, M. Zanatta, E. Gilioli, V. Milman, K. Refson, B. Wehinger, B. Winkler, A. Fontana, and G. Monaco, Emergence of crystal-like atomic dynamics in glasses at the nanometer scale, Phys. Rev. Lett. 110, 185503 (2013).
  * Baldi  _et al._ [2014] 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 (2014).
  * Wuttke  _et al._ [1995] J. Wuttke, W. Petry, G. Coddens, and F. Fujara, Fast dynamics of glass-forming glycerol, Phys. Rev. E 52, 4026 (1995).
  * Mizuno  _et al._ [2014] H. Mizuno, S. Mossa, and J.-L. Barrat, Acoustic excitations and elastic heterogeneities in disordered solids, Proceedings of the National Academy of Sciences 111, 11949–11954 (2014), https://www.pnas.org/doi/pdf/10.1073/pnas.1409490111 .
  * Gelin  _et al._ [2016] S. Gelin, H. Tanaka, and A. Lematre, Anomalous phonon scattering and elastic correlations in amorphous solids, Nature Materials 15, 1177 (2016).
  * Shimada  _et al._ [2018] M. Shimada, H. Mizuno, and A. Ikeda, Anomalous vibrational properties in the continuum limit of glasses, Physical Review E 97, 10.1103/physreve.97.022609 (2018).
  * Mizuno  _et al._ [2017] H. Mizuno, H. Shiba, and A. Ikeda, Continuum limit of the vibrational properties of amorphous solids, Proceedings of the National Academy of Sciences 114, 10.1073/pnas.1709015114 (2017).
  * Wang  _et al._ [2019a] L. Wang, L. Berthier, E. Flenner, P. Guan, and G. Szamel, Sound attenuation in stable glasses, Soft Matter 15, 7018–7025 (2019a).
  * Wang  _et al._ [2019b] L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel, and E. Flenner, Low-frequency vibrational modes of stable glasses, Nature Communications 10, 10.1038/s41467-018-07978-1 (2019b).
  * Saitoh and Mizuno [2021] K. Saitoh and H. Mizuno, Sound damping in soft particle packings: the interplay between configurational disorder and inelasticity, Soft Matter 17, 4204 (2021).
  * Hu and Tanaka [2022] Y.-C. Hu and H. Tanaka, Origin of the boson peak in amorphous solids, Nature Physics 18, 669 (2022).
  * Szamel and Flenner [2022] G. Szamel and E. Flenner, Microscopic analysis of sound attenuation in low-temperature amorphous solids reveals quantitative importance of non-affine effects, The Journal of Chemical Physics 156, 144502 (2022).
  * Mahajan and Ciamarra [2021] S. Mahajan and M. P. Ciamarra, Unifying description of the vibrational anomalies of amorphous materials, Phys. Rev. Lett. 127, 215504 (2021).
  * Horbach  _et al._ [2001] J. Horbach, W. Kob, and K. Binder, High frequency sound and the boson peak in amorphous silica, The European Physical Journal B - Condensed Matter and Complex Systems 19, 531 (2001).
  * Shimada  _et al._ [2019] M. Shimada, H. Mizuno, and A. Ikeda, Vibrational spectrum derived from local mechanical response in disordered solids., Soft matter (2019).
  * Lerner and Bouchbinder [2021] E. Lerner and E. Bouchbinder, Low-energy quasilocalized excitations in structural glasses, The Journal of Chemical Physics 155, 200901 (2021).
  * Mézard  _et al._ [1999] M. Mézard, G. Parisi, and A. Zee, Spectra of euclidean random matrices, Nuclear Physics B 559, 689 (1999).
  * Martin-Mayor  _et al._ [2001] V. Martin-Mayor, M. Mézard, G. Parisi, and P. Verrocchio, The dynamical structure factor in topologically disordered systems, The Journal of Chemical Physics 114, 8068 (2001).
  * Grigera  _et al._ [2001] T. S. Grigera, V. Martín-Mayor, G. Parisi, and P. Verrocchio, Vibrational spectrum of topologically disordered systems, Physical Review Letters 87, 10.1103/physrevlett.87.085502 (2001).
  * Ciliberti  _et al._ [2003] S. Ciliberti, T. S. Grigera, V. Martín-Mayor, G. Parisi, and P. Verrocchio, Brillouin and boson peaks in glasses from vector euclidean random matrix theory, The Journal of Chemical Physics 119, 8577 (2003).
  * Grigera  _et al._ [2011] T. S. Grigera, V. Martin-Mayor, G. Parisi, P. Urbani, and P. Verrocchio, On the high-density expansion for euclidean random matrices, Journal of Statistical Mechanics: Theory and Experiment 2011, P02015 (2011).
  * Ganter and Schirmacher [2010] C. Ganter and W. Schirmacher, Rayleigh scattering, long-time tails, and the harmonic spectrum of topologically disordered systems, Phys. Rev. B 82, 094205 (2010).
  * Ganter and Schirmacher [2011] C. Ganter and W. Schirmacher, Euclidean random matrix theory: low-frequency non-analyticities and rayleigh scattering, Philosophical Magazine 91, 1894 (2011).
  * Folli  _et al._ [2017] V. Folli, G. Ruocco, and W. Schirmacher, Moment-preserving theory of vibrational dynamics of topologically disordered systems, Frontiers in Physics 5 (2017).
  * Schirmacher  _et al._ [2019] W. Schirmacher, V. Folli, C. Ganter, and G. Ruocco, Self-consistent euclidean-random-matrix theory, Journal of Physics A: Mathematical and Theoretical 52, 464002 (2019).
  * Mehta [2014] M. L. Mehta,  _Random Matrices_ , 3rd ed., Pure and Applied Mathematics (Academic Press, 2014).
  * Goetschy and Skipetrov [2013] A. Goetschy and S. E. Skipetrov, Euclidean random matrices and their applications in physics (2013).
  * Conyuh and Beltukov [2021] D. A. Conyuh and Y. M. Beltukov, Random matrix approach to the boson peak and ioffe-regel criterion in amorphous solids, Phys. Rev. B 103, 104204 (2021).
  * Vogel and Fuchs [2023] F. Vogel and M. Fuchs, Vibrational phenomena in glasses at low temperatures captured by field theory of disordered harmonic oscillators, Physical Review Letters 130, 10.1103/physrevlett.130.236101 (2023).
  * Knyazev [2001] A. V. Knyazev, Toward the optimal preconditioned eigensolver: Locally optimal block preconditioned conjugate gradient method, SIAM Journal on Scientific Computing 23, 517–541 (2001), https://doi.org/10.1137/S1064827500366124 .
  * Ashcroft and Mermin [1976] N. W. Ashcroft and N. D. Mermin,  _Solid State Physics_ (Holt-Saunders, 1976).
  * Wang  _et al._ [2021] L. Wang, G. Szamel, and E. Flenner, Low-frequency excess vibrational modes in two-dimensional glasses, Physical Review Letters 127, 10.1103/physrevlett.127.248001 (2021).
  * Wang  _et al._ [2022] L. Wang, L. Fu, and Y. Nie, Density of states below the first sound mode in 3d glasses, The Journal of Chemical Physics 157, 074502 (2022).
  * Wang  _et al._ [2023] L. Wang, G. Szamel, and E. Flenner, Scaling of the non-phononic spectrum of two-dimensional glasses, The Journal of Chemical Physics 158, 126101 (2023).
  * Ramos [2022] M. A. Ramos,  _Low-Temperature Thermal and Vibrational Properties of Disordered Solids_ (WORLD SCIENTIFIC (EUROPE), 2022).
  * Shintani and Tanaka [2008] H. Shintani and H. Tanaka, Universal link between the boson peak and transverse phonons in glass, Nature Materials 7, 870 (2008).
  * Mizuno and Ikeda [2018] H. Mizuno and A. Ikeda, Phonon transport and vibrational excitations in amorphous solids, Physical Review E 98, 10.1103/physreve.98.062612 (2018).
  * Matsumoto and Nishimura [1998] M. Matsumoto and T. Nishimura, Mersenne twister, ACM Transactions on Modeling and Computer Simulation 8, 3 (1998).
  * [43] The slight shift of the maximum position in Fig. 2 arises from the spline interpolation through the discrete data set.
  * Tanguy [2023] A. Tanguy, Vibrations and heat transfer in glasses: the role played by disorder (2023), arXiv:2307.15038 [cond-mat.dis-nn] .
  * Maxwell [1864] J. C. Maxwell, L. on the calculation of the equilibrium and stiffness of frames, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 27, 294 (1864).
  * Marruzzo  _et al._ [2013] 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).
  * Akkermans and Maynard [1985] E. Akkermans and R. Maynard, Weak localization and anharmonicity of phonons, Phys. Rev. B 32, 7850 (1985).
  * Sarkar  _et al._ [2004] S. K. Sarkar, G. S. Matharoo, and A. Pandey, Universality in the vibrational spectra of single-component amorphous clusters, Phys. Rev. Lett. 92, 215503 (2004).
  * Schirmacher  _et al._ [1998] W. Schirmacher, G. Diezemann, and C. Ganter, Harmonic vibrational excitations in disordered solids and the “boson peak”, Physical Review Letters 81, 136 (1998).
  * Zeller and Pohl [1971] R. C. Zeller and R. O. Pohl, Thermal conductivity and specific heat of noncrystalline solids, Phys. Rev. B 4, 2029 (1971).
  * Anderson  _et al._ [1972] P. W. Anderson, B. I. Halperin, and C. M. Varma, Anomalous low-temperature thermal properties of glasses and spin glasses, Philosophical Magazine 25, 1 (1972).
  * Anderson  _et al._ [2010] G. W. Anderson, A. Guionnet, and O. Zeitouni,  _An introduction to random matrices_ , 118 (Cambridge university press, 2010).
  * Taraskin  _et al._ [2001] S. N. Taraskin, Y. L. Loh, G. Natarajan, and S. R. Elliott, Origin of the boson peak in systems with lattice disorder, Phys. Rev. Lett. 86, 1255 (2001).
  * Wu  _et al._ [2023] Z. W. Wu, Y. Chen, W.-H. Wang, W. Kob, and L. Xu, Topology of vibrational modes predicts plastic events in glasses, Nature Communications 14, 2955 (2023).