# Ioffe-Regel criterion and viscoelastic properties of amorphous solids

D. A. Conyuh  Y. M. Beltukov<EMAIL_ADDRESS>Ioffe Institute, 194021
Saint-Petersburg, Russia

###### Abstract

We show that viscoelastic effects play a crucial role in the damping of
vibrational modes in harmonic amorphous solids. The relaxation of a given
plane wave is described by a memory function of a semi-infinite one-dimensions
mass-spring chain. The initial vibrational energy spreads from the first site
of the chain to infinity. In the beginning of the chain, there is a barrier,
which significantly reduces the decay of vibrational energy below the Ioffe-
Regel frequency. To obtain the parameters of the chain, we present a
numerically stable method, based on the Chebyshev expansion of the local
vibrational density of states.

##  I Introduction

Damping of vibrational modes plays a crucial role in the thermal conductivity
of amorphous dielectrics (glasses). Low-frequency vibrations are well-defined
phonons with a long mean free path. However, in a wide range of temperatures,
the heat transfer in glasses is determined by another type of delocalized
vibrations, which are known as diffusons Allen and Feldman (1993); Allen  _et
al._ (1999). The crossover between low-frequency phonons and diffusons at
higher frequencies is known as the Ioffe-Regel crossover Allen  _et al._
(1999); Beltukov  _et al._ (2013).

In amorphous solids, the attenuation of plane elastic waves (sound) is
governed by multiple mechanisms: scattering on two-level systems Anderson  _et
al._ (1972); Phillips (1972); Jäckle (1972) and soft modes Buchenau  _et al._
(1992); Ji  _et al._ (2019), thermally activated relaxation processes Jäckle
_et al._ (1976); Tielbürger  _et al._ (1992), and scattering induced by
structural and elastic disorder Ruocco  _et al._ (1999); Dell’Anna  _et al._
(1998); Beltukov  _et al._ (2016, 2018); Gelin  _et al._ (2016). The last
contribution is temperature independent and dominates the attenuation in the
THz frequency range Ruocco  _et al._ (1999); Damart  _et al._ (2017).

In the low-frequency range, there are phonons with a well-defined dispersion
law $\omega(\mathbf{q})$ and a weak damping
$\Gamma(\mathbf{q})\ll\omega(\mathbf{q})$. In this case, the initial plane
wave with the wavevector $\mathbf{q}$ oscillates with the frequency
$\omega(\mathbf{q})$ with a slow exponential decay. This attenuation can be
described using the damped harmonic oscillator (DHO) model.

However, the damping increases rapidly with increasing the wavevector
$\mathbf{q}$. For some wavevector $|\mathbf{q}|=q_{c}$, the damping becomes
comparable to the frequency, $\Gamma(\mathbf{q})\sim\omega(\mathbf{q})$, which
corresponds to the Ioffe-Regel criterion. For $q\gtrsim q_{c}$, the notion of
the dispersion law $\omega(\mathbf{q})$ could not be applied. It was shown
that the DHO model can not be used for an accurate analysis of vibrational
properties for frequencies above the Ioffe-Regel crossover Baldi  _et al._
(2016); Buchenau (2014). Viscoelastic properties are important to study the
high-frequency vibrations of amorphous solids Luo  _et al._ (2020);
Ranganathan  _et al._ (2017); Lemaître and Maloney (2006).

In the theory of viscoelastic relaxation of liquids, it is known that memory
effects are important for the relaxation of density fluctuations Hansen and
McDonald (1990). These memory effects can be presented using Mori continued
fraction Mori (1965). In this paper we show that the vibrational relaxation in
harmonic amorphous solids can also be described using a general viscoelastic
model with some memory function $K(t)$. In terms of vibrations, the continued
fraction representation corresponds to a semi-infinite mass-spring chain,
which reproduces the same memory effects. We present a stable method to find
parameters of the arbitrary number of sites in the mass-spring chain.

Another powerful tool to analyze the general properties of disordered systems
is the random matrix theory (RMT). Depending on the inherent symmetry
properties of different disordered systems, various random matrix ensembles
are used Evers and Mirlin (2008). It was shown that the Wishart ensemble
naturally arises in the study of vibrational properties due to the requirement
of mechanical stability Beltukov  _et al._ (2013); Beltukov (2015); Conyuh and
Beltukov (2020). In this paper we apply the RMT to find the memory function
and the corresponding representation using the mass-spring chain.

This paper is organized as follows. In Section II we consider a general
viscoelastic relaxation of vibrations in harmonic amorphous solids. Section
III demonstrates that the relaxation is described by a Green function in the
form of continued fraction, which corresponds to a one-dimensional mass-spring
chain. In Section IV we obtain the parameters of the chain in the framework of
the RMT. Section V demonstrates that the same approach can be used to analyze
any given numerical dynamical matrix. In Section VI we discuss the properties
of obtained mass-spring chains and compare them to the Ioffe-Regel crossover.

##  II Viscoelastic damping and the memory function

The general equation of motion of a solid near equilibrium position can be
written as

$|\ddot{u}(t)\rangle=-\hat{M}|u(t)\rangle,$ (1)

where $\hat{M}$ is $N\times N$ dynamical matrix with $N$ being the number of
degrees of freedom. $N$-dimensional vector $|u(t)\rangle$ describes the
deviation of atoms from the equilibrium position at time $t$.

To study the relaxation of a plane wave with wavevector $\mathbf{q}$, we can
solve Eq. (1) with initial conditions

$|u(0)\rangle=0,\quad|\dot{u}(0)\rangle=|\mathbf{q}\rangle$ (2)

for any given dynamical matrix $\hat{M}$. The relaxation of the initial plane
wave is described by the projection

$u_{\mathbf{q}}(t)=\Big{\langle}\langle u(t)|\mathbf{q}\rangle\Big{\rangle},$
(3)

where the big angle brackets denote the averaging over different realizations
of the dynamical matrix $\hat{M}$. Assuming the normalization
$\langle\mathbf{q}|\mathbf{q}\rangle=1$, we obtain $u_{\mathbf{q}}(0)=0$ and
$\dot{u}_{\mathbf{q}}(0)=v_{0}$ from Eq. (2).

In a simplified model, the relaxation of $u_{\mathbf{q}}(t)$ can be described
using the DHO model

$m_{\mathbf{q}}\ddot{u}_{\mathbf{q}}(t)+\eta_{\mathbf{q}}\dot{u}_{\mathbf{q}}(t)+k_{\mathbf{q}}u_{\mathbf{q}}(t)=0,$
(4)

where the mass $m_{\mathbf{q}}$, the damping $\eta_{\mathbf{q}}$, and the
stiffness $k_{\mathbf{q}}$ may depend on the wavevector $\mathbf{q}$. In terms
of the DHO model, the Ioffe-Regel crossover separates weakly decaying long-
wave vibrational modes and overdamped short-wave vibrational modes. However,
it is important to take into account the frequency dependence of the damping,
which results in a nonlocal-in-time viscoelastic equation. We will also take
into account that the initial equation (1) is time-reversal and do not have
the inherent energy dissipation.

Figure 1: (Color online) a) Damped harmonic oscillator with the mass
$m_{\mathbf{q}}$ (the ball), the stiffness $k^{\rm inst}_{\mathbf{q}}$ (the
spring), and the damping $\eta_{\mathbf{q}}$ (the dashpot). b) Harmonic
oscillator with a general viscoelastic element defined by the memory function
$K_{\mathbf{q}}(t)$ (the rectangle).

To analyze the relaxation process, we consider the resolvent

$\hat{G}(z)=\left\langle\frac{1}{z-\hat{M}}\right\rangle,$ (5)

where $z$ is a complex parameter. The relaxation of a plane wave with initial
conditions ${u_{\mathbf{q}}(0)=0}$, ${\dot{u}_{\mathbf{q}}(0)=v_{0}}$ can be
written using the resolvent $\hat{G}(z)$ as (see Appendix A)

$u_{\mathbf{q}}(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\tilde{u}_{\mathbf{q}}(\omega)e^{i\omega
t}d\omega,$ (6)

where

$\displaystyle\tilde{u}_{\mathbf{q}}(\omega)=-v_{0}\mathcal{G}_{\mathbf{q}}\big{(}(\omega-i0)^{2}\big{)},$
(7)
$\displaystyle\mathcal{G}_{\mathbf{q}}(z)=\langle\mathbf{q}\big{|}\hat{G}(z)|\mathbf{q}\rangle.$
(8)

We can present the Green function $\mathcal{G}_{\mathbf{q}}(z)$ as the
Stieltjes transform of the spatial Fourier transform of eigenmodes:

$\displaystyle\mathcal{G}_{\mathbf{q}}(z)=\int_{0}^{\infty}\frac{{\cal
F}_{\mathbf{q}}(\omega)}{z-\omega^{2}}d\omega,$ (9) $\displaystyle{\cal
F}_{\mathbf{q}}(\omega)=\left\langle\sum_{n}\langle\mathbf{q}|n\rangle\langle
n|\mathbf{q}\rangle\delta(\omega-\omega_{n})\right\rangle,$ (10)

where $|n\rangle$ is $n$-th eigenmode. The Fourier transform of eigenmodes is
closely related to the structure factor, which is
$S_{\mathbf{q}}(\omega)=k_{B}Tq^{2}/(m\omega^{2})\mathcal{F}_{\mathbf{q}}(\omega)$
Shintani and Tanaka (2008); Beltukov  _et al._ (2016). In this paper we use
the scalar model to simplify the notations. All qualitative results that we
obtain can be applied to the vector model as well. It was shown that
vibrations in the scalar and vector models belong to the same class of
universality Skipetrov and Beltukov (2018).

In the DHO model with the natural frequency $\omega_{\mathbf{q}}$ and the
frequency-independent damping rate $\Gamma_{\mathbf{q}}$, we have

$\displaystyle\mathcal{F}^{\textsc{dho}}_{\mathbf{q}}(\omega)$
$\displaystyle=\frac{2}{\pi}\frac{\omega^{2}\Gamma_{\mathbf{q}}}{\big{(}\omega^{2}-\omega_{\mathbf{q}}^{2}\big{)}^{2}+\omega^{2}\Gamma_{\mathbf{q}}^{2}},$
(11) $\displaystyle\mathcal{G}^{\textsc{dho}}_{\mathbf{q}}(z)$
$\displaystyle=-\frac{1}{\omega_{\mathbf{q}}^{2}-z+\Gamma_{\mathbf{q}}\sqrt{-z}}.$
(12)

It corresponds to Eq. (4) with the stiffness
$k_{\mathbf{q}}=m_{\mathbf{q}}\omega_{\mathbf{q}}^{2}$ and the damping
$\eta_{\mathbf{q}}=m_{\mathbf{q}}\Gamma_{\mathbf{q}}$. In the DHO model, the
mass $m_{\mathbf{q}}$ can be chosen arbitrarily.

In the general case, we can present the Green function
$\mathcal{G}_{\mathbf{q}}(z)$ as

$\mathcal{G}_{\mathbf{q}}(z)=-\frac{m_{\mathbf{q}}}{k^{\rm
inst}_{\mathbf{q}}-m_{\mathbf{q}}z+\mathcal{G}_{1\mathbf{q}}(z)}.$ (13)

It corresponds to the viscoelastic equation of motion

$m_{\mathbf{q}}\ddot{u}_{\mathbf{q}}(t)+k^{\rm
inst}_{\mathbf{q}}u_{\mathbf{q}}(t)+\int_{-\infty}^{t}K_{\mathbf{q}}(t-t^{\prime})u_{\mathbf{q}}(t^{\prime})dt^{\prime}=0$
(14)

with the mass $m_{\mathbf{q}}$, the instantaneous stiffness $k^{\rm
inst}_{\mathbf{q}}$, and the memory function

$K_{\mathbf{q}}(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathcal{G}_{1\mathbf{q}}\big{(}(\omega-i0)^{2}\big{)}e^{i\omega
t}d\omega.$ (15)

In order to find $m_{\mathbf{q}}$, $k^{\rm inst}_{\mathbf{q}}$, and
$K_{\mathbf{q}}(t)$, we can present the Green function
$\mathcal{G}_{\mathbf{q}}(z)$ as a series

$\mathcal{G}_{\mathbf{q}}(z)=\sum_{k=0}^{\infty}\frac{{\cal
F}_{\mathbf{q}}^{(k)}}{z^{k+1}}$ (16)

with moments

$\mathcal{F}_{\mathbf{q}}^{(k)}=\int_{0}^{\infty}\omega^{2k}\mathcal{F}_{\mathbf{q}}(\omega)d\omega.$
(17)

For large $z$, $\mathcal{G}_{\mathbf{q}}(z)$ is $1/z$ due to the normalization
${\cal F}_{\mathbf{q}}^{\smash{(0)}}=\int_{0}^{\infty}{\cal
F}_{\mathbf{q}}(\omega)d\omega=1$. For any harmonic system described by Eq.
(1), all moments $\mathcal{F}_{\mathbf{q}}^{(k)}$ are finite. In this case, we
can assume that $\mathcal{G}_{1\mathbf{q}}(z)$ is also $1/z$ for large $z$. It
results in the following values

$\displaystyle k^{\rm inst}_{\mathbf{q}}$ $\displaystyle=m_{\mathbf{q}}{\cal
F}_{\mathbf{q}}^{(1)},$ (18) $\displaystyle m_{\mathbf{q}}$
$\displaystyle=\Big{[}\mathcal{F}_{\mathbf{q}}^{(2)}-\big{(}\mathcal{F}_{\mathbf{q}}^{(1)}\big{)}^{2}\Big{]}^{-1},$
(19) $\displaystyle\mathcal{G}_{1\mathbf{q}}(z)$
$\displaystyle=m_{\mathbf{q}}z-k^{\rm
inst}_{\mathbf{q}}-\frac{m_{\mathbf{q}}}{\mathcal{G}_{\mathbf{q}}(z)}.$ (20)

The decreasing of $\mathcal{G}_{1\mathbf{q}}(z)$ for large $z$ corresponds to
the absence of the instantaneous component in the memory function
$K_{\mathbf{q}}(t)$.

##  III Continued fraction and one-dimensional chain

We can repeatedly apply the same type of presentation for the Green function:

$\mathcal{G}_{n+1,\mathbf{q}}(z)=m_{n\mathbf{q}}z-a_{n\mathbf{q}}-\frac{b_{n\mathbf{q}}^{2}}{\mathcal{G}_{n\mathbf{q}}(z)}$
(21)

with starting Green function $\mathcal{G}_{1\mathbf{q}}(z)$ defined by Eq.
(20). As a result, we obtain the Mori continued fraction Mori (1965):

$\displaystyle\mathcal{G}_{\mathbf{q}}(z)=-\frac{m_{\mathbf{q}}}{k^{\rm
inst}_{\mathbf{q}}-m_{\mathbf{q}}z+\mathcal{G}_{1\mathbf{q}}(z)}$
$\displaystyle=-\frac{m_{\mathbf{q}}}{k^{\rm
inst}_{\mathbf{q}}-m_{\mathbf{q}}z-\frac{b_{1\mathbf{q}}^{2}}{a_{1\mathbf{q}}-m_{1\mathbf{q}}z+\mathcal{G}_{2\mathbf{q}}(z)}}$
$\displaystyle=-\frac{m_{\mathbf{q}}}{k^{\rm
inst}_{\mathbf{q}}-m_{\mathbf{q}}z-\frac{b_{1\mathbf{q}^{2}}}{a_{1\mathbf{q}}-m_{1\mathbf{q}}z-\frac{b_{2\mathbf{q}^{2}}}{a_{2\mathbf{q}}-m_{2\mathbf{q}}z+\mathcal{G}_{3\mathbf{q}}(z)}}}=\dots$
(22)

For each Green function $\mathcal{G}_{n\mathbf{q}}(z)$ we can find the
corresponding function $\mathcal{F}_{n\mathbf{q}}(\omega)$, the moments
$\mathcal{F}_{n\mathbf{q}}^{(k)}$, and the memory function
$K_{n\mathbf{q}}(t)$ using the same relations as in Eqs. (9), (15)–(17). A
comprehensive set of relations is given in Appendix B.

As before, we assume that each Green function $\mathcal{G}_{n\mathbf{q}}(z)$
is $1/z$ for large $z$, which determines the relation between coefficients
$a_{n}$, $b_{n}$, and $m_{n}$:

$\displaystyle a_{n\mathbf{q}}=m_{n\mathbf{q}}{\cal
F}_{n\mathbf{q}}^{(1)},\quad b_{n\mathbf{q}}^{2}=m_{n\mathbf{q}},$ (23)
$\displaystyle
m_{n\mathbf{q}}=\Big{[}\mathcal{F}_{n\mathbf{q}}^{(2)}-\big{(}\mathcal{F}_{n\mathbf{q}}^{(1)}\big{)}^{2}\Big{]}^{-1}.$
(24)

The recurrence relation (21) can be presented in the form

$\displaystyle\mathcal{F}_{n+1,\mathbf{q}}(\omega)=\frac{m_{n\mathbf{q}}\mathcal{F}_{n\mathbf{q}}(\omega)}{\left|\mathcal{G}_{n\mathbf{q}}\big{(}(\omega-i0)^{2}\big{)}\right|^{2}}.$
(25)

The mass $m_{n\mathbf{q}}$ can be defined by the normalization condition
${\cal
F}_{n+1,\mathbf{q}}^{(0)}=\int_{0}^{\infty}\mathcal{F}_{n+1,\mathbf{q}}(\omega)d\omega=1$,
which is equivalent to Eq. (24). The details of the numerical realization of
this recurrence procedure are discussed in Appendix D.

Figure 2: (Color online) a) A general case of the semi-infinite one-
dimensional mass-spring model. The displacements of masses
$m_{\mathbf{q}},m_{1\mathbf{q}},m_{2\mathbf{q}},\dots$ are denoted by
$u_{\mathbf{q}},u_{1\mathbf{q}},u_{2\mathbf{q}},\dots$ respectively. Each
spring is denoted by its stiffness. b) Finite representation of the same chain
with 3 sites. The tail of the chain is replaced by a general viscoelastic
element defined by the memory function $K_{3\mathbf{q}}(t)$ and denoted by the
rectangle.

The continued fraction (22) can be presented using a solution of the following
infinite system:

$\displaystyle m_{\mathbf{q}}\ddot{u}_{\mathbf{q}}(t)$ $\displaystyle=-k^{\rm
inst}_{\mathbf{q}}u_{\mathbf{q}}(t)+b_{1\mathbf{q}}u_{1\mathbf{q}}(t),$ (26)
$\displaystyle m_{n\mathbf{q}}\ddot{u}_{n\mathbf{q}}(t)$
$\displaystyle=-a_{n\mathbf{q}}u_{n\mathbf{q}}(t)+b_{n\mathbf{q}}u_{n-1,\mathbf{q}}(t)$
$\displaystyle\qquad\qquad+b_{n+1,\mathbf{q}}u_{n+1,\mathbf{q}}(t),\quad n\geq
1,$ (27)

where $u_{0\mathbf{q}}(t)\equiv u_{\mathbf{q}}(t)$ and

$\displaystyle u_{n\mathbf{q}}(t)$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}\tilde{u}_{n\mathbf{q}}(\omega)e^{i\omega
t}d\omega,$ (28) $\displaystyle\tilde{u}_{n\mathbf{q}}(\omega)$
$\displaystyle=-\frac{\tilde{u}_{n-1,\mathbf{q}}(\omega)}{b_{n\mathbf{q}}}\mathcal{G}_{n\mathbf{q}}\big{(}(\omega-i0)^{2}\big{)}.$
(29)

It corresponds to the dynamics of a semi-infinite one-dimensional mass-spring
chain (Fig. 2a) with masses $m_{\mathbf{q}}$ and $m_{n\mathbf{q}}$, horizontal
springs with stiffnesses $b_{n\mathbf{q}}$ between consequent masses, and
vertical springs with stiffnesses

$\displaystyle k_{\mathbf{q}}$ $\displaystyle=k^{\rm
inst}_{\mathbf{q}}-b_{1\mathbf{q}},$ (30) $\displaystyle k_{n\mathbf{q}}$
$\displaystyle=a_{n\mathbf{q}}-b_{n\mathbf{q}}-b_{n+1,\mathbf{q}}$ (31)

between masses and the common ground. The initial condition of the chain is
$u_{\mathbf{q}}(0)=0$, $\dot{u}_{\mathbf{q}}(0)=v_{0}$,
$u_{n\mathbf{q}}(0)=0$, $\dot{u}_{n\mathbf{q}}(0)=0$, which defines
$u_{\mathbf{q}}(t)$ for $t>0$.

Using Eq. (29), the last term in Eq. (27) can be replaced by the corresponding
memory function, which is known as the Mori-Zwanzig procedure Hansen and
McDonald (1990)

$\displaystyle
b_{n+1,\mathbf{q}}u_{n+1,\mathbf{q}}(t)=-\int_{-\infty}^{t}K_{n+1,\mathbf{q}}(t-t^{\prime})u_{n\mathbf{q}}(t^{\prime})dt^{\prime}.$
(32)

It corresponds to a finite chain shown in Fig. 2b. Therefore, the memory
function $K_{n\mathbf{q}}(t)$ describes a response of the tail of the mass-
spring chain starting from the site with the number $n$. Using the relations
given in Appendix B, the corresponding function
$\mathcal{F}_{n\mathbf{q}}(\omega)$ can be considered as a local vibrational
density of states (LVDOS) on site $n$ for the chain under constraint
$u_{k\mathbf{q}}(t)=0$ for $k<n$.

##  IV A random matrix approach

General vibrational properties could be studied using the random matrix
approach Conyuh and Beltukov (2020). This approach is based on two main
properties of amorphous solids: mechanical stability and the invariance under
continuous translation of an amorphous body. In the beginning of this section,
we briefly discuss the main points of the random matrix approach.

The mechanical stability of amorphous solids is equivalent to the positive
definiteness of the dynamical matrix $\hat{M}$. Any positive definite matrix
$\hat{M}$ can be written as $\hat{M}=\hat{A}\hat{A}^{T}$ and vice versa,
$\hat{A}\hat{A}^{T}$ is positive definite for any (not necessarily square)
matrix $\hat{A}$ Bhatia (2009). Therefore, we can consider a $N\times K$
random matrix $\hat{A}$ to obtain a mechanically stable system with the
dynamical matrix in the form of the Wishart ensemble
$\hat{M}=\hat{A}\hat{A}^{T}$.

Each column of the matrix $\hat{A}$ represents a bond with a positive
potential energy $U_{k}=\frac{1}{2}\big{(}\sum_{i}A_{ik}u_{i}\big{)}{}^{2}$
with $u_{i}$ being the displacement of $i$-th atom from the equilibrium
position Beltukov and Parshin (2016); Conyuh and Beltukov (2020). Each row of
the matrix $\hat{A}$ corresponds to some degree of freedom. In the random
matrix approach, the parameter

$\varkappa=\frac{K}{N}-1$ (33)

plays a crucial role. It is a relative difference between the number of bonds
$K$ and the number of degrees of freedom $N$. In a stable system with a finite
rigidity, the number of bonds should be larger than the number of degrees of
freedom, which is known as the Maxwell counting rule. For $0<\varkappa\lesssim
1$, the parameter $\varkappa$ has the same effect as the parameter $z-z_{c}$
in the jamming transition O’Hern  _et al._ (2003). In real amorphous solids,
one can estimate $\varkappa=0.3$ – 1 depending on the number and the type of
covalent bonds Conyuh and YM (2020).

The second important mechanical property of amorphous solid is the invariance
under continuous translation. It means that the bond energy $U_{k}$ should not
depend on the shift $u_{i}\to u_{i}+const$. Therefore, the matrix $\hat{A}$
obeys the _sum rule_ $\sum_{i}A_{ik}=0$. It means that the matrix elements
$A_{ij}$ are _correlated_. In the minimal model, we can assume that amorphous
solid consists of statistically equivalent random bonds. In this case the
pairwise correlations between matrix elements $A_{ij}$ can be written as

$\langle A_{ik}A_{jl}\rangle=\frac{1}{N}C_{ij}\delta_{kl},$ (34)

where $\hat{C}$ is some correlation matrix. One can see that the correlation
matrix $\hat{C}$ is proportional to the average dynamical matrix:
$\hat{C}=\tfrac{N}{K}\big{\langle}\hat{M}\big{\rangle}$. For a system with
statistically equivalent bonds, the correlation matrix $\smash{\hat{C}}$ is a
regular matrix, which describes a lattice with some dispersion law
$\omega_{\textrm{cor}}(\mathbf{q})$. We assume that there is only one branch
of $\omega_{\textrm{cor}}(\mathbf{q})$. In the general case, one can apply the
summation over different branches below.

Using the random matrix approach, it can be shown that statistical properties
of the random matrix $\hat{M}$ are related to the known correlation matrix
$\hat{C}$. In the thermodynamic limit $N\to\infty$, there is a fundamental
duality relation Burda  _et al._ (2004)

$z\hat{G}(z)=Z\hat{G}_{\textrm{cor}}(Z)$ (35)

between resolvents $\hat{G}(z)=\big{\langle}(z-\hat{M})^{-1}\big{\rangle}$ and
$\hat{G}_{\textrm{cor}}(Z)=(Z-\hat{C})^{-1}$ where complex parameters $z$ and
$Z$ are related by a conformal map $Z(z)$ defined by the equation

$Z\big{(}\varkappa+1+\mathcal{M}_{\textrm{cor}}(Z)\big{)}=z.$ (36)

The contour $Z(z)$ in the complex plane for $z=(\omega-i0)^{2}$ is known as
the critical horizon Burda  _et al._ (2006). The moment-generating function
$\mathcal{M}_{\textrm{cor}}(Z)$ for the correlation matrix $\hat{C}$ is
defined as

$\mathcal{M}_{\textrm{cor}}(Z)=\frac{Z}{N}\operatorname{Tr}\hat{G}_{\textrm{cor}}(Z)-1=\sum_{k=1}^{\infty}\frac{\mathcal{M}_{\rm
cor}^{(k)}}{Z^{k}}$ (37)

with moments

$\mathcal{M}_{\rm
cor}^{(k)}=\frac{1}{N}\operatorname{Tr}\hat{C}^{k}=\frac{1}{V_{\textsc{bz}}}\int_{\textsc{bz}}\omega_{\textrm{cor}}^{2k}(\mathbf{q})d\mathbf{q}.$
(38)

Here integration is performed over the first Brillouin zone, which has the
volume $V_{\textsc{bz}}$. The moment-generating function
$\mathcal{M}_{\textrm{cor}}(Z)$ can also be written as an integral over the
first Brillouin zone

$\displaystyle\mathcal{M}_{\textrm{cor}}(Z)=\frac{1}{V_{\textsc{bz}}}\int_{\textsc{bz}}\frac{\omega_{\textrm{cor}}^{2}(\mathbf{q})}{Z-\omega_{\textrm{cor}}^{2}(\mathbf{q})}d\mathbf{q}.$
(39)

###  IV.1 Relaxation of a plane wave

Since the correlation matrix $\hat{C}$ is regular, for a given wavevector
$\mathbf{q}$ we obtain

$\langle\mathbf{q}|\hat{G}_{\textrm{cor}}(Z)|\mathbf{q}\rangle=\frac{1}{Z-\omega_{\textrm{cor}}^{2}(\mathbf{q})}.$
(40)

Using the duality relation (35) and the conformal map $Z(z)$ defined by Eq.
(36), we can present the Green function $\mathcal{G}_{\mathbf{q}}(z)$ as

$\mathcal{G}_{\mathbf{q}}(z)=\frac{1}{z-\big{(}\varkappa+1+\mathcal{M}_{\textrm{cor}}(Z(z))\big{)}\omega_{\textrm{cor}}^{2}(\mathbf{q})}.$
(41)

For large $z$, we have $Z(z)=(1+\varkappa)z+O(1)$. Therefore, from Eq. (39),
we obtain the asymptotics

$\mathcal{M}_{\textrm{cor}}(Z(z))=\frac{(1+\varkappa)\mathcal{M}_{\rm
cor}^{(1)}}{z}+O\left(\frac{1}{z^{2}}\right).$ (42)

Therefore, the Green function $\mathcal{G}_{\mathbf{q}}(z)$ can be presented
in the form of Eq. (13) with

$\displaystyle k^{\rm inst}_{\mathbf{q}}=k_{\rm inst}=1/\mathcal{M}_{\rm
cor}^{(1)},$ (43) $\displaystyle
m_{\mathbf{q}}=\big{[}(\varkappa+1)\mathcal{M}_{\rm cor}^{(1)}\omega_{\rm
cor}^{2}(\mathbf{q})\big{]}^{-1},$ (44)
$\displaystyle\mathcal{G}_{1\mathbf{q}}(z)=\mathcal{G}_{1}(z)=\frac{\mathcal{M}_{\textrm{cor}}(Z(z))}{(1+\varkappa)\mathcal{M}_{\rm
cor}^{(1)}}.$ (45)

We omit the subscript $\mathbf{q}$ for the values, which do not depend on the
wavevector $\mathbf{q}$. In the framework of the RMT, only the first mass
$m_{\mathbf{q}}$ depends on the wavevector $\mathbf{q}$ in the one-dimensional
chain. All other parameters of the one-dimensional chain do not depend on
$\mathbf{q}$.

The Green function $\mathcal{G}_{1}(z)$ corresponds to the LVDOS

$\mathcal{F}_{1}(\omega)=\frac{\omega^{2}g(\omega)}{(1+\varkappa)\mathcal{M}_{\rm
cor}^{(1)}},$ (46)

where $g(\omega)=\frac{2\omega}{\pi
N}\operatorname{Im}\operatorname{Tr}\hat{G}\big{(}(\omega-i0)^{2}\big{)}$ is
the full vibrational density of states for a given $\varkappa$. The
denominator in Eq. (46) ensures the normalization of
$\mathcal{F}_{1}(\omega)$.

Figure 3: (Color online) The LVDOS $\mathcal{F}_{n}(\omega)$ shown as a ratio
$\mathcal{F}_{n}(\omega)/\mathcal{F}_{\infty}(\omega)$ for $\varkappa=1$ for
different values of $n$ (up to $n=100$).

Figure 4: (Color online) Parameters of the chain as a function of scaled site
number $n\omega_{c}/\omega_{\rm max}$ obtained using the RMT for different
values of $\varkappa$. For small values of $\varkappa$, dots merge into lines.
In panel (a), the dashed line marks the Ioffe-Regel criterion
$k_{n}/m_{n}=\omega_{c}^{2}$.

We can use the recurrence procedure (23)–(25) to obtain the LVDOS
$\mathcal{F}_{n}(\omega)$ and the parameters $a_{n}$, $b_{n}$, $m_{n}$ (all of
them do not depend on the wavevector $\mathbf{q}$). For $n\to\infty$, we
observe that $\mathcal{F}_{n}(\omega)$ converges to a stationary solution of
Eq. (25):

$\mathcal{F}_{\infty}(\omega)=\frac{16\omega^{2}}{\pi\omega_{\rm
max}^{3}}\sqrt{1-\frac{\omega^{2}}{\omega^{2}_{\rm max}}}.$ (47)

where $\omega_{\rm max}$ is the maximum frequency in the system. In this case

$\displaystyle b_{\infty}=\frac{a_{\infty}}{2}=\frac{4}{\omega_{\rm
max}^{2}},\quad m_{\infty}=\frac{16}{\omega_{\rm max}^{4}},\quad
k_{\infty}=0.$ (48)

For example, we consider the correlation matrix $\hat{C}$ as a regular matrix
on a simple cubic lattice with a unit lattice constant. The non-diagonal
elements are $C_{ij}=-1$ if atoms $i$ and $j$ are neighbors and $C_{ij}=0$
otherwise. The diagonal elements are $C_{ii}=6$. The corresponding dispersion
law is

$\omega_{\rm cor}^{2}({\bf
q})=4\Big{(}\sin^{2}\frac{q_{x}}{2}+\sin^{2}\frac{q_{y}}{2}+\sin^{2}\frac{q_{z}}{2}\Big{)}.$
(49)

In this case, the moment-generating function is $\mathcal{M}_{\rm
cor}(Z)=\frac{1}{2}W_{s}\big{(}\frac{Z}{2}-3\big{)}$, where $W_{s}$ is the
third Watson integral Zucker (2011). Using recurrence relation (25), we obtain
the LVDOS $\mathcal{F}_{n}(\omega)$ and the values of $a_{n}$, $b_{n}$,
$m_{n}$.

Figure 3 shows $\mathcal{F}_{n}(\omega)/\mathcal{F}_{\infty}(\omega)$ for
different values of $n$ for $\varkappa=1$. One can observe that
$\mathcal{F}_{n}(\omega)$ is a smooth function, which gradually approaches
$\mathcal{F}_{\infty}(\omega)$ for $n\to\infty$.

The chain parameters for different values of $\varkappa$ are presented in Fig.
4 as $k_{n}/(m_{n}\omega_{c}^{2})$, $m_{n}/m_{\infty}$, and
$b_{n}/b_{\infty}$. The characteristic frequency $\omega_{c}\sim\varkappa$ is
used for the scaling and will be discussed in the next section. One can
observe that $k_{n}/m_{n}\sim\omega_{c}^{2}$ for $n\lesssim\omega_{\rm
max}/\omega_{c}$ and $k_{n}/m_{n}\to 0$ for $n\to\infty$. For the better
understanding of the behaviour of the obtained mass-spring chain, we consider
the low-frequency approximation, which can be analyzed analytically.

###  IV.2 Low-frequency approximation

In order to study the low-frequency dynamics, we investigate the behavior of
the Green function $\mathcal{G}_{\mathbf{q}}(z)$ for small values of $z$,
which corresponds to small values of $Z$. In this case, the generating
function $\mathcal{M}_{\textrm{cor}}(Z)$ can be approximated as a linear
function of $Z$ for small $Z$:

$\mathcal{M}_{\textrm{cor}}(Z)=-1-\beta Z+o(Z),$ (50)

where

$\beta=\frac{1}{V_{\textsc{bz}}}\int_{\textsc{bz}}\frac{d\mathbf{q}}{\omega_{\textrm{cor}}^{2}(\mathbf{q})}.$
(51)

For example, for a simple cubic lattice with the dispersion (49), the constant
is $\beta=w_{s}/2$ where $w_{s}=0.505462$ is the third Watson constant Zucker
(2011). For large $Z$, any generating function of the form (37) have to decay
as $\mathcal{M}_{\textrm{cor}}(Z)\sim 1/Z$. Therefore, to study the low-
frequency behavior, we can consider the regularized form of Eq. (50):

$\mathcal{M}_{\textrm{cor}}(Z)=\frac{1}{\beta Z-1},$ (52)

In this case, Eqs. (43)–(45) becomes

$\displaystyle k_{\rm inst}=\beta,$ (53) $\displaystyle
m_{\mathbf{q}}=\frac{\beta}{(1+\varkappa)\omega_{\rm cor}^{2}(\mathbf{q})},$
(54)
$\displaystyle\mathcal{G}_{1}(z)=\frac{\beta}{1+\varkappa}\mathcal{M}_{\textrm{cor}}(Z(z)).$
(55)

Using Eq. (36), we can present $\mathcal{M}_{\textrm{cor}}(Z(z))$ in the
recurrence form:

$\mathcal{M}_{\textrm{cor}}(Z(z))=-\frac{\varkappa+1}{\varkappa+2-\beta
z+\mathcal{M}_{\textrm{cor}}(Z(z))}.$ (56)

Thus, we automatically obtain the continued fraction:

$\mathcal{G}(z)=-\frac{m_{\mathbf{q}}}{k_{\rm
inst}-m_{\mathbf{q}}z-\frac{b^{2}}{a-mz-\frac{b^{2}}{a-mz-\ldots}}},$ (57)

where

$\displaystyle a=\beta\frac{2+\varkappa}{1+\varkappa},\quad
b=\frac{\beta}{\sqrt{1+\varkappa}},\quad m=\frac{\beta^{2}}{1+\varkappa}.$
(58)

The corresponding mass-spring chain is homogeneous starting from $n=1$, which
is shown in Fig. (5). Each mass in the semi-infinite chain (except the first
one) is connected to the ground by a spring with the stiffness $k=a-2b$. It
corresponds to the natural frequency

$\omega_{c}=\sqrt{\frac{k}{m}}=\frac{\sqrt{1+\varkappa}-1}{\sqrt{\beta}},$
(59)

which is the minimal frequency of vibrations, which can propagate along the
chain. For small $\varkappa$, we have
$\omega_{c}\approx\varkappa/(2\sqrt{\beta})$. This frequency is used for the
scaling in Fig. 4. The dashed line in Fig. 4a corresponds to the homogeneous
chain obtained in the low-frequency approximation.

Figure 5: (Color online) A regular one-dimensional mass-spring model, which
corresponds to the low-frequency approximation.

Above the frequency $\omega_{c}$, the chain effectively absorbs the
vibrational energy of the first site. Below $\omega_{c}$, the initial
excitation distributes over several sites in the beginning of the chain
without any further damping. Indeed, the frequency $\omega_{c}$ coincides with
the Ioffe-Regel frequency in the random matrix approach Conyuh and Beltukov
(2020). In the low-frequency approximation (52) used in this section, no
damping of vibrations below the Ioffe-Regel frequency can be observed. In
other words, all frequencies $\omega<\omega_{c}$ are localized in the
beginning of the one-dimensional chain. However, a precise evaluation of the
chain parameters shows that $k_{n}$ gradually decreases to zero for large $n$
(see Fig. 4). It corresponds to a finite relaxation time for frequencies
$\omega<\omega_{c}$, which will be discussed in Section VI.

From Eq. (56), we can find $\mathcal{M}_{\textrm{cor}}(Z(z))$ explicitly:

$\mathcal{M}_{\textrm{cor}}(Z(z))=\frac{1}{2}\big{(}mz-2b-k\\\
-\sqrt{mz-k}\sqrt{mz-4b-k}\Big{)}.$ (60)

Using Eqs. (45) and (15), we can find the corresponding memory function $K(t)$
for $k\ll b$ and $\varkappa\ll 1$ (see Appendix C). In this case the equation
of motion (14) becomes

$m_{\mathbf{q}}\ddot{u}_{\mathbf{q}}(t)+\gamma_{s}\dot{u}_{\mathbf{q}}(t)+k_{s}u_{\mathbf{q}}(t)\\\
+\int_{-\infty}^{t}K_{l}(t-t^{\prime})u_{\mathbf{q}}(t^{\prime})dt^{\prime}=0$
(61)

with the short-term damping $\gamma_{s}=\sqrt{mb}$, the short-term stiffness
$k_{s}=k_{\rm inst}-b-k/2=\sqrt{kb}$ and the long-term Bessel memory function

$K_{l}(t)=\frac{\sqrt{bk}}{t}J_{1}\left(t\sqrt{\frac{k}{m}}\right),\quad t>0.$
(62)

Equation (61) defines the natural frequency of the first site in the chain

$\omega_{0}(\mathbf{q})=\sqrt{\frac{k_{s}}{m_{\mathbf{q}}}}=\sqrt{\frac{\varkappa}{2}}\omega_{\rm
cor}(\mathbf{q}).$ (63)

If this frequency is above the Ioffe-Regel frequency $\omega_{c}$, a strong
damping is observed. In terms of the wavevector, it corresponds to $q>q_{c}$,
where the Ioffe-Regel wavenumber is

$q_{c}=\frac{1}{v_{\rm cor}}\sqrt{\frac{\varkappa}{2\beta}}.$ (64)

Here we use $\omega_{\rm cor}(\mathbf{q})=v_{\rm cor}q$ for small $q$. For
$q<q_{c}$, one can find the resonance frequency in Eq. (61)

$\omega(\mathbf{q})=v_{\rm cor}^{2}q\sqrt{\beta(2q_{c}^{2}-q^{2})}.$ (65)

In this case, the memory function $K_{l}(t)$ cancels out the damping term
$\gamma_{s}\dot{u}_{\mathbf{q}}(t)$ and modifies the resonance frequency:
$\omega(\mathbf{q})>\omega_{0}(\mathbf{q})$. For $q\ll q_{c}$, we have
$\omega(\mathbf{q})=\sqrt{2}\omega_{0}(\mathbf{q})$.

Figure 6: (Color online) Relaxation of plane waves in the framework of the RMT
for $\varkappa=0.01$ (a) and $\varkappa=0$ (b) for the same set of wavenumbers
$q$. Solid lines show the precise result, which corresponds to the general
chain (Fig. 2). Dashed lines represent a relaxation in the low-frequency
approximation, which corresponds to the regular chain (Fig. 5).

Figure 6 shows the relaxation of $u_{\mathbf{q}}(t)$ for a simple cubic
lattice with the dispersion $\omega_{\rm cor}(\mathbf{q})$ defined by Eq.
(49). For $\varkappa=0.01$, the Ioffe-Regel wavenumber is $q_{c}\approx 0.14$.
Figure 6a shows the relaxation for two wavenumbers below $q_{c}$ and two
wavenumbers above $q_{c}$. Solid lines show the exact solution while dashed
lines show the low-frequency approximation (61). One can see a good agreement
between them. However, for large time $t$ and $q<q_{c}$, one can observe the
slow relaxation of the exact solution while the low-frequency approximation
has stationary oscillations.

For $\varkappa=0$, the short-term stiffness $k_{s}$ and the long-term memory
function $K_{l}(t)$ vanish, and we obtain viscous damping without returning
force

$m_{\mathbf{q}}\ddot{u}_{\mathbf{q}}(t)+\gamma_{s}\dot{u}_{\mathbf{q}}(t)=0.$
(66)

Figure 6b shows this behavior both for the exact solution and the low-
frequency approximation.

##  V A numerical analysis of dynamical matrices

The same approach can be used to obtain the chain representation for a given
dynamical matrix or a given ensemble of dynamical matrices. These dynamical
matrices can be obtained by using molecular dynamics simulations or using
numerical random matrix models.

In the previous section, we consider the thermodynamic limit $N\to\infty$ in
the framework of the RMT. In this case, there are no fluctuations of
vibrational properties. However, in a finite system, fluctuations may be
important, especially for the stability of the parameters of the mass-spring
chain. In this section, we demonstrate that the proposed recurrence algorithm
may be used for a finite numerical system as well.

We use the numerical random matrix model in the form of the correlated Wishart
ensemble $\hat{M}=\hat{A}\hat{A}^{T}$ which can be controlled by the same
parameter $\varkappa$. For simplicity, we consider a simple cubic lattice with
random bonds and unit lattice constant as in the RMT. However, in contrast to
the RMT, the numerical random matrix model has the finite interaction radius.

For $\varkappa=0$ the matrix $\hat{A}$ is square and the number of bonds $K$
is equal to the number of degrees of freedom $N$. We can consider the
following structure of the non-diagonal elements of the matrix $\hat{A}$
Beltukov  _et al._ (2013); Conyuh and Beltukov (2020):

$A_{ij}=\left\\{\begin{array}[]{ll}\frac{1}{2}\xi_{ij}&\text{if $i$ and $j$
are neighbors},\\\ 0&\text{otherwise},\end{array}\right.$ (67)

where $\xi_{ij}$ are independent Gaussian random numbers with zero mean and
unit variance. The diagonal elements are defined using the sum rule
$A_{ii}=-\sum_{j\neq i}A_{ji}$. This procedure results in the same correlation
matrix $\hat{C}$ that was used in the framework of the RMT in the previous
Section.

For $\varkappa>0$ we can use two realizations of square random matrices
defined by Eq. (67): $\hat{A}^{(0)}$ and $\hat{A}^{(1)}$. The resulting
rectangular matrix $\hat{A}$ can be obtained by inserting $\varkappa N$
randomly chosen columns of the matrix $\hat{A}^{(1)}$ into the matrix
$\hat{A}^{(0)}$. This random insertion of the new columns corresponds to the
random addition of new bonds to the vibrational system.

We use the Kernel Polynomial Method (KPM) Beltukov  _et al._ (2016) to obtain
the initial function $\mathcal{F}_{\mathbf{q}}(\omega)$ for the recurrence
relation (25). The initial function $\mathcal{F}_{\mathbf{q}}(\omega)$ was
calculated for a system with $N=100^{3}$ atoms and different values of the
parameter $\varkappa$ and the wavenumber $q$. The averaging over $10^{2}$ –
$10^{3}$ realizations was applied. An example for $\varkappa=1$ and $q=1$ is
shown in Fig. 7.

Figure 7: (Color online) The function $\mathcal{F}_{\mathbf{q}}(\omega)$
calculated using the KPM for a numerical random matrix with $N=100^{3}$ atoms,
the parameter $\varkappa=1$, and the wave number $q=1$. Averaging over 100
realizations was applied. Vertical dashed line shows the chosen position of
$\omega_{\rm max}$.

We use the Chebyshev expansion of $\mathcal{F}_{\mathbf{q}}(\omega)$ in the
frequency range $0\leq\omega\leq\omega_{\rm max}$ to evaluate the recurrence
relation (25) (see Appendix D). The choice of $\omega_{\rm max}$ is important
because there is an exponential tail in the high-frequency vibrational density
of states without any specific maximum frequency. For any finite system size,
there are big relative fluctuations in the high-frequency tail due to a small
number of vibrations there (see Fig. 7). It may lead to additional
fluctuations of the parameters $a_{n\mathbf{q}}$, $b_{n\mathbf{q}}$,
$m_{n\mathbf{q}}$ of the obtained chain. Thus, we leave a small number of
vibrational modes above $\omega_{\rm max}$ and do not use them in the
Chebyshev expansion, which significantly reduces the fluctuations of the
obtained parameters.

The KPM is also based on the Chebyshev expansion Weiße  _et al._ (2006).
However, the maximum frequency in the KPM, $\omega_{\rm max}^{\textsc{kpm}}$,
should be larger than any frequency in the system for stability purposes.
Therefore, we remap the obtained Chebyshev expansion to another one with
slightly smaller maximum frequency $\omega_{\rm max}$ to drop a small number
of high-frequency modes.

Figure 8: (Color online) Parameters of the chain as a function of scaled site
number $\varkappa n/\omega_{\rm max}$ obtained using the numerical random
matrix model (67) for different values of the parameter $\varkappa$ and the
wavenumber $q$.

The resulting parameters $k_{\mathbf{q}n}$, $b_{\mathbf{q}n}$,
$m_{\mathbf{q}n}$ of the mass-spring chain are presented in Fig. 8 for
different values of the parameter $\varkappa$ and the wavenumber $q$. The
results are similar to those obtained in the framework of the RMT (Fig. 4). In
the given scale, all calculated values $k_{n}/(\varkappa^{2}m_{n})$ almost
coincide in Fig. 8a. The masses $m_{\mathbf{q}n}$ and stiffnesses
$b_{n\mathbf{q}}$ are close to their stationary values $m_{\infty}$ and
$b_{\infty}$ respectively. For different wavenumbers $q$, we observe similar
parameters of the mass-spring chain.

##  VI Discussion

We have shown that the relaxation of a plane wave with the wavevector
$\mathbf{q}$ coincides with the relaxation of the first site in the semi-
infinite one-dimensional mass-spring chain (Fig. 2). The parameters of the
chain form the continued fraction representation of the Green function
$\mathcal{G}_{\mathbf{q}}(z)$.

Each site in the chain is connected to the ground by a spring with a stiffness
$k_{n\mathbf{q}}$. We observe the same behavior of $k_{n\mathbf{q}}$ in both
RMT and numerical analysis of finite random matrices: in the beginning of the
chain $k_{n\mathbf{q}}\sim k>0$ (denoted as a barrier in Fig. 4a) and then it
gradually goes to zero.

For small values of $\varkappa$, we observe decaying oscillations of
$k_{n\mathbf{q}}$. It is known that oscillations of the coefficients in the
continued fraction are related to singularities in the density of states
Gaspard and Cyrot-Lackmann (1973); Hodges (1977). In our case, the vibrational
density of states has a steep behavior near the Ioffe-Regel frequency for
small $\varkappa$ Conyuh and Beltukov (2020); Conyuh  _et al._ (2019).
However, these oscillations are not important for qualitative analysis.

For large enough site number $n$, the parameters
$a_{n\mathbf{q}},b_{n\mathbf{q}},m_{n\mathbf{q}},k_{n\mathbf{q}}$ converge to
their stationary values $a_{\infty},b_{\infty},m_{\infty},k_{\infty}$, which
depend only on the maximum frequency in the system $\omega_{\rm max}$. For
electronic systems, it is known, that the coefficients of the continued
fraction have stationary values, which depend on the width of the energy band
Gaspard and Cyrot-Lackmann (1973). Since $k_{\infty}=0$, the tail of the mass-
spring chain is free and homogeneous. This tail can be considered as a simple
thermal bath, which finally absorbs the initial vibrational energy Adelman and
Garrison (1976); Li and Weinan (2007). At the same time, no explicit damping
is introduced in the mass-spring chain, which corresponds to the absence of
damping in the initial equation of motion (1).

In the studied models, all masses $m_{n\mathbf{q}}$ are close to the
stationary value $m_{\infty}$ (except the first one). Thus, the chain can
easily absorb the vibrational energy above the Ioffe-Regel frequency
$\omega_{c}=\sqrt{k/m_{\infty}}$. For frequencies below $\omega_{c}$, the
absorption is much smaller because these frequencies are “forbidden” in the
beginning of the chain. Depending on the form and the width of the barrier,
there is a relatively small absorption below $\omega_{c}$. The dynamics of the
chain can be mapped to a discrete version of a one-dimensional Schrödinger
equation with the potential energy $k_{n\mathbf{q}}/m_{n\mathbf{q}}$ and the
energy $\omega^{2}$. Therefore, the region with
$k_{n\mathbf{q}}/m_{n\mathbf{q}}>\omega^{2}$ acts as the tunneling barrier (or
a high-pass filter).

The first mass in the chain strongly depends on the wavevector $\mathbf{q}$:
$m_{\mathbf{q}}\sim\omega_{\rm cor}^{-2}(\mathbf{q})$. It results in a natural
frequency $\omega_{0}(\mathbf{q})=\sqrt{\varkappa/2}\omega_{\rm
cor}(\mathbf{q})$ of the first site. If this frequency is smaller than the
Ioffe-Regel frequency $\omega_{c}$ for a given wavevector $\mathbf{q}$, then
the damping of this vibrational mode is relatively slow. It corresponds to the
notion of phonons with well-defined dispersion $\omega(\mathbf{q})$. If
$\omega_{0}(\mathbf{q})>\omega_{c}$, a strong damping is observed, which
corresponds to the notion of diffusons above the Ioffe-Regel frequency.

One can note that $\omega_{\textrm{cor}}^{2}(\mathbf{q})$ is proportional to
the Laplacian on the corresponding lattice. Therefore, in the framework of the
RMT, the equation of motion (14) can be rewritten in the real space as

$\rho\ddot{u}(\mathbf{r},t)+k_{\rm inst}\Delta
u(\mathbf{r},t)+\int_{-\infty}^{t}K(t-t^{\prime})\Delta
u(\mathbf{r},t^{\prime})dt^{\prime}=0.$ (68)

This viscoelastic equation is not local in time, but local in space. In a
general case, the instantaneous stiffness $k^{\rm inst}_{\mathbf{q}}$ and the
memory function $K_{\mathbf{q}}(t)$ may depend on the wavevector $\mathbf{q}$
which results in additional spatial convolution in Eq. (68). In the low-
frequency approximation, from Eq. (61) we obtain

$\rho\ddot{u}(\mathbf{r},t)+\gamma_{s}\Delta\dot{u}(\mathbf{r},t)+k_{s}\Delta
u(\mathbf{r},t)\\\ +\int_{-\infty}^{t}K_{l}(t-t^{\prime})\Delta
u(\mathbf{r},t^{\prime})dt^{\prime}=0$ (69)

with long-term Bessel memory function $K_{l}(t)$ defined by Eq. (62). For
$\varkappa=0$, from Eq. (66) we obtain a viscous equation without any
returning force

$\rho\ddot{u}(\mathbf{r},t)+\gamma_{s}\Delta\dot{u}(\mathbf{r},t)=0.$ (70)

The case $\varkappa=0$ is known as the isostatic state in the jamming
transition Wyart  _et al._ (2005). In this case $\omega_{c}=0$ and the entire
low-frequency range is occupied by diffusons. Equation (70) is consistent with
a model of random walks of atomic displacements for isostatic case Beltukov
_et al._ (2013).

The above equation was obtained in the scalar model. In general case,
$\mathbf{u}(\mathbf{r},t)$ is a vector and all constants and memory functions
in the above equations become tensors.

The continued fraction presentation (22) is well-known in the classical theory
of moments Akhiezer (1965). However, the direct evaluation of moments leads to
stability and performance issues for $n\gtrsim 30$. The proposed method is
based on the Chebyshev expansion (Appendix D) and shows the numerical
stability both for the RMT (Section IV) and numerical random matrices (Section
V). It takes about one hour on a modern computer to find up to $10^{5}$
coefficients in the continued fraction from the known function
$\mathcal{F}_{\mathbf{q}}(\omega)$. Usually, the most time-consuming part is
the calculation of $\mathcal{F}_{\mathbf{q}}(\omega)$ and its averaging over
different realizations of the dynamical matrix $\hat{M}$.

For a given dynamical matrix $\hat{M}$, the continued fraction (22) can be
also obtained using the Lanczos method using $|\mathbf{q}\rangle$ as a
starting vector Dagotto (1994). However, there are two important issues
concerning the Lanczos method. It is known that the Lanczos method is unstable
due to the loss of orthogonality. To stabilize the algorithm, an additional
reorthogonalization is required, which decreases the performance Saad (1992).
The second issue is that the Lanczos method is highly sensitive to a small
number of high-frequency eigenmodes. There is no simple way to discard these
eigenvalues or average the resulting parameters
$a_{n\mathbf{q}},b_{n\mathbf{q}},m_{n\mathbf{q}}$ over different realizations
of the dynamical matrix $\hat{M}$ in the framework of the Lanczos method. In
the proposed algorithm, the averaging is performed directly on
$\mathcal{F}_{\mathbf{q}}(\omega)$ and a small number of high-frequency modes
can be easily removed as was discussed in Section V.

##  VII Conclusion

We have shown that viscoelastic relaxation of plane waves in amorphous solids
can be considered as dynamics of the one-dimensional semi-infinite chain. The
first site in this chain represents the initial plane wave while the rest of
the chain represents the memory effects. The initial vibrational energy
gradually spreads along the chain, which results in the vibrational relaxation
of the initial plane wave.

In the beginning of the chain, there is a natural barrier for frequencies
$\omega<\omega_{c}$, which corresponds to the Ioffe-Regel crossover. In the
framework of the RMT, the memory function does not depend on the wavevector
$\mathbf{q}$, which results in the viscoelastic equation, which is local in
space but not local in time. In the low-frequency approximation, the long-term
memory function can be described by the Bessel function.

The proposed method demonstrates numerical stability for studied theoretical
and numerical random matrices.

##  VIII Acknowledgments

We wish to acknowledge D. A. Parshin and A. V. Shumilin for valuable
discussions. The authors thank the Russian Foundation for Basic Research
(project no. 19-02-00184) for the financial support.

##  Appendix A Evolution of plane waves

For initial conditions $|u(0)\rangle=0$ and
$|\dot{u}(0)\rangle=|\mathbf{q}\rangle$, the solution of (1) is

$|u(t)\rangle=v_{0}\sum_{n}|n\rangle\frac{\sin(\omega_{n}t)}{\omega_{n}}\langle
n|\mathbf{q}\rangle,$ (71)

where $|n\rangle$ is $n$-th eigenvector of the matrix $\hat{M}$ and
$\omega_{n}$ is the corresponding eigenfrequency. Therefore, the projection to
the plane wave is

$u_{\mathbf{q}}(t)=\left\langle
v_{0}\sum_{n}\langle\mathbf{q}|n\rangle\frac{\sin(\omega_{n}t)}{\omega_{n}}\langle
n|\mathbf{q}\rangle\right\rangle.$ (72)

Using the Fourier transform, we can write

$u_{\mathbf{q}}(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\tilde{u}_{\mathbf{q}}(\omega)e^{i\omega
t}d\omega,$ (73)

where

$\tilde{u}_{\mathbf{q}}(\omega)=\lim_{\epsilon\to
0^{+}}\int_{0}^{\infty}u_{\mathbf{q}}(t)e^{-i\omega t-\epsilon t}dt\\\
=\left\langle v_{0}\lim_{\epsilon\to
0^{+}}\sum_{n}\langle\mathbf{q}|n\rangle\langle
n|\mathbf{q}\rangle\int_{0}^{\infty}\frac{\sin(\omega_{n}t)}{\omega_{n}}e^{-i\omega
t-\epsilon t}dt\right\rangle\\\ =\left\langle v_{0}\lim_{\epsilon\to
0^{+}}\sum_{n}\langle\mathbf{q}|n\rangle\frac{1}{\omega_{n}^{2}-(\omega-i\epsilon)^{2}}\langle
n|\mathbf{q}\rangle\right\rangle\\\
=-v_{0}\big{\langle}\mathbf{q}\big{|}\hat{G}\big{(}(\omega-i0)^{2}\big{)}\big{|}\mathbf{q}\big{\rangle}.$
(74)

It is worth to note that $u_{\mathbf{q}}(t)$ defined by Eq. (73) is zero for
$t<0$.

##  Appendix B Viscoelastic relations

In this Appendix we provide the most important relation between different
viscoelastic functions. Relations between the Green function
$\mathcal{G}_{n\mathbf{q}}(z)$, the LVDOS $\mathcal{F}_{n\mathbf{q}}(\omega)$,
and the moments $\mathcal{F}_{n\mathbf{q}}^{(k)}$ are:

$\displaystyle\mathcal{F}_{n\mathbf{q}}(\omega)$
$\displaystyle=\frac{2\omega}{\pi}\operatorname{Im}\mathcal{G}_{n\mathbf{q}}\big{(}(\omega-i0)^{2}\big{)},$
(75) $\displaystyle\mathcal{G}_{n\mathbf{q}}(z)$
$\displaystyle=\int_{0}^{\infty}\frac{{\cal
F}_{n\mathbf{q}}(\omega)}{z-\omega^{2}}d\omega=\sum_{k=0}^{\infty}\frac{{\cal
F}_{n\mathbf{q}}^{(k)}}{z^{k+1}},$ (76)
$\displaystyle\mathcal{F}_{n\mathbf{q}}^{(k)}$
$\displaystyle=\int_{0}^{\infty}\omega^{2k}\mathcal{F}_{n\mathbf{q}}(\omega)d\omega.$
(77) Relations between the Green function $\mathcal{G}_{n\mathbf{q}}(z)$ and
the memory function $K_{n\mathbf{q}}(t)$ are: $\displaystyle
K_{n\mathbf{q}}(t)$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathcal{G}_{n\mathbf{q}}\big{(}(\omega-i0)^{2}\big{)}e^{i\omega
t}d\omega,$ (78) $\displaystyle\mathcal{G}_{n\mathbf{q}}(z)$
$\displaystyle=\int_{0}^{\infty}K_{n\mathbf{q}}(t)e^{-t\sqrt{-z}}d\omega.$
(79) Relations between the LVDOS $\mathcal{F}_{n\mathbf{q}}(\omega)$ and the
memory function $K_{n\mathbf{q}}(t)$ are: $\displaystyle K_{n\mathbf{q}}(t)$
$\displaystyle=-\theta(t)\int_{0}^{\infty}\frac{\mathcal{F}_{n\mathbf{q}}(\omega)}{\omega}\sin(\omega
t)d\omega,$ (80) $\displaystyle\mathcal{F}_{n\mathbf{q}}(\omega)$
$\displaystyle=-\frac{2\omega}{\pi}\int_{0}^{\infty}K_{n\mathbf{q}}(t)\sin(\omega
t)dt,$ (81)

where $\theta(t)$ is the Heaviside step function.

##  Appendix C Low-frequency memory function

In the case $k\ll b$, we can write the generating function (60) in the form

$\mathcal{M}_{0}(Z)=\mathcal{M}_{s}(Z)+\mathcal{M}_{p}(Z)+\mathcal{M}_{l}(Z)+o\left(\frac{k}{b}\right),$
(82)

where the main term and two kinds of perturbations have the following form

$\displaystyle\mathcal{M}_{s}(Z)$
$\displaystyle=\frac{1}{2}\left(mz-2b-\sqrt{mz}\sqrt{mz-4b}\right),$ (83)
$\displaystyle\mathcal{M}_{p}(Z)$
$\displaystyle=\frac{k}{2}\left(\frac{\sqrt{-bmz}}{mz}+\frac{mz-2b}{\sqrt{mz}\sqrt{mz-4b}}\right),$
(84) $\displaystyle\mathcal{M}_{l}(Z)$
$\displaystyle=\sqrt{b(k-mz)}-\sqrt{-bmz}.$ (85)

Using Eq. (15) with $z=(\omega-i0)^{2}$, we obtain the corresponding memory
functions

$\displaystyle K_{s}(t)$
$\displaystyle=-\frac{2b}{t}J_{2}\big{(}\tilde{t}\big{)}\theta(t),$ (86)
$\displaystyle K_{p}(t)$
$\displaystyle=\frac{k}{2}\sqrt{\frac{b}{m}}\Bigg{(}J_{0}\big{(}\tilde{t}\big{)}\left[\tilde{t}-\frac{\pi}{2}\tilde{t}H_{1}\big{(}\tilde{t}\big{)}\right]$
$\displaystyle\qquad-
J_{1}\big{(}\tilde{t}\big{)}\left[2-\frac{\pi}{2}\tilde{t}H_{0}\big{(}\tilde{t}\big{)}\right]-1\Bigg{)}\theta(t),$
(87) $\displaystyle K_{l}(t)$
$\displaystyle=\frac{\sqrt{bk}}{t}J_{1}\left(t\sqrt{\frac{k}{m}}\right)\theta(t),$
(88)

where $\tilde{t}=2t\sqrt{b/m}$ is the scaled time, $J_{n}$ is the Bessel
function, $H_{n}$ is the Struve function, and $\theta(t)$ is the Heaviside
step function.

The memory function $K_{s}(t)$ is the main short-term memory function, which
is not zero for $k=0$. It coincides with the memory function for the semi-
infinite free mass-spring chain Li and Weinan (2007); Kwidzinski and Bulla
(2016). The memory function $K_{p}(t)$ is a perturbation of the short-term
memory function for nonzero $k$. The memory function $K_{l}(t)$ is a long-term
memory perturbation since it depends on another scaled time $t\sqrt{k/m}$,
which scales with $k$.

In the low-frequency approximation, the memory functions $K_{s}(t)$ and
$K_{p}(t)$ can be considered as an instantaneous response:

$\int_{-\infty}^{t}K_{s}(t-t^{\prime})u_{\mathbf{q}}(t^{\prime})dt^{\prime}\\\
\approx\int_{-\infty}^{t}K_{s}(t-t^{\prime})\big{[}u_{\mathbf{q}}(t)+(t^{\prime}-t)\dot{u}_{\mathbf{q}}(t))\big{]}dt^{\prime}\\\
=-bu_{\mathbf{q}}(t)+\sqrt{mb}\dot{u}_{\mathbf{q}}(t),$ (89)
$\int_{-\infty}^{t}K_{p}(t-t^{\prime})u_{\mathbf{q}}(t^{\prime})dt^{\prime}\\\
\approx\int_{-\infty}^{t}K_{p}(t-t^{\prime}))\big{[}u_{\mathbf{q}}(t)+(t^{\prime}-t)\dot{u}_{\mathbf{q}}(t))\big{]}dt^{\prime}=-\frac{k}{2}u_{\mathbf{q}}(t).$
(90)

##  Appendix D Chebyshev expansion

In order to calculate the next LVDOS $\mathcal{F}_{n+1}(\omega)$ using the
recurrence relation (25), one should calculate
${\mathcal{G}_{n}\big{(}(\omega-i0)^{2}\big{)}}$ from the known LVDOS
$\mathcal{F}_{n}(\omega)$. It can be done using the Chebyshev expansion of the
form

$\mathcal{F}_{n}(\omega)=\frac{4\omega^{2}}{\omega_{\rm
max}^{3}}\sqrt{1-\frac{\omega^{2}}{\omega_{\rm
max}^{2}}}\sum_{k}c_{n,k}U_{k}\left(\frac{2\omega^{2}}{\omega_{\rm
max}^{2}}-1\right),$ (91)

where $U_{k}$ is the Chebyshev polynomial of the second kind. The coefficients
$c_{n,k}$ can be obtained using the orthogonal relation:

$c_{n,k}=\frac{4}{\pi\omega_{\rm max}^{2}}\int_{0}^{\omega_{\rm
max}}\\!\mathcal{F}_{n}(\omega)U_{k}\left(\frac{2\omega^{2}}{\omega_{\rm
max}^{2}}-1\right)d\omega.$ (92)

First several moments of $\mathcal{F}_{n}(\omega)$ are

$\displaystyle\mathcal{F}_{n}^{(0)}$ $\displaystyle=\pi
c_{n,0}\frac{\omega_{\rm max}^{2}}{4},$ (93)
$\displaystyle\mathcal{F}_{n}^{(1)}$
$\displaystyle=\pi\left(2c_{n,0}+c_{n,1}\right)\frac{\omega_{\rm
max}^{4}}{16},$ (94) $\displaystyle\mathcal{F}_{n}^{(2)}$
$\displaystyle=\pi\left(5c_{n,0}+4c_{n,1}+c_{n,2}\right)\frac{\omega_{\rm
max}^{6}}{64}.$ (95)

The normalization condition implies $\mathcal{F}_{n}^{(0)}=1$, which
corresponds to $c_{n,0}=4/(\pi\omega_{\rm max}^{2})$. Using Chebyshev
expansion (91), we can evaluate the Green function:

$\mathcal{G}_{n}\big{(}(\omega-i0)^{2}\big{)}=\int_{0}^{\infty}\frac{{\cal
F}_{n}(\omega_{1})}{(\omega-i0)^{2}-\omega_{1}^{2}}d\omega_{1}\\\
=\pi\sum_{k}c_{n,k}T_{k+1}\left(\frac{2\omega^{2}}{\omega_{\rm
max}^{2}}-1\right)+\frac{i\pi}{2\omega}\mathcal{F}_{n}(\omega),$ (96)

where $T_{k}(x)=(U_{k}(x)-U_{k-2}(x))/2$ is the Chebyshev polynomial of the
first kind.

For the stationary LVDOS $\mathcal{F}_{\infty}(\omega)$ defined in Eq. (47),
all coefficients $c_{n,k}$ except $c_{n,0}$ are zero. In this case, the Green
function has constant absolute value
$\big{|}\mathcal{G}_{\infty}{\big{(}(\omega-i0)^{2}\big{)}}\big{|}=4/\omega_{\rm
max}^{2}$ for $0\leq\omega\leq\omega_{\rm max}$. Thus, for
$\mathcal{F}_{n}(\omega)=\mathcal{F}_{\infty}(\omega)$ and $m_{n}=m_{\infty}$,
the recurrence relation (25) gives the same LVDOS
$\mathcal{F}_{n+1}(\omega)=\mathcal{F}_{\infty}(\omega)$.

For numerical purposes, we can calculate $\mathcal{F}_{n}(\omega)$ on a finite
number of points, which are known as Chebyshev nodes:

$\omega_{j}=\omega_{\rm max}\sin\left(\pi\frac{2j+1}{4N_{\rm pts}}\right),$
(97)

where $N_{\rm pts}$ is the number of points. In this case, the transformation
between $c_{n,k}$ and $\mathcal{F}_{n}(\omega_{j})$ can be performed using the
Fast Fourier Transform. The number $N_{\rm pts}$ corresponds to the maximum
degree of Chebyshev polynomials that we would take into account. The number
$N_{\rm pts}$ should be several times bigger than the maximum number $n$ used
in the recurrence relation.

Using the values on the Chebyshev nodes, the recurrence relation (25) can be
calculated directly

$\displaystyle\mathcal{F}^{\prime}_{n+1}(\omega_{j})$
$\displaystyle=\frac{\mathcal{F}_{n}(\omega_{j})}{\left|\mathcal{G}_{n}\big{(}(\omega_{j}-i0)^{2}\big{)}\right|^{2}},$
(98) $\displaystyle\mathcal{F}_{n+1}(\omega_{j})$
$\displaystyle=\frac{\mathcal{F}^{\prime}_{n+1}(\omega_{j})}{\mathcal{F}_{n+1}^{\prime(0)}},$
(99)

producing the following parameters:

$m_{n}=\Big{[}\mathcal{F}_{n+1}^{\prime(0)}\Big{]}^{-1},\quad
a_{n}=m_{n}\mathcal{F}_{n}^{(1)},\quad b_{n}=m_{n}^{2}.$ (100)

Then we can use again the Chebyshev expansion (91) for
$\mathcal{F}_{n+1}(\omega)$ and Eq. (96) to obtain
$\mathcal{G}_{n+1}\big{(}(\omega-i0)^{2}\big{)}$.

## References

  * Allen and Feldman (1993) P. B. Allen and J. L. Feldman, Phys. Rev. B 48, 12581 (1993).
  * Allen  _et al._ (1999) P. B. Allen, J. L. Feldman, J. Fabian, and F. Wooten, Phil. Mag. B 79, 1715 (1999).
  * Beltukov  _et al._ (2013) Y. Beltukov, V. Kozub, and D. Parshin, Physical Review B 87, 134203 (2013).
  * Anderson  _et al._ (1972) P. W. Anderson, B. Halperin, and C. M. Varma, Philosophical Magazine 25, 1 (1972).
  * Phillips (1972) W. Phillips, Journal of Low Temperature Physics 7, 351 (1972).
  * Jäckle (1972) J. Jäckle, Zeitschrift für Physik A Hadrons and nuclei 257, 212 (1972).
  * Buchenau  _et al._ (1992) U. Buchenau, Y. M. Galperin, V. Gurevich, D. Parshin, M. Ramos, and H. Schober, Physical Review B 46, 2798 (1992).
  * Ji  _et al._ (2019) W. Ji, M. Popović, T. W. de Geus, E. Lerner, and M. Wyart, Physical Review E 99, 023003 (2019).
  * Jäckle  _et al._ (1976) J. Jäckle, L. Piché, W. Arnold, and S. Hunklinger, Journal of Non-Crystalline Solids 20, 365 (1976).
  * Tielbürger  _et al._ (1992) D. Tielbürger, R. Merz, R. Ehrenfels, and S. Hunklinger, Physical review B 45, 2750 (1992).
  * Ruocco  _et al._ (1999) G. Ruocco, F. Sette, R. Di Leonardo, D. Fioretto, M. Krisch, M. Lorenzen, C. Masciovecchio, G. Monaco, F. Pignon, and T. Scopigno, Physical Review Letters 83, 5583 (1999).
  * Dell’Anna  _et al._ (1998) R. Dell’Anna, G. Ruocco, M. Sampoli, and G. Viliani, Physical review letters 80, 1236 (1998).
  * Beltukov  _et al._ (2016) Y. Beltukov, C. Fusco, D. Parshin, and A. Tanguy, Physical Review E 93, 023006 (2016).
  * Beltukov  _et al._ (2018) Y. Beltukov, D. Parshin, V. Giordano, and A. Tanguy, Physical Review E 98, 023005 (2018).
  * Gelin  _et al._ (2016) S. Gelin, H. Tanaka, and A. Lemaître, Nature materials 15, 1177 (2016).
  * Damart  _et al._ (2017) T. Damart, A. Tanguy, and D. Rodney, Physical Review B 95, 054203 (2017).
  * Baldi  _et al._ (2016) G. Baldi, V. M. Giordano, B. Ruta, and G. Monaco, Physical Review B 93, 144204 (2016).
  * Buchenau (2014) U. Buchenau, Physical Review E 90, 062319 (2014).
  * Luo  _et al._ (2020) H. Luo, A. Gravouil, V. Giordano, W. Schirmacher, and A. Tanguy, Physical Review E 102, 033003 (2020).
  * Ranganathan  _et al._ (2017) R. Ranganathan, Y. Shi, and P. Keblinski, Journal of Applied Physics 122, 145103 (2017).
  * Lemaître and Maloney (2006) A. Lemaître and C. Maloney, Journal of statistical physics 123, 415 (2006).
  * Hansen and McDonald (1990) J.-P. Hansen and I. R. McDonald,  _Theory of simple liquids_ (Elsevier, 1990).
  * Mori (1965) H. Mori, Progress of Theoretical Physics 34, 399 (1965).
  * Evers and Mirlin (2008) F. Evers and A. D. Mirlin, Reviews of Modern Physics 80, 1355 (2008).
  * Beltukov (2015) Y. M. Beltukov, JETP Letters 101, 345 (2015).
  * Conyuh and Beltukov (2020) D. A. Conyuh and Y. M. Beltukov, “A random matrix approach to the boson peak and ioffe-regel criterion in amorphous solids,” (2020), arXiv:2007.12288 .
  * Shintani and Tanaka (2008) H. Shintani and H. Tanaka, Nature materials 7, 870 (2008).
  * Skipetrov and Beltukov (2018) S. Skipetrov and Y. Beltukov, Physical Review B 98, 064206 (2018).
  * Bhatia (2009) R. Bhatia,  _Positive definite matrices_ , Vol. 24 (Princeton University Press, 2009).
  * Beltukov and Parshin (2016) Y. M. Beltukov and D. A. Parshin, JETP letters 104, 552 (2016).
  * O’Hern  _et al._ (2003) C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E 68, 011306 (2003).
  * Conyuh and YM (2020) D. Conyuh and B. YM, JETP Letters 112, 513 (2020).
  * Burda  _et al._ (2004) Z. Burda, A. Görlich, A. Jarosz, and J. Jurkiewicz, Physica A 343, 295 (2004).
  * Burda  _et al._ (2006) Z. Burda, A. Görlich, J. Jurkiewicz, and B. Wacław, The European Physical Journal B-Condensed Matter and Complex Systems 49, 319 (2006).
  * Zucker (2011) I. Zucker, Journal of Statistical Physics 145, 591 (2011).
  * Weiße  _et al._ (2006) A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, Reviews of modern physics 78, 275 (2006).
  * Gaspard and Cyrot-Lackmann (1973) J. Gaspard and F. Cyrot-Lackmann, Journal of Physics C: Solid State Physics 6, 3077 (1973).
  * Hodges (1977) C. Hodges, Journal de Physique Lettres 38, 187 (1977).
  * Conyuh  _et al._ (2019) D. Conyuh, Y. Beltukov, and D. Parshin, Journal of Physics: Conference Series 1400, 044026 (2019).
  * Adelman and Garrison (1976) S. Adelman and B. Garrison, The Journal of Chemical Physics 65, 3751 (1976).
  * Li and Weinan (2007) X. Li and E. Weinan, Physical Review B 76, 104107 (2007).
  * Wyart  _et al._ (2005) M. Wyart, L. E. Silbert, S. R. Nagel, and T. A. Witten, Physical Review E 72, 051306 (2005).
  * Akhiezer (1965) N. Akhiezer,  _The Classical Moment Problem and Some Related Questions in Analysis_ , University mathematical monographs (Oliver & Boyd, 1965).
  * Dagotto (1994) E. Dagotto, Reviews of Modern Physics 66, 763 (1994).
  * Saad (1992) Y. Saad,  _Numerical methods for large eigenvalue problems_ (Manchester University Press, 1992).
  * Kwidzinski and Bulla (2016) N. Kwidzinski and R. Bulla, arXiv preprint arXiv:1608.00616 (2016).