{'fraction_non_alphanumeric': 0.033380059575959346, 'fraction_numerical': 0.014280707902575783, 'mean_word_length': 4.227965713537918, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Ioffe-Regel criterion and diffusion of vibrations in random lattices

We consider diffusion of vibrations in 3d harmonic lattices with strong force-constant disorder. Above some frequency w_IR, corresponding to the Ioffe-Regel crossover, notion of phonons becomes ill defined. They cannot propagate through the system and transfer energy. Nevertheless most of the vibrations in this range are not localized. We show that they are similar to diffusons introduced by Allen, Feldman et al., Phil. Mag. B 79, 1715 (1999) to describe heat transport in glasses. The crossover frequency w_IR is close to the position of the boson peak. Changing strength of disorder we can vary w_IR from zero value (when rigidity is zero and there are no phonons in the lattice) up to a typical frequency in the system. Above w_IR the energy in the lattice is transferred by means of diffusion of vibrational excitations. We calculated the diffusivity of the modes D(w) using both the direct numerical solution of Newton equations and the formula of Edwards and Thouless. It is nearly a constant above w_IR and goes to zero at the localization threshold. We show that apart from the diffusion of energy, the diffusion of particle displacements in the lattice takes place as well. Above w_IR a displacement structure factor S(q,w) coincides well with a structure factor of random walk on the lattice. As a result the vibrational line width Gamma(q)=D_u q^2 where D_u is a diffusion coefficient of particle displacements. Our findings may have important consequence for the interpretation of experimental data on inelastic x-ray scattering and mechanisms of heat transfer in glasses.

We consider diffusion of vibrations in 3d harmonic lattices with strong force-constant disorder. Above some frequency ωIR, corresponding to the Ioffe-Regel crossover, notion of phonons becomes ill defined. They cannot propagate through the system and transfer energy. Nevertheless most of the vibrations in this range are not localized. We show that they are similar to diffusons introduced by Allen, Feldman et al., Phil. Mag. B 79, 1715 (1999) to describe heat transport in glasses. The crossover frequency ωIR is close to the position of the boson peak. Changing strength of disorder we can vary ωIR from zero value (when rigidity is zero and there are no phonons in the lattice) up to a typical frequency in the system. Above ωIR the energy in the lattice is transferred by means of diffusion of vibrational excitations. We calculated the diffusivity of the modes D(ω) using both the direct numerical solution of Newton equations and the formula of Edwards and Thouless. It is nearly a constant above ωIR and goes to zero at the localization threshold. We show that apart from the diffusion of energy, the diffusion of particle displacements in the lattice takes place as well. Above ωIR a displacement structure factor S(q, ω) coincides well with a structure factor of random walk on the lattice. As a result the vibrational line width Γ(q) = Duq 2 where Du is a diffusion coefficient of particle displacements. Our findings may have important consequence for the interpretation of experimental data on inelastic x-ray scattering and mechanisms of heat transfer in glasses.

I. INTRODUCTION
Propagation of vibrational excitations in disordered systems is one of the advanced problems in condensed matter physics. In particular, transport mediated by these excitations is responsible for the thermal conductivity of amorphous dielectrics (glasses). However mechanisms of heat transfer in glasses above the plateau region are still poorly understood.
At low temperatures below 1 K the low frequency plane long wave acoustical phonons are well defined excitations which transfer the heat in glasses. At these temperatures the thermal conductivity κ(T ) ∝ T 2 and is controlled by a resonant scattering of phonons on two-level systems (TLS) 1,2 . Between 4 K and 20 K the thermal conductivity κ(T ) saturates and displays a well known plateau 3 . As was shown in 4 it can be explained by resonant scattering of phonons by quasilocal vibrations (QLV). The QLV, together with TLS and phonons are vibrational excitations responsible for many universal properties of glasses 5 .
However, above approximately 20 K the thermal conductivity rises again and finally saturates on the level of one order of magnitude higher, at temperatures about several hundreds Kelvin 6 . As generally believed, the origin of this second rise of the thermal conductivity (above the plateau) is not related to phonons. It was established long ago [7][8][9] , that in this temperature (frequency) range the mean free path of phonons l becomes of the order of their wave length λ (or even smaller, of the order of interatomic distance). Correspondingly, the Ioffe-Regel criterion for phonons 10 becomes violated. The existence of such crossover was confirmed by molecular dynamics calculations for some real and model glasses 11,12 and disordered lattices 13,14 .
In the regime of such strong scattering a standard concept of plane waves (phonons) with a well defined wave vector q becomes inapplicable. The question then arises: what physical mechanism is responsible for the heat transfer in glasses in this temperature range? The numerical simulations show that majority of the vibrational modes in the corresponding frequency range are not localized [15][16][17] .
As was shown in 18-20 , a lower limit of the thermal conductivity of amorphous solids above 30 K can be correctly estimated within the framework of the Einstein's model 21 . It was assumed that the mechanism of heat transport above the plateau is a random walk of thermal energy between clusters of neighboring atoms vibrating with random phases. In fact, a diffusion mechanism for the heat transfer in this temperature range was proposed.
At the same time, delocalized vibrations in glasses of a new type, different from phonons, were introduced. They were called diffusons [22][23][24][25][26] . These are vibrations spreading through the system not ballistically, as phonons (on distances of the order of mean free path) but by means of diffusion. It is an important class of excitations which occupy in glasses the dominant part of the spectrum 26 . In these papers it was put forward the hypothesis that the boundary between phonons and diffusons is determined by the Ioffe-Regel criterion for phonons. Since diffusons are delocalized excitations, they may be responsible for the thermal conductivity of glasses above the plateau.
The similar conclusion was made by the authors of 27,28 . They considered the case of strong scattering of phonons in disordered lattices with a significant fraction of randomly located missing sites, but which is still far from the percolation threshold. It was shown that, in contrast to the electronic case, the Ioffe-Regel criterion is inaccurate in the prediction of phonon localization. Instead of localization, the vibrational transport above the Ioffe-Regel threshold becomes diffusive with approximately constant energy diffusivity D(ω). The diffusivity was calculated by numerical solution of the Newton equations for particle displacements. Similar calculations but for real glasses were done in the papers 29,30 using molecular dynamics methods.
The diffusons above the Ioffe-Regel crossover were identified also in granular jammed systems with repulsive forces between the particles 31,32 . They also have diffusivity which is independent of frequency ω. It was calculated making use of the Kubo-Greenwood formula for the thermal conductivity derived in 23 . In jammed systems the Ioffe-Regel crossover frequency ω IR can vary. It is shifted to zero when the system approaches the jamming transition point and rigidity goes to zero.
Therefore, as we believe, it is important to study properties of diffusons systematically in systems where they exist. They bring a new physics to our understanding of vibrational properties in strongly disordered systems and energy/heat transfer in glasses. To study these properties, we should have a model being sufficiently simple but still allowing to describe all of them.
Since we consider harmonic models, the simplest but still rather general, are (scalar or vector) models where particles, placed in equilibrium positions, are connected by random elastic springs. The equilibrium positions can be taken on a lattice 13,14,33,34 , or randomly distributed in space 35,36 . With some exceptions, there is no principal difference between these two cases because equilibrium positions do not enter to the dynamical matrix. The only important features are the type of disorder in elastic spring constants and the topology of the bonds. If all spring constants are positive, one can study different situations taking different distributions of random springs to explain existing experimental data 33 .
However a problem appears when some spring constants take negative values 13,34,36 . Even if the number of negative springs and their absolute values are relatively small, in the absence of any correlation with positive springs they produce an inevitable mechanical instability in all infinite systems. Therefore, a situation of strong force-constant disorder when appreciable concentration of springs can take negative values comparable with positive ones, is not possible in these models. On the other hand, an inclusion of negative spring constants can considerably improve an agreement between theory and existing experimental data 13,34,36 .
One can solve this problem mathematically, using a stable random matrix approach [37][38][39][40] . In this approach positive and negative springs are tangly correlated with each other. These correlations automatically guaranty the mechanical stability of the system independently of character and strength of the force constant disorder. In the present paper we are going to use this approach to investigate diffusion of vibrational excitations in disordered lattices with strong force-constant disorder. Some of our preliminary results have been presented in short form elsewhere 40 .
The paper is organized as follows. In Section II for the sake of clarity of further consideration, we outline the main properties of the model. We consider disordered lattices with strong force-constant disorder, described by stable positive definite random dynamical matrix AA T having positive eigenvalues only. Matrix A is a random matrix (not necessary symmetric) built on simple cubic lattice, with statistically independent matrix elements between the nearest neighbors A i =j , having zero mean We show that the density of states g(ω) is not zero at ω = 0 and phonons cannot propagate through the lattice. Similarly to systems at jamming transition point, the rigidity of the lattice is also zero. However the physical reason is different. In our case it is due to high concentration of negative springs (about 45%) in the system what makes it extremely soft. The participation ratio P (ω) indicates that all modes with exception of high frequency part are delocalized. As it is shown by further investigation, all of them are diffusons. In Section III we consider slightly additively deformed dynamical matrix AA T + µM 0 which has phonon-like excitations at small frequencies. Here positive definite matrix M 0 (random or non-random) is independent of A. µ is a parameter of the model which can vary in the interval 0 µ < ∞. Analyzing properties of this matrix, we calculate the Young modulus E, sound velocity, the density of states and participation ratio, the dynamical structure factor S(q, ω), the phonon dispersion law ω q , and also their mean free path l(ω). Comparison of the later with phonon wave length λ determines the Ioffe-Regel crossover frequency ω IR . It goes to zero when µ → 0. We show that above ω IR , phonons cease to exist. They are transformed to diffusons. In Section IV we consider properties of diffusons. The number of diffusing physical quantities coincides with the number of integrals of motion. In a closed free mechanical system there are two integrals of motion, momentum and energy. In Section IV A we investigate diffusion of momentum. We show that for all masses equal, it is equivalent to the diffusion of particle displacements since center of inertia is conserved. The displacement structure factor S(q, ω) coincides well with the structure factor S rw (q, ω) of the random walk on a lattice. We introduce new additional diffusion coefficient D u which is diffusivity of particle displacements. It is different from the energy diffusivity D(ω) investigated in [22][23][24][25][26] . We calculate the correlation function of particle displacements C(r, ω) and radius of diffuson. In Section IV B, we investigate the diffusion of energy D(ω) using two different approaches. The first approach is a direct numerical solution of Newton's equations. In the second approach the diffusivity is calculated by means of Edwards and Thouless formula which relates the energy diffusivity with an infinitesimal change of boundary conditions. Both approaches give similar results. We show that diffusivity D(ω) is independent of frequency in the diffuson range. In Section V we discuss scaling properties of the model (their dependence on parameters V and µ). We show that they are similar with systems near jamming transition point. In Section VI we discuss the obtained results and compare them with experiment.

II. A RANDOM MATRIX APPROACH
In harmonic approximation vibrational properties of a mechanical system of N particles are determined by the dynamical matrix the force constant matrix and m i are the particle masses.
The matrices M and Φ are real, symmetric and positively definite matrices N × N ( for simplicity we will consider a scalar model). The condition of positive definiteness is important. It ensures mechanical stability of the system. One can always present every real, symmetric and positive definite matrix M in the following form 37,38 Here A is some real matrix of a general form (not necessarily symmetric). And, vice versa, for every real matrix A the product AA T is always a positively definite symmetric matrix. Matrices AA T belong to Wishart ensemble 41 . The eigenvalue distribution for such kind of large random matrices was firstly investigated in 42 . For a free mechanical system it is necessary to satisfy also conditions 43 (for simplicity we consider below all masses m i = 1). It ensures that the potential energy of the system and forces between the particles depend only on the differences of particle displacements u i − u j . As a result, the potential energy and forces are not changed under any translation of the system as a whole. These conditions are necessary (but not sufficient) for existence of low frequency acoustic phonon-like modes in the system. If conditions (2) are violated, we have spatially pinned system where propagation of Goldstone modes (phonons) is not possible. In structural glasses in many cases (as, for example, in vitreous silica or amorphous silicon) a mass disorder is not important and we usually deal with the force constant disorder. It is related to fluctuations of valence bond lengths and valence bond angles because of an absence of crystalline ordering. Since valence forces depend exponentially on the distances between the atoms, they can experience strong fluctuations. Due to positional disorder there are also fluctuations of long distance Coulomb forces in non covalent materials. Thus the force-constant disorder plays an essential role in glassy dynamics. Therefore, one may expect that some important properties of glasses can be reproduced if we take matrix A as a random one.
As soon as the random elastic spring constants −M ij connecting the particles are fixed, the exact equilibrium particle positions are no longer important for particle dynamic on a long length scales much bigger than the interatomic distances. They do not enter to the dynamical matrix M . Therefore, it is reasonable to consider harmonic lattice models involving only force constant disorder.
If disorder is sufficiently strong, then it automatically includes the coordination number disorder since any weak bonding is equivalent to a negligibly small interaction between the two neighbors. Also, if the average coordination number is sufficiently big, its exact value and, therefore, the type of the lattice is of no importance as well. The similar considerations concern polarization of the modes. A vectorial character of vibrations in real glasses makes the issue to be more complicated. Many universal properties of glasses (as, for example, the thermal conductivity) are not related to polarization of the modes. Therefore, for better understanding of the physics involved it is often instructive to exploit scalar models. In these models the problem of zero frequency modes existing in vectorial isostatic lattices 44 does not exist 13 . For these reasons, different scalar models were successfully used in glassy physics in the past 13,33,35,36 .
Due to the reasons mentioned above we, as in 40 , consider the case of a simple cubic lattice with N particles and lattice constant a 0 = 1. Each particle has its unique integer index i which takes values from 1 to N . We construct the random matrix A as follows. The nondiagonal elements A ij (for i = j) we take as independent random numbers from Gaussian distribution with zero mean A ij = 0 and unit variance A 2 ij = V 2 = 1 if i-th and j-th particles are nearest neighbors. For each particle in a simple cubic lattice there are six nearest neighbors. As a result, for given i we have 6 non zero non-diagonal elements A ij for matrix A. Non-diagonal elements A ij and A ji are statistically independent from each other (matrix A is non-symmetric). All other nondiagonal elements (for non-nearest neighbors) A ij = 0. To ensure the property (2) the diagonal elements A ii are calculated as follows Then, according to Eq. (1), the Eq. (2) will be also met. The matrix M is then constructed according to Eq. (1). As was shown in 39 , in a 3d simple cubic lattice each particle is connected by elastic springs with 24 neighbors. The elastic spring constants are random and can be either positive or negative. A negative spring by definition is a spring which expands after initial stretching and shrinks after initial contraction. The effect of negative spring constants on atomic vibrations was discussed in different papers 13,34,35,[45][46][47][48] .
To elucidate where this coordination number 24 in 3d case comes from, let us first consider as example a 2d simple square lattice. In this case each particle interacts with 12 neighbors shown on Fig. 1. In accordance to Eq. (1), matrix elements connecting central black particle with its 4 nearest neighbors are of the type From Eq. (4) it follows that diagonal elements meet the following relations One has to insert them in Eq. (5) Since averaged values A i =j = 0 and different nondiagonal matrix elements A ij are statistically independent from each other, the average value M 01 is determined by the first two quadratic terms in Eq. (8). As a result, it is non-zero and negative. It corresponds to positive average elastic spring k 01 = −M 01 between particles 0 and 1 Though, according to Gaussian distribution of A i =k , the spring constant k 01 can take negative values as well. All the aforesaid is valid for other nearest neighbor matrix elements M 02 , M 03 and M 04 . The next nearest neighbor matrix elements are given by It is easy to see that In 3d case for simple cubic lattice there are 6 springs of the type M 01 , 12 springs of the type M 05 and 6 springs of the type M 09 . As a result all together we have 24 particles interacting with the central black particle. All these 24 spring constants can be either positive or negative but to ensure the mechanical stability of the whole system they are correlated with each other in a rather tangly way.
Distributions of different spring constants are shown on Fig. 2. The distribution of k 01 is asymmetric with positive mean value. The distributions of k 05 and k 09 are even (with zero average value) and for k 09 are given by zeroth-order Macdonald function 39 which logarithmically diverges at k = 0. The resulting distribution of all spring constants was calculated numerically in 40 . The number of negative springs was found to be about 45%. One can find a similarity between our spring constant distributions and dynamical matrix element distributions obtained in 14 for IC-glass, in 49 for simple fluid with shortranged interactions (see Fig. 1 in these papers), and in 50 for realistic model of amorphous silicon (see Figs 2.12,2.13). Though it is difficult to compare our scalar model with vector models analyzed in 14,49,50 .
Concluding this part, we can easily include into consideration the next neighbor shell for matrix A. Then, in addition to the previous case, the matrix elements of the type A 05 should be taken into account. As a result the coordination number for matrix M in simple cubic lattice increases up to 90. Just opposite, applying some additional constraints, we can reduce the coordination number from 24 to smaller numbers or make it fluctuating quantity, etc. We have checked that all these modifications can lead to quantitative changes but do not change qualitatively the main results of the paper. Therefore, we will restrict our consideration by the simplest case outlined before.  As follows from the figure, the spectrum is gapless i.e. g(ω) is nonzero at ω → 0. In spite of the fact that the conditions (2) are fulfilled, we do not see the expected phonon modes with their DOS g ph (ω) ∝ ω 2 for ω → 0. It means that phonons as plane wave excitations cannot propagate through the lattice. This result is not changed qualitatively 40 neither by including next neighbor shells to build matrix A nor by switch to vector model.

P(ω)
As was shown in 40 , such a behavior of DOS at ω → 0 is related to the fact that the affine assumptions are violated and the macroscopic elasticity theory becomes inapplicable in this case. The average value of the static Young modulus of the lattice E ∝ 1/N . Therefore, in the thermodynamic limit (N → ∞) E → 0. As a result, the rigidity of the lattice and sound velocity are also tend to zero. This unusual behavior is due to a presence of high concentration of negative springs (45%) in the lattice which makes it to be extremely soft.
To determine whether vibrational modes are localized or delocalized, we have calculated the participation ratio Here e i (ω) is i-th particle projection of the normalized eigenvector with frequency ω. As one can see from the inset of Fig. 3, all modes with exception of small high frequency part are delocalized. They have P (ω) ≈ 0.2 which is independent of the system size. This value is close to the theoretical value 1/3 for Porter-Thomas distribution of e 2 i (ω) 39,51 . We have verified also that the level spacing distribution obeys the Wigner-Dyson statistics 39 . It also indicates the mode delocalization. As we will show in Section IV, all these delocalized gapless vibrational modes can be identified as diffusons. They spread in the lattice by means of diffusion.
To elucidate a spacial structure of the eigenmodes for matrix M = AA T we considered as an example a two dimensional square lattice with N = 400 × 400 particles and calculated eigenvector e i (ω min ) (i = 1, 2, ..., N ) for the lowest frequency ω min in the system. The result is shown on Fig. 4. Particles with positive and negative displacements are shown by white and black dots correspondingly. As one can see from the figure, the mode is delocalized. Its spatial structure is random (fractal) and has nothing to do with a plane wave. Similar picture takes place in a 3d case.

III. PHONONS
To introduce phonons into the picture we should have finite rigidity of the lattice. The rigidity can be introduced by different means. Since a sum of positive definite matrices is a positive definite matrix, then simplest possi-bility is to add to the random matrix AA T a "crystalline part" 40 Here A is the same random matrix built on a 3d simple cubic lattice with a 0 = 1 as in the previous Section. Matrix M 0 is a positively definite crystal dynamical matrix for the same lattice with unit masses, and all spring constants (between the nearest neighbors) equal to unity. As was shown in 40 the tune parameter µ 0 controls the rigidity of the lattice.
Adding the regular part µM 0 , changes the distribution of spring constants k 01 between the nearest neighbors, as shown on Fig. 2. The average value is equal to k 01 = 2+µ. At small values of µ 1 the change is negligible. The distribution mainly consists from strongly fluctuating part AA T (compare the distributions of k 01 for µ = 0 and µ = 0.1). Therefore, it is not obvious at all that such small perturbation is able to introduce a finite rigidity and phonons into the system. A strong scattering of phonons by diffusons may leave the diffuson spectrum unchanged.
Also we can consider elastic springs in matrix µM 0 to be fluctuating quantities, somehow distributed in the closed interval [0, µ]. Otherwise we can cut out a big amount of springs µ from the lattice, so the phonons cease to exist in the term µM 0 at all (see Section VIII A). Another nontrivial possibility is shown in Section VIII B. In the paper we limit ourselves by the most simple case described by Eq. (13). To find the rigidity (as a function of µ), we calculated numerically the Young modulus E of the lattice with dynamical matrix given by Eq. (13) for µ = 0. In modeling of amorphous solids, the standard method to do that is to use Irving-Kirkwood stress tensor formula 52 . However, it is difficult to implement this procedure in our case of strong local fluctuations of elastic springs when microscopic displacement field u(r) is not a differentiable function of atomic positions 40 .
Therefore, to avoid these difficulties we, as in 40 , used a direct numerical method. We took a very big cubic sample with N = L × L × L = 10 6 particles and side L − 1 to reduce fluctuations and possible non-affine response. According to the standard textbook formula of the macroscopic elasticity theory (see Eq. 5.2 in 53 ), Young modulus is given by E = σ zz /u zz . Here σ zz is the stress, and u zz is the strain. The component u zz gives the relative lengthening of the sample ∆L/(L−1). Then we fixed the strain and calculated the stress σ zz .
For that we fixed particles on the left hand side of our cubic sample and displaced all particles on the opposite (right hand) side by the unit distance ∆L = 1. Since Newton equations are linear, the final result is independent of the value of the step strain used. In other two directions we used the periodic boundary conditions. Then, solving the system of linear Newton equations, we found the new equilibrium positions of all other particles in the sample and calculated restoring forces f i acting on the displaced particles on the right boundary. Due to randomness of the elastic bonds, the restoring forces are also random. Letf be the average restoring force. Then, by definition, the stress σ zz = i f i /L 2 =f and the Young modulus E can be calculated as follows To avoid confusion, we remind that we are using here a scalar version of the elasticity theory. Therefore, all forces in the lattice are parallel (or antiparallel) to the particle displacements.
The results of these calculations are shown on Fig. 5 for cubic sample with N = 10 6 particles 54 . As we can see from the fit, the Young modulus has a following dependence on µ: As a result, for µ 1 we have a usual crystal, where disorder is relatively small and relation (15) is obvious. For µ 1 the force constant disorder is strong. The fluctuations of the nondiagonal matrix elements M i =j are much bigger than the averaged values 39,40 . In this case Young modulus E ∝ √ µ. It is much bigger than the crystal result (15). Strong fluctuations of the positive and negative elastic springs which in average almost compensate each other make the lattice much more rigid than in the case of crystal. Therefore for µ 1 one can not consider our lattice as a simple superposition of two systems AA T and µM 0 . The origin of this behavior E ∝ √ µ is unclear and it should be elucidated in further work (see also Section V). But below we will support our numerical findings by calculation of the sound velocity and of the phonon density of states (for small ω) and by a comparison of the latter with total DOS calculated numerically for matrix (13). Below in this paper we will consider the case of strong and moderate force constant disorder when 0 µ 1. Filled and open diamonds correspond to phonon contribution to the DOS below and above the Ioffe-Regel crossover frequency ωIR correspondingly (see further text for details). Inset: dependence To calculate the phonon contribution to the DOS at small ω, we need to know the sound velocity v at zero frequency. It is related to the Young modulus in a standard way: (since all particle masses m i = 1 and lattice constant a 0 = 1). Then for the phonon DOS (in the scalar model) we have The total DOS g(ω), normalized to unity and calculated numerically by the kernel polynomial method (KPM) 55,56 for dynamical matrix (13) and different values of µ, is shown on Fig. 6. We see from the figure that for µ = 0 the DOS at low enough frequencies is proportional to ω 2 which corresponds to acoustical phonon excitations. Thus, introducing finite values of µ, we open up a soft phonon gap in the gapless diffuson spectrum, existing at µ = 0. The DOS in the gap, as we will show in the paper, is built by acoustic phonon-like modes and at low frequencies goes to zero as g(ω) ∝ ω 2 . The term phonon gap is motivated since, if conditions (2) are violated, then addition µM 0 to random matrix AA T opens a hard gap in the gapless vibrational spectrum (see Fig. 8 below). Just above this gap the DOS has a sharp maximum at frequency ω max which we will identify with the width of the gap. As follows from the figure, the maximum frequency for µ 1 increases as ω max ∝ √ µ. Above the maximum the vibrational excitations remain to be diffusons (see Section IV).
One can try to explain the dependence ω max ∝ √ µ as follows. In the absence of random part AA T the dynamical matrix M is determined by the crystalline part µM 0 only. Then (for simple cubic lattice) we have well defined phonon modes with dispersion law The maximum frequency in this case is equal to However the sound velocity in this pure crystallyne lattice case v cryst = √ µ. Though according to Eqs. (17, It means that simple superposition approach does not work in this case and physical picture is more complicated. As we will show in Section V the Young modulus E depends also on the amplitude of the random part AA T . Since the DOS g(ω) is normalized to unity for all values of µ, we conclude from Fig. 6 (comparing the DOS for µ = 0 with DOS for µ = 0) that vibrations corresponding to the maximum for µ = 0 were pushed out from the region of small frequencies ω < ω max for µ = 0. We see also from the figure that, after initial ω 2 dependence, the DOS for µ = 0 increases much faster than ω 2 . It is a clear signature of the presence of the boson peak in our disordered lattice. As we will show further (see Table I), the frequency ω max is correlated with position of the boson peak ω b (the maximum in the reduced DOS g(ω)/ω 2 ). Therefore appearance of the boson peak in disordered systems is not necessarily related to the acoustic van Hove singularity in crystals as was proposed recently 13,34,57 .
The straight lines on the Fig. 6 correspond to the phonon DOS g ph (ω) determined by Eq. (18) with the sound velocity v = √ E and E calculated from Fig. 5. One can see a good agreement of the total g(ω) at low frequencies with the phonon contribution g ph (ω). From that we can conclude that at least the low frequency excitations in the phonon gap are the usual long-wave acoustical phonons. However, actually, as we will show further, nearly all excitations in the gap up to the frequencies close to ω max correspond to phonons, but with a nonlinear dispersion law. This conclusion is supported by calculations of the participation ratio P (ω). It is shown in Fig. 7 for various values of µ. For µ = 0, one can clearly distinguish in the function P (ω) a presence of the two different frequency regions. As follows from Fig. 6, the low frequency part (below ω max ) corresponds to the phonons. In this range the participation ratio increases with decreasing frequency. It is related to increase of the phonon mean free path l(ω) as ω → 0 (see Fig. 12). In the high frequency part (above ω max ) P (ω) is approximately inde- pendent of the frequency and coincides with participation ratio for µ = 0. As we will show in Section IV this range corresponds to diffusons. A similar rise of the participation ratio with decreasing frequency was found recently in 2d Lennard-Jones glasses 58 (see Fig. 1b of this paper). It is important to emphasize that for existence of the acoustical phonon excitations the conditions (2) are crucial. If they are not obeyed, then, instead of soft phonon gap in the vibrational spectrum shown on Fig. 6, we have a hard gap shown on Fig. 8. Inside the hard gap there are no vibrations at all. The dynamical matrix M in this case was taken in the same form (13). But diagonal elements A ii of the matrix A were taken as independent Gaussian random variables with average A ii = 0 and unite variance A 2 ii = 1. As a result the condition (4) (and therefore (2)) was violated and we have got a spatially pinned lattice where low frequency acoustical phonon modes cannot exist. However, the width of the hard gap in this case has the same µ dependence as the width of the phonon gap, ω max ∝ √ µ.
To find the phonon dispersion curve (dependence of the phonon frequency ω on the wave vector q) and phonon mean free path l(ω) we should calculate space and time Fourier transform of the particle displacement field u(r, t). For that we ascribed to all the particles at the initial moment t = 0 random displacements u(r, 0) (from Gaussian distribution with zero mean and unit variance) and zero velocities. Then, numerically solving Newton equations (with all masses m i = 1) we analyzed the particle dynamics at t = 0. In calculations we used Runge-Kutta-4 method with sufficiently small time step ∆t = 0.01. We have checked that in this case the total energy of the system is conserved over the whole investigated time interval T with relative precision 10 −7 without use of any damping. The calculated values of particle displacements also have relative precision higher then 10 −7 (we compared the results with time step ∆t = 0.01 with results obtained with two times smaller time step ∆t = 0.005).
Let u(r i , t) be the i-th particle displacement as a function of particle coordinate r i and time t. We define the displacement structure factor (DSF) of the displacement field as follows For better frequency resolution, the upper time limit T was taken sufficiently large (T = 3000), while the integration time step was chosen as ∆t = 0.01. Since vectors r i in a cubic lattice are discrete, the wave vectors q ≡ q n are also discrete and are defined on the corresponding reciprocal lattice. For example, for cubic sample L × L × L and q 100 direction we have q n = 2πn/L where integer numbers n are −L/2 n L/2. One can show (see Section VIII C) that definition (20) is equivalent to the usual expression Here e i (ω j ) -is i-th component of the eigenvector of the dynamical matrix M corresponding to i-th particle and eigenfrequency ω j 43 . The normalized density of states is related to the structure factor by the sum rule According to definition (20) S(0, ω) = 0 since the position of center of inertia is conserved and i u(r i , t) = 0.
To analyze phonon excitations, we have found the maximum of S(q, ω) as a function of ω for each discrete value of q n , for several values of µ. As an example, the results for µ = 0.1 and one q direction are shown on Fig. 9. For the fitting curves we used the Lorentz distribution From this fit we can find both the phonon frequency ω q and the phonon line width ∆ω. The results for ω q are shown on Fig. 10 for three values of µ and q 100 . For sufficiently small values of wave vector q we see a nice linear dispersion curve ω q = vq, with the sound velocity v given by Eq. (17). It is independent of the q direction (i.e. the sound velocity is isotropic). With increase of q, the frequency ω q shows a pronounced negative dispersion of the group velocity v g = dω q /dq and approaches the maximum frequency ω max where the dependence ω q saturates. In this q region we observed a weak anisotropy of the dispersion curves for µ = 1. At smaller values of µ the dependence ω q is isotropic. Since ω max ∝ √ µ, the vertical axis on Fig. 10 scales approximately as √ µ and the horizontal axis scales as µ 1/4 (sound velocity v ∝ √ E ∝ µ 1/4 , and q max ≈ ω max /v ∝ µ 1/4 as well).
The strong negative dispersion of the group velocity v g for big q values can be explained by avoided crossing principle (or level repulsion effect) due to the coupling of phonons to quasilocal vibrations near frequency ω max , corresponding to sharp maximum in DOS g(ω) (see Fig. 6). Similar phenomenon exists in polariton physics 59 . The dip in the participation ratio P (ω) for µ = 0.001, µ = 0.01 and µ = 0.1 at ω ≈ ω max (see Fig. 7) evidences in favor of this idea. The vibrations inside the dip correspond to frequencies near ω max and have smaller participation ratio than the others. Therefore they can be referred to as quasilocal vibrations. In the following we will see that this strong scattering is also responsible for the deep minimum in the diffusivity D(ω) at ω ≈ ω max (see Fig. 20).
The negative dispersion of the group velosity v g is responsible also for the pronounced rise of the phonon DOS above the ω 2 dependence, given by Eq. (18). It is clearly seen on the Fig. 6. Indeed, taking the dispersion into account and disregarding weak anisotropy (taking place only for µ = 1) we can write instead of Eq. (18) Here v g (ω) = dω/dq is the group velocity shown in Insets on Fig. 10. Taking for q(ω) and v g (ω) the data from Fig. 10 we obtain the points (filled and open diamonds) shown on Fig. 6. Since they perfectly coincide with numerical data for g(ω) below ω max , we conclude that all the excitations in the phonon gap belong to phonons (with nonlinear dispersion at higher values of q).
The phonon line width ∆ω can be also found from fits similar to those shown on Fig. 9. It is related to the phonon life time τ = 1/2∆ω. The factor 2 takes into account that ∆ω corresponds to decay of the amplitude of the vibration. The results are shown on Fig. 11. As follows from this figure, ∆ω ∝ ω 4 and does not depend on the direction of q. We think that this frequency dependence is not due to Rayleigh scattering of phonons on a static disorder. In such a case ∆ω would be proportional to q 4 . Due to nonlinear dispersion in ω q , these dependencies do not correspond to each other. More likely, the phonon line width is due to strong resonant scattering of phonons by quasilocal vibrations responsible for the sharp peak in the DOS, similar to those introduced in 4 . The deep minimum in the diffusivity D(ω) around frequency ω max also supports this idea (see Fig. 20). We hope to investigate this important question in future work.
With known value of ∆ω, the phonon mean free path l(ω) can be calculated as follows The phonons are well defined excitations if their mean free path l(ω) exceeds the phonon wave length λ = 2π/q (Ioffe-Regel criterium for phonons). As we will see in the next Section, phonons transform to diffusons when l(ω) ≈ λ/2. We will call the corresponding crossover frequency as ω IR . Fig. 12 shows the ratio l(ω)/λ as a function of ω for several values of µ and different directions of the wave vector q. The boundary between filled and open symbols (the full horizontal line) corresponds to frequency ω IR . Thus filled and open symbols on Figs. 6, 10, 11, 12 belong to phonons with frequencies below and above the Ioffe-Regel crossover frequency correspondingly. Usually in glasses the Ioffe-Regel crossover frequency ω IR is correlated with position of the boson peak ω b , see [60][61][62][63][64] and references therein. It is the frequency where the reduced DOS g(ω)/ω 2 has a maximum. We also have a rather sharp boson peak in our disordered lattices 40 . As follows from Fig. 6   from phonons having negative dispersion of the group velocity dω q /dq. Similar conclusion was made recently for 2d and 3d Lennard-Jones glasses 58,65,66 . The right side of the boson peak consists from diffuson modes shifted from the region of small frequencies 0 < ω < ω max by additional term µM 0 and further modified by interaction with phonons. But more work is necessary to elucidate the precise structure of these modes.
The frequencies ω max , ω IR , and ω b are collected in Table I for different µ. As we can see from the table, ω IR is close to the frequency ω max and to the position of the boson peak ω b . Above ω IR phonons cease to exist as well defined excitations. They are smoothly transformed to diffusons which we will consider in the next Section. The relative number of phonons in the lattice can be estimated as follows These values are also given in the Table I. We see that for all investigated values of µ the relative number of phonons in the lattice is small. It is in agreement with similar estimates for amorphous silicon 26 .   The frequency of maximum in DOS ωmax, the frequency of the Ioffe-Regel crossover ωIR and the boson peak frequency ω b for various µ. Star * means that ωIR was found for q 100 direction. N ph is a relative number of phonons in the lattice.

IV. DIFFUSONS
In this section we are going to consider properties of diffusons. As is well known, the diffusion phenomenon usually takes place for physical quantities which are conserved. In a free closed mechanical system we have two integrals of motion, momentum and energy. Therefore one should discriminate between diffusion of momentum and energy.

A. Diffusion of momentum
First let us consider diffusion of momentum. Usually the diffusion of momentum is related to viscosity in the system. When all particle masses being equal (m i = 1), the diffusion of momentum is equivalent to the diffusion of particle displacements. It is because in our system the position of the center of inertia is conserved and we can put it at the origin of the coordinate system. Then the sum of all particle displacements vanishes i.e. it is an integral of motion. The diffusion of displacements in this case looks like a diffusion of "particles" in a lattice where the total number of particles is conserved. By analogy with diffusion of "particles" the information about diffusivity of displacements is absorbed in the displacement structure factor S(q, ω) (20). We remind that to calculate this structure factor we ascribed at the initial moment t = 0 the random displacements to all the particles with Gaussian distribution (with zero mean and unit variance) and velocities equal to zero. So the condition (27) at t = 0 was satisfied. Therefore let us analyze now this structure factor in the diffuson frequency range.
Consider first the case of µ = 0 when phonons are absent and only diffusons are present in the lattice. Fig. 13 shows the structure factor S(q, ω) as a function of wave vector q for three different directions in q space (symbols) and for three different frequencies ω. Let us compare this displacement structure factor with structure factor of the random walk S rw (q, ω) on the lattice.
As was shown in 67 for the case of the random walk on a lattice, S rw (q, ω) is given by expression It is a Lorentzian, with a width Γ(q) given by where D rw is a diffusion constant of the random walk. In a simple cubic lattice (with lattice constant a 0 = 1) the function Q(q) reads Q(q) = 2 sin 2 q x 2 + sin 2 q y 2 + sin 2 q z 2 .
For small values of q 1, Q(q) = q and in the continuum limit we have the well known result for the diffusion structure factor S rw (q, ω) = 2D rw q 2 D 2 rw q 4 + ω 2 .
Let us note that the structure factor (28) has a maximum at q values obeying the condition We can specify it as a dispersion law for diffusons. The width of the maximum is Γ(q). For q 1, Γ(q) = D rw q 2 . A comparison of the displacement structure factor S(q, ω), (20), and the structure factor of the random walk S rw (q, ω), (28), is shown on Fig. 13. One fitting parameter was the diffusion coefficient D rw in Eq. (29). From comparison of these data we obtained D rw ≈ 0.7. It means that the diffusion coefficient of particle displacements D u ≈ 0.7 (see Section VI). Another fitting parameter was a height h(ω) of the random walk structure factor in the maximum. According to Eq. (28), in the maximum Γ(q) = ω and h(ω) = 1/ω, but to fit the data points on Fig. 13 we used slightly higher values of h(ω).
The small difference between h(ω) and 1/ω can be explained by different frequency dependencies of the density of states g(ω) for vibrations and for the random walk (following from the sum rule similar to Eq. (22)). As we can see from the figure, for the investigated frequencies the fit is perfect. With increasing frequency above ω ≈ 2 − 3, the fitting becomes more and more poor since we approach the localization threshold at ω loc ≈ 5.5±0.5 (see below) which is not described well by a simple model of Markovian random walk. Now let us consider a behavior of a correlation function. The correlation function of particle displacements at some frequency ω, expressed through eigenvectors e r (ω) of the dynamical matrix M , reads C(r, ω) = r e r +r (ω)e r (ω).
Let us compare this correlation function with correlation function of the random walk. For distances bigger than the period of the lattice (a 0 = 1) we can make use of the limit of small q 1 and integrate Eq. (31) for the random walk structure factor taken in approximation of continuous media. As a result, we derive Fig. 14 shows a good agreement of our correlation function (34) with the correlation function of the random walk (35). For all investigated frequencies the numerical data collapse together and become indistinguishable from the theoretical prediction (35). We can see also on this figure the anticorrelation phenomenon (the region of negative values of the correlation function). As follows from Eq. (35), the correlation function of the random walk changes its sign for the first time at It is also in a good agreement with our numerical results. Therefore we can call a corresponding value of r found from Eq. (36) as a radius of diffuson. It is a typical size of the regions vibrating with frequency ω and having the same sign of all particle displacements. According to (36), the radius of diffuson is given by At ω = 0 the correlation function (35) decays slowly as 1/r. In disordered systems at critical point the correlation function decays as C(r) ∝ 1/r d−D2 where d is the space dimension and D 2 is a correlation dimension. From this we conclude that in our case D 2 = 2 what corresponds to diffusion. Now let us analyze the displacement structure factor S(q, ω) for µ = 0. For better visual effect we will show a map of the function S(q, ω) on the plane (ω, q) for different directions in q space. To do that, for each frequency ω we have found the maximum S(q, ω) as a function of q along some directions in q space. Then we normalized function S(q, ω) along this line ω=const to the magnitude of this maximum.
The results are shown on Fig. 15 for four different values of µ and two directions in q space. The white color corresponds to the maximum when normalized structure factor S n (q, ω) = 1 while the black color to the case where S n (q, ω) = 0. For µ = 0 we can see clearly two types of excitations in the lattice. At low enough frequencies, below ω IR , we see phonons with well defined dispersion law ω q , the same as in the previous Section. At the Ioffe-Regel crossover frequency ω IR , the structure factor strongly broadens and phonon dispersion line disappears. Above ω IR the displacement structure factor coincides well with the structure factor for µ = 0 case shown on Fig. 15a, which corresponds to diffusons. The maximum of the normalized structure factor S n (q, ω) (white regions) agrees well with Eq. (32) (with the same diffusion coefficient D rw ) giving the maximum of the random walk structure factor S rw (q, ω) (black line). It means that diffusion coefficient of particle displacements is independent of µ. Deviations from S rw (q, ω) take place at high frequencies near the localization threshold.
For µ = 0 the radius of diffuson (37) takes a maximum value at ω ≈ ω IR . At smaller frequencies we have well defined phonons. Since ω IR ∝ √ µ and D rw ≈ 1 we can write for 0 < µ 1 The value r c plays a role of correlation length in our lattice. It diverges when µ → 0. The physical meaning of this length is that it by the order of the value coincides with the Ioffe-Regel wave length λ IR = 2π/q IR corresponding to frequency ω IR (see Section V). Samples with size smaller than r c have no phonon-like modes at al.
To compare phonon and diffuson structure factors, a cross section of the structure factor S n (q, ω) in q space for q z = 0 and frequency ω = 0.5 is shown on Fig. 16 for µ = 0 and µ = 0.1. At the left side (a) of this figure we see the structure factor of diffuson. On the right side we see the structure factor of phonon (b). As compared with phonon structure factor, the diffuson structure factor is much more broadened.

B. Diffusion of energy
Now let us consider the diffusion of energy. The diffusion of energy is different from diffusion of particle displacements (see Section VI). The first approach to calculate the diffusivity of energy D(ω) for vibrations with frequency ω is a direct numerical solution of Newton's equations. For that we have used the Runge-Kutta-4 method with time step ∆t = 0.01 applied to a cubic sample with N = L × L × L particles (lattice constant a 0 = 1) and with free boundary conditions along the x direction. Along other two directions we take the periodic boundary conditions. Assuming zero initial conditions for displacements and velocities of all the particles, let us apply external forces with frequency ω and random phases ϕ i to all the parti- where ωT 1. The right and the left sides of the sample have coordinates x r,l = ±L/2. In such a way we excite vibrations with frequencies near frequency ω distributed in a small frequency interval (ω −1/T, ω +1/T ). In calculations we used T = 5 for all frequencies ω. We started our calculations at time t 0 = −5T when the external force is still negligible.
After applying the force to the central layer x = 0, vibrations will spread to the left and to the right ends of the sample. The average squared distance to the energy diffusion front we define as usual Here x i is the x coordinate of the i-th particle, E i (t) is the energy of i-th particle and sum is taken over all particles in the sample. E tot = i E i (t) is the total energy of the system. It is independent of time after the external force f ext i (t) becomes negligibly small (i.e. for t > 5T ). The energy of i-th particle E i (t) we define as a sum of the kinetic energy and a half of the potential energy of connected bonds (m i = 1) Here v i (t) =u i (t) is a particle velocity. Summation over all particles in Eq. (40) we can divide in two steps. First we sum over all particles in the layer x and then we sum over all layers. Let E(x, t) be a total energy confined to the layer x at time t. Having in mind that in our case we have lattice constant a 0 = 1 and sample size L 1, we can change summation over different layers to integration over coordinate x for times where R(t) 1. We will apply this method to the case of µ = 0 (i.e. for the lattice without phonons). The results are shown on Fig. 17. As we can see from the figure for small and middle frequencies, R 2 (t) ∝ t. Therefore for these frequencies vibrations indeed spread along the x axis by means of diffusion. The slope of the lines decreases with frequency ω. For calculating the slope, we take the time interval ∆t where, on the one hand t > 5T , and on the other hand, R L/2. From the slope of R 2 (t) we can calculate the diffusivity of modes D(ω) using one dimensional formula

D(ω)
This diffusivity is shown by black dots on Fig. 18. At small frequencies it is approximately constant, then it decreases with frequency approaching zero at the localization threshold, ω loc ≈ 5.5 ± 0.5. At higher frequencies above ω loc the dependence R 2 (t) saturates with increasing t. This indicates localization of the vibrational modes.
The difference between delocalized and localized modes is clearly seen if we examine the dependence E(x, t) as a function of coordinate x at some moment t for two different frequencies below and above the localization threshold. These two points for investigation are shown on Fig. 17. Black diamond corresponds to delocalized mode with frequency ω = 2 and has coordinates t = 234 and R 2 = 166. The distribution of energy E(x, t) over the x , a.u. sample calculated numerically at this moment is shown by black diamonds on Fig. 19. The data are perfectly fitted by solid line drawn according to the solution of diffusion equation in 1d case with value of R 2 = 166. Black triangle on Fig. 17 corresponds to localized mode with frequency ω = 6 and has coordinates t = 900 and R 2 = 22. The distribution of energy E(x, t) over the sample calculated numerically at this moment is shown by black triangles on Fig. 19. This distribution is drastically different from the previous case. For localized modes we expect the usual exponential decay The fit of the numerical data with this function and R 2 = 22 is shown on Fig. 19. The fit is perfect except for the central point at x = 0 which lies noticeably above prediction of Eq. (44). The coefficients in Eqs. (43,44) were taken to satisfy the obvious rules (45) To find the diffusivity D(ω) for µ = 0, the method of numerical solution of Newton's equations is not appropriate, because in this case we have phonons in the lattice with long mean free paths. Correspondingly samples with much bigger sizes are necessary to use this approach. Therefore for µ = 0 we used a second approach.
In this approach, the diffusivity D(ω i ) at eigenfrequency ω i was calculated by means of the formula of Edwards and Thouless 69 where L is the length of the sample and ∆ω i is sensitivity of the eigenfrequency ω i to a twist of boundary conditions. More precisely, we defined the diffusivity as follows: where ϕ is the angle of twisting, and c is some constant of the order of unity. It will be determined from comparison with the Newton method. The averaging in Eq. 47 is performed over frequencies ω in the small interval (ω − δω, ω +δω) with δω = 0.25 and/or over several thousands realizations.

D(ω)
For all other bonds M jk = M jk . So ∆ω i is the difference between i-th eigenfrequencies of matrices M and M Twisting of boundary conditions was performed for x direction only. For others two directions the periodic boundary conditions were used.
For µ = 0 the results for D(ω) are shown on Fig. 18 for three different cubic samples (full lines). We compared these results with numerical solution of Newton equations for µ = 0 (black dots) and get for the constant c ≈ 1. Then we used this c value for µ = 0. The results are shown on Fig. 20. For µ = 0 we see clearly two different frequency regions in the function D(ω).
At low frequencies, diffusivity increases with decreasing of ω. This range corresponds to the phonons. Indeed, the diffusivity of phonons D(ω) can be calculated as follows Open symbols on Fig. 20 show contribution calculated from this equation (just below Ioffe-Regel threshold). We see a good agreement with Edwards and Thouless formula. After a deep minimum at frequency ω ≈ ω max the diffusivity D(ω) saturates at a constant level (independent of µ) coinciding with D(ω) for µ = 0. The diffusivity in this range corresponds to diffusons. Similar behavior of D(ω) was found recently in jammed systems 31,32 . The deep minimum in the diffusivity at ω ≈ ω max corresponds to strong scattering of phonons by the quasilocal vibrations near the sharp peaks in the DOS g(ω) (see Fig. 6).

V. SCALING RELATIONS
Finally, the concept of diffusons allows us to establish useful scaling relations between observable values and important parameters of our model. One parameter is µ. It has a dimensionality of frequency squared. The second important parameter of the model is the variance of nondiagonal elements A ij of the random matrix A which we hitherto considered to be equal to unity The parameter V has dimension of frequency and assigns the scale of typical frequencies in the system. In particular, the normalized density of states g(ω) for µ = 0 shown on Fig 3 has the following scaling relation Since for ω below ω IR in our disordered lattice we have phonons with ω = vq (here v is sound velocity) and above ω IR we have diffusons with ω = Dq 2 (here D = D rw ) we can write at the Ioffe-Regel threshold the order of the magnitude estimates From these equations it follows that Since, according to Eq. (17) v = √ E (the units of mass and length we put equal to unity, m = a 0 = 1), we find for the Young modulus a useful relation Because, as we have shown in Section IV, the diffusion coefficient D is independent of µ, the Young modulus has the same µ dependence as ω IR √ µ (see inset on Fig. 6).
It is in a full agreement with Fig. 5 for µ 1. From the dimensionality considerations (since V has a dimension of frequency) we have for the diffusivity D the following estimate It is quite natural since the typical diffusion jump length is of the order of lattice constant a 0 = 1 and typical jump frequency is of the order of typical frequency in the system V . Therefore, the diffusivity D V a 2 0 . Taking this into account, for the Young modulus and sound velocity at small frequencies we have for i.e. the Young modulus is proportional to the characteristic frequency in the system V . The correlation length (38) (58) Though our paper is not aimed at jamming transition and we consider completely different model, it is interesting to note that these scaling relations are identical to those found in jamming transition 32 . Authors 32 study a model of amorphous packing of frictionless spheres interacting via the repulsive pair potential where the distance between the centers of particles i and j is denoted by r ij and the sum of their radii by σ ij . This model system, irrespective of the value of α, exhibits a jamming/unjamming transition at T = 0 at a packing fraction φ = φ c at which the particles are just touching each other and there is no overlap 70 . At densities lower than φ c particles are free to rearrange while above φ c at ∆φ ≡ φ − φ c , the system behaves as a weakly connected amorphous solid with an average coordination number that scales as a power law with an exponent where z c = 2d, with d being the space dimension. It was found that different quantities exhibit scaling behavior near the jamming point. According to 32 the Ioffe-Regel crossover frequency ω * and the shear modulus G behave as (we use below the notation of the paper 32 ) The transverse sound velocity v t and the diffusivity in The applied pressure p and the plateau in the density of states D 0 depend on the packing fraction as follows 70 Thus if we put then the crossover frequency ω IR , the Young modulus E, sound velocity v, the diffusivity at the plateau D, and the density of states g(ω) in our model have the same scaling as the crossover frequency ω * , the shear modulus G, transverse sound velocity v t , the diffusivity in the plateau d 0 , and the density of states D 0 in the jamming transition model respectively. In particular, the parameters µ and V in our model are equivalent to pressure p and inverse density of states 1/D 0 in the jamming transition model correspondingly.
In the paper we mainly considered a case of strong disorder, µ V 2 . Taking into account Eq. (64) we find that the small parameter of our model coincides with the small parameter ∆φ in the jamming transition model. The mean free path at the crossover as follows from (58) and (64) is given by what also coincides with 32 . It would be very interesting to investigate physical reasons for this striking "mapping" of two models to each other in more details in a future work.

VI. DISCUSSION
We have developed a stable random matrix approach to describe vibrations in strongly disordered systems, which have properties similar to what one observes in granular matter at the jamming transition point, in jammed systems and, finally, in real glasses. This approach has one important advantage in comparison to other models. It describes mechanical systems which are always stable independently of the degree of disorder. Previous random matrix models 13,34,36 suffer from an inherent mechanical instability that occurs at some critical amount of disorder. As a result they are limited by consideration of "relatively weak" or "moderate" disorder. We use scalar model and take the dynamical matrix in the form M = AA T + µM 0 . Here A is a random matrix N × N built on a simple cubic lattice with N particles and interaction between nearest neighbors only. The only non zero non-diagonal matrix elements A ij between the nearest neighbors are taken as independent random numbers from Gaussian distribution with zero mean A ij = 0 and unit variance A 2 ij = V 2 = 1. The variance controls the degree of disorder in the lattice. To ensure the important property (2) the diagonal elements are calculated as a minus sum of non-diagonal elements A ii = − j =i A ji . M 0 is a crystalline dynamical matrix with unit springs between the nearest neighbors. As a result each particle in this lattice is connected by random elastic springs with 24 surrounding particles. Since matrix AA T is always positive definite, such form of the dynamical matrix guarantees the mechanical stability of the system for any positive value of µ.
If the first term AA T is responsible for the disorder in the system, the second term µM 0 describes the ordered part of the Hamiltonian. The parameter µ controls the relative amplitude of this part and the rigidity of the lattice. It can vary in the interval 0 µ < ∞, changing the rigidity and relative amount of disorder. In this paper we have mainly considered the case of strong and moderate disorder when 0 µ V 2 and fluctuating part of the dynamical matrix is bigger then the ordered part. In this case the Young modulus of the lattice E ∝ V √ µ. The parameter µ plays the same role as pressure in jammed systems.
We have found that the delocalized vibrational excitations in this disordered lattice are of two types. At low frequencies below the Ioffe-Regel crossover, ω < ω IR , they are the usual phonons (plane waves) which can be characterized by frequency ω and wave vector q. However, with increasing of ω, due to the disorder-induced scattering, the phonon line width ∆ω increases rapidly as ∆ω ∝ ω 4 and at some frequency ω ≈ ω IR the phonon mean free path l becomes of the order of the wave length λ. Though this crossover is not sharp and has no critical behavior at ω = ω IR , the structure of the eigenmodes at higher frequencies quite soon become very different from the plane waves.
As a result, at higher frequencies the original notion of phonons is lost and delocalized vibrational modes have a diffusive nature. They are similar to diffusons introduced by Allen and Feldman, et al. 26 . The diffusons again can be characterized by frequency ω, but have no well defined wave vector q. Above ω ≈ ω IR the structure factor of particle displacements S(q, ω) becomes very similar to the structure factor S rw (q, ω) of a random walk on the lattice. The former has a broad maximum as a function of q at q = ω/D u , where D u V is a diffusion coefficient of the particle displacements.
The displacement structure factor S(q, t) in the diffuson range, for small q 1/a 0 , decays as following, S(q, t) ∝ exp(−D u q 2 t). As a result the vibrational line width Γ(q) = D u q 2 . Such quadratic dependence of Γ(q) was found in many glasses in the experiments on inelastic x-ray scattering, see for example 71,72 and references therein. It was also found in molecular dynamic simulation of amorphous silicon 73 . However in these and other papers this line width was attributed to phonons without discussion of its physical origin. We guess that the observed q 2 dependence of Γ(q) has nothing to do with phonons and is in fact related to diffusons. However, a more detailed investigation is necessary for a definite conclusion.
The crossover between phonons and diffusons takes place at the Ioffe-Regel crossover frequency ω IR which is close to the position of the boson peak. Since for phonons ∆ω ∝ ω 4 and for diffusons Γ(q) = D u q 2 , there should exist a crossover from ω 4 to q 2 dependence of the line width. Such a crossover was indeed found recently in inelastic x-ray scattering in lithium diborate glass 62 , densified vitreous silica 74 , vitreous silica [75][76][77] , glassy sorbitol 78 and glycerol glass 79 . The crossover frequency was found to be close to the BP position.
As a result, if our guess is true, we can calculate the diffusion coefficient of particle displacements, D u = Γ(q)/q 2 , from the experimental line width Γ(q) in the range, where it is proportional to q 2 . Taking into account that D u ≈ a 2 0 /τ where a 0 is the lattice constant and τ is an average time for a jump, we come to the order of the value estimate D u ≈ 1 mm 2 /sec for a 0 ≈ 2Å and τ ≈ 0.4 × 10 −13 sec. Let us compare this value with experimental data.
In the paper 76 it was found that in vitreous silica Γ/( ω) 2 = 0.07 meV −1 for q ≥ 2 nm −1 . Taking the sound velocity v L = 5250 m sec −1 for q = 2 nm −1 we get for diffusion coefficient D u = 1.3 mm 2 /sec. Let us compare this value with the diffusivity of energy D(ω) for small ω in the same glass. We expect that both coefficients should be of the same order of magnitude. The diffusivity of energy D(ω) in vitreous silica was calculated in the paper 29 . It was obtained that D(0) = 1.4 mm 2 /sec. A close estimate D(0) = 1.1 mm 2 /sec was given in 30 . As one can see the agreement between D u and D(0) is unexpectedly good. In glycerol glass 80 we found the diffusivity about factor of two smaller, D u = 0.46 mm 2 /sec. For amorphous silicon from molecular dynamic calculations 73 we get D ul = 3.2 mm 2 /sec for longitudinal vibrations, and D ut = 1.2 mm 2 /sec for transverse vibrations. For the diffusivity of energy we have in this glass the estimate 26 D(0) = 0.6 mm 2 /sec.
Since ω IR ∝ √ µ (and independent of the strength of disorder V ), we can vary the Ioffe-Regel crossover frequency and, therefore, the relative number of phonons N ph in the system, changing the parameter µ. It is zero when µ = 0 and there are no phonons in the lattice.
In this case all delocalized vibrations are diffusons. If 0 < µ 1 we have phonons, but their relative number is small. One can show that in this case N ph ∝ µ 3/4 . In the opposite case, µ 1, the disorder is relatively small and nearly all vibrations in the lattice are well defined plane waves, i.e. phonons.
In amorphous silicon the relative number of phonons (plane waves) was estimated to be only 4% from all of the vibrational modes in the system 26 . The estimates show that in our model we have such a small amount of propagating modes, as in a-Si, for µ ≈ 0.1. In the silica glass we can estimate the relative number of phonons from the data 11 . Taking into account that Ioffe-Regel crossover frequency in amorphous silica was estimated to be 11 ν IR = 1 THz, and integrating density of states 11 up to this frequency we come to the relative number N ph = 0.002 ± 0.0005. As a result in the typical glass such as amorphous silica only 0.2 % of all modes are phonons. As follows from Table I it corresponds to very small values of µ < 0.01. It means that small amount of phonons in disordered systems is a signature of strong disorder.
Usually the phenomenon of diffusion takes place for conserved quantities. In our system we have two integrals of motion. They are the momentum and the energy of the lattice. Therefore, first of all, one has to discriminate the diffusion of particle momentums (or particle displacements) from the diffusion of energy. Conservation of displacement is related to conservation of the center of inertia in the system. As a result, the diffusion of particle displacements has the same diffusion coefficient as the diffusion of particle momentums.
The diffusion coefficient of displacements/momentums D u/v is hidden in the displacement structure factor S(q, ω) (20). Comparing this structure factor with the structure factor of the random walk on the lattice, we found that for the case of µ = 0 the diffusion coefficient D u/v = D rw = 0.7. We can check that it is indeed the diffusion coefficient of particle displacements/momentums in a similar way we used for finding the diffusivity of energy D(ω) in Section IV B.
Let us consider a cubic random lattice L × L × L with µ = 0 and unit masses m i = 1 with periodic boundary conditions. At initial moment t = 0 let us displace all particles in a thin layer around the central layer (with coordinate x = 0) according to Gaussian distribution Here the thickness of the layer x 0 should be small enough in comparison to the sample size L, i.e. x 0 L/2. Initial velocitiesu(0) of all the particles are equal to zero.
After initial displacements in the thin central layer, the particle displacements will diffuse to the left and to the right ends of the sample. Solving numerically the Newton equations, we find the average squared distance to the displacement diffusion front, similar to Eq. (40) Since the center of inertia does not move, the total displacement of all particles u tot is independent of time and equal to the total displacement at t = 0. From the slope of R 2 u (t) we can calculate the diffusion coefficient of the displacements D u as follows similar to Eq. (42).
In the same way we can calculate the diffusion of momentum. For that at the moment t = 0 initial displacements of all the particles we put equal to zero. However initial velocities v =u(0) in the thin central layer we take distributed similar to Eq. (67) Then, as in the previous case, solving numerically the Newton equations we find Since the total momentum is conserved, v tot is also independent of time and equal to its initial value at t = 0. From the slope of R 2 v (t) we can calculate the diffusion coefficient of the momentum D v using one dimensional equation similar to Eq. (69).
In both cases we have obtained for diffusion coefficients D u and D v the same value as was derived from the structure factor, D u ≈ D v ≈ D rw = 0.7. It confirms our statement that the displacement structure factor S(q, ω) gives us the information about diffusion of particle displacements (or momentums). The diffusion of momentum is usually related to viscosity η of the medium. Therefore in the case of µ = 0 our lattice has no rigidity but has a finite value of viscosity.
In disordered lattices the diffusion of energy is different from the diffusion of particle displacements (momentums). In the harmonic approximation the eigenmodes with different frequencies do not interact with each other. Therefore the energy cannot be transferred from one eigenmode to other eigenmodes. It means that energy of every eigenmode E(ω i ) is conserved (with time). The total energy E tot is just a sum of these eigenmode contributions As a result, instead of one integral of motion (the total energy E tot ), in a scalar harmonic system with N particles we have N integrals of motion E(ω i ). And for each frequency ω i we have its own unique energy diffusivity D(ω i ). At this point our model decidedly confirms the physical picture suggested in papers 22-26 for amorphous silicon. We believe that it can be applied to some other glasses as well.
Usually this diffusivity is hidden in a displacement/momentum structure factor of the 4-th order. However, we calculated the diffusivity of energy D(ω) in a different way using two different approaches as it was discussed in Section IV B. The first approach is based on the direct solution of Newton equations. In the second approach we calculated the diffusivity using Edwards and Thouless formula 69 . Both approaches give the same result.
In the first approach we used a short external force pulse ∆t exciting vibrations in a small space region of the lattice and in a small frequency interval ∆ω ≈ 1/∆t near frequency ω. Then on a time scale t ∆t the energy diffused through the lattice. Using Newton equations of motion we calculated this diffusion directly. It was supposed that the interval ∆ω is much bigger than the interlevel spacing δω and therefore the former consists of many eigenmodes. In the thermodynamic limit δω ∝ 1/N → 0 if N → ∞. Therefore in an infinite system we can take the interval ∆ω arbitrary small. The energy diffusion coefficient D(ω) in this case is a function of frequency ω. Approaching the localization threshold ω loc the diffusivity D(ω) should go to zero.
We applied this method for µ = 0, when there are no phonons in the lattice. In this case we obtained for diffusivity at zero frequency D(0) ≈ 0.4, i.e. the value about factor of two smaller than for diffusivity of displacements, D u . However this approach is rather difficult to implement for computer simulations in the case when µ = 0. In this case we have phonons in the lattice with long mean free paths. And samples with much bigger sizes are necessary.
Therefore, to calculate the diffusivity D(ω) for arbitrary value of µ (including the case of µ = 0), we used another approach. In this approach Edwards and Thouless formula 69 , D(ω i ) = cL 2 |∆ω i |, was used. It relates the diffusivity D(ω i ) with shift of the eigenfrequencies ∆ω i due to change of the boundary conditions in one direction. The proportionality coefficient c we found from the comparison with the Newton method for µ = 0. In this case both methods result in the same frequency dependence of D(ω).
Here g(ω) is density of states and C(ω, T ) is specific heat of harmonic oscillator Localized modes have D(ω i ) = 0 and make no contribution to κ(T ). If functions g(ω) and D(ω) are approximately constant in some frequency interval (the case that we have, for example, in our picture for ω > ω IR ), then we find from Eq. 74 that approximately κ(T ) ∝ T in the corresponding temperature range 31 . It explains a quasilinear temperature dependence of the thermal conductivity above the plateau observed in glasses 6 . With increasing frequency the functions g(ω) and D(ω) finally drop to zero and thermal conductivity saturates at some constant level independent of temperature. Thus the conception of diffusons gives clear explanation for the temperature dependence of the thermal conductivity of glasses and other disordered systems.
Summarizing, using a stable random matrix approach we have presented a consequent theory of vibrational properties in strongly disordered systems. In these systems a relative amount of phonons is small and almost all delocalized vibrations are diffusons. The diffusons play an important role and are responsible for the transport properties of glasses at higher temperatures. Presumably they are also accounted for the mysterious q 2 dependence of the vibrational line width Γ(q) observed in many experiments on inelastic x-ray scattering in glasses. Therefore we think that it is necessary to take them into account in interpretation of experimental data.

VII. ACKNOWLEDGMENTS
We are very grateful to V. L. Gurevich and Anne Tanguy for many stimulating discussions and gratefully acknowledge interesting discussions with B. Rufflé and E. Courtens as well. One of the authors (DAP) thanks the University Lyon 1 for hospitality. This work was supported by St. Petersburg Government (diploma project no. 2.4/29-06/143C), Dynasty Foundation, RF President Grant "Leading Scientific Schools" NSh-5442.2012.2 and Russian Ministry of Education and Science (contract N 14.740.11.0892). Consider here the case when some part of springs µ are cut out from the matrix µM 0 in dynamical matrix (13). The value of parameter µ = 1 we will keep fixed. Let parameter p gives the percentage of remaining springs. The percolation threshold in the simple cubic lattice for bond percolation problem is at p c ≈ 25% 81 . If p < p c , then there is no infinite cluster of connected springs and therefore matrix µM 0 with cut out springs itself has no acoustical phonon-like modes at all. Nevertheless, the full dynamical matrix (13) still has well defined phonon modes with density of states ∝ ω 2 for all positive values of p even below the percolation threshold. The normalized density of states g(ω) for µ = 1 and different values of p is shown on Fig. 21. The straight lines show the phonon contribution to the DOS calculated from Eq. (18) with sound velocity given by Eq. (17). The Young modulus E was calculated numerically using Eq. (14) for the lattice with N = 10 6 particles (one realisation) in the same way as it was done in Section III. The details of these calculations will be published elsewhere.

B. Superposition of two random matrices
Another (less obvious) possibility to get phonons is to add to the random dynamical matrix AA T a random matrix βBB T . Here β is a parameter and the random matrix B is build in the same way as random matrix A but they are statistically independent from each other. Though both terms AA T and βBB T taken separately have zero rigidity (and do not have phonons) their superposition introduces a finite rigidity E to the system. The rigidity changes when we vary parameter β as E ∝ √ β and goes to zero when β → 0. So the scaling relations in this case for β V 2 are the same as in Section V with replacement of µ by β. The preliminary results obtained within this approach are shown on Fig. 22. Further details will be published elsewhere.
These two examples show clearly, that appearance of phonons in the system is not related to the crystalline order in the term µM 0 . The issue is more complicated. We are going to discuss this problem in more details elsewhere.

C. Displacement structure factor
Let us consider the displacement structure factor given by Eq. (20) We will assume that initial velocities of all particles at t = 0 are zero. Then the displacement of i-th particle u(r i , t) as a function of time can be written in the form a j e i (ω j ) cos(ω j t).
Here e i (ω j ) -is eigenvector of the dynamical matrix M corresponding to i-th particle and eigenfrequency ω j . The eigenvectors satisfy equations N j=1 M ij e j (ω k ) = ω 2 k e i (ω k ).
Using (79), one can write the coefficients a j in (77) in terms of the particle displacements for t = 0 The initial displacements u(r i , 0) are independent Gaussian random variables with zero mean and unit variance u(r i , 0) = 0, u(r i , 0)u(r j , 0) = δ ij .
Basing on this equation and making use of (80) and of (79) one can prove that the coefficients a j are also independent random Gaussian variables a j = 0, a i a j = δ ij .
Using this property, one can evaluate the average (76) as