Ioffe-Regel criterion and viscoelastic properties of amorphous solids
IIoffe-Regel criterion and viscoelastic properties of amorphous solids
D. A. Conyuh and Y. M. Beltukov ∗ Ioffe Institute, 194021 Saint-Petersburg, Russia (Dated: January 1, 2021)We show that viscoelastic effects play a crucial role in the damping of vibrational modes inharmonic amorphous solids. The relaxation of a given plane wave is described by a memory functionof a semi-infinite one-dimensions mass-spring chain. The initial vibrational energy spreads from thefirst site of the chain to infinity. In the beginning of the chain, there is a barrier, which significantlyreduces the decay of vibrational energy below the Ioffe-Regel frequency. To obtain the parametersof the chain, we present a numerically stable method, based on the Chebyshev expansion of thelocal vibrational density of states.
I. INTRODUCTION
Damping of vibrational modes plays a crucial rolein the thermal conductivity of amorphous dielectrics(glasses). Low-frequency vibrations are well-definedphonons with a long mean free path. However, in a widerange of temperatures, the heat transfer in glasses is de-termined by another type of delocalized vibrations, whichare known as diffusons [1, 2]. The crossover between low-frequency phonons and diffusons at higher frequencies isknown as the Ioffe-Regel crossover [2, 3].In amorphous solids, the attenuation of plane elasticwaves (sound) is governed by multiple mechanisms: scat-tering on two-level systems [4–6] and soft modes [7, 8],thermally activated relaxation processes [9, 10], and scat-tering induced by structural and elastic disorder [11–15]. The last contribution is temperature independentand dominates the attenuation in the THz frequencyrange [11, 16].In the low-frequency range, there are phonons witha well-defined dispersion law ω ( q ) and a weak damp-ing Γ( q ) (cid:28) ω ( q ). In this case, the initial plane wavewith the wavevector q oscillates with the frequency ω ( q )with a slow exponential decay. This attenuation can bedescribed using the damped harmonic oscillator (DHO)model.However, the damping increases rapidly with increas-ing the wavevector q . For some wavevector | q | = q c , thedamping becomes comparable to the frequency, Γ( q ) ∼ ω ( q ), which corresponds to the Ioffe-Regel criterion. For q (cid:38) q c , the notion of the dispersion law ω ( q ) could notbe applied. It was shown that the DHO model can notbe used for an accurate analysis of vibrational propertiesfor frequencies above the Ioffe-Regel crossover [17, 18].Viscoelastic properties are important to study the high-frequency vibrations of amorphous solids [19–21].In the theory of viscoelastic relaxation of liquids, itis known that memory effects are important for the re-laxation of density fluctuations [22]. These memory ef-fects can be presented using Mori continued fraction [23]. ∗ Electronic address: [email protected]
In this paper we show that the vibrational relaxationin harmonic amorphous solids can also be described us-ing a general viscoelastic model with some memory func-tion K ( t ). In terms of vibrations, the continued fractionrepresentation corresponds to a semi-infinite mass-springchain, which reproduces the same memory effects. Wepresent a stable method to find parameters of the arbi-trary number of sites in the mass-spring chain.Another powerful tool to analyze the general proper-ties of disordered systems is the random matrix theory(RMT). Depending on the inherent symmetry proper-ties of different disordered systems, various random ma-trix ensembles are used [24]. It was shown that theWishart ensemble naturally arises in the study of vibra-tional properties due to the requirement of mechanicalstability [3, 25, 26]. In this paper we apply the RMT tofind the memory function and the corresponding repre-sentation using the mass-spring chain.This paper is organized as follows. In Section II weconsider a general viscoelastic relaxation of vibrationsin harmonic amorphous solids. Section III demonstratesthat the relaxation is described by a Green function inthe form of continued fraction, which corresponds to aone-dimensional mass-spring chain. In Section IV we ob-tain the parameters of the chain in the framework of theRMT. Section V demonstrates that the same approachcan be used to analyze any given numerical dynamicalmatrix. In Section VI we discuss the properties of ob-tained mass-spring chains and compare them to the Ioffe-Regel crossover. II. VISCOELASTIC DAMPING AND THEMEMORY FUNCTION
The general equation of motion of a solid near equilib-rium position can be written as | ¨ u ( t ) (cid:105) = − ˆ M | u ( t ) (cid:105) , (1)where ˆ M is N × N dynamical matrix with N being thenumber of degrees of freedom. N -dimensional vector | u ( t ) (cid:105) describes the deviation of atoms from the equi-librium position at time t . a r X i v : . [ c ond - m a t . d i s - nn ] D ec k q inst u q m q u q m q K q ( t ) k q inst η q a b FIG. 1: (Color online) a) Damped harmonic oscillator withthe mass m q (the ball), the stiffness k inst q (the spring), andthe damping η q (the dashpot). b) Harmonic oscillator with ageneral viscoelastic element defined by the memory function K q ( t ) (the rectangle). To study the relaxation of a plane wave with wavevec-tor q , we can solve Eq. (1) with initial conditions | u (0) (cid:105) = 0 , | ˙ u (0) (cid:105) = | q (cid:105) (2)for any given dynamical matrix ˆ M . The relaxation of theinitial plane wave is described by the projection u q ( t ) = (cid:68) (cid:104) u ( t ) | q (cid:105) (cid:69) , (3)where the big angle brackets denote the averaging overdifferent realizations of the dynamical matrix ˆ M . Assum-ing the normalization (cid:104) q | q (cid:105) = 1, we obtain u q (0) = 0 and˙ u q (0) = v from Eq. (2).In a simplified model, the relaxation of u q ( t ) can bedescribed using the DHO model m q ¨ u q ( t ) + η q ˙ u q ( t ) + k q u q ( t ) = 0 , (4)where the mass m q , the damping η q , and the stiffness k q may depend on the wavevector q . In terms of the DHOmodel, the Ioffe-Regel crossover separates weakly decay-ing long-wave vibrational modes and overdamped short-wave vibrational modes. However, it is important to takeinto 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 energydissipation.To analyze the relaxation process, we consider the re-solvent ˆ G ( z ) = (cid:28) z − ˆ M (cid:29) , (5)where z is a complex parameter. The relaxation of aplane wave with initial conditions u q (0) = 0, ˙ u q (0) = v can be written using the resolvent ˆ G ( z ) as (see Ap-pendix A) u q ( t ) = 12 π (cid:90) ∞−∞ ˜ u q ( ω ) e iωt dω, (6)where ˜ u q ( ω ) = − v G q (cid:0) ( ω − i (cid:1) , (7) G q ( z ) = (cid:104) q (cid:12)(cid:12) ˆ G ( z ) | q (cid:105) . (8) We can present the Green function G q ( z ) as the Stieltjestransform of the spatial Fourier transform of eigenmodes: G q ( z ) = (cid:90) ∞ F q ( ω ) z − ω dω, (9) F q ( ω ) = (cid:42)(cid:88) n (cid:104) q | n (cid:105)(cid:104) n | q (cid:105) δ ( ω − ω n ) (cid:43) , (10)where | n (cid:105) is n -th eigenmode. The Fourier transformof eigenmodes is closely related to the structure factor,which is S q ( ω ) = k B T q / ( mω ) F q ( ω ) [13, 27]. In thispaper we use the scalar model to simplify the notations.All qualitative results that we obtain can be applied tothe vector model as well. It was shown that vibrationsin the scalar and vector models belong to the same classof universality [28].In the DHO model with the natural frequency ω q andthe frequency-independent damping rate Γ q , we have F dho q ( ω ) = 2 π ω Γ q (cid:0) ω − ω q (cid:1) + ω Γ q , (11) G dho q ( z ) = − ω q − z + Γ q √− z . (12)It corresponds to Eq. (4) with the stiffness k q = m q ω q and the damping η q = m q Γ q . In the DHO model, themass m q can be chosen arbitrarily.In the general case, we can present the Green function G q ( z ) as G q ( z ) = − m q k inst q − m q z + G q ( z ) . (13)It corresponds to the viscoelastic equation of motion m q ¨ u q ( t ) + k inst q u q ( t ) + (cid:90) t −∞ K q ( t − t (cid:48) ) u q ( t (cid:48) ) dt (cid:48) = 0 (14)with the mass m q , the instantaneous stiffness k inst q , andthe memory function K q ( t ) = 12 π (cid:90) ∞−∞ G q (cid:0) ( ω − i (cid:1) e iωt dω. (15)In order to find m q , k inst q , and K q ( t ), we can present theGreen function G q ( z ) as a series G q ( z ) = ∞ (cid:88) k =0 F ( k ) q z k +1 (16)with moments F ( k ) q = (cid:90) ∞ ω k F q ( ω ) dω. (17)For large z , G q ( z ) is 1 /z due to the normalization F (0) q = (cid:82) ∞ F q ( ω ) dω = 1. For any harmonic system describedby Eq. (1), all moments F ( k ) q are finite. In this case, wecan assume that G q ( z ) is also 1 /z for large z . It resultsin the following values k inst q = m q F (1) q , (18) m q = (cid:104) F (2) q − (cid:0) F (1) q (cid:1) (cid:105) − , (19) G q ( z ) = m q z − k inst q − m q G q ( z ) . (20)The decreasing of G q ( z ) for large z corresponds to theabsence of the instantaneous component in the memoryfunction K q ( t ). III. CONTINUED FRACTION ANDONE-DIMENSIONAL CHAIN
We can repeatedly apply the same type of presentationfor the Green function: G n +1 , q ( z ) = m n q z − a n q − b n q G n q ( z ) (21)with starting Green function G q ( z ) defined by Eq. (20).As a result, we obtain the Mori continued fraction [23]: G q ( z ) = − m q k inst q − m q z + G q ( z )= − m q k inst q − m q z − b q a q − m q z + G q ( z ) = − m q k inst q − m q z − b q a q − m q z − b q a q − m q z + G q ( z ) = . . . (22)For each Green function G n q ( z ) we can find the corre-sponding function F n q ( ω ), the moments F ( k ) n q , and thememory function K n q ( t ) using the same relations as inEqs. (9), (15)–(17). A comprehensive set of relations isgiven in Appendix B.As before, we assume that each Green function G n q ( z )is 1 /z for large z , which determines the relation betweencoefficients a n , b n , and m n : a n q = m n q F (1) n q , b n q = m n q , (23) m n q = (cid:104) F (2) n q − (cid:0) F (1) n q (cid:1) (cid:105) − . (24)The recurrence relation (21) can be presented in the form F n +1 , q ( ω ) = m n q F n q ( ω ) (cid:12)(cid:12) G n q (cid:0) ( ω − i (cid:1)(cid:12)(cid:12) . (25)The mass m n q can be defined by the normalization condi-tion F (0) n +1 , q = (cid:82) ∞ F n +1 , q ( ω ) dω = 1, which is equivalentto Eq. (24). The details of the numerical realization ofthis recurrence procedure are discussed in Appendix D. b q b q k q k q k q b q k q inst – b q u q u q u q u q m q m q m q m q b q + b q k q k q b q k q inst – b q u q u q u q m q m q m q ab K q ( t ) FIG. 2: (Color online) a) A general case of the semi-infiniteone-dimensional mass-spring model. The displacements ofmasses m q , m q , m q , . . . are denoted by u q , u q , u q , . . . re-spectively. Each spring is denoted by its stiffness. b) Finiterepresentation of the same chain with 3 sites. The tail of thechain is replaced by a general viscoelastic element defined bythe memory function K q ( t ) and denoted by the rectangle. The continued fraction (22) can be presented using asolution of the following infinite system: m q ¨ u q ( t ) = − k inst q u q ( t ) + b q u q ( t ) , (26) m n q ¨ u n q ( t ) = − a n q u n q ( t ) + b n q u n − , q ( t )+ b n +1 , q u n +1 , q ( t ) , n ≥ , (27)where u q ( t ) ≡ u q ( t ) and u n q ( t ) = 12 π (cid:90) ∞−∞ ˜ u n q ( ω ) e iωt dω, (28)˜ u n q ( ω ) = − ˜ u n − , q ( ω ) b n q G n q (cid:0) ( ω − i (cid:1) . (29)It corresponds to the dynamics of a semi-infinite one-dimensional mass-spring chain (Fig. 2a) with masses m q and m n q , horizontal springs with stiffnesses b n q betweenconsequent masses, and vertical springs with stiffnesses k q = k inst q − b q , (30) k n q = a n q − b n q − b n +1 , q (31)between masses and the common ground. The initialcondition of the chain is u q (0) = 0, ˙ u q (0) = v , u n q (0) =0, ˙ u n q (0) = 0, which defines u q ( t ) for t > b n +1 , q u n +1 , q ( t ) = − (cid:90) t −∞ K n +1 , q ( t − t (cid:48) ) u n q ( t (cid:48) ) dt (cid:48) . (32)It corresponds to a finite chain shown in Fig. 2b. There-fore, the memory function K n q ( t ) describes a responseof the tail of the mass-spring chain starting from the sitewith the number n . Using the relations given in Ap-pendix B, the corresponding function F n q ( ω ) can be con-sidered as a local vibrational density of states (LVDOS)on site n for the chain under constraint u k q ( t ) = 0 for k < n . IV. A RANDOM MATRIX APPROACH
General vibrational properties could be studied usingthe random matrix approach [26]. This approach is basedon two main properties of amorphous solids: mechanicalstability and the invariance under continuous translationof an amorphous body. In the beginning of this section,we briefly discuss the main points of the random matrixapproach.The mechanical stability of amorphous solids is equiva-lent to the positive definiteness of the dynamical matrixˆ M . Any positive definite matrix ˆ M can be written asˆ M = ˆ A ˆ A T and vice versa, ˆ A ˆ A T is positive definite forany (not necessarily square) matrix ˆ A [29]. Therefore,we can consider a N × K random matrix ˆ A to obtain amechanically stable system with the dynamical matrix inthe form of the Wishart ensemble ˆ M = ˆ A ˆ A T .Each column of the matrix ˆ A represents a bond witha positive potential energy U k = (cid:0) (cid:80) i A ik u i (cid:1) with u i being the displacement of i -th atom from the equilibriumposition [26, 30]. Each row of the matrix ˆ A correspondsto some degree of freedom. In the random matrix ap-proach, the parameter κ = KN − K and the number of degrees of freedom N . In a stable system with a finite rigidity, the numberof bonds should be larger than the number of degrees offreedom, which is known as the Maxwell counting rule.For 0 < κ (cid:46)
1, the parameter κ has the same effect as theparameter z − z c in the jamming transition [31]. In realamorphous solids, one can estimate κ = 0 . U k should not de-pend on the shift u i → u i + const . Therefore, the matrixˆ A obeys the sum rule (cid:80) i A ik = 0. It means that the ma-trix elements A ij are correlated . In the minimal model,we can assume that amorphous solid consists of statisti-cally equivalent random bonds. In this case the pairwisecorrelations between matrix elements A ij can be writtenas (cid:104) A ik A jl (cid:105) = 1 N C ij δ kl , (34)where ˆ C is some correlation matrix. One can see thatthe correlation matrix ˆ C is proportional to the average dynamical matrix: ˆ C = NK (cid:10) ˆ M (cid:11) . For a system with sta-tistically equivalent bonds, the correlation matrix ˆ C is aregular matrix, which describes a lattice with some dis-persion law ω cor ( q ). We assume that there is only onebranch of ω cor ( q ). In the general case, one can apply thesummation over different branches below.Using the random matrix approach, it can be shownthat statistical properties of the random matrix ˆ M arerelated to the known correlation matrix ˆ C . In the ther-modynamic limit N → ∞ , there is a fundamental dualityrelation [33] z ˆ G ( z ) = Z ˆ G cor ( Z ) (35)between resolvents ˆ G ( z ) = (cid:10) ( z − ˆ M ) − (cid:11) and ˆ G cor ( Z ) =( Z − ˆ C ) − where complex parameters z and Z are relatedby a conformal map Z ( z ) defined by the equation Z (cid:0) κ + 1 + M cor ( Z ) (cid:1) = z. (36)The contour Z ( z ) in the complex plane for z = ( ω − i is known as the critical horizon [34]. The moment-generating function M cor ( Z ) for the correlation matrixˆ C is defined as M cor ( Z ) = ZN Tr ˆ G cor ( Z ) − ∞ (cid:88) k =1 M ( k )cor Z k (37)with moments M ( k )cor = 1 N Tr ˆ C k = 1 V bz (cid:90) bz ω k cor ( q ) d q . (38)Here integration is performed over the first Brillouinzone, which has the volume V bz . The moment-generatingfunction M cor ( Z ) can also be written as an integral overthe first Brillouin zone M cor ( Z ) = 1 V bz (cid:90) bz ω ( q ) Z − ω ( q ) d q . (39) A. Relaxation of a plane wave
Since the correlation matrix ˆ C is regular, for a givenwavevector q we obtain (cid:104) q | ˆ G cor ( Z ) | q (cid:105) = 1 Z − ω ( q ) . (40)Using the duality relation (35) and the conformal map Z ( z ) defined by Eq. (36), we can present the Green func-tion G q ( z ) as G q ( z ) = 1 z − (cid:0) κ + 1 + M cor ( Z ( z )) (cid:1) ω ( q ) . (41)For large z , we have Z ( z ) = (1 + κ ) z + O (1). Therefore,from Eq. (39), we obtain the asymptotics M cor ( Z ( z )) = (1 + κ ) M (1)cor z + O (cid:18) z (cid:19) . (42) / max n ( )( ) = 1 n = n = n = FIG. 3: (Color online) The LVDOS F n ( ω ) shown as a ratio F n ( ω ) / F ∞ ( ω ) for κ = 1 for different values of n (up to n =100). Therefore, the Green function G q ( z ) can be presented inthe form of Eq. (13) with k inst q = k inst = 1 / M (1)cor , (43) m q = (cid:2) ( κ + 1) M (1)cor ω ( q ) (cid:3) − , (44) G q ( z ) = G ( z ) = M cor ( Z ( z ))(1 + κ ) M (1)cor . (45)We omit the subscript q for the values, which do notdepend on the wavevector q . In the framework of theRMT, only the first mass m q depends on the wavevector q in the one-dimensional chain. All other parameters ofthe one-dimensional chain do not depend on q .The Green function G ( z ) corresponds to the LVDOS F ( ω ) = ω g ( ω )(1 + κ ) M (1)cor , (46)where g ( ω ) = ωπN Im Tr ˆ G (cid:0) ( ω − i (cid:1) is the full vibra-tional density of states for a given κ . The denominatorin Eq. (46) ensures the normalization of F ( ω ).We can use the recurrence procedure (23)–(25) to ob-tain the LVDOS F n ( ω ) and the parameters a n , b n , m n (all of them do not depend on the wavevector q ). For n → ∞ , we observe that F n ( ω ) converges to a stationarysolution of Eq. (25): F ∞ ( ω ) = 16 ω πω (cid:115) − ω ω . (47)where ω max is the maximum frequency in the system. Inthis case b ∞ = a ∞ ω , m ∞ = 16 ω , k ∞ = 0 . (48) k n m n c a Barrier = 1= 0.5= 0.2= 0.1 = 0.05= 0.02= 0.01= 0.0051.001.011.02 b n b b n c / max m n m c FIG. 4: (Color online) Parameters of the chain as a functionof scaled site number nω c /ω max obtained using the RMT fordifferent values of κ . For small values of κ , dots merge intolines. In panel (a), the dashed line marks the Ioffe-Regelcriterion k n /m n = ω c . For example, we consider the correlation matrix ˆ C asa regular matrix on a simple cubic lattice with a unitlattice constant. The non-diagonal elements are C ij = − i and j are neighbors and C ij = 0 otherwise.The diagonal elements are C ii = 6. The correspondingdispersion law is ω ( q ) = 4 (cid:16) sin q x q y q z (cid:17) . (49)In this case, the moment-generating function is M cor ( Z ) = W s (cid:0) Z − (cid:1) , where W s is the third Watsonintegral [35]. Using recurrence relation (25), we obtainthe LVDOS F n ( ω ) and the values of a n , b n , m n .Figure 3 shows F n ( ω ) / F ∞ ( ω ) for different values of n for κ = 1. One can observe that F n ( ω ) is a smoothfunction, which gradually approaches F ∞ ( ω ) for n → ∞ .The chain parameters for different values of κ are pre-sented in Fig. 4 as k n / ( m n ω c ), m n /m ∞ , and b n /b ∞ . Thecharacteristic frequency ω c ∼ κ is used for the scalingand will be discussed in the next section. One can ob-serve that k n /m n ∼ ω c for n (cid:46) ω max /ω c and k n /m n → n → ∞ . For the better understanding of the be-haviour of the obtained mass-spring chain, we considerthe low-frequency approximation, which can be analyzedanalytically. B. Low-frequency approximation
In order to study the low-frequency dynamics, we in-vestigate the behavior of the Green function G q ( z ) forsmall values of z , which corresponds to small values of Z . In this case, the generating function M cor ( Z ) can beapproximated as a linear function of Z for small Z : M cor ( Z ) = − − βZ + o ( Z ) , (50)where β = 1 V bz (cid:90) bz d q ω ( q ) . (51)For example, for a simple cubic lattice with the dispersion(49), the constant is β = w s / w s = 0 . Z , any generatingfunction of the form (37) have to decay as M cor ( Z ) ∼ /Z . Therefore, to study the low-frequency behavior, wecan consider the regularized form of Eq. (50): M cor ( Z ) = 1 βZ − , (52)In this case, Eqs. (43)–(45) becomes k inst = β, (53) m q = β (1 + κ ) ω ( q ) , (54) G ( z ) = β κ M cor ( Z ( z )) . (55)Using Eq. (36), we can present M cor ( Z ( z )) in the recur-rence form: M cor ( Z ( z )) = − κ + 1 κ + 2 − βz + M cor ( Z ( z )) . (56)Thus, we automatically obtain the continued fraction: G ( z ) = − m q k inst − m q z − b a − mz − b a − mz − ... , (57)where a = β κ κ , b = β √ κ , m = β κ . (58)The corresponding mass-spring chain is homogeneousstarting from n = 1, which is shown in Fig. (5). Eachmass in the semi-infinite chain (except the first one) isconnected to the ground by a spring with the stiffness k = a − b . It corresponds to the natural frequency ω c = (cid:114) km = √ κ − √ β , (59)which is the minimal frequency of vibrations, whichcan propagate along the chain. For small κ , we have ω c ≈ κ / (2 √ β ). This frequency is used for the scaling in b bk k kbk inst – bu q u q m q m m mu q u q FIG. 5: (Color online) A regular one-dimensional mass-springmodel, which corresponds to the low-frequency approxima-tion.
Fig. 4. The dashed line in Fig. 4a corresponds to the ho-mogeneous chain obtained in the low-frequency approxi-mation.Above the frequency ω c , the chain effectively absorbsthe vibrational energy of the first site. Below ω c , the ini-tial excitation distributes over several sites in the begin-ning of the chain without any further damping. Indeed,the frequency ω c coincides with the Ioffe-Regel frequencyin the random matrix approach [26]. In the low-frequencyapproximation (52) used in this section, no damping of vi-brations below the Ioffe-Regel frequency can be observed.In other words, all frequencies ω < ω c are localized inthe 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 ω < ω c , which will be discussed in Section VI.From Eq. (56), we can find M cor ( Z ( z )) explicitly: M cor ( Z ( z )) = 12 (cid:0) mz − b − k − √ mz − k √ mz − b − k (cid:17) . (60)Using Eqs. (45) and (15), we can find the correspond-ing memory function K ( t ) for k (cid:28) b and κ (cid:28) m q ¨ u q ( t ) + γ s ˙ u q ( t ) + k s u q ( t )+ (cid:90) t −∞ K l ( t − t (cid:48) ) u q ( t (cid:48) ) dt (cid:48) = 0 (61)with the short-term damping γ s = √ mb , the short-termstiffness k s = k inst − b − k/ √ kb and the long-termBessel memory function K l ( t ) = √ bkt J (cid:32) t (cid:114) km (cid:33) , t > . (62)Equation (61) defines the natural frequency of the firstsite in the chain ω ( q ) = (cid:115) k s m q = (cid:114) κ ω cor ( q ) . (63) u q ( t ) = 0.01 a t u q ( t ) = 0 b q = 0.1 q = 0.12 q = 0.16 q = 0.28 FIG. 6: (Color online) Relaxation of plane waves in the frame-work of the RMT for κ = 0 .
01 (a) and κ = 0 (b) for the sameset of wavenumbers q . Solid lines show the precise result,which corresponds to the general chain (Fig. 2). Dashed linesrepresent a relaxation in the low-frequency approximation,which corresponds to the regular chain (Fig. 5). If this frequency is above the Ioffe-Regel frequency ω c , astrong damping is observed. In terms of the wavevector,it corresponds to q > q c , where the Ioffe-Regel wavenum-ber is q c = 1 v cor (cid:114) κ β . (64)Here we use ω cor ( q ) = v cor q for small q . For q < q c , onecan find the resonance frequency in Eq. (61) ω ( q ) = v q (cid:112) β (2 q c − q ) . (65)In this case, the memory function K l ( t ) cancels out thedamping term γ s ˙ u q ( t ) and modifies the resonance fre-quency: ω ( q ) > ω ( q ). For q (cid:28) q c , we have ω ( q ) = √ ω ( q ).Figure 6 shows the relaxation of u q ( t ) for a simple cu-bic lattice with the dispersion ω cor ( q ) defined by Eq. (49).For κ = 0 .
01, the Ioffe-Regel wavenumber is q c ≈ . q c and two wavenumbers above q c . Solid linesshow the exact solution while dashed lines show the low-frequency approximation (61). One can see a good agree-ment between them. However, for large time t and q < q c , one can observe the slow relaxation of the exact solutionwhile the low-frequency approximation has stationary os-cillations.For κ = 0, the short-term stiffness k s and the long-term memory function K l ( t ) vanish, and we obtain vis-cous damping without returning force m q ¨ u q ( t ) + γ s ˙ u q ( t ) = 0 . (66)Figure 6b shows this behavior both for the exact solutionand the low-frequency approximation. V. A NUMERICAL ANALYSIS OFDYNAMICAL MATRICES
The same approach can be used to obtain the chainrepresentation for a given dynamical matrix or a givenensemble of dynamical matrices. These dynamical ma-trices can be obtained by using molecular dynamics sim-ulations or using numerical random matrix models.In the previous section, we consider the thermody-namic limit N → ∞ in the framework of the RMT. Inthis case, there are no fluctuations of vibrational prop-erties. However, in a finite system, fluctuations may beimportant, especially for the stability of the parametersof the mass-spring chain. In this section, we demonstratethat the proposed recurrence algorithm may be used fora finite numerical system as well.We use the numerical random matrix model in the formof the correlated Wishart ensemble ˆ M = ˆ A ˆ A T which canbe controlled by the same parameter κ . For simplicity,we consider a simple cubic lattice with random bonds andunit lattice constant as in the RMT. However, in contrastto the RMT, the numerical random matrix model has thefinite interaction radius.For κ = 0 the matrix ˆ A is square and the number ofbonds 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 ˆ A [3, 26]: A ij = (cid:26) ξ ij if i and j are neighbors , , (67)where ξ ij are independent Gaussian random numberswith zero mean and unit variance. The diagonal elementsare defined using the sum rule A ii = − (cid:80) j (cid:54) = i A ji . Thisprocedure results in the same correlation matrix ˆ C thatwas used in the framework of the RMT in the previousSection.For κ > A (0) and ˆ A (1) . The result-ing rectangular matrix ˆ A can be obtained by inserting κ N randomly chosen columns of the matrix ˆ A (1) into thematrix ˆ A (0) . This random insertion of the new columnscorresponds to the random addition of new bonds to thevibrational system.We use the Kernel Polynomial Method (KPM) [13] toobtain the initial function F q ( ω ) for the recurrence rela-tion (25). The initial function F q ( ω ) was calculated for q ( ) = 1, q = 1 max FIG. 7: (Color online) The function F q ( ω ) calculated usingthe KPM for a numerical random matrix with N = 100 atoms, the parameter κ = 1, and the wave number q = 1.Averaging over 100 realizations was applied. Vertical dashedline shows the chosen position of ω max . a system with N = 100 atoms and different values ofthe parameter κ and the wavenumber q . The averagingover 10 – 10 realizations was applied. An example for κ = 1 and q = 1 is shown in Fig. 7.We use the Chebyshev expansion of F q ( ω ) in the fre-quency range 0 ≤ ω ≤ ω max to evaluate the recurrencerelation (25) (see Appendix D). The choice of ω max isimportant because there is an exponential tail in thehigh-frequency vibrational density of states without anyspecific maximum frequency. For any finite system size,there are big relative fluctuations in the high-frequencytail due to a small number of vibrations there (see Fig. 7).It may lead to additional fluctuations of the parameters a n q , b n q , m n q of the obtained chain. Thus, we leave asmall number of vibrational modes above ω max and donot use them in the Chebyshev expansion, which signif-icantly reduces the fluctuations of the obtained parame-ters.The KPM is also based on the Chebyshev expan-sion [36]. However, the maximum frequency in the KPM, ω kpm max , should be larger than any frequency in the sys-tem for stability purposes. Therefore, we remap the ob-tained Chebyshev expansion to another one with slightlysmaller maximum frequency ω max to drop a small num-ber of high-frequency modes.The resulting parameters k q n , b q n , m q n of the mass-spring chain are presented in Fig. 8 for different valuesof the parameter κ and the wavenumber q . The resultsare similar to those obtained in the framework of theRMT (Fig. 4). In the given scale, all calculated values k n / ( κ m n ) almost coincide in Fig. 8a. The masses m q n and stiffnesses b n q are close to their stationary values m ∞ and b ∞ respectively. For different wavenumbers q ,we observe similar parameters of the mass-spring chain. k n q m n q a = 1, q = 0.5= 1, q = 1= 0.5, q = 0.5= 0.5, q = 1= 0.25, q = 0.5= 0.25, q = 10.981.001.02 b n q b b n / max m n q m c FIG. 8: (Color online) Parameters of the chain as a functionof scaled site number κ n/ω max obtained using the numericalrandom matrix model (67) for different values of the param-eter κ and the wavenumber q . VI. DISCUSSION
We have shown that the relaxation of a plane wavewith the wavevector q coincides with the relaxation of thefirst site in the semi-infinite one-dimensional mass-springchain (Fig. 2). The parameters of the chain form thecontinued fraction representation of the Green function G q ( z ).Each site in the chain is connected to the ground by aspring with a stiffness k n q . We observe the same behaviorof k n q in both RMT and numerical analysis of finite ran-dom matrices: in the beginning of the chain k n q ∼ k > κ , we observe decaying oscillationsof k n q . It is known that oscillations of the coefficientsin the continued fraction are related to singularities inthe density of states [37, 38]. In our case, the vibrationaldensity of states has a steep behavior near the Ioffe-Regelfrequency for small κ [26, 39]. However, these oscillationsare not important for qualitative analysis.For large enough site number n , the parameters a n q , b n q , m n q , k n q converge to their stationary values a ∞ , b ∞ , m ∞ , k ∞ , which depend only on the maximumfrequency in the system ω max . For electronic systems,it is known, that the coefficients of the continued frac-tion have stationary values, which depend on the widthof the energy band [37]. Since k ∞ = 0, the tail of themass-spring chain is free and homogeneous. This tailcan be considered as a simple thermal bath, which fi-nally absorbs the initial vibrational energy [40, 41]. Atthe same time, no explicit damping is introduced in themass-spring chain, which corresponds to the absence ofdamping in the initial equation of motion (1).In the studied models, all masses m n q are close to thestationary value m ∞ (except the first one). Thus, thechain can easily absorb the vibrational energy above theIoffe-Regel frequency ω c = (cid:112) k/m ∞ . For frequencies be-low ω c , the absorption is much smaller because these fre-quencies are “forbidden” in the beginning of the chain.Depending on the form and the width of the barrier, thereis a relatively small absorption below ω c . The dynamicsof the chain can be mapped to a discrete version of aone-dimensional Schr¨odinger equation with the potentialenergy k n q /m n q and the energy ω . Therefore, the re-gion with k n q /m n q > ω acts as the tunneling barrier(or a high-pass filter).The first mass in the chain strongly depends on thewavevector q : m q ∼ ω − ( q ). It results in a natural fre-quency ω ( q ) = (cid:112) κ / ω cor ( q ) of the first site. If thisfrequency is smaller than the Ioffe-Regel frequency ω c for a given wavevector q , then the damping of this vi-brational mode is relatively slow. It corresponds to thenotion of phonons with well-defined dispersion ω ( q ). If ω ( q ) > ω c , a strong damping is observed, which corre-sponds to the notion of diffusons above the Ioffe-Regelfrequency.One can note that ω ( q ) is proportional to the Lapla-cian on the corresponding lattice. Therefore, in theframework of the RMT, the equation of motion (14) canbe rewritten in the real space as ρ ¨ u ( r , t ) + k inst ∆ u ( r , t ) + (cid:90) t −∞ K ( t − t (cid:48) )∆ u ( r , t (cid:48) ) dt (cid:48) = 0 . (68)This viscoelastic equation is not local in time, but localin space. In a general case, the instantaneous stiffness k inst q and the memory function K q ( t ) may depend on thewavevector q which results in additional spatial convo-lution in Eq. (68). In the low-frequency approximation,from Eq. (61) we obtain ρ ¨ u ( r , t ) + γ s ∆ ˙ u ( r , t ) + k s ∆ u ( r , t )+ (cid:90) t −∞ K l ( t − t (cid:48) )∆ u ( r , t (cid:48) ) dt (cid:48) = 0 (69)with long-term Bessel memory function K l ( t ) defined byEq. (62). For κ = 0, from Eq. (66) we obtain a viscousequation without any returning force ρ ¨ u ( r , t ) + γ s ∆ ˙ u ( r , t ) = 0 . (70)The case κ = 0 is known as the isostatic state in thejamming transition [42]. In this case ω c = 0 and the en- tire low-frequency range is occupied by diffusons. Equa-tion (70) is consistent with a model of random walks ofatomic displacements for isostatic case [3].The above equation was obtained in the scalar model.In general case, u ( r , t ) is a vector and all constants andmemory functions in the above equations become tensors.The continued fraction presentation (22) is well-knownin the classical theory of moments [43]. However, thedirect evaluation of moments leads to stability and per-formance issues for n (cid:38)
30. The proposed method isbased on the Chebyshev expansion (Appendix D) andshows the numerical stability both for the RMT (Sec-tion IV) and numerical random matrices (Section V). Ittakes about one hour on a modern computer to find up to10 coefficients in the continued fraction from the knownfunction F q ( ω ). Usually, the most time-consuming partis the calculation of F q ( ω ) and its averaging over differentrealizations of the dynamical matrix ˆ M .For a given dynamical matrix ˆ M , the continued frac-tion (22) can be also obtained using the Lanczos methodusing | q (cid:105) as a starting vector [44]. However, there aretwo important issues concerning the Lanczos method.It is known that the Lanczos method is unstable dueto the loss of orthogonality. To stabilize the algorithm,an additional reorthogonalization is required, which de-creases the performance [45]. The second issue is thatthe Lanczos method is highly sensitive to a small num-ber of high-frequency eigenmodes. There is no simpleway to discard these eigenvalues or average the result-ing parameters a n q , b n q , m n q over different realizations ofthe dynamical matrix ˆ M in the framework of the Lanc-zos method. In the proposed algorithm, the averaging isperformed directly on F q ( ω ) and a small number of high-frequency modes can be easily removed as was discussedin Section V. VII. CONCLUSION
We have shown that viscoelastic relaxation of planewaves in amorphous solids can be considered as dynamicsof the one-dimensional semi-infinite chain. The first sitein this chain represents the initial plane wave while therest of the chain represents the memory effects. The ini-tial vibrational energy gradually spreads along the chain,which results in the vibrational relaxation of the initialplane wave.In the beginning of the chain, there is a natural bar-rier for frequencies ω < ω c , which corresponds to theIoffe-Regel crossover. In the framework of the RMT, thememory function does not depend on the wavevector q ,which results in the viscoelastic equation, which is lo-cal in space but not local in time. In the low-frequencyapproximation, the long-term memory function can bedescribed by the Bessel function.The proposed method demonstrates numerical stabil-ity for studied theoretical and numerical random matri-ces.0 VIII. ACKNOWLEDGMENTS
We wish to acknowledge D. A. Parshin and A. V. Shu-milin for valuable discussions. The authors thank theRussian Foundation for Basic Research (project no. 19-02-00184) for the financial support.
Appendix A: Evolution of plane waves
For initial conditions | u (0) (cid:105) = 0 and | ˙ u (0) (cid:105) = | q (cid:105) , thesolution of (1) is | u ( t ) (cid:105) = v (cid:88) n | n (cid:105) sin( ω n t ) ω n (cid:104) n | q (cid:105) , (A1)where | n (cid:105) is n -th eigenvector of the matrix ˆ M and ω n isthe corresponding eigenfrequency. Therefore, the projec-tion to the plane wave is u q ( t ) = (cid:42) v (cid:88) n (cid:104) q | n (cid:105) sin( ω n t ) ω n (cid:104) n | q (cid:105) (cid:43) . (A2)Using the Fourier transform, we can write u q ( t ) = 12 π (cid:90) ∞−∞ ˜ u q ( ω ) e iωt dω, (A3)where˜ u q ( ω ) = lim (cid:15) → + (cid:90) ∞ u q ( t ) e − iωt − (cid:15)t dt = (cid:42) v lim (cid:15) → + (cid:88) n (cid:104) q | n (cid:105)(cid:104) n | q (cid:105) (cid:90) ∞ sin( ω n t ) ω n e − iωt − (cid:15)t dt (cid:43) = (cid:42) v lim (cid:15) → + (cid:88) n (cid:104) q | n (cid:105) ω n − ( ω − i(cid:15) ) (cid:104) n | q (cid:105) (cid:43) = − v (cid:10) q (cid:12)(cid:12) ˆ G (cid:0) ( ω − i (cid:1)(cid:12)(cid:12) q (cid:11) . (A4)It is worth to note that u q ( t ) defined by Eq. (A3) is zerofor t < Appendix B: Viscoelastic relations
In this Appendix we provide the most important re-lation between different viscoelastic functions. Relationsbetween the Green function G n q ( z ), the LVDOS F n q ( ω ),and the moments F ( k ) n q are: F n q ( ω ) = 2 ωπ Im G n q (cid:0) ( ω − i (cid:1) , (B1) G n q ( z ) = (cid:90) ∞ F n q ( ω ) z − ω dω = ∞ (cid:88) k =0 F ( k ) n q z k +1 , (B2) F ( k ) n q = (cid:90) ∞ ω k F n q ( ω ) dω. (B3) Relations between the Green function G n q ( z ) and thememory function K n q ( t ) are: K n q ( t ) = 12 π (cid:90) ∞−∞ G n q (cid:0) ( ω − i (cid:1) e iωt dω, (B4) G n q ( z ) = (cid:90) ∞ K n q ( t ) e − t √− z dω. (B5)Relations between the LVDOS F n q ( ω ) and the memoryfunction K n q ( t ) are: K n q ( t ) = − θ ( t ) (cid:90) ∞ F n q ( ω ) ω sin( ωt ) dω, (B6) F n q ( ω ) = − ωπ (cid:90) ∞ K n q ( t ) sin( ωt ) dt, (B7)where θ ( t ) is the Heaviside step function. Appendix C: Low-frequency memory function
In the case k (cid:28) b , we can write the generating function(60) in the form M ( Z ) = M s ( Z ) + M p ( Z ) + M l ( Z ) + o (cid:18) kb (cid:19) , (C1)where the main term and two kinds of perturbations havethe following form M s ( Z ) = 12 (cid:16) mz − b − √ mz √ mz − b (cid:17) , (C2) M p ( Z ) = k (cid:18) √− bmzmz + mz − b √ mz √ mz − b (cid:19) , (C3) M l ( Z ) = (cid:112) b ( k − mz ) − √− bmz. (C4)Using Eq. (15) with z = ( ω − i , we obtain the corre-sponding memory functions K s ( t ) = − bt J (cid:0) ˜ t (cid:1) θ ( t ) , (C5) K p ( t ) = k (cid:114) bm (cid:32) J (cid:0) ˜ t (cid:1) (cid:104) ˜ t − π tH (cid:0) ˜ t (cid:1)(cid:105) − J (cid:0) ˜ t (cid:1) (cid:104) − π tH (cid:0) ˜ t (cid:1)(cid:105) − (cid:33) θ ( t ) , (C6) K l ( t ) = √ bkt J (cid:32) t (cid:114) km (cid:33) θ ( t ) , (C7)where ˜ t = 2 t (cid:112) b/m is the scaled time, J n is the Besselfunction, H n is the Struve function, and θ ( t ) is the Heav-iside step function.The memory function K s ( t ) is the main short-termmemory function, which is not zero for k = 0. It coin-cides with the memory function for the semi-infinite freemass-spring chain [41, 46]. The memory function K p ( t )1is a perturbation of the short-term memory function fornonzero k . The memory function K l ( t ) is a long-termmemory perturbation since it depends on another scaledtime t (cid:112) k/m , which scales with k .In the low-frequency approximation, the memory func-tions K s ( t ) and K p ( t ) can be considered as an instanta-neous response: (cid:90) t −∞ K s ( t − t (cid:48) ) u q ( t (cid:48) ) dt (cid:48) ≈ (cid:90) t −∞ K s ( t − t (cid:48) ) (cid:2) u q ( t ) + ( t (cid:48) − t ) ˙ u q ( t )) (cid:3) dt (cid:48) = − bu q ( t ) + √ mb ˙ u q ( t ) , (C8) (cid:90) t −∞ K p ( t − t (cid:48) ) u q ( t (cid:48) ) dt (cid:48) ≈ (cid:90) t −∞ K p ( t − t (cid:48) )) (cid:2) u q ( t )+( t (cid:48) − t ) ˙ u q ( t )) (cid:3) dt (cid:48) = − k u q ( t ) . (C9) Appendix D: Chebyshev expansion