Quasi-localized vibrational modes, Boson peak and sound attenuation in model
mass-spring networks

Shivam Mahajan1 and Massimo Pica Ciamarra1,2,3$\star$

1 Division of Physics and Applied Physics, School of Physical and Mathematical
Sciences, Nanyang Technological University, Singapore

2 CNRS@CREATE LTD, 1 Create Way, #08-01 CREATE Tower, Singapore 138602

3 CNR–SPIN, Dipartimento di Scienze Fisiche, Università di Napoli Federico II,
I-80126, Napoli, Italy

⋆<EMAIL_ADDRESS>
## Abstract

We introduce an algorithm that constructs disordered mass-spring networks
whose elastic properties mimic that of glasses by tuning the fluctuations of
the local elastic properties, keeping fixed connectivity and controlling the
prestress. In two dimensions, the algorithm reproduces the dependence of
gasses’ vibrational properties, such as quasi-localised vibrational modes and
Boson peak, on the degree of stability. The sound attenuation displays
Rayleigh scattering and disorder-broadening regimes at different frequencies,
and the attenuation rate increases with increased stability. Our results
establish a strong connection between the vibrational features of disordered
solids and the fluctuations of the local elastic properties and provide a new
approach to investigating glasses’ vibrational anomalies.

###### Contents

  1. 1 Introduction
  2. 2 Numerical model and protocols
  3. 3 $f$ dependence of mechanical and geometrical properties
  4. 4 Disorder Parameter and Fluctuations of Elastic Properties
  5. 5 Vibrational Spectra
    1. 5.1 Boson Peak
    2. 5.2 Quasi-Localized Modes
  6. 6 Phonon Attenuation
  7. 7 Conclusions
  8. A Local elasticity

##  1 Introduction

The density of vibrational states (vDOS) controls solids’ specific heat and
transport properties [1, 2]. The vDOS of amorphous solids differs
qualitatively from that of crystals. In crystals, the low frequency
vibrational excitations are plane waves (phonons) distributed in frequency
according to Debye’s law, $D_{p}(\omega)\propto\omega^{d-1}$ in $d$ spatial
dimensions. On the contrary, amorphous materials have an excess of low-
frequency modes over Debye’s prediction that induces a peak in the reduced
density of states $D(\omega)/\omega^{d-1}$ at the Boson peak frequency,
$\omega_{\rm bp}$, in the terahertz regime for molecular solids. Previous
works have attributed the Boson peak to elastic disorder [3, 4, 5], localized
harmonic/anharmonic vibrations [6, 7, 8, 9], broadening of van Hove
singularities [10, 11] (but see [12]). In addition, in amorphous solids the
low-frequency excitations comprises both phononic-like modes and additional
quasi-localized vibrational modes (QLMs) that appears to be universally
distributed in frequency as $D_{\rm loc}(\omega)=A_{4}\omega^{4}$ [13, 14, 15,
16], with $A_{4}$ decreasing as the stability of the material increases [17].
The low-frequency $\omega^{4}$ scaling of the density of soft-localized modes
may originate from general considerations on the properties of localized
excitations in disordered systems [18] (see, however [19]). The universality
of the scaling law [15, 14] and its dimensionality independence suggest that
this scaling has a mean-field origin [20, 21, 22]. Finally, in amorphous
materials the extended low-frequency modes are not phonons: even in the
absence of temperature induced anharmonic effects [23], phonons of wave vector
$\kappa$ attenuate with a rate $\Gamma(\kappa)$ exhibiting a crossover from a
Rayleigh scattering [24] regime, $\Gamma\propto\kappa^{d+1}$, to a disordered-
broadening regime, $\Gamma\propto\kappa^{2}$ [25, 26, 27, 28, 29], as $\kappa$
increases.

The squared vibrational eigenfrequencies $\omega^{2}$ are the eigenvalues of
the matrix of the second derivatives of the energy with respect to the
particle positions or Hessian matrix. As such, the vibrational anomalies of
amorphous materials may possibly be rationalized within random matrices [30,
31, 32, 33, 34, 35]. Previous works primarily focused on the eigenvalues of
Wishart matrices, which are positively defined and hence may model stable
systems. A mean field [36] random-matrix approach suggests that the Boson peak
may originate from the reduction in coordination number driving the system
toward isostaticity [37, 31] and from hierarchical energy landscape. These two
scenarios are possibly relevant in colloidal hard-sphere-like glasses and
highly connected molecular systems. The random matrix approach may also be
used to investigate QLMs. In this case, the issue is determining the random
matrix ensemble reproducing the $D_{\rm loc}(\omega)$ distribution
characterizing amorphous solids or, equivalently, the correlations to be
enforced on the matrix. Research in this direction [38] succeeded in
reproducing a pseudogap, $D(\omega)\propto\omega^{\alpha}$, with an exponent
$\alpha<4$. In this research direction, the issue is integrating the two
approaches to random matrices that reproduce at the same time Boson peak and
quasi-localized modes, as well as their correlations.

Other approaches recovered the $\omega^{4}$ distribution by describing an
amorphous material as an elastic continuum punctuated by defects, possibly
anharmonic or interacting [39, 21, 22, 40]. Localized vibrations may thus
cause all vibrational anomalies of glasses, considering that they may induce
the Boson peak [6, 7, 8, 9] and control sound attenuation in Rayleigh’s theory
[24]. Recent numerical results supported this scenario in three dimensions by
demonstrating a relation between QLMs’ frequency distribution and Boson peak
frequency [17, 41], $A_{4}\propto\omega_{\rm bp}^{-5}$. In two-dimensions,
vibrations with frequencies close to the Boson peak consist of phonons
hybridized with QLMs [42].

In this manuscript, we investigate the physical origin of the vibrational
anomalies of amorphous solids by creating mass-spring networks, or
equivalently Hessian matrices, that reproduce them and their relationships.
Rather than looking for the ensemble of random matrices exhibiting the
anomalies of interest, we study how to vary an amorphous solid’s mass-spring
network to modulate them. Similar approaches have been introduced to induce a
Boson peak in unstressed networks [43] or suppress $A_{4}$ by artificially
reducing the prestress in stressed ones [44]. Here, we introduce an algorithm
to generate systems with varying degrees of elastic disorder and prestress. We
find that the boson peak frequency $A_{4}$ satisfies the relation
$A_{4}\propto\omega_{\rm bp}^{-5}$ and that $A_{4}$ decreases as the degree of
elastic disorder is suppressed. Our networks also reproduce sound’s
attenuation crossover from a Rayleigh scattering to a disordered broadening
regime observed in amorphous solids. In the Rayleigh scattering regime, the
attenuation rate relates to the material properties as predicted by
fluctuating elasticity theory [3, 4], consistently with our observation of no
varying correlations in the elastic properties. Our results support a deep
connection between the vibrational anomalies and glasses and clarify their
relationship with the fluctuations of the local elastic properties.

The paper is organised as follows. We introduce our approach to construct
mass-spring networks in Sec. 2. We illustrate how it influences geometrical
and macroscopic elastic properties in Sec. 3, and the local elastic ones in
Sec. 4. Sec. 5 demonstrates that the vibrational spectrum of our networks
reproduces the vibrational anomalies of amorphous materials and their
correlations [41]. Finally, Sec. 6 shows that sound attenuation in our
networks crossovers from Rayleigh scattering regime $\Gamma\sim\omega^{3}$ to
disorder-broadening regime $\Gamma\sim\omega^{2}$, as expected for two-
dimensional solids. We summarise our results and discuss future research
directions in the conclusions.

##  2 Numerical model and protocols

Our mass-spring network generating algorithm takes as input the disordered
mass-spring network associated with the linear response regime of an amorphous
solid, which generally has bond-depending elastic constants and rest lengths.
We transform this original network by swapping the attributes of randomly
selected bond pairs, i.e., by exchanging their spring constants and rest
lengths, as schematically illustrated in Fig. 1. Henceforth, the algorithm
does not vary the connectivity. We define the fraction of swapped bonds as
$f=2N_{\rm swap}/N_{b}$, where $N_{\rm swap}$ is the number of swap moves,
$N_{b}$ the number of bonds in the network, and the factor $2$ accounts for
the fact that each swapping event involves two bonds. Hence, for $f=0$ we
retain our original network, while for $f=1$ each bond has been swapped once
on average. After the bond swapping, we minimize the energy of the new network
bringing it into a mechanically stable configuration and study its vibrational
properties.

We remark that, on increasing $f$, the bond randomization procedure destroys
the correlations in the local elastic properties of the initial network more
effectively. However, regardless of the $f$ value, the elastic properties of
the final network might or might not exhibit correlations that build up during
the final minimization procedure. The exact relation between $f$ and
correlation in the elastic properties, if any, needs to be determined a
posteriori.

The swapping procedure may change the network’s prestress. To ascertain if the
changes in the prestress correlate with changes in the vibrational properties,
for selected $f$ values we also generate a set of networks at fixed prestress,
which we tune by varying the density after the bond-swapping procedure.

Figure 1:  Schematic of the bond-swapping algorithm we use to tune the elastic
properties of a disordered mass-spring network. In each bond swapping event,
we randomly select two elastic springs and swap the values of their elastic
constants and rest lengths. We obtain different network by varying the
fraction of swapped bonds $f=2N_{\rm swap}/N_{b}$, where $N_{\rm swap}$ is the
number of swap moves and $N_{b}$ the number of bonds in the network, and then
minimizing the energy of the resulting network.

We have implemented our bond-swapping procedure in two dimensions. To generate
our initial network, we consider systems of particles interacting via the
Weeks-Chandler Andersen (WCA) potential

$U_{ij}(r_{ij})=4\epsilon\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right]+\epsilon;~{}~{}r_{ij}\leq
2^{1/6}r_{ij},$ (1)

with $r_{ij}$ the distance between interacting particles,
$\sigma_{ij}=(\sigma_{i}+\sigma_{j})/2.0$ where $\sigma_{i}$ is a particle
diameter drawn from a uniform random distribution in the range [0.8:1.2]. We
equilibrate systems at fixed number density $\rho=1.2$ at high temperature
$T=4\epsilon$, and then instantaneously quench them into amorphous solid
configurations by minimizing the energy via the conjugate-gradient algorithm
[45].

We generated initial mass-spring networks with $N=1024$ to $360000$ particles.
Each network is fed to our algorithm to create different networks by swapping
a fraction $f$ of the bonds. For each $N$ and $f$, we average our data over
$200$ independent initial networks unless otherwise specified.

##  3 $f$ dependence of mechanical and geometrical properties

Figure 2:  Panel (a) illustrates the radial distribution function of $N=40000$
particle systems. Panel (b) shows the distribution of the magnitude of the
spring forces $F$ normalized by their average value. Data are averaged over
200 configurations. Panel (c) illustrates the dependence of shear modulus
$\mu_{0}$ on the swapping fraction $f$. Panel (d) shows the variation of pre-
stress $e$ with $f$. Data (b)-(d) are for $N=1024$. Full symbols refer to
different $f$ values, e.g., as in (b). Open symbols for $f\geq 0.6$ represent
systems compressed to impose a fixed pre-stress value.

We investigate the influence of the bond-swapping procedure on the geometrical
and mechanical properties of the elastic network in Fig. 2. Panel a shows that
bond swapping does not affect two-point correlations as the radial
distribution function is de-facto $f$-independent. In panel b, we study the
$f$ dependence of the distribution of the interparticle forces $\rm{F}$. We
obtain our initial $f=0$ network by minimising the energy of a system of
particles interacting via a repulsive potential. Consequently, for $f=0$ all
interparticle forces are positive, i.e., repulsive. The swapping protocol
changes the interparticle forces’ distribution by inducing tensile forces.
These changes influence the network’s mechanical properties by increasing the
shear modulus and reducing the pressure, as in Fig. 2c. The reduction in
pressure results from the emergence of tensile forces. To rationalise the
origin of the shear modulus’ change, we decompose $\mu_{0}=\mu_{\rm
a}-\mu_{\rm n.a.}$ in its affine and non-affine contribution. We find
$\mu_{0}$ increases with $f$ as the affine contribution grows while the non-
affine contribution decreases.

As the swapping fraction $f$ increase, the prestress [46, 47]
$e=(d-1)\langle(-U^{\prime}(r_{ij})/r_{ij}U^{\prime\prime}(r_{ij})\rangle_{ij}$
monotonically decreases, as illustrated in Fig. 2d. The prestress change is
small and comparable to that reported in previous works [44, 48], and mostly
occurs on moving from $f=0$ to $f>0$, i.e., as tensile forces appear in the
system. To assess if this small change in prestress could sensibly affect the
vibrational properties, at selected $f$ values we create additional networks
by changing the density to impose a fixed prestress value. In Fig. 2d, we use
open symbols to indicate the pre-stress value of these configurations.

##  4 Disorder Parameter and Fluctuations of Elastic Properties

Figure 3:  (a) Scaled distribution of the single-particle shear modulus. (b)
Correlation function of the single-particle shear modulus defined in Eq. 3.
Symbols represent different $f$ values as in Fig. 2d.

Our investigation of the local elastic properties starts from the measure of
particle-based elastic ones [49], which we obtain by studying how a per-
particle defined stress tensor varies under imposed external deformations
[41]. We provide details in Appendix A. The distribution of the single-
particle shear modulus $\mu_{i}=\mu_{xy,i}$ [49] has a central Gaussian peak
and long, asymmetric tails, as illustrated in Fig. 3a, where the moduli are
scaled by their averaged value $\mu_{0}$ (equal to the macroscopic shear
modulus). We previously observed analogous distributions in three dimensions
[41]. Fig. 3a clarifies that the scaled variance
$\gamma_{1}=\sigma^{2}_{\mu_{i}}/\mu_{0}^{2}$ of the single-particle shear
modulus decreases as $f$ increases.

The spatial heterogeneities of the local elastic constant are commonly
assessed via the investigation of the disorder parameter $\gamma$ introduced
in Schirmacher’s fluctuating elasticity [3, 50]. In this theory, the elastic
constants are assumed to be short range correlated. Therefore, according to
central limit theory, the fluctuations ${\sigma}_{\it{w}}$ of the shear
modulus coarse-grained over a large length scale ${\it w}$ scale
asymptotically as the inverse of the number of particles in that region,
$N_{\it w}\propto{\it w^{d}}$ in $d$ spatial dimensions. Therefore, the
normalised fluctuations of the shear modulus asymptotically scale as

$\frac{{\sigma}_{\it{w}}^{2}}{{\mu}_{0}^{2}}\propto\frac{\gamma}{N_{\it{w}}},$
(2)

where $\gamma$ is proportional to the fluctuations of the single-particle
modulus and the correlation volume.

The assumption of spatially uncorrelated elastic constants does hold for
amorphous solids. Indeed, stability induces long-range shear correlations
($\propto r^{-d}$ in $d$ spatial dimensions) with quadrupolar spatial symmetry
in the shear stress [51] and the shear modulus [49, 52]. These correlations
have been numerically observed in model systems [51, 49, 52]via the study of a
radial correlation function that occurs for the quadrupolar symmetry [49, 52],

$C(r)=\frac{\langle\mu(0)\mu(r)\rangle-\langle\mu\rangle^{2}}{\langle\mu^{2}\rangle-\langle\mu\rangle^{2}}\cos(4\theta),$
(3)

where $\theta$ is the angle between ${\bf r}$ and the $x$-axis (as
$\mu=\mu_{xy}$). Our model amorphous networks reproduce these correlations, as
illustrated in Fig. 3b. The figure demonstrates that curves corresponding to
different $f$ values collapse onto each other. This observation indicates that
the bond-swapping does not alter the elastic constant’s spatial correlations,
as conversely the asymptotic $r^{-2}$ decay would have set in after an
$f$-dependent length.

While the local shear modulus has long-range correlations, the anisotropy of
these correlations ensures that $\langle\mu(0)\mu(r)\rangle=0$ at all $r$.
Equivalently, the radial correlation function of the local shear modulus
appears $\delta$-correlated if the quadrupolar symmetry is not taken into
account. Because of this, the fluctuations of the shear modulus of $N_{w}$
particles enclosed in a compact volume are insensitive to anisotropic
correlations and scale as if there were no correlations, Eq. 2. Indeed, many
previous works verified Eq. 2, e.g. [4, 53, 23, 41, 44, 54].

The coarse-grained normalised fluctuations of the single-particle shear
modulus of our model system also satisfy Eq. 2, as we demonstrate in Fig. 4a.
The asymptotic value of the scaled fluctuations define the disorder parameter
$\gamma$, we find to decrease as the degree of swapping $f$ increases, as in
Fig. 4b. The disorder parameter depend on the fluctuations of the single-
particle shear modulus and on its spatial correlations. Since Fig. 3b proves
that spatial correlation are $f$-independent, we expect $\gamma$ to be
proportional to the scaled fluctuations of the single-particle shear modulus,
$\gamma\propto\gamma_{1}$, as we observe in the inset of Fig. 4b.

Figure 4:  Panel a illustrates the dependence of fluctuations of the local
shear modulus on the coarse-graining length scale $w$. Panel (b) shows the
asymptotic value of Schirmacher’s disorder parameter $\gamma$ as a function of
swapping fraction $f$. The disorder parameter is proportional to the
normalised fluctuations of the single-particle shear modulus. The inset
illustrates that $\gamma$ is to a good approximation proportional to the
normalised fluctuations of the single-particle shear modulus $\gamma_{1}$.

These investigations demonstrate that the local elastic properties of our
model systems have the same correlations observed in amorphous materials,
provided the elastic correlation length stays constant. The dependence of the
disorder parameter $\gamma$ on $f$ indicates that swapping leads to the
networks resembling those of glasses with increased stability.

##  5 Vibrational Spectra

We now show that the bond-swapping algorithm leads to elastic networks whose
vibrational properties exhibit a Boson peak and QLMs and investigate how these
vibrational anomalies relate to the disorder parameter $\gamma$.

Figure 5:  Panel (a) illustrates the density of states of system with
$N=160000$ as a function of $\omega/\omega_{0}$ for different swapping
fractions $f$. The Boson peak frequency shifts towards higher frequency on
increasing $f$. Panel (b) shows the reduced $D(\omega)$, as a function of
$\omega/\omega_{0}$. The Boson peak strength $(D(\omega_{bp})/omega_{bp}$
increases with the disorder parameter, as in the inset.

###  5.1 Boson Peak

We determine the vibrational density of states of large $N=160000$ systems by
Fourier transforming the velocity auto-correlation function. Fig. 5(a)
illustrates the vDOS for different $f$ values, upon scaling the frequency by
$\omega_{0}=c_{s}/a_{0}$, with $c_{s}=\sqrt{\mu_{0}/\rho}$ the shear-wave
speed and $a_{0}=\rho^{-1/2}$ the interparticle spacing. The reduced vDOS
$D(\omega)/\omega$ exhibits a boson peak at characteristic frequency which
increases with $f$. Fig. 5(b) illustrates the reduced density of states, to
highlight that the boson peak strength decreases as $f$ increases, or
equivalently, the disordered parameter decreases, as illustrated in the inset.
These results further confirm that swapping leads to networks resembling those
of stable glasses, which have a reduced Boson peak anomaly [50, 3, 7, 17, 41].

###  5.2 Quasi-Localized Modes

We investigate the vibrational spectrum’ low-frequency end via the direct
diagonalization of the Hessian matrix. We focus on small $N=1024$ systems to
shift the lowest phonon frequency ($\omega_{min}\propto c_{s}/L$) upwards,
exposing the QLMs, and perform averages over $50000$ realizations for each $f$
value. For all swapping probabilities, $f$, QLMs are distributed in frequency
as $A_{4}\omega^{4}$, as illustrated in Fig. 6a. The amplitude $A_{4}$
decreases on increasing $f$, as in the inset, again suggesting that networks
with larger $f$ mimic those of glasses with increased stability.

Our swapping procedure leads to a small change in pre-stress $e$, which could
be considered responsible [44] for the observed change in $A_{4}$. To assess
this possibility, for value of $f>0$ we have prepared additional networks by
slightly changing the density, and re-minimizing the energy, after the
swapping procedure. We considered density changes that induce a net zero-
change in the pre-stress, as shown with open symbols in Fig. 2d. Fig. 6
demonstrates that the low-density end of the vibrational density of states of
these additional networks match those of the original ones, clarifying that
small changes in the prestress do not affect the vibrational properties.

Figure 6:  (a) $D_{\rm loc}=A_{4}\omega^{4}$ scaling of low-frequency density
of states. Different closed symbols indicate different swapping fraction $f$.
The inset illustrates the $f$ dependence of the prefactor $A_{4}$. Open
symbols for $f\geq 0.6$ represent compressed systems of constant prestress in
both the main panel and the inset. (b) A plot of $D_{\rm loc}\omega_{\rm
bp}^{5}/n$ vs $\omega$, with $n$ a number density that weakly depends on $f$
as in the inset, leads to a data collapse. For each $f$, data are averaged
over $50,000$ configurations of $N=1024$ particles.

The amplitude $A_{4}$ has the units of a density of modes over a frequency to
the power $5$. If $\omega_{\rm bp}$ is the QLMs’ characteristic frequency,
then $A_{4}=\omega_{\rm bp}^{5}/n$ with $n$ the number density. In Fig. 6, we
find that if $D(\omega)\omega_{\rm bp}^{5}/n$ is plotted versus $\omega^{4}$
data for different $f$ collapse on a master curve, $A_{4}\omega_{\rm bp}^{5}$
having a week $f$ dependence, particularly for $f>0$, as in the inset. This
result establishes a close correspondence between Boson peak and QLMs. The
estimated values of $n$ are essentially constant. A similar result holds in
three-dimensional glasses [17, 41]. In that case, however, the density of
modes $n$ resulted smaller by a factor of ten. However, Ref. [16] found more
stable glasses to have a smaller $n$.

We estimate the typical size of a mode as $Np$, with $p$ the mode-
participation ratio. Fig. 7 shows a crossover from extended modes, where $Np$
is of the order of $N=1024$, to more localised ones as the frequency
decreases. Interestingly, the localised modes involve around $Np\simeq 100$
particles, regardless of the swapping fraction $f$. This result is consistent
with previous ones suggesting that the size of the localised modes relates to
the correlation length of the elastic properties [41], as we do not see this
length changing with $f$ (see Fig. 3b). Besides, the large asymptotic $Np$
value indicates that the low-frequency modes might not be truly localised.
Indeed, previous results have shown that the lowest frequency modes have an
$N$ independent $Np$ value in three and four spatial dimensions [55],
qualifying them as localised. Conversely, $Np$ grows with $N$ in two-spatial
dimensions, indicating that two-dimensional low-frequency modes are not truly
localised.

Figure 7:  The volume associated with the localized modes ratio of modes in
$N=1024$ as a function of $\omega/\omega_{0}$. Data are averaged over 50000
realizations. For clarity, we show results for selected $f$ values.

##  6 Phonon Attenuation

We now discuss how swapping influences phonon attenuation. To evaluate the
phonon attenuation rate, $\Gamma$, we excite [56, 57] a transverse acoustic
wave by giving each particle a velocity
$\bold{v}_{i}^{0}=\bold{A}_{T}\cos(\boldsymbol{\kappa}{\bf r}_{i}^{0})$, where
$\bf{A}_{T}\boldsymbol{\kappa}=0$, considering $\boldsymbol{\kappa}$ in which
one among $\kappa_{x}$ and $\kappa_{y}$ is zero, and evaluate the velocity
auto-correlation function:

$C(t)=\frac{\sum_{i=1}^{N}{\bold{v}_{i}(0).\bold{v}_{i}(t)}}{\sum_{i=1}^{N}{\bold{v}_{i}(0).\bold{v}_{i}(0)}}.$
(4)

We remind we work by definition in the linear response regime as we are
considering the response of a system of masses and springs. For each
$\kappa=|\boldsymbol{\kappa}|$, we average this correlation function over 30
phonons from independent samples for $N\leq 360000$. Finally, we extract
attenuation rate $\Gamma$ and frequency $\omega$ as a function of wave-vector
$\kappa$ by fitting the velocity autocorrelation function to a damped
oscillation, $\cos(\omega t)e^{-\Gamma t/2}$. As an example of this procedure,
we show in the inset of Fig 8 the velocity autocorrelation function for
$\kappa=\frac{2\pi}{L}(3,0,0)$ in a $N=160000$ system and its damped
exponential fit.

Figure 8:  Dependence of the phonons’ attenuation rate on the frequency
normalized by $\omega_{0}$. Colors correspond to different swapping fractions
$f$, as in Fig. 6b. Different symbols correspond to different system sizes,
$N$=40000 ($+$), 90000 ($\bigtriangledown$), 160000 ($\star$), 250000
($\square$) and 360000 ($\circ$). The inset shows the velocity auto-
correlation function $C(t)$ of transverse phonons of wavevector $k$ excited at
$t=0$, with superimposed damped cosine wave fit, for $N=160000$, $f=0.0$, and
$\kappa^{2}=9$.

Fig 8 illustrates the dependence of the attenuation parameter on
$\omega/\omega_{0}$. At all $f$ values, we observe the crossover from strong
Rayleigh scattering $\Gamma\sim\omega^{3}$ to disorder-broadening regime
$\Gamma\sim\omega^{2}$ with increasing frequency $\omega$ as found in glasses.
At fixed $\omega/\omega_{0}$, the attenuation parameter $\Gamma$ reduces with
increasing $f$. We remark that these results do not suffer from size effects,
as we explicitly show combining data for different N values.

Rayleigh’s original model [24] explains the $\omega^{4}$ scaling of the
attenuation rate by describing the elastic medium as an elastic continuum
punctuated by isolated defects. In this model, the attenuation rate depends on
the number density of defects $n$, their size $\xi$, and the deviation of
their shear elastic properties from the background. In terms of the disorder
parameter $\gamma$, Rayleigh’s prediction results [41]

$\Gamma\frac{\xi}{c_{s}}\propto\gamma\left(\frac{\omega\xi}{c_{s}}\right)^{4}.$
(5)

The same prediction is recovered by a version of fluctuating elasticity theory
that considers the elastic properties to be correlated over a length scale
$\xi$ [50, 41]. In the present system two-dimensional system, at variance with
the three-dimensional case, the typical size of the localised modes is
constant. This suggests that $\xi$ is constant so that Rayleigh’s prediction
simplifies to $\Gamma\frac{a_{0}}{c_{s}}\propto\gamma\left(\frac{\omega
a_{0}^{4}}{c_{s}}\right)^{4}$. This prediction is obtained by fluctuating
elasticity theory at constant elastic correlation length [3]. Eq. 5 with $\xi$
constant has also been recovered by studying the disordered induced broadening
of the width of phonon bands [58].

Fig. 9 illustrates the dependence of $\Gamma\omega_{0}^{2}/\omega^{d+1}$ on
$\omega/\omega_{0}$, and demonstrates that the low-frequency attenuation rate
is well described by the theoretical prediction of Eq. 5 with $\xi$ constant.
This result confirms previous investigation of sound attenuation in two-
dimensional systems [44]. In three dimensions, Eq. 5 also holds, but the
correlation length $\xi$, identified with the size of the soft modes,
decreases as the glass becomes more stable [41].

Eq. 5 predicts the dependence of the attenuation rate on the disorder
parameter and the elastic correlation length up to a constant factor.
Estimating this factor suggests that the fluctuating elasticity framework
quantitatively underestimates the damping coefficient [59]. Possibly, this
occurs as the theory uses linear elasticity to account for the scattering by
sources of length scale $\xi$, without considering that elasticity theory only
holds on larger length scales [16]. In addition, fluctuating elasticity theory
assumes the local elastic constant to be short-range correlated, while
conversely, they are long-range correlated, as we illustrated in Fig. 3b.

Figure 9:  Scaled attenuation rate as a function of the frequency. The
asymptotic value in the $\omega\to 0$ limit is compared to the predictions of
classical FET. This prediction corresponds to Eq. 5 with $\xi$ constant, as
FET assumes the local elastic properties to be $\delta$-correlated. Symbols
identify the system size as in Fig. 8.

##  7 Conclusions

We have introduced an algorithm that modulates the elastic disorder of a mass-
spring network without affecting its connectivity by swapping the properties
of randomly selected bonds. We have generated a series of disordered networks
and investigated their elastic properties using this algorithm and an initial
network derived from a two-dimensional model system. Increasing the fraction
of swapped bonds suppresses the fluctuation of the shear modulus, reduces the
Boson peak’s amplitude and the Boson peaks’ frequency over the natural
frequency $\omega_{0}$, and increases the typical frequency of the low-
frequency modes. As such, bond-swapping leads to networks whose vibrational
properties increasingly resemble those of stable glasses. The observed changes
in vibrational properties do not relate to variation in the connectivity
occurring close to the jamming point [60], or analogous, to the emergence of
many weak contacts leading to a small effective connectivity [61]. Similarly,
they do not originate from changes in the prestress.

The generated networks reproduce the sound wave attenuation rate’s crossover
from a low-frequency Rayleigh scattering to a disordered broadening regime
characterizing amorphous solids in the harmonic approximation. The attenuation
rate in the Rayleigh scattering regime scales with the material properties as
predicted by fluctuating elasticity theory, which works under the constant
shear modulus correlation length assumption.

Our investigation of the bond-swapping algorithm to a two-dimensional model
network reproduces vibrational anomalies and sound attenuation observed in
two-dimensional amorphous solids. Soft modes are dimensionality dependent as
in two-dimensions they are not truly localised [62, 14], a feature possibly
connected to the dimensionality dependence of the glasses’ relaxation dynamics
[63, 64, 65, 66]. It would be interesting to assess if the bond-swapping
algorithm reproduces three-dimensional systems’ vibrational properties. More
generally, the influence of dimensionality on sound attenuation requires
further investigations.

#### Acknowledgments

We thank E. Lerner and E. Bouchbinder for comments on a earlier version of
this manuscript.

#### Funding information

We acknowledge support from the Singapore Ministry of Education through the
Academic Research Fund Tier 1 (RG56/21) and Tier 2 (MOE-T2EP50221-0016) and
are grateful to the National Supercomputing Centre (NSCC) of Singapore for
providing the computational resources.

##  Appendix A Local elasticity

Within two-dimensional linear elasticity the macroscopic stress and the
macroscopic strain are related by
$\sigma_{\alpha}=c_{\alpha\beta}\epsilon_{\beta}$, with
$\alpha,\beta\in\\{xx,yy,xy\\}$ and $c_{\alpha\beta}$ the stiffness tensor.
Here, we define the local stiffness matrix $c^{\rm
cg}_{\alpha\beta}=\frac{d\sigma^{\rm cg}_{\alpha}}{d\epsilon_{\beta}}$ as the
ratio between a locally defined stress and the macroscopic strain [41, 52]. We
define the coarse grained stress as $\sigma^{\rm
cg}_{\alpha}(w)=\langle\sigma_{\alpha}^{(i)}\rangle$, where the average is
over all particles $i$ in the coarse-graining volume, and
$\sigma_{\alpha}^{(i)}$ is a per-particle stress [49].

We note that other approaches could be used to define local elastic properties
[67, 68], e.g., by introducing a locally defined stresses and strains. These
diverse definitions converge for large coarse-gaining lengths but differ at
finite $w$. The definition we have adopted here recover self-averaging. In
addition, with this definition the statistics of the elastic properties coarse
grained over a sub-region containing $N_{0}$ particles of a $N\gg N_{0}$
system match those of a $N_{0}$ particle system [41].

Practically, we evaluate $c_{\alpha\beta}$ by monitoring the change in the
stresses of the particles in response to small deformation followed by energy
minimization, to capture the non-affine contribution to the elasticity, making
sure we work in the linear response regime. We coarse grain the single-
particle elastic properties over square regions of side length $w$.

## References

  * [1] C. Kittel, P. McEuen and P. McEuen,  _Introduction to solid state physics_ , vol. 8,  Wiley New York (1996).
  * [2] R. C. Zeller and R. O. Pohl,  _Thermal conductivity and specific heat of noncrystalline solids_ ,  Phys. Rev. B 4, 2029 (1971),  10.1103/PhysRevB.4.2029.
  * [3] W. Schirmacher,  _Thermal conductivity of glassy materials and the “boson peak"_ ,  Europhysics Letters (EPL) 73(6), 892 (2006),  10.1209/epl/i2005-10471-9.
  * [4] 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),  10.1038/srep01407.
  * [5] S. R. Elliott,  _A Unified Model for the Low-Energy Vibrational Behaviour of Amorphous Solids_ ,  Europhys. Lett 19(3), 201 (1992),  10.1209/0295-5075/19/3/009.
  * [6] W. A. Phillips and A. Anderson, eds.,  _Amorphous solids: low-temperature properties_ ,  Springer,  ISBN 0387103309 (1981).
  * [7] E. Maurer and W. Schirmacher,  _Local oscillators vs elastic disorder: A comparison of two models for the boson peak_ ,  Journal of Low Temperature Physics 137(3-4), 453 (2004),  10.1023/B:JOLT.0000049065.04709.3e.
  * [8] U. Buchenau, Y. M. Galperin, V. L. Gurevich, D. A. Parshint, M. A. Ramost and H. R. Schober,  _Interaction of soft modes and sound waves in glasses_ ,  Phys. Rev. B 46(5), 2798 (1992),  10.1103/PhysRevB.46.2798.
  * [9] B. Rufflé, D. A. Parshin, E. Courtens and R. Vacher,  _Boson Peak and its Relation to Acoustic Attenuation in Glasses_ ,  Phys. Rev. Lett 100, 015501 (2008),  10.1103/PhysRevLett.100.015501.
  * [10] 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),  10.1103/PhysRevLett.86.1255.
  * [11] A. I. Chumakov, G. Monaco, A. Monaco, W. A. Crichton, A. Bosak, R. Rüffer, A. Meyer, F. Kargl, L. Comez, D. Fioretto, H. Giefers, S. Roitsch _et al._ ,  _Equivalence of the Boson Peak in Glasses to the Transverse Acoustic van Hove Singularity in Crystals_ ,  Phys. Rev. Lett. 106, 225501 (2011),  10.1103/PhysRevLett.106.225501.
  * [12] L. Zhang, Y. Wang, Y. Chen, J. Shang, A. Sun, X. Sun, S. Yu, J. Zheng, Y. Wang, W. Schirmacher and J. Zhang,  _Disorder-induced vibrational anomalies from crystalline to amorphous solids_ ,  Phys. Rev. Research 3, L032067 (2021),  10.1103/PhysRevResearch.3.L032067.
  * [13] H. Mizuno, S. Mossa and J.-L. Barrat,  _Acoustic excitations and elastic heterogeneities in disordered solids_ ,  Proceedings of the National Academy of Sciences 111(33), 11949 (2014),  10.1073/pnas.1409490111.
  * [14] G. Kapteijns, E. Bouchbinder and E. Lerner,  _Universal nonphononic density of states in 2d, 3d, and 4d glasses_ ,  Physical Review Letters 121(5), 055501 (2018),  10.1103/PhysRevLett.121.055501.
  * [15] D. Richard, K. González-López, G. Kapteijns, R. Pater, T. Vaknin, E. Bouchbinder and E. Lerner,  _Universality of the nonphononic vibrational spectrum across different classes of computer glasses_ ,  Physical Review Letters 125(8), 085502 (2020),  10.1103/PhysRevLett.125.085502.
  * [16] C. Rainone, E. Bouchbinder and E. Lerner,  _Pinching a glass reveals key properties of its soft spots_ ,  Proceedings of the National Academy of Sciences 117(10), 5228 (2020),  10.1073/pnas.1919958117.
  * [17] L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel and E. Flenner,  _Low-frequency vibrational modes of stable glasses_ ,  Nature Communications 10(1) (2019),  10.1038/s41467-018-07978-1.
  * [18] V. Gurarie and J. T. Chalker,  _Bosonic excitations in random media_ ,  Physical Review B 68, 134207 (2003),  10.1103/PhysRevB.68.134207.
  * [19] A. Kumar, I. Procaccia and M. Singh,  _Density of quasi-localized modes in athermal glasses_ ,  Europhysics Letters 135, 66001 (2021),  10.1209/0295-5075/ac27e6.
  * [20] H. Ikeda,  _Universal non-mean-field scaling in the density of states of amorphous solids_ ,  Phys. Rev. E 99, 050901 (2019),  10.1103/PhysRevE.99.050901.
  * [21] E. Bouchbinder, E. Lerner, C. Rainone, P. Urbani and F. Zamponi,  _Low-frequency vibrational spectrum of mean-field disordered systems_ ,  Phys. Rev. B 103, 174202 (2021),  10.1103/PhysRevB.103.174202.
  * [22] C. Rainone, P. Urbani, F. Zamponi, E. Lerner, and E. Bouchbinder,  _Mean-field model of interacting quasilocalized excitations in glasses_ ,  SciPost Phys. Core 4, 008 (2021),  10.21468/SciPostPhysCore.4.2.008.
  * [23] H. Mizuno, G. Ruocco and S. Mossa,  _Sound damping in glasses: Interplay between anharmonicities and elastic heterogeneities_ ,  Physical Review B 101, 174206 (2020),  10.1103/PhysRevB.101.174206.
  * [24] J. W. Strutt (Lord Rayleigh),  _On the Transmission of Light through an Atmosphere containing Small Particles in Suspension, and on the Origin of the Blue of the Sky_ ,  Philos. Magazine 47, 375 (1903),  10.1017/cbo9780511703997.052.
  * [25] C. Masciovecchio, G. Baldi, S. Caponi, L. Comez, S. Di Fonzo, D. Fioretto, A. Fontana, A. Gessini, S. C. Santucci, F. Sette, G. Viliani, P. Vilmercati _et al._ ,  _Evidence for a Crossover in the Frequency Dependence of the Acoustic Attenuation in Vitreous Silica_ ,  Phys. Rev. Lett. 97, 035501 (2006),  10.1103/PhysRevLett.97.035501.
  * [26] 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),  10.1103/PhysRevLett.104.195501.
  * [27] G. Monaco and V. M. Giordano,  _Breakdown of the debye approximation for the acoustic modes with nanometric wavelengths in glasses_ ,  Proceedings of the national Academy of Sciences 106(10), 3659 (2009),  10.1073/pnas.0808965106.
  * [28] B. Ruta, G. Baldi, F. Scarponi, D. Fioretto, V. M. Giordano and G. Monaco,  _Acoustic excitations in glassy sorbitol and their relation with the fragility and the boson peak_ ,  Journal of Chemical Physics 137(21), 214502 (2012),  10.1063/1.4768955.
  * [29] 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),  10.1103/PhysRevLett.110.185503.
  * [30] T. S. Grigera, V. Martín-Mayor, G. Parisi and P. Verrocchio,  _Vibrations in glasses and euclidean random matrix theory_ ,  Journal of Physics: Condensed Matter 14, 2167 (2002),  10.1088/0953-8984/14/9/306.
  * [31] M. L. Manning and A. J. Liu,  _A random matrix definition of the boson peak_ ,  EPL (Europhysics Letters) 109, 36002 (2015),  10.1209/0295-5075/109/36002.
  * [32] M. Baggioli and A. Zaccone,  _Universal Origin of Boson Peak Vibrational Anomalies in Ordered Crystals and in Amorphous Materials_ ,  Physical Review Letters 122 (2019),  10.1103/PhysRevLett.122.145501.
  * [33] M. Pernici,  _Mean-field density of states of a small-world model and a jammed soft spheres model_ ,  10.48550/ARXIV.2001.02622 (2020).
  * [34] 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),  10.1103/PhysRevB.103.104204.
  * [35] H. Ikeda and M. Shimada,  _Vibrational density of states of jammed packing at high dimensions: Mean-field theory_ ,  Phys. Rev. E 106, 024904 (2022),  10.1103/PhysRevE.106.024904.
  * [36] S. Franz, G. Parisi, P. Urbani and F. Zamponi,  _Universal spectrum of normal modes in low-temperature glasses_ ,  Proceedings of the National Academy of Sciences 112(47), 14539 (2015),  https://doi.org/10.1073/pnas.1511134112.
  * [37] M. Wyart, S. R. Nagel and T. A. Witten,  _Geometric origin of excess low-frequency vibrational modes in weakly connected amorphous solids_ ,  Europhysics Letters (EPL) 72, 486 (2005),  10.1209/epl/i2005-10245-5.
  * [38] E. Stanifer, P. K. Morse, A. A. Middleton and M. L. Manning,  _Simple random matrix model for the vibrational spectrum of structural glasses_ ,  Phys. Rev. E 98, 042908 (2018),  10.1103/PhysRevE.98.042908.
  * [39] A. Moriel,  _Internally stressed and positionally disordered minimal complexes yield glasslike nonphononic excitations_ ,  Phys. Rev. Lett. 126, 088004 (2021),  10.1103/PhysRevLett.126.088004.
  * [40] W. Ji, M. Popović, T. W. J. de Geus, E. Lerner and M. Wyart,  _Theory for the density of interacting quasilocalized modes in amorphous solids_ ,  Phys. Rev. E 99, 023003 (2019),  10.1103/PhysRevE.99.023003.
  * [41] S. Mahajan and M. P. Ciamarra,  _Unifying description of the vibrational anomalies of amorphous materials_ ,  Phys. Rev. Lett. 127, 215504 (2021),  10.1103/PhysRevLett.127.215504.
  * [42] E. Lerner and E. Bouchbinder,  _Boson-peak vibrational modes in glasses feature hybridized phononic and quasilocalized excitations_ ,  10.48550/ARXIV.2210.10326 (2022).
  * [43] Y. Nie, H. Tong, J. Liu, M. Zu and N. Xu,  _Role of disorder in determining the vibrational properties of mass-spring networks_ ,  Frontiers of Physics 12(3), 1 (2017),  https://doi.org/10.1007/s11467-017-0668-8.
  * [44] G. Kapteijns, D. Richard, E. Bouchbinder and E. Lerner,  _Elastic moduli fluctuations predict wave attenuation rates in glasses_ ,  The Journal of Chemical Physics 154(8), 081101 (2021),  https://doi.org/10.1063/5.0038710.
  * [45] R. Fletcher,  _Conjugate gradient methods for indefinite systems_ ,  Springer Berlin Heidelberg,  ISBN 978-3-540-38129-7,  https://doi.org/10.1007/BFb0080116 (1976).
  * [46] M. Shimada, H. Mizuno, L. Berthier and A. Ikeda,  _Low-frequency vibrations of jammed packings in large spatial dimensions_ ,  Physical Review E 101(5), 052906 (2020),  https://doi.org/10.1103/PhysRevE.101.052906.
  * [47] E. DeGiuli, A. Laversanne-Finot, G. Düring, E. Lerner and M. Wyart,  _Effects of coordination and pressure on sound attenuation, boson peak and elasticity in amorphous solids_ ,  Soft Matter 10(30), 5628 (2014),  https://doi.org/10.1039/C4SM00561A.
  * [48] K. González-López, E. Bouchbinder and E. Lerner,  _Variability of mesoscopic mechanical disorder in disordered solids_ ,  Journal of Non-Crystalline Solids 604, 122137 (2023),  https://doi.org/10.1016/j.jnoncrysol.2023.122137.
  * [49] H. Tong, S. Sengupta and H. Tanaka,  _Emergent solidity of amorphous materials as a consequence of mechanical self-organisation_ ,  Nature Communications 11(1), 1 (2020),  10.1038/s41467-020-18663-7.
  * [50] W. Schirmacher,  _Some comments on fluctuating-elasticity and local oscillator models for anomalous vibrational excitations in glasses_ ,  Journal of Non-Crystalline Solids 357, 518 (2011),  10.1016/j.jnoncrysol.2010.07.052.
  * [51] A. Lemaître,  _Stress correlations in glasses_ ,  Journal of Chemical Physics 149(10), 104107 (2018),  10.1063/1.5041461.
  * [52] S. Mahajan, J. Chattoraj and M. P. Ciamarra,  _Emergence of linear isotropic elasticity in amorphous and polycrystalline materials_ ,  Physical Review E 103, 52606 (2021),  10.1103/PhysRevE.103.052606.
  * [53] M. Tsamados, A. Tanguy, C. Goldenberg and J. L. Barrat,  _Local elasticity map and plasticity in a model Lennard-Jones glass_ ,  Physical Review E 80(2), 026112 (2009),  10.1103/PhysRevE.80.026112.
  * [54] A. Shakerpoor, E. Flenner and G. Szamel,  _Stability dependence of local structural heterogeneities of stable amorphous solids_ ,  Soft Matter 16, 914 (2020),  10.1039/c9sm02022e.
  * [55] G. Kapteijns, E. Bouchbinder and E. Lerner,  _Universal nonphononic density of states in 2d, 3d, and 4d glasses_ ,  Phys. Rev. Lett. 121, 055501 (2018),  10.1103/PhysRevLett.121.055501.
  * [56] S. Gelin, H. Tanaka and A. Lematre,  _Anomalous phonon scattering and elastic correlations in amorphous solids_ ,  Nature Materials 15(11), 1177 (2016),  10.1038/nmat4736.
  * [57] H. Mizuno and A. Ikeda,  _Phonon transport and vibrational excitations in amorphous solids_ ,  Phys. Rev. E 98 (2018),  10.1103/PhysRevE.98.062612.
  * [58] E. Bouchbinder and E. Lerner,  _Universal disorder-induced broadening of phonon bands: From disordered lattices to glasses_ ,  New Journal of Physics 20(7), 073022 (2018),  10.1088/1367-2630/aacef4.
  * [59] 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(14), 144502 (2022),  https://doi.org/10.1063/5.0085199.
  * [60] M. Wyart, L. E. Silbert, S. R. Nagel and T. A. Witten,  _Effects of compression on the vibrational modes of marginally jammed solids_ ,  Phys. Rev. E 72, 051306 (2005),  10.1103/PhysRevE.72.051306.
  * [61] N. Xu, M. Wyart, A. J. Liu and S. R. Nagel,  _Excess vibrational modes and the boson peak in model glasses_ ,  Phys. Rev. Lett. 98, 175502 (2007),  10.1103/PhysRevLett.98.175502.
  * [62] H. Mizuno, H. Shiba and A. Ikeda,  _Continuum limit of the vibrational properties of amorphous solids._ ,  Proceedings of the National Academy of Sciences of the United States of America 114(46), E9767 (2017),  10.1073/pnas.1709015114.
  * [63] E. Flenner and G. Szamel,  _Fundamental differences between glassy dynamics in two and three dimensions_ ,  Nature Communications 6(1), 7392 (2015),  10.1038/ncomms8392.
  * [64] H. Shiba, Y. Yamada, T. Kawasaki and K. Kim,  _Unveiling dimensionality dependence of glassy dynamics: 2d infinite fluctuation eclipses inherent structural relaxation_ ,  Phys. Rev. Lett. 117, 245701 (2016),  10.1103/PhysRevLett.117.245701.
  * [65] S. Vivek, C. P. Kelleher, P. M. Chaikin and E. R. Weeks,  _Long-wavelength fluctuations and the glass transition in two dimensions and three dimensions_ ,  Proc. Natl. Acad. Sci. U.S.A. 114(8), 1850 (2017),  10.1073/pnas.1607226113.
  * [66] Y.-W. Li, C. K. Mishra, Z.-Y. Sun, K. Zhao, T. G. Mason, R. Ganapathy and M. P. Ciamarra,  _Long-wavelength fluctuations and anomalous dynamics in 2-dimensional liquids_ ,  Proceedings of the National Academy of Sciences 116(46), 22977 (2019),  10.1073/PNAS.1909319116.
  * [67] H. Mizuno, S. Mossa and J.-L. Barrat,  _Measuring spatial distribution of the local elastic modulus in glasses_ ,  Physical Review E 87(4), 042306 (2013),  https://doi.org/10.1103/PhysRevE.87.042306.
  * [68] A. Cakir and M. Pica Ciamarra,  _Emergence of linear elasticity from the atomistic description of matter_ ,  Journal of Chemical Physics 145(5), 054507 (2016),  10.1063/1.4960184.


In Fig. <ref> we present the HSE mass profiles resulting from each of the fits.

Posterior distributions of the parameters obtained in the fit of the NIKA2 and R18 pressure bins for the gNFW pressure model.

Reduced $\chi^2$ of the gNFW (solid) and NFW (dashed) model fits to the NIKA2 and R18 pressure bins. The different colors show the NIKA2 bins that have been used.

Posterior distributions of the parameters obtained in the fit of the NIKA2 and R18 pressure bins, combined with the XMM-Newton electron density, for the NFW pressure model. The values of $P_{zero}$ have been multiplied by $10^3$.

HSE mass profiles obtainted with the gNFW (left) and NFW (right) fit methods. Each color represents the profile obtained using the corresponding NIKA2 pressure bins.


#  Physical Wireless Resource Virtualization for Software-Defined Whole-Stack
Slicing ††thanks: This work is supported in part by the NSF awards 1827211 and
1821962.

Matthias Sander-Frigau, Tianyi Zhang, Hongwei Zhang, Ahmed E. Kamal, Arun K.
Somani  Department of Electrical and Computer Engineering, Iowa State
University
<EMAIL_ADDRESS>
###### Abstract

Radio access network (RAN) virtualization is gaining more and more ground and
expected to re-architect the next-generation cellular networks. Existing RAN
virtualization studies and solutions have mostly focused on sharing
communication capacity and tend to require the use of the same PHY and MAC
layers across network slices. This approach has not considered the scenarios
where different slices require different PHY and MAC layers, for instance, for
radically different services and for whole-stack research in wireless living
labs where novel PHY and MAC layers need to be deployed concurrently with
existing ones on the same physical infrastructure. To enable whole-stack
slicing where different PHY and MAC layers may be deployed in different
slices, we develop _PV-RAN_ , the first open-source virtual RAN platform that
enables the sharing of the same SDR physical resources across multiple slices.
Through API Remoting, PV-RAN enables running paravirtualized instances of
OpenAirInterface (OAI) at different slices without requiring modifying OAI
source code. PV-RAN effectively leverages the inter-domain communication
mechanisms of Xen to transport time-sensitive I/Q samples via shared memory,
making the virtualization overhead in communication almost negligible. We
conduct detailed performance benchmarking of PV-RAN and demonstrate its low
overhead and high efficiency. We also integrate PV-RAN with the CyNet wireless
living lab for smart agriculture and transportation.

##  I Introduction

Recent advances in 5G virtualization [1], through the use of Network Function
Virtualization (NFV) and Software Defined Networking (SDN), have led to the
development of the concept of Network Slicing [2]. Although there is no
consensus on the strict definition of a Network Slice, it generally refers to
an end-to-end logical portion of the network resources. IETF defines Network
Slicing as the collection of a set of technologies to create specialized and
dedicated logical networks as a service (NaaS) in support of network service
differentiation and meeting the diversified requirements from vertical
industries [3]. Thus, network slicing has the capability for the network to
provide a multitude of services that demand very specific requirements in
terms of latency, reliability, and bandwidth. According to telecommunications
Standard Developing Organizations (SDOs), network slicing is a key enabler for
accommodating diverse 5G services such as ultra-reliable and low-latency
communications (URLLC), extreme mobile broadband (xMBB), and massive machine-
type communication (mMTC).

Critical to network virtualization and slicing is the virtualization and
slicing of radio access networks (RANs) in addition to core networks and
edge/cloud resources. Despite much progress in RAN virtualization and slicing,
existing work has mostly focused on sharing communication capacity, and the
different slices tend to use the same PHY and MAC layers [4, 5, 6].
Nonetheless, the heterogeneous requirements of 5G and beyond wireless services
may well demand different PHY and MAC layers. For instance, optimal URLLC and
mMTC solutions may well differ in their PHY and MAC layers due to the
radically different requirements that URLLC and mMTC pose to communication
timeliness, reliability, throughput, and energy efficiency. Therefore, there
exists the unmet need for RAN virtualization solutions that enable different
slices to run different PHY and MAC layers. For convenience, we call this type
of RAN virtualization _whole-stack slicing_. Whole-stack slicing is not only
important for production 5G-and-beyond systems, it is also important for
enabling whole-stack research and innovation in wireless living labs where
novel solutions in the PHY and MAC layers can be experimented together with
existing solutions [7, 8].

To enable the use of different PHY and MAC layers at different RAN slices, we
need to enable the sharing of the physical wireless radios, a.k.a. remote
radio heads (RRHs), across slices. In this study, we use USRP software-
defined-radios (SDRs) [9] as the physical radios, and we propose the _PV-RAN_
platform that virtualizes the RAN SDRs. PV-RAN uses the Xen hypervisor [10] as
the base virtualization platform for the virtual RANs, and its current
implementation assumes that the open-source wireless software platform
OpenAirInterface (OAI) [11] is used as the protocol implementation platform.
As shown in Figure 1,

Figure 1: (a) Traditional C-RAN where the fronthaul connects slices (with
virtualized network stacks) to their dedicated remote-radio-heads (RRHs). (b)
PV-RAN where the fronthaul connects a shared RRH to multiple slices.

PV-RAN spans two main components of each RAN:

  * •

_PV-Back_ in the Dom0 of the base station computer: it interacts with the SDR
and implements the logic necessary to virtualize and multiplex SDR radio
resources.

  * •

_PV-Front_ in each DomU of the base station computer: it runs the full network
protocol stack (e.g., from PHY to PDCP in LTE) with a modified version of the
USRP low-level interface of OAI that accounts for SDR radio resource
virtualization.

Compared to the C-RAN architecture that makes use of at least one radio per
slice, the PV-RAN platform can support multiple slices using the same radio.

The contributions of this work are as follows:

  * •

We present the first open-source design and implementation of the
virtualization of physical SDR resources to enable whole-stack slicing.
Besides supporting whole-stack slicing in production systems, PV-RAN allows
researchers to prototype new cellular network PHY and MAC layers using
exclusively open-source software: PV-Back and PV-Front use Xen
paravirtualization and its inter-domain communication mechanism to transport
I/Q samples, and the virtualized slices use OAI as the wireless network
protocol stack.

  * •

PV-RAN enables running paravirtualized instances of OAI without requiring
modifying OAI source code. This is accomplished through a novel API Remoting
method that intercepts UHD library function calls from paravirtualized OAIs in
DomUs and redirect them to the privileged domain Dom0 where the library call
is executed and forwarded to the SDRs.

  * •

PV-RAN ensures slice isolation in shared SDR access, and it features novel and
effective techniques such as lightweight, in-band timestamping synchronization
between PV-Back and PV-Front.

  * •

We perform detailed performance benchmarking of PV-RAN, and demonstrate that
the inter-domain communication channel of Xen is a fast and efficient
interface for the transport of time-sensitive I/Q samples. In comparison to
OAI running in bare-metal computers without virtualization, the overhead due
to virtualization is insignificant.

  * •

We have integrated PV-RAN with the software-defined cyberinfrastructure CyNet
[7], the first field-deployment of open-source, virtualized cellular platforms
for smart agriculture and transportation.

The rest of this paper is organized as follows. We discuss related work in
Section II. in Section III, we present the system model and introduce the Xen
hypervisor and CyNet wireless living lab. We present the design and
implementation of PV-RAN in Section IV. We evaluate the performance of PV-RAN
in Section V. We share concluding remarks in Section VI.

##  II Related Work

Virtualized RAN (vRAN) centralizes and virtualizes a baseband-unit pool (vBBU)
at a strategic location (e.g., central office), where the vBBU can be deployed
on commercial-off-the-shelf (COTS) servers rather than proprietary hardware.
vRAN allows operators to rapidly expand the capacity and coverage of the
network, thus minimizing capital expenditure (CAPEX) and operating expenditure
(OPEX) [12]. Garcia et al. [13] proposed an analytical framework for vRAN
called FluidRAN that optimizes the placement of vRAN functions jointly with
the routing policy according to network and computing resources. Nikaein et
al. [14] have considered potential vRAN architectures that could meet real-
time deadlines and intensive computational and I/O requirements.

With the availability of SDR platforms and OAI and the growing interest toward
virtualization, researchers can prototype low-cost cellular networks and
evaluate the performance of different virtualized RAN architectures [15]. In
[16], Tran et al. evaluated a C-RAN testbed based on the OAI RRH/BBU split and
virtualized the OAI BBU-pool in VMWare virtual machines. Their study provides
interesting insight into the computational requirements of vBBU under
different radio resource configurations. The H2020 5G-PPP SliceNet project
[17] introduced a slice-friendly virtualized 5G RAN where the Central Unit
(CU) and Distributed Unit (DU) run in LXC containers. In [18], Trindade et al.
containerized OAI in Docker and used Kubernetes to orchestrate the RRH and BBU
containers. Nikaein et al. [19] investigated the performance of virtualized
BBUs under containers (Docker, LXC) and KVM virtual machines.

Esmaeily et al. [5] recently introduced a testbed 5GIIK capable of performing
E2E network slicing. They rely on the ETSI NFV MANO framework to manage and
orchestrate the resources required for creating, managing, and delivering
services through different slices. The FlexRAN [4] platform allows the
programmable control of the underlying RAN infrastructure via virtual control
functions, and it adds a virtualization layer over the RAN infrastructure to
enable communication capacity allocation across slices.

The aforementioned studies have touched on various aspects of RAN
virtualization, but they have not considered virtualizing the physical
wireless resources for whole-stack slicing. The work closest to PV-RAN is that
of the Wireless Spectrum Hypervisor [20]. It considered multiplexing several
OFDM signal streams through a shared SDR. It cannot support slices running
non-OFDM PHY layers; the spectrum hypervisor requires moving the OFDM
modulation component from the user slices to the hypervisor layer, thus not
being transparent to the slice users; it also uses ZeroMQ (a.k.a. ZMQ) bus for
communication between user slices and the hypervisor, which introduces much
higher overhead than the lightweight inter-domain communication mechanism
adopted by PV-RAN, as we show in Section V.

##  III Preliminaries

System model. As shown in Figure 2,

Figure 2: System model

the system consists of an orchestrator and a set of RANs. Each RAN has a base
station (BS) and a set of User Equipment (UEs), with each BS consisting of a
SDR and a computer. The SDR serves as the remote radio head (RRH), and the
computer runs the cellular network stacks realized as software. A network
slice consists of a RAN slice at one or more RANs, and each RAN slice runs a
cellular network stack implemented using OpenAirInterface (OAI) [11] with
potentially different PHY and MAC layers. We consider Frequency-Division-
Multiplexing (FDM) in RAN slicing such that different slices of a RAN operate
in non-overlapping frequency bands. The specific frequency division strategies
at different RANs could be different or the same, and this is coordinated
through the central orchestrator and defined by software. Within this context,
we study the systems issues of how to virtualize the BSes so that a BS SDR can
be shared among multiple slices of the associated RAN to support whole-stack
slicing (i.e., with each slice running potentially different PHY and MAC
layers).

For simplicity of exposition, here we assume that each RAN only has one BS
SDR, but the virtualization solution PV-RAN is readily extensible to the case
where each RAN may have multiple SDRs. As we will explain shortly, a key task
of PV-RAN is multiplexing/demultiplexing I/Q samples between the BS SDR and
the network stacks of different RAN slices. Therefore, the architecture
considered here resembles the RAN functional split Option 8 as defined by 3GPP
[21]. The PV-RAN design, however, is readily extensible to other functional
split options such as Options 7-1 and 7-2 where frequency-domain I/Q samples
need to be multiplex-/demultiplexed between the BS SDR and different network
stacks running on the BS computer. Similarly, the PV-RAN framework can be
extended to network slicing strategies involving Time Division Multiplexing
(TDM), and the network stacks could be based on other open-source cellular
software platforms such as srsLTE. These are all interesting topics for future
exploration, but their detailed studies are beyond the scope of this work.

Xen hypervisor. Considering the security benefits of virtual machines (VMs) as
compared with containers [22], as well as the facts that VMs can be made
lightweight [22] and that paravirtualization enables efficient I/O operations,
we use the Xen hypervisor [10] as the base virtualization platform in our PV-
RAN implementation.111Other paravirtualization platforms such as KVM may also
be used, but detailed study is beyond the scope of this work. In addition,
thanks to its microkernel architecture, Xen provides a greater degree of
isolation between the hypervisor and domains than KVM does. In particular, we
use Xen to manage the BS computer VMs used by individual slices of a RAN. Xen
hypervisor runs directly on top of the hardware and provides abstraction and
isolation required for the operation of multiple guest OSes hosted on the same
physical computer. A Xen-based virtualization system consists of the
hypervisor and the privileged domain Dom0 that has direct hardware access and
can manage one or more unprivileged guest domains (DomUs).

In Xen, Dom0 has privileged access to hardware (e.g., SDR) while DomUs are
typically not allowed to access hardware directly. However, Xen hypervisor
supports an I/O virtualization mechanism, known as the split-driver model,
allowing DomUs to access hardware through the use of virtual device drivers.
In particular, a front-end device driver lies in a DomU and a back-end device
driver resides in Dom0 and talks to the native device driver. The front-end
and back-end drivers use I/O ring buffers to communicate requests/replies. I/O
ring buffers rely on the grant table mechanism to share memory pages between
domains, and they use _event channels_ , an asynchronous notification
mechanism, to notify a domain when there are waiting data in the I/O ring
buffers. Thus, grant tables and event channels are convenient Xen facilities
that allow to create I/O channels to establish bi-directional communication
between domains. As we will explain in detail in Section IV, PV-RAN adopts the
split-driver model for efficient I/Q sample streaming between Dom0 and DomUs.

CyNet wireless living lab. CyNet [7] is a field-deployed wireless living lab,
consisting of two cellular RANs deployed at the Iowa State University Curtiss
Research Farm and Research Park in the City of Ames, Iowa, USA. Each RAN has
2-3 USRP X310s as the BS SDRs and several USRP B210s as UEs. This work is
motivated by the need for CyNet to partition its RANs into different slices
that employ different PHY and MAC layers to support research and education in
wireless, agriculture, and transportation respectively. PV-RAN has been tested
and deployed with the CyNet infrastructure.

##  IV PV-RAN Design and Implementation

In our implementation of PV-RAN, the BS SDR is the USRP X310, and the UE SDR
is the USRP B210. For ease of discussion, we will refer to these specific SDRs
and related software modules (e.g., the SDR driver UHD) in this section, but
the design and implementation strategies of PV-RAN are readily
applicable/extensible to other SDR platforms.

###  IV-A Architecture

Figure 3

Figure 3: Architecture of PV-RAN

shows the architecture of the PV-RAN. The software of the full OAI network
stacks of user slices are deployed in the DomUs of the Xen-virtualized BS
computer, and they communicate with the physical SDR via the PV-Front in the
individual DomUs; the virtual channel carries the I/Q samples through the Xen
libvchan interface; PV-Back is deployed in Dom0 to perform all the baseband
I/O operations from/to the USRP X310 via Dom0’s kernel native network driver;
USRP X310 integrates the antennas and acts as the RRH, and it communicates
with Dom0 using the USRP Hardware Driver (UHD). The PV-RAN controller in Dom0
acts as a middleware layer between the WebUI Slice Orchestrator and the Xen
toolstack xl by enforcing security and control. It enables the WebUI Slice
Orchestrator, a Web portal for software-defined orchestration and management
of DomUs, to send orchestration messages to the Xen toolstack that controls
DomUs. When an end-user requests the creation of a new DomU through the Web
portal, an orchestration message is sent to the PV-RAN controller which will
convert it into an xl command, which in turn will create the new DomU user VM.

In order to run multiple heterogeneous paravirtualized OAI protocol stacks in
different DomUs (i.e., user slices), each isolated OAI instance must appear to
have dedicated access to the SDR. Since only one process can access the USRP
Hardware Driver (UHD) at a time, we must decouple RF I/O (a.k.a. baseband I/O)
from the PHY layer to allow multiple concurrent OAI instances to
simultaneously access the USRP device. This is achieved by virtualizing the
SDR hardware resource and delegating all RF I/O operations to a shim layer
located in the Xen Dom0 and between the PHY layer of each paravirtualized OAI
instance and the actual UHD. This shim layer is called the _PV virtualization
layer_. The virtualization layer provides the following capabilities:

  * •

Slice isolation: the virtualization layer provides a dedicated virtual channel
to each paravirtualized OAI instance so that the OAI instances in different
user slices are isolated from one another. The virtual channels are created
using the libvchan interface, a Xen interface for inter-domain communication
(IDC) that mainly relies on gntalloc, a userspace grant allocation driver that
allocates shared pages, and gntdev that maps the granted pages into the
address space of a userland process.

  * •

Virtualization transparency: the OAI source code in user slices (i.e., DomUs)
does not need to be changed in order to use PV-RAN. This is accomplished
through _API Remoting_ , a software-based paravirtualization technique, that
virtualizes the SDR at the programming API level. In particular, the
virtualization layer intercepts SDR function calls at run-time from the OAI
instances running in DomUs, replaces them with API stubs, and forwards them to
Dom0 which in turn interacts with the physical SDR to actually execute the
function calls.

In what follows, we elaborate on the PV-RAN virtualization layer.

###  IV-B PV virtualization layer

Figure 4

Figure 4: PV-RAN virtualization layer

shows the internals of the PV virtualization layer. It consists of several
cooperating modules located in DomU and Dom0. In what follows, we first
present API Remoting for transparent virtualization, and then we present PV-
Back for SDR virtualization and inter-domain streaming of I/Q samples.
(Information on more detailed implementation strategies can be found in the
technical report [23].)

####  IV-B1 Transparent virtualization via API Remoting

OAI provides low-level interfaces to different types of SDR platforms such as
HackRF, BladeRF, and USRP. In this paper we focus the discussion on a commonly
used base station SDR — USRP X310, but the design and implementation
strategies can be applied to other SDRs. The USRP interface (usrp_lib.cpp) is
a software module that sits between the PHY layer and the SDR RF front-end. It
implements functions that rely on UHD API functions to configure the USRP X310
(channel, bandwidth, gains, etc.) and perform RF I/O operations (i.e., send
and receive time-domain I/Q samples). Since Dom0 interacts directly with the
USRP X310 via UHD, all the UHD API calls in DomU must be offloaded to Dom0.
One naive approach to addressing this challenge is to directly modify the
source code of OAI in each DomU. While OAI is an open-source software allowing
for modifications, a software-based mechanism that enables the sharing of an
SDR device among DomU guests without requiring OAI modification would be much
more desirable. To this end, we develop the API Remoting method that comprises
the following modules:

  * •

UHD API Remoting Frontend: this module is a dynamic library that relies on
library interposition. It intercepts at run-time each UHD library function
call and substitutes it with a stub function or wrapper that will execute
another piece of code that forwards the call through the UHD API Transport to
the UHD API Remoting Backend. The dynamic library is loaded during execution
into the memory space associated with OAI lte-softmodem process. This is
achieved through the use of the LD_PRELOAD environment variable that instructs
the loader to load the dynamic library. Each stub makes a call to dlsym() with
RTLD_NEXT to locate the address of the related symbol in memory. Since OAI
uses C++ bindings of the UHD API, we have to use mangled symbol names. Each
UHD API C++ function has a mangled name that can be found in the dynamic
symbol table. While the mangled symbol name is directly available for non-
virtual C++ functions, the procedure to retrieve the mangled symbol name of
virtual functions (e.g., UHD API setters and getters) is not trivial. It
involves complex lookup in the virtual table which we will cover in our future
work.

  * •

UHD API Transport: this transport component is a Xen event channel that allows
the UHD API Remoting Frontend in PV-Front to send short control messages to
the UHD API Remoting Backend in PV-Back. Each control message represents a UHD
library function that must be executed by the UHD in Dom0. For instance, upon
start-up, the OAI instance in DomU (in particular, usrp_lib.cpp) will first
issue the uhd::device::find() call to find the USRP device connected. This
call is then intercepted by the UHD API Remoting Frontend, and, upon
interception, the stub calls another function that sends an INIT request to
Dom0’s UHD API Remoting Backend, which in turn 1) informs PV-Back in Dom0 to
start streaming for this OAI slice and 2) returns to the UHD API Remoting
Frontend in DomU the type of USRP device connected.

  * •

UHD API Remoting Backend: this module handles the call requests received over
the UHD API Transport and executes them on the SDR device using the actual UHD
API. In its initial state, it waits for an incoming INIT request from a DomU’s
UHD API Remoting Frontend. This message is transmitted upon the execution of
lte-softmodem in a DomU. The UHD API Remoting Backend processes the message
and spawns RX/TX streamer threads which will be discussed shortly.

####  IV-B2 SDR virtualization through PV-Back in Dom0

Based on API remoting, PV-Back in Dom0 interacts with PV-Front in DomUs to
perform RF I/O operations on behalf of paravirtualized OAI instances in
different DomUs (i.e., user slices). It is designed to (i) initialize the USRP
device, configure the RX/TX frequencies, gains, bandwidths, and rates of the
communication channel according to OAI base station (e.g., eNodeB)
configuration in DomU, and map UHD RX/TX streamers to a dedicated SDR radio
channel; (ii) create bi-directional virtual channels and event channels for
I/Q sample streaming and UHD API transport respectively; and (iii) run high-
priority RX/TX streamer threads for the continuous I/Q sample streaming
between the USRP device and paravirtualized OAI instances in DomUs through the
stream-based libvchan communication interface (to be discussed shortly).

The RX/TX streamer threads are in charge of streaming I/Q samples. They are
created on-demand upon the reception of an INIT request triggered when lte-
softmodem is started in a DomU. Each paravirtualized OAI instance triggers the
creation of two streamer objects: one RX streamer object that is the Dom0
interface to receive I/Q samples, and one TX streamer object that is the Dom0
interface to transmit I/Q samples. PV-Back maps both streamers to one SDR
radio channel per DomU and spawns two threads for the reception/transmission
of I/Q samples. Besides continuously receiving and transmitting I/Q samples,
these two threads enable _lightweight, in-band timestamping synchronization_
between PV-Back in Dom0 and and PV-Front in DomU as follows:

  * •

RX streamer thread: it continuously reads nsamples (e.g., 7680 I/Q samples for
25-PRB channels) from the USRP device and then writes those samples to the
virtual channel. During the first run, this thread generates a timestamp
rx_timestamp and transmits it to PV-Front. This timestamp is then used by
OAI’s trx_usrp_read function as the time at which the first sample was
received. (This timestamp is also used to compute the tx_timestamp as we
discuss shortly.) After the first run, each time the thread performs an
iteration of USRP device read, it uses the previous rx_timestamp value and
increases it by a fixed number depending on channel bandwidth (e.g., 7680 for
25-PRB channels); so does that OAI instance (in particular, the trx_usrp_read
function) in the corresponding DomU. This way, the OAI instance in DomU and
the TX streamer thread in Dom0 are synchronized in rx_timestamp without any
explicit coordination/messaging after the first run.

  * •

TX streamer thread: it continuously reads nbytes (e.g., 30720 bytes for the
7680 I/Q samples for 25-PRB channels) from the virtual channel and then
transmits those samples to the USRP device. During the very first run, this
thread is locked and waits for the RX streamer thread to compute the
rx_timestamp. Once computed, the thread execution is resumed and it computes
the appropriate tx_timestamp, which is used by PV-Back to inform the UHD
send() function at what time the first sample must be sent. The tx_timestamp
is computed considering the channel bandwidth. For instance, for 25-PRB
channels, tx_timestamp is computed as rx_timestamp + 30640. The OAI instance
in the corresponding DomU (in particular, functions in lte-ru.c) also
calculates, in the same way as the TX streamer thread in PV-Back, the
tx_timestamp based on the rx_timestamp it receives from the RX streamer
thread, and it uses it in the OAI code execution in DomU. This way, the TX
streamer thread in Dom0 and the OAI instance in DomU are synchronized in
tx_timestamp without any explicit coordination/messaging. Then, for each
iteration of the virtual channel read, the thread calculates the tx_timestamp
again by using the previous tx_timestamp value and increasing it by a fixed
number (e.g., 7680 for 25-PRB channels).

####  IV-B3 Inter-domain streaming of I/Q samples

PV-Back in Dom0 and PV-Front in DomUs create shared memory pages to
communicate with each other. This functionality is available through the use
of the grant table mechanism that allows kernel memory pages to be shared
between domains. Each domain possesses a grant table which is a set of pages
shared between itself and the Xen hypervisor. The domains involved in this
mechanism are referred to as the client and the server, where the server
offers the memory used for communication and transmits its credentials (grant
reference, event channel ID) to the client in a Xen directory service known as
XenStore. The XenStore is a centralized database that stores key-value pairs
and relies on the Linux kernel interface XenBus for messaging. To get access
to the server’s shared memory, the client needs the server’s domain ID and the
XenStore path in which the server offered its credentials. This phase is known
as the rendezvous procedure.

Now, let’s observe the mechanism in more detail by considering a DomU $A$, the
client, wishing to communicate with another DomU $B$, the server:

  * •

DomU $B$ creates a page entry in its grant table and then advertises the index
of this entry (grant reference) to DomU $A$ through a dedicated Xen event
channel. The kernel driver gntalloc is the one in charge of creating grant
references.

  * •

Upon reception of the grant reference, DomU $A$ validates the grant and maps
the page via the gntdev device into the address space of the application
running in DomU $A$.

Once the rendezvous procedure has been completed, Xen makes use of shared data
structures called ring buffers for bulk data transfer. Historically those ring
buffers have been used by split drivers to communicate I/O requests-responses.

In this paper we use Xen virtual channel library libvchan, a library that
implements a datagram-based (packet-based and stream-based) interface on top
the standard Xen ring buffers. Libvchan allows to specify the size of the
rings and whether or not to perform blocking I/O operations on the virtual
channel. As OAI deals with streams of I/Q samples, we use the stream-based
communication interface that is composed of libxenvchan_read and
libxenvchan_write functions to read I/Q samples from a buffer and write I/Q
samples to a buffer respectively. The blocking I/O flag is set on the bi-
directional virtual channel to ensure that I/O operations are performed
sequentially.

##  V Measurement Evaluation

We have implemented the PA-RAN design [23] and integrated it with the CyNet
wireless living lab [7]. A recorded demo of the operational PV-RAN can be
found at [24]. Here we conduct detailed measurement benchmarking to
characterize the performance of the PV-RAN platform.

###  V-A PV-RAN testbed

Figure 5: PV-RAN testbed at ISU Research Park

We set up a testbed in the basement of the Iowa State University (ISU)
Research Park to conduct our measurement evaluation (see Figure 5), with the
basement having significantly less interference from surrounding wireless
equipment. The testbed is composed of:

  * •

One PV-RAN platform running two slices: the base station computer hosting PV-
RAN under the Xen hypervisor is a Dell Precision 3630 with 64GB of RAM and 8
CPU cores operating at 3.1GHz. Each DomU slice has been set up as a
paravirtualized guest under Ubuntu 16.04 and is configured with 4 vCPUs and
6GB of RAM. Each DomU slice runs OpenAirInterface as the LTE protocol stack.
Bi-directional virtual channels transport I/Q samples between each DomU PV-
Front and Dom0 PV-Back. Dom0 runs under Ubuntu 18.04 and has direct I/O access
to the USRP X310 which is connected to the base station computer through OS2
single-mode fiber.

  * •

Two UEs: each UE runs on an Intel NUC7i7BNB under Ubuntu 16.04 with a low-
latency kernel. Each NUC is connected via USB 3.0 to a USRP B210 equipped with
two VERT900 antennas that can operate in both the Educational Broadband
Service (EBS) and TV White Space (TVWS) frequencies.

###  V-B Network performance tuning

The USRP X310 is connected via single-mode fiber to one of the SFP+ network
interface of the base station computer which provides a 10Gbps network speed.
Table I summarizes the parameters that must be fine-tuned to allow the best
network performance possible. We used ifconfig and ethtool to configure the
network interface. It is recommended to use jumbo frames (for larger-packet
transmission support), to increase the size of RX/TX kernel buffers and
Ethernet buffer rings, to disable Ethernet frame pauses (RX/TX autoneg off),
and to minimize the raise of interrupts to the CPU during packet reception to
improve CPU utilization.

Network interface parameters
---
MTU size |  8000 bytes
Queue length |  1000 bytes
RX/TX autoneg |  off
Ethernet buffer rings |  4096 bytes
Interrupt moderation |  3 $\mu$secs
TABLE I: Hardware and software setup

###  V-C Throughput and Latency

To assess the smooth operation of PV-RAN, we compare the throughput of one PV-
RAN slice running OAI against a non-virtualized instance of OAI running in a
bare-metal server. The bare-metal server running OAI has the same
specification as the base station computer used for PV-RAN. Both the PV-RAN
OAI and Bare-metal OAI use the same LTE bandwidth of 25 PRBs (5MHz) and
operate at the LTE downlink frequency and uplink frequency of 2.685GHz and
2,565GHz respectively. The main eNB parameters are shown in Table II along
with UE parameters.

eNB configuration
---
LTE bandwidth | 25 PRBs
Downlink frequency | 2.685 GHz
Uplink frequency | 2.565 GHz
PDSCH reference signal power | -16 dB
UE configuration
RX gain | 120 dB
TX gain | 0 dB
Max UE power | 0 dB
TABLE II: eNB and UE configuration

We attach one UE to each eNB and measure the UDP throughput for 180 seconds
using the iperf tool. The throughput in the PV-RAN OAI and bare-metal OAI are
similar. For instance, Figure 6

Figure 6: Throughput comparison between PV-RAN OAI and Bare-metal OAI

depicts the throughput between each eNB and their associated UE in one single
run. The mean throughput in the PV-RAN OAI and bare-metal OAI are 6.24Mbps and
5.66Mbps, respectively; the difference is mostly due to environmental and
wireless channel dynamics.

Next, we evaluate the performance when two OAI slices are running on PV-RAN.
In this experiment, we measure the RTT and throughput between the OAI eNB and
UE in each slice. Each OAI eNB uses a channel bandwidth of 25 PRBs. Both
slices operate in TVWS bands: slice 1 runs at 595MHz (downlink) and 499MHz
(uplink), and slice 2 at 580MHz (downlink) and 484MHz (uplink). Each slice
occupies a different radio channel on the USRP X310. We connect one UE to each
slice and measure the RTT with the ping command and the UDP throughput with
the iperf command. The ICMP echo requests have been sent from each UE to their
respective eNB. Figure 7

Figure 7: RTT for each PV-RAN OAI slice

illustrates the RTT variation over time for each UE-eNB pair. The average RTT
is of 25ms for slice 1 and 29ms for slice 2. We witness similar trends with
OAI running on bare-metal servers. The throughput of each OAI slice is
depicted in Figure 8.

Figure 8: UDP Throughput for each PV-RAN OAI slice

It can be seen that slice 1 achieves a stable throughput of $\sim$6Mbps on
average while slice 2 achieves $\sim$4Mbps on average, which are also
comparable to throughput achieved with bare-metal OAI. The reason why slice 1
has slightly lower RTT and higher throughput than slice 2 is because the UE in
slice 1 is closer to its eNB and thus have better channel quality. (The slice
throughput is less stable at the beginning phase when UEs and eNBs dynamically
tune their communication parameters such as transmission power and MCS
schemes.)

###  V-D CPU utilization

It is of paramount importance to understand the CPU consumption of our PV-RAN
platform. For this study, all CPU cores ran at their maximum frequency in
Dom0. Each slice spawns two threads in PV-Back, namely stream-rx that receives
I/Q samples from the USRP device and writes them to the virtual channel, and
stream-tx that reads I/Q samples from the virtual channel and sends them to
the USRP device. We pinned each thread onto a single CPU core using
pthread_setaffinity_np to have a better understanding of CPU usage at a thread
level. In this experiment, the useful Linux task monitoring tools top and htop
have been used to collect CPU consumption of individual threads.

The CPU core usage per thread is reported in Figure 9

Figure 9: Core usage per thread for Dom0 PV-Back

for a single OAI slice running on PV-RAN. The main thread listens for incoming
PV-Front requests and is responsible for the creation of USRP device object
and virtual channel data structures. It has a mean CPU core consumption of
33.3%. It can be observed that an additional thread named zero-copy-recv is
spawned by the call to the UHD recv() function in the stream_rx thread. This
thread is one of the most influential thread in terms of CPU consumption for
OAI eNBs. It helps to minimize unnecessary buffer copy operations between the
NIC and userspace by bypassing buffer copy operations from the NIC to the
kernel space and from the kernel space to the userspace. stream-rx and zero-
copy-recv consume 20.7% and 20% respectively. These tasks play a major role in
CPU consumption, while in contrast the thread stream-tx has a mean consumption
of only 8% of its CPU core.

We also performed a comparative analysis of the total CPU consumption (over
all CPU cores) between OAI running in bare-metal and OAI running as a PV-RAN
slice. The CPU usage over all cores is depicted in Figure 10.

Figure 10: Mean CPU usage over all cores for different bandwidths

It can be seen that, for a LTE bandwidth of 25 PRBs, OAI eNB in bare-metal
consumes about 2.6% of CPU, and it consumes 3.2% of CPU when running as a PV-
RAN slice, a slight increase over bare-metal operation. When the LTE channel
bandwidth is doubled (i.e., with 50 PRBs), the CPU usage for OAI in bare-metal
increases by 0.85% while the CPU usage for OAI in PV-RAN increases by only
0.4%. So the CPU usage difference between bare-metal and PV-RAN operations
actually decreases with increasing channel bandwidth. Therefore, the CPU usage
overhead due to virtualization in PV-RAN is small, and tends to be negligible
as channel bandwidth increases.

###  V-E Latency overhead

To complete our measurement study, we focus on latency overhead which, if not
controlled to be low, can be detrimental to I/Q sample streaming and cellular
communication performance. Since we add an extra layer between OAI USRP
interface and UHD I/O functions to transfer I/Q samples, we need to evaluate
the latency overhead introduced by the PV virtualization layer. For this case
study, in particular, we add timestamps in PV-Front and PV-Back to measure the
latency overhead caused by the virtual channels. Figure 11

Figure 11: Latency overhead of inter-domain I/Q sample streaming

shows the latency overhead in PV-RAN. Surprisingly, for a channel bandwidth of
25 PRBs, the mean overhead introduces by the virtual channel is small and only
5$\mu$s on average. For a channel bandwidth of 50 PRBs which doubles the ring
buffer size in the virtual channel and tends to introduce higher latency, the
latency overhead is also only 6$\mu$s on average. This low latency overhead
and high efficiency of PV-RAN explain the good throughput performance of PV-
RAN discussed earlier.

Based on the above latency results, we decide to evaluate the theoretical
upper-bound capacity of our PV-RAN platform. More specifically, we verify how
the inter-domain I/Q sample streaming scales with respect to the number of OAI
slices for different LTE bandwidths. In addition, the Wireless Spectrum
Hypervisor [20] uses a ZMQ socket with the PUB-SUB (Publish-Subscribe) pattern
for I/Q sample streaming. Therefore, we also compare the PV-RAN method of
shared-memory inter-domain I/Q sample streaming against the ZMQ PUB-SUB I/Q
sample streaming method. To this end, we create a variant of PV-RAN that uses
a ZMQ socket to publish I/Q samples between PV-Front and PV-Back.

The latency overhead observed for ZMQ PUB-SUB varies between 65$\mu$s and
90$\mu$s, significantly higher than that with the PV-RAN method. From these
results we can infer the maximum capacity for ZMQ PUB-SUB between Dom0 and
DomU too. Figure 12

Figure 12: OAI slice scalability for different LTE bandwidths

shows for different LTE bandwidths that shared-memory clearly outperforms ZMQ
PUB-SUB, with shared-memory inter-domain I/Q sample streaming being 7 times
faster than ZMQ PUB-SUB inter-domain I/Q sample streaming on average. While
ZMQ PUB-SUB remains a viable solution for general communications between
processes on the same host, it is not a good candidate for intensive real-time
inter-domain I/Q sample streaming.

##  VI Concluding Remarks

As the first open-source platform for physical wireless resource
virtualization, PV-RAN enables whole-stack slicing where different PHY and MAC
layers may be adopted for diverse communication services and wireless living
lab innovations. Using efficient techniques such as API Remoting, shared-
memory I/Q sample streaming, and light-weight in-band timestamping
synchronization, PV-RAN enables transparent, lightweight virtualization of
SDRs, and it serves a solid foundation for further development. For instance,
the PV-RAN platform can be extended to support SDR virtualization using TDM
and hybrid-TDM-FDM, and the PV-RAN implementation can be ported to other open-
source platforms such as KVM and srsLTE to expand community access to PV-RAN
capabilities.

## References

  * [1] F. Alvarez, D. Breitgand _et al._ , “An edge-to-cloud virtualized multimedia service platform for 5G networks,” _IEEE Transactions on Broadcasting_ , vol. 65, no. 2, pp. 369–380, 2019.
  * [2] H. Zhang, N. Liu, X. Chu, K. Long, A.-H. Aghvami, and V. C. Leung, “Network slicing based 5g and future mobile networks: mobility, resource management, and challenges,” _IEEE Communications Magazine_ , vol. 55, no. 8, pp. 138–145, 2017.
  * [3] S. Peng, R. Chen, and G. Mirsky, “Packet network slicing using segment routing,” Tech. Rep. Draft-penglsr-network-slicing-00. IETF, Tech. Rep., 2019\.
  * [4] “flexran: A flexible and programmable platform for software-defined radio access networks.”
  * [5] A. Esmaeily, K. Kralevska, and D. Gligoroski, “A Cloud-based SDN/NFV Testbed for End-to-End Network Slicing in 4G/5G,” _arXiv:2004.10455_ , 2020.
  * [6] M. Yang, Y. Li, D. Jin, L. Zeng, X. Wu, and A. V. Vasilakos, “Software-Defined and Virtualized Future Mobile and Wireless Networks: A Survey,” _Mobile Networks and Applications_ , vol. 20, 2015.
  * [7] “CyNet: End-to-End Software-Defined Cyberinfrastruture for Smart Agriculture and Transportation,” https://www.ece.iastate.edu/~hongwei/group/projects/CyNet.html.
  * [8] Y. Xie, H. Zhang, and P. Ren, “Unified scheduling for predictable communication reliability in industrial cellular networks,” in _IEEE ICII_ , 2018.
  * [9] E. Research, “USRP Software-Defined Ratios,” https://www.ettus.com/product.
  * [10] P. Barham, B. Dragovic, K. Fraser, S. Hand, T. Harris, A. Ho, R. Neugebauer, I. Pratt, and A. Warfield, “Xen and the Art of Virtualization,” in _ACM SOSP_ , 2003.
  * [11] N. Nikaein, M. K. Marina, S. Manickam, A. Dawson, R. Knopp, and C. Bonnet, “OpenAirInterface: A Flexible Platform for 5G Research,” _ACM SIGCOMM Computer Communication Review_ , vol. 44, no. 5, pp. 33–38, 2014.
  * [12] V. Suryaprakash, P. Rost, and G. Fettweis, “Are heterogeneous cloud-based radio access networks cost effective?” _IEEE Journal on Selected Areas in Communications_ , vol. 33, no. 10, pp. 2239–2251, 2015.
  * [13] A. Garcia-Saavedra, X. Costa-Perez, D. J. Leith, and G. Iosifidis, “FluidRAN: Optimized vran/mec orchestration,” in _IEEE INFOCOM_ , 2018.
  * [14] N. Nikaein, E. Schiller, R. Favraud, R. Knopp, I. Alyafawi, and T. Braun, “Towards a cloud-native radio access network,” in _Advances in mobile cloud computing and big data in the 5G era_. Springer, 2017, pp. 171–202.
  * [15] F. Kaltenberger, A. P. Silva, A. Gosain, L. Wang, and T.-T. Nguyen, “Openairinterface: Democratizing innovation in the 5g era,” _Computer Networks_ , p. 107284, 2020.
  * [16] T. X. Tran, A. Younis, and D. Pompili, “Understanding the computational requirements of virtualized baseband units using a programmable cloud radio access network testbed,” in _IEEE ICAC_ , 2017.
  * [17] I. Sanchez-Navarro, A. S. Mamolar, Q. Wang, and J. M. A. Calero, “5gtoponet: Real-time topology discovery and management on 5g multi-tenant networks,” _Future Generation Computer Systems_ , vol. 114, pp. 435–447, 2020.
  * [18] I. Trindade, C. Nahum, C. Novaes, D. Cederholm, G. Patra, and A. Klautau, “C-ran virtualization with openairinterface,” _arXiv preprint arXiv:1908.07503_ , 2019.
  * [19] N. Nikaein, “Processing radio access network functions in the cloud: Critical issues and modeling,” in _ACM MCS_ , 2015.
  * [20] F. A. de Figueiredo, R. Mennes, I. Jabandžić, X. Jiao, and I. Moerman, “A baseband wireless spectrum hypervisor for multiplexing concurrent ofdm signals,” _Sensors_ , vol. 20, no. 4, p. 1101, 2020.
  * [21] “Study on new radio access technology: Radio access architecture and interfaces,” _3GPP TR38.801 Release 14_.
  * [22] F. Manco, J. Mendes, K. Yasukata, C. Lupu, S. Kuenzer, C. Raiciu, F. Schmidt, S. Sati, and F. Huici, “My VM is Lighter (and Safer) than your Container,” in _ACM SOSP_ , 2017.
  * [23] M. Sander-Frigau, T. Zhang, H. Zhang, A. E. Kamal, and A. K. Somani, “Physical wireless resource virtualization for software-defined whole-stack slicing,” Iowa State University, Tech. Rep. ISU-DNC-TR-20-02 (https://www.ece.iastate.edu/~hongwei/group/publications/PV-RAN-TR.pdf), 2020\.
  * [24] “PV-RAN Demo,” https://www.ece.iastate.edu/~hongwei/group/projects/CyNet/PV-RAN-CyNet-Demo.mp4.