Présentée Par

Beltiukov Iaroslav

M M M M M M Mme

Jean-Louis Barrat

Sergey Skipetrov

Anne Tanguy

Benoit Rufflé

Vladimir Lorman

Veniamin Kozub

Michel Dyakonov

Examinateur Rapporteur Rapporteur

Matrices aléatoires et propriétés vibrationnelles de solides amorphes dans le domaine terahertz

Keywords: amorphous solids, vibrations, random matrices

It is well known that various amorphous solids have many universal properties. One of them is the temperature dependence of the thermal conductivity. However, the microscopic mechanism of the heat transfer above 20 K is still poorly understood. Recent numerical simulations of amorphous silicon and silica show that vibrational modes in the corresponding frequency range (called "diffusons") are delocalized, however they are completely different from low-frequency acoustic phonons.

In this work we present a stable random matrix model of an amorphous solid. In this model one can vary the strength of disorder going from a perfect crystal to extremely disordered soft medium without macroscopic rigidity. We show that real amorphous solids are close to the second limiting case, and that diffusons occupy the dominant part of the vibrational spectrum. The crossover frequency between acoustic phonons and diffusons is determined by the Ioffe-Regel criterion. Interestingly, this crossover frequency practically coincides with the boson peak position. We also show that, as a function of frequency, the diffusivity and the vibrational density of states of diffusons are practically constant. As a result, the thermal conductivity is a linear function of temperature up to rather high temperatures and then saturates. This conclusion is in agreement with numerous experimental data.

Further, we consider a numerical model of amorphous silicon-like materials and investigate the role of disorder for longitudinal and transverse vibrations. We also show that the random matrix theory can be successfully applied to estimate the vibrational density of states of granular jammed systems.

Résumé

Il est bien connu que divers solides amorphes ont de nombreuses propriétés universelles. L'une d'entre elles est la variation de la conductivité thermique en fonction de la température. Cependant, le mécanisme microscopique du transfert de chaleur dans le domaine de température supérieure à 20 K est encore mal compris. Simulations numériques récentes du silicium et de la silice amorphes montrent que les modes de vibration dans la gamme de fréquences correspondante (dits diffusons ) sont délocalisés. En même temps ils sont complètement différents des phonons acoustiques de basse fréquence.

Dans ce travail, nous présentons un modèle stable de matrice aléatoire d'un solide amorphe. Dans ce modèle, on peut faire varier le degré de désordre allant du cristal parfait jusqu'au milieu mou extrêmement désordonné sans rigidité macroscopique. Nous montrons que les solides amorphes réels sont proches du deuxième cas limite, et que les diffusons occupent la partie dominante du spectre de vibration. La fréquence de transition entre les phonons acoustiques et diffusons est déterminée par le critère Ioffe-Regel. Fait intéressant, cette fréquence de transition coïncide pratiquement avec la position du pic boson. Nous montrons également que la diffusivité et la densité d'états de vibration de diffusons sont pratiquement constantes en fonction de la fréquence. Par conséquent, la conductivité thermique est une fonction linéaire de la température dans le domaine allant à des températures relativement élevées, puis elle sature. Cette conclusion est en accord avec de nombreuses données expérimentales.

En outre, nous considérons un modèle numérique de matériaux de type de silicium amorphe et étudions le rôle du désordre pour les vibrations longitudinales et transverses. Nous montrons aussi que la théorie des matrices aléatoires peut être appliquée avec succès pour estimer la densité d'états vibrationnels des systèmes granulaires bloqués.

Mots clés solides amorphes, vibrations, matrices aléatoires.

iii Establishing the general properties of vibrations in amorphous solids (glasses) is one of the key problems in the physics of disordered systems. Amorphous solids have approximately the same local order of atoms as in crystals, but there is no long-range order (Fig. 1). The disorder in atomic positions strongly modifies the macroscopic properties such as the thermal conductivity [START_REF] Eucken | Über die Temperaturabhängigkeit der Wärmeleitfähigkeit fester Nichtmetalle[END_REF][START_REF] Berman | Thermal conductivity of glasses at low temperatures[END_REF]. Fig. 2 shows the thermal conductivity of amorphous and crystalline SiO 2 .

The thermal conductivity of dielectrics is completely determined by transport properties of vibrations. The thermal conductivity of a phonon gas was first obtained by Debye [START_REF] Debye | Zustandsgleichung und Quantenhypothese mit einem Anhang über Wärmeleitung. Vorträge über die kinetische Theorie der Materie und der Elektrizität[END_REF][START_REF] Kittel | Interpretation of the thermal conductivity of glasses[END_REF]

κ = 1 3 Cvl, ( 1
)
where C is the specific heat, v is the sound velocity and l is the mean free path of phonons. The thermal conductivity of nearly perfect crystals at high temperatures Filled symbols are the first results of Eucken [START_REF] Eucken | Über die Temperaturabhängigkeit der Wärmeleitfähigkeit fester Nichtmetalle[END_REF]]. Open symbols are more recent Pohl's results [START_REF] Zeller | Thermal conductivity and specific heat of noncrystalline solids[END_REF][START_REF] Raychaudhuri | Thermal conductivity of neutronirradiated silica[END_REF][START_REF] Vandersande | Simple apparatus for the measurement of thermal diffusivity between 80-500 K using the modified Ångström method[END_REF][START_REF] Cahill | Thermal conductivity of amorphous solids above the plateau[END_REF][START_REF] Cahill | Lattice vibrations and heat transport in crystals and glasses[END_REF]. The black bottom line is the existing theory of the thermal conductivity of amorphous solids [START_REF] Buchenau | Interaction of soft modes and sound waves in glasses[END_REF]]. The straight red line shows the linear dependence (it will be obtained in Section 2.4).

is κ ∝ 1/T (Fig. 2), which is related to umklapp processes (resulting from a small anharmonicity). At small temperatures (much less than the Debye temperature) the probability of umklapp processes is exponentially small. In this case, the mean free path l is limited by the sample size L. According to the Debye law, the lowtemperature specific heat is C ∝ T 3 [START_REF] Debye | Zur Theorie der spezifischen Wärmen[END_REF]]. Therefore, the low-temperature thermal conductivity of crystals is κ ∝ T 3 (Fig. 2).

Amorphous solids have bond lengths and bond angles, which differ by ∼10% from the crystalline values. However, the temperature dependence of the thermal conductivity is completely different, and the difference can be more than four orders of magnitude (Fig. 2). Other amorphous dielectrics have qualitatively the same temperature dependence of the thermal conductivity (Fig. 3).

At low temperatures below 1 K, the low-frequency long-wave acoustical phonons are well-defined excitations which transfer the heat in glasses. At these temperatures Comparison of the thermal conductivity of different amorphous dielectrics: GeO 2 , Se, SiO 2 , solid varnish GE 7031, polymer PMMA [START_REF] Zeller | Thermal conductivity and specific heat of noncrystalline solids[END_REF].

the thermal conductivity κ(T ) ∝ T 2 is controlled by a resonant scattering of phonons by two-level systems (TLS) [START_REF] Hunklinger | Thermal and elastic anomalies in glasses at low temperatures[END_REF][START_REF] Phillips | Two-level states in glasses[END_REF]. Between 4 K and 20 K the thermal conductivity κ(T ) saturates and displays a wellknown plateau [START_REF] Zeller | Thermal conductivity and specific heat of noncrystalline solids[END_REF]. As was shown in [START_REF] Buchenau | Interaction of soft modes and sound waves in glasses[END_REF]], it can be explained by a 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 [START_REF] Parshin | Soft-potential model and universal properties of glasses (review)[END_REF]]. Above approximately 20 K the thermal conductivity rises again (approximately linearly, κ ∝ T ) and finally saturates at the level of one order of magnitude higher, at temperatures about several hundred Kelvin [START_REF] Cahill | Thermal conductivity of amorphous solids above the plateau[END_REF].

A microscopic mechanism of the heat transfer in glasses in the temperature range above the plateau (from 20 K up to the glass transition temperature T g ∼ 1400 K [START_REF] Brückner | Properties and structure of vitreous silica. I[END_REF]) is still poorly understood. As generally believed, the origin of this second rise of the thermal conductivity (above the plateau) is not related to acoustic phonons. However, the existing models (Einstein model and the model of the minimum thermal conductivity) do not follow from the microscopic description of atomic vibrations, and they have no κ ∝ T dependence [START_REF] Cahill | Lattice vibrations and heat transport in crystals and glasses[END_REF].

It was established long ago [START_REF] Birch | The thermal conductivity of rocks and its dependence upon temperature and composition[END_REF][START_REF] Kittel | Interpretation of the thermal conductivity of glasses[END_REF][START_REF] Graebner | Phonon localization in glasses[END_REF], that in the temperature (frequency) range under consideration the mean free path of acoustic phonons l becomes of the order of their wavelength λ (or Figure 4. The numerical calculations of the VDOS of the amorphous SiO 2 using molecular dynamics methods [START_REF] Jin | Dynamic structure factor and vibrational properties of SiO 2 glass[END_REF].

even smaller, of the order of interatomic distance). Correspondingly, the Ioffe-Regel criterion for phonons [START_REF] Ioffe | Non-crystalline, amorphous and liquid electronic semiconductors[END_REF] becomes violated. Molecular dynamics calculations confirmed the existence of such crossover for some real and model glasses [START_REF] Taraskin | Ioffe-Regel crossover for plane-wave vibrational excitations in vitreous silica[END_REF][START_REF] Schober | Vibrations and relaxations in a soft sphere glass: boson peak and structure factors[END_REF]] and disordered lattices [START_REF] Schirmacher | Harmonic vibrational excitations in disordered solids and the "boson peak[END_REF]Taraskin and Elliott 2002b].

In the regime of such strong scattering, the concept of plane waves (phonons) 1with 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 [START_REF] Jin | Dynamic structure factor and vibrational properties of SiO 2 glass[END_REF][START_REF] Oligschleger | Dynamics of SiO 2 glasses[END_REF][START_REF] Taraskin | Nature of vibrational excitations in vitreous silica[END_REF]. The vibrational density of states (VDOS) of amorphous SiO 2 in the corresponding frequency range is approximately constant (Fig. 4). The same behavior for g(ω) was found in the soft-sphere glass [START_REF] Schober | Low-frequency vibrations and relaxations in glasses[END_REF][START_REF] Schober | Low-frequency vibrations in a model glass[END_REF] and amorphous Se [START_REF] Oligschleger | Dynamics of Se glasses[END_REF][START_REF] Oligschleger | Calculation of vibrational properties of selenium[END_REF]. In other glasses, the VDOS has a broad maximum and then decays to zero [START_REF] Hafner | Propagating and localized vibrational modes in Ni-Zr glasses[END_REF][START_REF] Meshkov | Low-frequency dynamics of Lennard-Jones glasses[END_REF]; [START_REF] Ballone | Embedded-atom model of glass-forming Si-metal alloys[END_REF][START_REF] Abraham | Vibrational dynamics and boson peak in a supercooled polydisperse liquid[END_REF].

At the same time, delocalized vibrations in glasses of a new type, different from plane wave-like phonons, were introduced. They were called diffusons [START_REF] Allen | Thermal conductivity of glasses: theory and application to amorphous Si[END_REF][START_REF] Allen | Thermal conductivity of disordered harmonic solids[END_REF][START_REF] Feldman | Thermal conductivity and localization in glasses: Numerical study of a model of amorphous silicon[END_REF][START_REF] Feldman | Numerical study of lowfrequency vibrations in amorphous silicon[END_REF][START_REF] Allen | Diffusons, locons and propagons: character of atomic vibrations in amorphous Si[END_REF]]. 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 [START_REF] Allen | Diffusons, locons and propagons: character of atomic vibrations in amorphous Si[END_REF]. In these papers, the hypothesis was put forward 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.

A similar conclusion was made by the authors of [START_REF] Sheng | Heat conductivity of amorphous solids: simulation results on model structures[END_REF][START_REF] Sheng | Phonon transport in strong-scattering media[END_REF]. 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 [START_REF] Feldman | Realistic model calculations based on the Kubo theory for the thermal conductivity of amorphous insulators[END_REF][START_REF] Yu | Thermal conductivity computed for vitreous silica and methyl-doped silica above the plateau[END_REF] using molecular dynamics methods. The difference between electron and phonon localization is related to a very strong scattering of electrons by perturbations in electrostatic potentials associated, say, with structural disorder [START_REF] Taraskin | Low-frequency vibrational excitations in vitreous silica: the Ioffe-Regel limit[END_REF].

frequency ω ir can vary. It is shifted to zero when the system approaches a so-called jamming transition point, and rigidity goes to zero.

Therefore, as we believe, it is important to study properties of diffusons and the Ioffe-Regel criterion systematically. 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. This work is organized as follows. In Chapter 1 we present a stable random matrix model of an amorphous solid. Starting from general properties of the dynamical matrix, we obtain a disordered medium with random bond strength between atoms on a simple cubic lattice. In this model, we can vary the strength of disorder between two limiting cases. One limit corresponds to extremely disordered soft medium without macroscopic rigidity. In the other limit one has a perfect crystal. We also compare our results with classical results of the random matrix theory.

In Chapter 2 we analyze transport properties of vibrations in the suggested model amorphous solid. We numerically calculate different physical quantities such as the dynamical structure factor, the Ioffe-Regel crossover and the diffusivity of the vibrational energy. As a result, we obtain the thermal conductivity, which is in agreement with the experimental data. We also obtain general scaling relations and compare them with ones for jammed granular systems.

In Chapter 3 we consider a Stillinger-Weber model of an amorphous silicon. We show that transverse and longitudinal vibrations in silicon-like amorphous materials have sufficiently different properties (the VDOS, the dynamical structure factor, the diffusivity and other related quantities). To emphasize this difference we study the effect of the local bending rigidity, which mostly affect the transverse vibrations.

In Chapter 4 we show that the random matrix theory can be successfully applied to estimate the vibrational density of states. As an example, we consider jammed granular systems near the jamming transition.

In the Appendices, we present some important technical information. In Appendix A we present an original derivation of some classical results of the random matrix theory, including the density of states and the level statistics of the Wishart ensemble. In Appendix B we consider the kernel polynomial method. It is a very efficient numerical method for eigenvalue analysis, which can be applied to dynamical matrices. In Appendix C we present a general method of the decomposition of the VDOS to its transverse and longitudinal components.

Chapter 1

The random matrix approach

In this Chapter we consider disordered lattices with a strong force-constant disorder, described by a stable positive definite random dynamical matrix AA T having positive eigenvalues only. We show that the vibrational density of states g(ω)

is not zero at ω = 0 and phonons cannot propagate through the lattice. We explain this by the fact that the system is extremely soft, and the macroscopic rigidity is zero. The participation ratio P (ω) indicates that all modes with an exception of a high-frequency part are delocalized. Further investigation shows that all of them are diffusons. In Section 1.8 we introduce slightly additively deformed dynamical matrix AA T + µM 0 which has phonon-like excitations at small frequencies. Here the positive definite matrix M 0 (random or non-random) is independent of A and µ is a parameter of the model which can vary in the interval 0 µ < ∞.

Introduction

Masses of atomic nuclei are much larger than electron masses, so we can separate the electron motion and the motion of nuclei. Therefore, motion of atoms can be described by the classical Newton's equation of motion

m i riα = - ∂U ∂r iα (1.1)
where m i is the mass of ith atom, α indicates spatial direction and U is the total potential energy, which depends on the atomic positions r 1 , . . . , r N . In the solid state each atom vibrates around a certain equilibrium position R i . In this case, we can linearize the equation of motion (1.1)

üiα = - jβ M iα,jβ u jβ (1.2)
with atomic displacements u i = √ m i (r i -R i ) and the dynamical matrix

M iα,jβ = 1 √ m i m j ∂ 2 U ∂r iα ∂r jβ
.

(1.

3)

The linear equation (1.2) corresponds to the eigenvalue problem

ω 2 u iα = jβ
M iα,jβ u jβ (1.4) where ε = ω 2 are eigenvalues of the dynamical matrix M .

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 a 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 can say that dynamical matrix of an amorphous system has random elements with some non-trivial correlations between them. The main restriction is the mechanical stability of the whole system. All eigenfrequencies ω in such a system are real, and the dynamical matrix M is positive definite (it is also real and symmetric by the definition (1.4)). However, not every random symmetrical dynamical matrix is positive definite. We will solve this problem using a following mathematical approach.

Every real symmetric and positive definite matrix M can be presented in the form [START_REF] Bhatia | Positive definite matrices[END_REF]]

M = AA T (1.5)
where A is some real (not necessarily symmetric nor square) matrix. And, vice versa, for every real matrix A the product AA T is a positive definite symmetric matrix [START_REF] Bhatia | Positive definite matrices[END_REF]]. One may assume that in amorphous solids the dynamical matrix has a form M = AA T , where correlations between elements of the matrix A are less important than correlations between elements of the matrix M .

A vectorial character of vibrations in real glasses makes the issue be more complicated, so in this and in the next Chapters we will use a so-called scalar model, where we for simplicity omit indices α and β. One can imagine that atoms can vibrate along x direction only, and polarization of the modes is of no importance. Different scalar models were successfully used in glassy physics in the past [START_REF] Schirmacher | Harmonic vibrational excitations in disordered solids and the "boson peak[END_REF][START_REF] Martín-Mayor | Dynamical structure factor in disordered systems[END_REF][START_REF] Grigera | Vibrations in glasses and Euclidean random matrix theory[END_REF][START_REF] Kantelhardt | Excess modes in the vibrational spectrum of disordered systems and the boson peak[END_REF]]. In this model, the dynamical matrix M is N × N matrix and the matrix A is N × K rectangular matrix.

Let us start the analysis from a simple case where all elements of the N × K matrix A (they have dimension of s -1 ) are independent and identically distributed random numbers with

A ij = 0 and A 2 ij = Ω 2 . (1.6)
It is so-called Wishart random matrix ensemble [START_REF] Wishart | The generalized product moment distribution in samples from a normal multivariate population[END_REF]]. In the limit

N, K → ∞, K N = const (1.7)
eigenvalues ε i of the matrix M = AA T have the Marĉhenko-Pastur distribution [START_REF] Marĉhenko | Distribution of eigenvalues for some sets of random matrices[END_REF] (see also Appendix A). The Wishart ensemble was investigated in the theory of financial markets [START_REF] Plerou | Random matrix approach to cross correlations in financial data[END_REF], complex networks [START_REF] Barthélemy | Large scale cross-correlations in Internet traffic[END_REF], and wireless communications [START_REF] Tulino | Random matrix theory and wireless communications[END_REF]. As far as we know, for vibrations in disordered solids this approach was not used so far (as an exception see the paper [START_REF] Gurarie | Bosonic excitations in random media[END_REF]). The Marĉhenko-Pastur distribution corresponds to the following vibrational density of states (VDOS)

g(ω) = 1 πω 2 0 ω (ω 2 -ω 2 -)(ω 2 + -ω 2 ), ω -< ω < ω + (1.8)
where

ω ± = ω 0 K/N ± 1 , ω 0 = √ N Ω.
(1.9)

If K < N there is N -K zero eigenvalues and the Marchenko-Pastur law (1.8) contains a delta-function term (1 -K/N )δ(ω). In Appendix A we present a rather simple derivation of Eq. (1.8). If K = N then the VDOS has a quarter circle form

g(ω) = 1 πω 2 0 4ω 2 0 -ω 2 , 0 < ω < 2ω 0 , (1.10)
which is similar to the well known Wigner semicircle [START_REF] Wigner | Characteristic vectors of bordered matrices with infinite dimensions[END_REF]]. This VDOS is approximately constant in a wide frequency range. Real and model amorphous systems also have a constant VDOS in some frequency range [START_REF] Jin | Dynamic structure factor and vibrational properties of SiO 2 glass[END_REF][START_REF] Schober | Low-frequency vibrations in a model glass[END_REF][START_REF] Oligschleger | Calculation of vibrational properties of selenium[END_REF][START_REF] Hafner | Propagating and localized vibrational modes in Ni-Zr glasses[END_REF][START_REF] Meshkov | Low-frequency dynamics of Lennard-Jones glasses[END_REF][START_REF] Ballone | Embedded-atom model of glass-forming Si-metal alloys[END_REF][START_REF] Abraham | Vibrational dynamics and boson peak in a supercooled polydisperse liquid[END_REF]. The corresponding distribution of eigenvalues ε = ω 2 is p(ε) ∝ 1/ √ ε. This singular-like behavior was observed in [Taraskin and Elliott 2002a] and [START_REF] Huang | Localization-delocalization transition in Hessian matrices of topologically disordered systems[END_REF] (see Fig. 2 of this paper). Therefore, we will consider the case K = N below. The general case K = N will be discussed in Chapter 4.

Sparse matrices

In the random medium model described by the Wishart ensemble, each of the elements M ij of the dynamical matrix M is not zero

M ij = k A ik A jk .
( 1.11) Obviously, this corresponds to the case of long-range interaction when each atom is connected by random forces with every atom in the system. However, this model is not justified from the physical point of view. In amorphous materials, only closely spaced atoms are bonded by elastic forces. Therefore, the more real case is when the number of nonzero elements m in each row of the matrix M is small as compared to N and does not depend on N . As a result, the matrix M is sparse. Usually, such sparse matrices arise in computer calculations of atomic vibrations in amorphous solids (and liquids). For example, in the case of the short-range order for a simple cubic lattice with the interaction only between nearest neighbors and the vector character of vibrations (in three-dimensional space), we have m = 18 + 3 = 21. For other lattices we obtain m = 24+18+3 = 45 for a bcc lattice and m = 36+18+3 = 57 for an fcc lattice. In the last two cases, we took into account all interactions in the first and second coordination shells.

Therefore, we get a more real case if we consider a sparse matrix A where each row contains only n nonzero matrix elements (with n ≪ N ). Then, each row of the matrix M = AA T will have approximately m = n 2 nonzero elements. At n 2 ≪ N , this corresponds to the case of a sparse matrix M .

If the nonzero elements of the matrix A are chosen randomly and n ≫ 1 then the VDOS is also described by the Eqs. (1.8) -(1.10) with

ω 0 = √ nΩ.
(1.12)

If n ≫ 1 we can use the obtained quarter-circle VDOS (1.10) by substituting the variance nΩ 2 /N instead of the variance Ω 2 ; that is

g(ω) = N πnΩ 2 √ 4nΩ 2 -ω 2 , 0 < ω < 2Ω √ n. (1.13)
At the same time we can take n ≪ N . This means that, at N ≫ n ≫ 1 the quarter-circle distribution for the density of states g(ω) is still valid even in the case when nonzero elements occupy only a small part of the matrix A. Such form of the distribution in our model is a universal law and does not depend on the distribution density ρ 0 (a ij ), the size of the system N , and the number of nonzero elements n for sufficiently large values of n. A numerical analysis confirms that with an increase of n the density of states g(ω) actually approaches quarter-circle distribution (for n ≫ 1) and, in this case, the inequality n ≪ N is possible (see Fig. 1.1). Already for values of n of the order of 10 and large values of N , we obtain the density of vibrational states that only slightly differs from the quarter-circle distribution. In this case, the vibrational spectrum (normalized to unity) does not depend on the size of the system N .

The symmetric sparse random matrix M = AA T considered in this Section is topologically equivalent to a tree (closed to itself on the system size) so that the number m = n 2 specifies the order of branching or the coordination number of this tree. However, the random bond structure in amorphous systems (glasses) more likely corresponds to the short-range order in the atomic arrangement, topologically similar to the bond structure existing in the corresponding crystals. It is clear that topologically the crystal structure differs fundamentally from the tree structure. In the tree structure, there are no small closed loops which are present in the lattice.

In the conclusion of this section, it should be noted that the singularity in the density of states g(ω) at ω → 0 manifests itself for small values of n (see Fig. 1.1 for n = 5). Similar singularity also exists in the density of eigenvalues the sparse random Hamiltonian H [START_REF] Dyson | The dynamics of a disordered linear chain[END_REF][START_REF] Rodgers | Density of states of sparse random matrices[END_REF][START_REF] Evangelou | 1/f noise and spectral singularities in strongly disordered electronic systems[END_REF][START_REF] Evangelou | A numerical study of sparse random matrices[END_REF]]. Taking into account that this singularity was first discovered in the density of vibrational states of a disordered one-dimensional chain by Dyson [START_REF] Dyson | The dynamics of a disordered linear chain[END_REF]], it sometimes is called the Dyson singularity. It has been believed that this singularity is an indication of strong fluctuations in a random medium and related quasi-localization of modes [START_REF] Evangelou | A numerical study of sparse random matrices[END_REF]]. A simple cubic lattice.

Cubic lattice with random bonds

Our purpose in this section is to construct a simple random matrix model of an amorphous system with certain physical properties: this structure should have a given topology of bonds, and the total potential energy U should not depend on the translation of the whole system. The latter property is necessary (but not enough as we will see below) to the propagation of low-frequency acoustic phonons. It corresponds to the sum rule in the dynamical matrix (here and below we assume that all masses m i = m are equal) [Taraskin and Elliott 2002b;[START_REF] Kühn | Translational invariance in models for lowtemperature properties of glasses[END_REF] i

M ij = j M ij = 0.
(1.14) Indeed, in this case the potential energy is .15) As soon as the dynamical matrix M is fixed, the exact atomic equilibrium positions are no longer important for dynamics on a long length scales much bigger than the interatomic distances. They do not enter the dynamical matrix M . Therefore, it is reasonable to consider harmonic lattice models involving only force-constant disorder.

U = m 2 ij M ij u i u j = - m 2 i,j<i M ij (u i -u j ) 2 . ( 1
As an example, we consider a simple cubic lattice with N = L 3 atoms, a lattice constant a 0 , and random bond strength between neighbor atoms. In all other aspects our system remains random without any periodicity (except for the topology of the bonds). The atoms have coordinates (a 0 i x , a 0 i y , a 0 i z ) and each index i α can take on values from 1 to L. Let us introduce the integer index i = i x + L(i y -1) + L 2 (i z -1). Each atom in the lattice is characterized by its unique index i running from 1 to N . Let us now construct a corresponding random matrix A. The element A ij is random and nonzero if the ith and jth atoms are nearest neighbors. Non-diagonal elements A ij and A ji are statistically independent of each other (matrix A is nonsymmetric). All other elements A ij are equal to zero except the diagonal element

A ii = - j =i A ji .
(1.16)

Then according to Eq. (1.11), the Eq. ( 1.14) will also be satisfied. As a result, for the simple cubic lattice, we obtain n = 7 nonzero elements in each row and in each column of the matrix A. In the dynamical matrix M = AA T the element M ij will be nonzero if ith atom will be nearest or next nearest to the jth atom (or it is the same atom for i = j). Fig. 1.4 shows the atoms interacting with the central (black) atom.

We calculated the density of vibrational states g(ω) for this cubic lattice. The mean value of nonzero elements of the matrix A was taken equal to zero A ij = 0 and the variance was taken to be A 2 ij = Ω 2 (Gaussian distribution). The results of the numerical analysis are shown by solid line in Fig. 1.5. According to the Debye law, the VDOS of acoustic phonons in the limit ω → 0 is where v is the sound velocity. However, Fig. 1.5 shows a weak (logarithmic) singularity in the VDOS at ω → 0. So we can suppose that the system under consideration has no acoustic phonons despite the sum rule (1.14). To elucidate a spatial structure of the eigenmodes for the matrix M = AA T , we consider as an example a two-dimensional square lattice with N = 400 × 400 particles and calculate eigenvector e i (ω min ) (i = 1, 2, ..., N ) for the lowest nonzero frequency ω min in the system. The result is shown in Fig. 1.6. 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. A similar picture takes place in the 3d case. In Section 2.3.1 we will show that in the limit ω → 0 the structure factor has power-law behavior, which confirms the fractal structure of a such low-frequency mode.

g ph (ω) ∝ ω d-1 v d , (1.17)

Participation ratio

One of the most important problems in the disordered systems is the problem of modes localization. As is well known from the seminal paper of Anderson [START_REF] Anderson | Absence of diffusion in certain random lattices[END_REF]], the sufficiently strong disorder in a system leads to localization of elementary excitations. To estimate the inverse strength of mode localization one usually introduce the participation ratio. As a rule, the participation ratio is defined by the expression

P (ω) = 1 N N i=1 e 4 i (ω) , (1.18)
where e i (ω) is the ith coordinate of the eigenvector corresponding to the eigenvalue ω 2 of the dynamical matrix M . In the case of completely localized mode

|e 1 | = 1, e 2 = e 3 = ... = e N = 0, P ∼ 1 N (1.19)
the participation ratio P (ω) decreases with the increase of the system size. In the case of completely delocalized mode

|e 1 | = |e 2 | = ... = |e N | = 1 √ N , P ∼ 1, (1.20)
the participation ratio does not depend on the system size N and is of the order of unity.

We performed the numerical calculations of the participation ratio for vibrational excitations in the cubic lattice with random bonds, which was introduced in the previous Section. Fig. 1.7 shows the result for lattices in different dimensions. Threedimensional cubic lattice has P (ω) ≈ 0.2, which does not depend on N . Therefore, almost all vibrations (except a small high-frequency range) are delocalized, but none of them are plane waves. According to the terminology proposed by [START_REF] Allen | Diffusons, locons and propagons: character of atomic vibrations in amorphous Si[END_REF]] they are referred to as the diffusons. In Section 2.3 we will show that these vibrations spread by mean of diffusion. 126 2 ,25 3 ,5 6 ). The black horizontal line corresponds to the theoretical value P (ω) = 1/3. The frequency has units of Ω.

Qualitatively our plot P (ω) coincides well with the results of the numerical calculations of the participation ratio for amorphous SiO 2 using molecular dynamics methods [START_REF] Jin | Dynamic structure factor and vibrational properties of SiO 2 glass[END_REF] in the frequency range 0 < ω < 120 meV (Fig. 1.8). It is interesting that in the two-dimensional (square) and one-dimensional lattices with random bonds (constructed in a similar manner) the participation ratio is one order of magnitude smaller than that in the cubic lattice (Fig. 1.7). By analogy with disordered electronic systems [START_REF] Abrahams | Scaling theory of localization: absence of quantum diffusion in two dimensions[END_REF], this can indicate to the localization of vibrational modes in these low-dimensional structures.

The numerical values of the participation ratio of vibrational modes P (ω) in various glasses according to the data obtained by molecular dynamics methods, as a rule, are in the range 0.2 P (ω) 0.6 [START_REF] Jin | Dynamic structure factor and vibrational properties of SiO 2 glass[END_REF][START_REF] Schober | Low-frequency vibrations in a model glass[END_REF][START_REF] Schober | Low-frequency vibrations and relaxations in glasses[END_REF][START_REF] Oligschleger | Dynamics of Se glasses[END_REF][START_REF] Oligschleger | Calculation of vibrational properties of selenium[END_REF][START_REF] Hafner | Propagating and localized vibrational modes in Ni-Zr glasses[END_REF][START_REF] Meshkov | Low-frequency dynamics of Lennard-Jones glasses[END_REF][START_REF] Ballone | Embedded-atom model of glass-forming Si-metal alloys[END_REF][START_REF] Abraham | Vibrational dynamics and boson peak in a supercooled polydisperse liquid[END_REF]. This is in good agreement with the results of the random matrix theory. For example, by assuming that the eigenvectors e i (ω) (i = 1, 2, ..., N ) of the random matrix M = AA T are unit vectors in the N -dimensional space

N i=1 e 2 i (ω) = 1, (1.21)
which are isotropically oriented in all possible directions, the quantity r = e 2 i (ω) will be distributed according to the Porter-Thomas law [START_REF] Haake | Quantum signatures of chaos[END_REF]]

p(r) = N 2πr exp - N r 2 .
(1.22) The numerical calculations of the participation ratio of the amorphous SiO 2 using molecular dynamics methods [START_REF] Jin | Dynamic structure factor and vibrational properties of SiO 2 glass[END_REF]].

As a result, we have

e 2 i (ω) = r = 1 N , e 4 i (ω) = r 2 = 3 N 2 (1.23)
and according to Eq. (1.18) the participation ratio is

P (ω) = 1/3.
(1.24) Fig. 1.7 shows that the participation ratio approaches to the theoretical value 1/3 with the increase of the space dimensionality.

In the literature [START_REF] Schober | Low-frequency vibrations in a model glass[END_REF], there is one more definition of the participation ratio for the vector model (in contrast to the above scalar model)

P 3 (ω) = 1 N N i=1 3 α=1 e 2 iα (ω) 2 , (1.25)
where index i indicates the index number of the atom (i = 1, 2, ..., N ) and the index α stands for the Cartesian projection of the displacement of the ith atom onto the α (α = x, y, z) axis. In this case, under the assumption that the unit vectors e iα (ω) are isotropically distributed in the 3N -dimensional space we have in analogy with Eqs. (1.23)

e 2 iα (ω) = r = 1 3N , e 4 iα (ω) = r 2 = 3 (3N ) 2 = 1 3N 2 .
( 1.26) we find that the participation ratio P 3 (ω) (1.25) is equal P 3 (ω) = 3/5 = 0.6.

(1.28)

In this sense, the values of P (ω) = 1/3 for the scalar model and P 3 (ω) = 0.6 for the vector model are equivalent to each other from the standpoint of the random matrix theory. The value of P 3 (ω) ≈ 0.6 was obtained in numerical calculations of the participation ratio of vibrational modes for the soft-sphere glass [START_REF] Schober | Localized low-frequency vibrational modes in glasses[END_REF]. Finally, making use of Eqs. (1.26) it can be shown that the participation ratios P O ≈ P Si ≈ 0.3 calculated numerically by Jin et al. [START_REF] Jin | Dynamic structure factor and vibrational properties of SiO 2 glass[END_REF] for amorphous SiO 2 are also in good agreement with the theoretical values P O = P Si = 1/3 that follows from formula (18) of this work. Summarizing we can conclude that participation ratio calculated in different papers for different glasses is in good agreement with predictions of the random matrix theory.

Levels statistic

The level statistics is another powerful criterion that makes it possible to judge about localization or delocalization of vibrational modes. If the modes are localized, Figure 1.10. The molecular dynamics calculations of the statistics of vibrational levels for amorphous clusters [START_REF] Sarkar | Universality in the vibrational spectra of single-component amorphous clusters[END_REF]]. The line is the theoretical prediction (1.32).

their frequencies are randomly distributed over the frequency axis according to the Poisson distribution law without any correlation with each other. For quantitative description let us introduce the normalized difference between the vibrational eigenfrequencies

s = ∆ω ∆ω , (1.29)
where ∆ω is the distance between the two neighboring frequencies that corresponds to the frequency ω and ∆ω is the mean distance between these frequencies. Then for localized modes the distribution function can be presented in the form Z(s) = exp(-s).

(1.30)

When the modes are delocalized the repulsion effect takes place and Z(s) ∝ s for small values of s ≪ 1. The level statistics of the Wishart ensemble (and the Wigner ensemble as well) is approximately

Z W (s) ≈ π 2 s exp - π 4 s 2 , (1.31)
It is so called Wigner surmise [START_REF] Haake | Quantum signatures of chaos[END_REF]]. The precise form of Z(s) can be written as [Forrester and Witte 2000] 1.32) where σ(t) is the solution of the nonlinear differential equation with an asymptotic behavior

Z(s) = d 2 ds 2 exp - (πs/2) 2 0 σ(t) t dt (
(tσ ′′ ) 2 + σ ′ (σ -tσ ′ )(4σ ′ -1) - 1 4 (σ ′ ) 2 = 0 (1.
σ(t → 0) = 1 π t 1/2 + 2 π 2 t + 4 π 3 - 1 3π t 3/2 .
(1.34)

We have calculated numerically the level statistics for the dynamical matrix of the simple cubic lattice with random bonds. The level statistics for different frequencies is presented in Fig. 1.9. In a wide frequency range we have Z(s) ∝ s for small s. In the frequency range 0 < ω < 4.7, the result is indistinguishable from Eq. (1.32). As was mentioned above, this statistics corresponds to the case of delocalized modes. Therefore, we conclude that the majority of vibrational modes in our system are delocalized. This is in a good agreement with the data presented in Fig. 1.7 for the participation ratio P (ω). In conclusion, we note that our results agree well with the molecular dynamics calculations of the statistics of vibrational levels for amorphous clusters (Fig. 1.10).

Young modulus and the absence of acoustic phonons

Let us consider the difference between the VDOS (Fig. 1.5) and the Debye prediction (1.17). In the three-dimensional case and the scalar model we have

g ph (ω) = a 3 0 2π 2 ω 2 v 3 (1.35)
where a 0 is the lattice constant, v is the sound velocity. The VDOS is normalized to unity, and the coefficient a 3 0 indicates a volume associated with each degree of freedom. According to the standard textbook formula of the macroscopic elasticity theory for isotropic medium (see Eq. 5.2 in [START_REF] Landau | Theory of Elasticity[END_REF]), Young Figure 1.12. The Young modulus (in units of mΩ 2 /a 0 ) averaged over 1000 realizations for different system sizes.

modulus is given by E = σ xx /u xx . Here σ xx is the stress, and u xx is the strain. The component u xx gives the relative lengthening of the sample. The sound velocity depends on the Young modulus as

v = E ρ 0 (1.36)
where the mass density is ρ 0 = m/a 3 0 . In modeling of amorphous solids, the standard method to calculate the Young modulus is to use Irving-Kirkwood stress tensor formula [START_REF] Irving | The Statistical Mechanical Theory of Transport Processes. IV. The Equations of Hydrodynamics[END_REF]. 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 equilibrium particle positions. Therefore, to avoid these difficulties we used a direct numerical method. We apply the force f 0 to each atom on the right side and the force -f 0 to each atom on the left side (Fig. 1.11). Since Newton equations are linear, the final result is independent of the value of the force f 0 . In other two directions we use the periodic boundary conditions. In this case, the stress is .37) and the strain is

σ xx = f 0 a 2 0 , ( 1
ξ = 1 a 0 L 2 (L -1) r u r + l u l = ūr -ūl a 0 (L -1)
, (1.38) where u l and u r are the displacements of the particles on the left and the right side. Their mean values are denoted by ūr and ūl respectively. The resulting Young modulus is .39) 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.

E = σ xx u xx = f 0 (L -1) a 0 (ū r -ūl ) . ( 1
Numerical calculations show that the Young modulus strongly fluctuates: its value is sufficiently different for different realizations of the system. The mean value of the Young modulus decreases with the system size as (Fig. 1.12)

E ∝ 1 N mΩ 2 a 0 .
(1.40)

In the thermodynamical limit (N → ∞) the Young modulus is 0, i.e. the macroscopic rigidity is 0. Therefore, the sound velocity v is also 0 (Eq. (1.36)), acoustic phonons cannot propagate, and the Debye prediction (1.35) is not applicable.

The distribution of the measured Young moduli is shown in Fig. 1.13. It has an exponential form

ρ(E) = αN a 0 mΩ 2 exp - αN a 0 mΩ 2 E , α = 0.39. (1.41)
The same distribution holds not only for a cubic sample but also for any rectangular parallelepiped sample. The result for 20 × 15 × 10 sample is presented in Fig. 1.13. The relative fluctuations of the Young modulus do not decrease with the system a p(u) size

N E -E 2 E 2 = 1. (1.42)
The fluctuations of displacements of the edge atoms are much larger than the mean displacement, and about 35% of atoms have displacements, which are opposite to the force (Fig. 1.14). Therefore, we cannot apply the elasticity theory to the system under consideration because we cannot introduce a strain as a smooth function of the coordinate.

Distribution of elements of the dynamical matrix

Let us consider why the Young modulus is zero. The macroscopic rigidity depends on the microscopical atomic interactions. We have shown that each atom interacts with 24 surrounding atoms in the three-dimensional case (Fig. 1.4). For clarity, we now consider the two-dimensional analog, where each atom interacts with 12 surrounding atoms (Fig. 1.15).

In Section 1.3 we defined matrix elements of the dynamical matrix M = AA T . Matrix elements between the central atom in Fig. 1.15 and its nearest neighbors have a form

M 01 = k A 0k A 1k = A 00 A 10 + A 01 A 11 . (1.43)
The diagonal elements A are defined by the sum rule (1.14)

A 00 = -(A 10 + A 20 + A 30 + A 40 ), (1.44)
A 11 = -(A 91 + A 51 + A 61 + A 01 ). (1.45)
From Eq. (1.43), we get .46) Since averaged values A ij = 0 and different non-diagonal matrix elements A ij are statistically independent of each other, the average value M 01 is determined by the first two quadratic terms in Eq. (1.46). 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 .47) According to Gaussian distribution of A i =k , there is a probability (1/4 in two dimensions and (1 -1/ √ 6)/2 ≈ 30% in three dimensions) that the spring constant is negative k 01 < 0. All the aforesaid is valid for other nearest neighbor matrix elements M 02 , M 03 , and M 04 .

M 01 = -A 2 10 -A 2 01 -A 10 (A 20 + A 30 + A 40 ) -A 01 (A 91 + A 51 + A 61 ). ( 1
k 01 = - m a 0 M 01 = m a 0 A 2 10 + m a 0 A 2 01 = 2mΩ 2 a 0 . ( 1
The negative spring constant ("negative spring") means unstable spring between two atoms: it expands after initial stretching and shrinks after initial contraction. However, the whole system is stable because the dynamical matrix M = AA T is Figure 1.16. Distributions of random elastic spring constants (in units of mΩ 2 /a 0 ) in 3d simple cubic lattice. The parameter µ will be defined in Section 1.8. Spring constants k 05 and k 09 do not depend on µ.

positive definite and there are always positive springs, which stabilize the system. The effect of negative spring constants on atomic vibrations was discussed in different papers [START_REF] Schirmacher | Harmonic vibrational excitations in disordered solids and the "boson peak[END_REF][START_REF] Dederichs | Resonance modes of interstitial atoms in fcc metals[END_REF][START_REF] Oshima | Surface phonons of the superconducting materials NbC(100) and TaC(100)[END_REF][START_REF] Erwin | Frenkel-Kontorova model of vacancy-Line interactions on Ga/Si(112)[END_REF][START_REF] Rösch | Electron-phonon interaction in the t-J model[END_REF][START_REF] Martín-Mayor | Dynamical structure factor in disordered systems[END_REF][START_REF] Taraskin | Origin of the boson peak in systems with lattice disorder[END_REF].

The matrix elements between next nearest neighbors are

M 05 = k A 0k A 5k = A 01 A 51 + A 04 A 54 ,
(1.48)

M 09 = k A 0k A 9k = A 01 A 91 . (1.49)
The distributions of these elements are symmetrical and in average we have M 05 = 0 and M 09 = 0. Therefore, the average spring constant of these bonds are zero, and there is 50% probability of negative springs. The same is valid for 6 other next nearest neighbor matrix elements M 06 , M 07 , M 08 and M 0,10 , M 0,11 , M 0,12 .

In the 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 in Fig. 1.16. 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 zerothorder Macdonald function which logarithmically diverges at k = 0. One can find a Figure 1.17. The probability to obtain a stable system if we change a randomly chosen spring constant by δk (in units of mΩ 2 /a 0 ). similarity between our spring constant distributions and dynamical matrix element distributions obtained in [Taraskin and Elliott 2002b] for IC-glass, in [START_REF] Huang | Localization-delocalization transition in Hessian matrices of topologically disordered systems[END_REF] for simple fluid with short-ranged interactions (see Fig. 1 in these papers), and in [START_REF] Christie | Modelling the structural and vibrational properties of amorphous materials[END_REF]] 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 [Taraskin and Elliott 2002b;[START_REF] Huang | Localization-delocalization transition in Hessian matrices of topologically disordered systems[END_REF][START_REF] Christie | Modelling the structural and vibrational properties of amorphous materials[END_REF]].

In the three-dimensional case, we have almost half ((4-1/ √ 6)/8 ≈ 45%) negative springs. The dynamical matrix M = AA T ensures the stability of the lattice, however, the high concentration of negative springs makes it extremely soft. It is a subtle equilibrium: if we decrease the spring constant of a random spring, the lattice certainly becomes unstable (Fig. 1.17). Therefore, in the thermodynamical limit N → ∞, we have a critical system.

Concluding this part, we can easily include into consideration the next neighbor shell for the 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 the 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 qualitatively change the main results. Therefore, we will restrict our consideration to the simplest case outlined before. Young modulus E (in units of mΩ 2 /a 0 ) as a function of µ for dynamical matrix M = AA T + µM 0 built on a simple cubic lattice with N = 100 × 100 × 100 particles (one realization). Black dots are calculated values, and the line is the best least-square fit.

The macroscopic rigidity

Real glasses have nonzero macroscopic rigidity, and the low-frequency modes are acoustic phonons. For example, numerical calculations show that only 0.2% of vibrational modes of amorphous SiO 2 is acoustic modes [START_REF] Taraskin | Ioffe-Regel crossover for plane-wave vibrational excitations in vitreous silica[END_REF]. To include phonons in the picture, we should have a 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 possibility is to add to the random matrix AA T a "crystalline part" (other possibilities are discussed in Section 1.9) M = AA T + µM 0 .

(1.50)

Here A is the same random matrix built on a 3d simple cubic lattice as in the Section 1.3. Matrix M 0 is a positive definite crystal dynamical matrix for the same lattice with nonzero matrix elements between the nearest neighbors M 0ij = Ω 2 (the same units as AA T term). We will see below that the dimensionless 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 in Fig. 1.16. The average value is equal to k 01 = mΩ 2 (2 + µ)/a 0 . At small values of µ ≪ 1 the change is negligible. The distribution mainly consists of strongly fluctuating part AA T (compare the distributions of k 01 for µ = 0 and µ = 0.1). Therefore, it is not obvious that such small perturbation is able to introduce a finite rigidity and phonons into the system.

The macroscopic rigidity can be also introduced by a non-regular addition (see Section 1.9). However, in this work we limit ourselves to the most simple case described by Eq. (1.50). The advantage of this form is that it comprises both cases of small (µ ≫ 1) and strong (µ ≪ 1) disorder.

To find the rigidity (as a function of µ), we calculated numerically the Young modulus E of the lattice with dynamical matrix given by Eq. ( 1.50) for µ = 0. We took a very big cubic sample with N = L × L × L = 10 6 particles to reduce fluctuations and possible non-affine response.

The results of these calculations are shown in Fig. 1.18 for cubic sample with N = 10 6 particles (we have checked that for µ > 10 -4 for such a big sample the fluctuations of Young modulus from sample to sample are small so we can use one realization only. It is different from the case µ = 0 where the relative fluctuations of the Young modulus are of the order of unity). As we can see from the fit, the Young modulus has the following dependence on µ:

E = mΩ 2 a 0 µ, µ ≫ 1, (1.51) E = 1.5 mΩ 2 a 0 √ µ, µ ≪ 1. (1.52)
As a result, for µ ≫ 1, we have a usual crystal , where the disorder is relatively small, and relation (1.51) is obvious. For µ ≪ 1 the force-constant disorder is strong.

The fluctuations of the non-diagonal matrix elements M i =j are much bigger than the averaged values. In this case Young modulus E ∝ √ µ. It is much bigger than the crystal result (1.51). 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 cannot 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 future work (see also Section 2.5).

But 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 VDOS calculated numerically for the matrix (1.50). Below we will consider the case of strong and moderate force-constant disorder when 0 µ 1.

The total VDOS g(ω), normalized to unity and calculated numerically by the kernel polynomial method (KPM, see Appendix B) for dynamical matrix (1.50) and different values of µ, is shown in Fig. 1.19. We see from the figure that the VDOS for µ = 0 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 VDOS in the gap, as we will show in the next Chapter, is built by acoustic phonon-like modes and at low frequencies goes to zero as g(ω) ∝ ω 2 . The term phonon gap is is used because if conditions (1.14) are violated, then addition µM 0 to random matrix AA T opens a hard gap in the gapless vibrational spectrum (see Fig. 1.21 below).

Just above this gap, the VDOS has a sharp maximum at frequency ω max , which we will identify with the width of the gap. As follows from Fig. 1.19, the maximum frequency for µ ≪ 1 increases as ω max ∝ √ µ. Above the maximum the vibrational excitations remain to be diffusons (see Section 2.3).

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 a simple cubic lattice) we have well-defined phonon modes with the dispersion law

ω 2 cryst = 4µΩ 2 sin 2 q x a 0 2 + sin 2 q y a 0 2 + sin 2 q z a 0 2 . (1.53)
The maximum frequency, in this case, is equal to ω max, cryst = 2 √ 3µΩ ∝ √ µ which qualitatively (but not quantitatively) explains the aforesaid dependence ∝ √ µ.

However the sound velocity in this pure crystalline lattice case v cryst = a 0 Ω √ µ.

Though according to Eqs. (1.36, 1.52

) v ∝ µ 1/4 for M = AA T + µM 0 what is much bigger than √ µ for small values of µ ≪ 1.
It means that simple superposition approach does not work in this case, and a physical picture is more complicated, and the Young modulus E also depends on the amplitude of the random part AA T .

Since the VDOS g(ω) is normalized to unity for all values of µ, we conclude from Fig. 1.19 (comparing the VDOS for µ = 0 with VDOS 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 VDOS 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 2.1), the frequency ω max is correlated with the position of the boson peak ω b (the maximum in the reduced VDOS g(ω)/ω 2 ). Therefore, the appearance of the boson peak in disordered systems is not necessarily related to the acoustic van Hove singularity in crystals as was proposed recently [START_REF] Schirmacher | Harmonic vibrational excitations in disordered solids and the "boson peak[END_REF][START_REF] Taraskin | Origin of the boson peak in systems with lattice disorder[END_REF][START_REF] Chumakov | Equivalence of the Boson Peak in Glasses to the Transverse Acoustic van Hove Singularity in Crystals[END_REF]].

The straight lines in Fig. 1.19 correspond to the phonon VDOS g ph (ω) determined by Eq. (1.35) with the sound velocity v = E/ρ 0 and E calculated from Fig. 1.18. 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. 1.20 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. 1.19, the low-frequency part (below ω max ) corresponds to 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. 2.4). In the highfrequency part (above ω max ) P (ω) is approximately independent of the frequency and coincides with participation ratio for µ = 0. As we will show in Section 2.3 this range corresponds to diffusons. A similar rise in the participation ratio with It is important to emphasize that the sum rule (1.14) is crucial for the existence of the acoustical phonon excitations. If they are not obeyed, then, instead of soft phonon gap in the vibrational spectrum shown in Fig. 1.19, we have a hard gap shown in Fig. 1.21. Inside the hard gap there are no vibrations at all. The dynamical matrix M , in this case, was taken in the same form (1.50). But diagonal elements A ii of the matrix A were taken as independent Gaussian random variables with average A ii = 0 and variance A 2 ii = Ω 2 . As a result, the condition (1.16) (and therefore (1.14)) 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 ∝ √ µ.

Non-crystalline origin of acoustical phonons

In this Section we demonstrate, that appearance of acoustical phonons (and macroscopic rigidity) in the system is not related to the crystalline order in the term µM 0 .

Lattices with cut out bonds

Let us consider the case when some part of springs are cut out from the matrix µM 0 in the dynamical matrix (1.50). 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% [START_REF] Stauffer | Introduction To Percolation Theory[END_REF]. If p < p c , then there is no infinite cluster of connected springs and, therefore, the matrix µM 0 with cut out springs has itself no acoustical phonon-like modes at all. Nevertheless, the full dynamical matrix (1.50) still has well-defined phonon modes with the VDOS g(ω) ∝ ω 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 in Fig. 1.22. The straight lines show the phonon contribution to the VDOS calculated from Eq. (1.35) with sound velocity given by Eq. (1.36). The Young modulus E was calculated numerically using Eq. (1.39) for the lattice with N = 10 6 atoms (one realization) in the same way, as it was done in Section 1.8.

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 of 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. In Section 2.5 we will show different scaling relations for the model M = AA T + µM 0 . The model

Conclusion

In this Chapter we demonstrated that dynamical sparse random matrices of the general form M = AA T + µM 0 with nonnegative eigenvalues ε = ω 2 can be successfully used to describe sufficiently general properties of vibrational spectra of amorphous solids. We have shown that the system with µ = 0 is a critical system with zero macroscopic rigidity in which acoustic phonons cannot propagate.

Compared to the currently used molecular dynamics methods, they have the advantage that the construction of the random dynamical matrix corresponding to a stable system requires significantly fewer efforts than the similar calculations for real glasses with specific interatomic interaction potentials. In many cases, the results appear to be very similar.

The study of the problem of the localization of these vibrational modes in the three-dimensional system led us to the conclusion that, despite a high degree of disorder, the majority of these modes are delocalized harmonic excitations. This is evidenced by the values of the participation ratio and the statistics of vibrational levels where the level repulsion effect clearly manifests itself. Our results are in good agreement with the corresponding results obtained for real glasses by molecular dynamics methods.

Chapter 2 Diffusion of vibrations

In this Chapter we consider transport properties of phonons in a model amorphous solid with the dynamical matrix M = AA T + µM 0 . Analyzing properties of this matrix, we calculate the dynamical structure factor S(q, ω), the phonon dispersion law ω q , and phonon mean-free path l(ω). Comparison of the later with phonon wavelength λ determines the Ioffe-Regel crossover frequency ω ir which goes to zero when µ → 0. We show that above ω ir , phonons cease to exist and they are transformed to diffusons. In Section 2.3 we consider properties of diffusons and show that they spread by means of diffusion. In Section 2.4 we find an asymptotic behavior of the thermal conductivity and compare it with the experimental data. In Section 2.5 we discuss scaling properties of the model (their dependence on parameters Ω and µ). We show that they map directly onto the scaling observed in systems near jamming transition point.

Dynamical structure factor

The dynamical structure factor is a self-correlation function of the normal modes [START_REF] Shintani | Universal link between the boson peak and transverse phonons in glass[END_REF]. In the scalar model (modes have no polarization, see previous Chapter) it can be defined as

S(q, ω) = π N N j=1 N i=1 u i (ω j )e iqR i 2 δ(ω -ω j ) (2.1)
where u i (ω j ) is the displacement of the ith atom for jth eigenmode. The prefactor π/N was chosen for convenience.

To find the phonon dispersion curve (dependence of the phonon frequency ω on the wavevector q) and phonon mean-free path l(ω) we calculate the dynamical structure factor (2.1) as a space and time Fourier transform of the atomic displacements u(R i , t) ≡ u i (t)

S(q, ω) = 2 N T max N i=1 e -iqR i Tmax 0 u(R i , t)e iωt dt 2 .
(2.2)

For that we ascribed to all the atoms at the initial moment t = 0 random displacements u(R i , 0) (from Gaussian distribution with zero mean and unit variance) and zero velocities. Then, numerically solving Newton equations (with all masses m i = m) we analyzed the atomic dynamics at t = 0. Let u(R i , t) be the i-th atomic displacement as a function of atomic coordinate R i and time t. One can show that the definitions (2.1) and ( 2.2) are equivalent. Another possibility to calculate S(q, ω) is the kernel polynomial method (KPM, see Appendix B). It gives the same results, it is more effective from the computational point of view and it will be used in the next Chapter.

Since equilibrium atom positions R i are discrete and form a cubic lattice, the wave vectors q 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/a 0 L where integer numbers n are -L/2 n L/2. In this case, the normalized density of states is related to the structure factor by the sum over the first Brillouin zone g(ω) = 1 π q S(q, ω).

(2.3)

Acoustic phonons

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 in Fig. 2.1. For the fitting curves we used the Lorentz distribution (fitting to the DHO model gives slightly better results, see the next Chapter)

S(q, ω) ∝ 1 (ω -ω q ) 2 + (∆ω) 2 .
(2.4)

From the fit, we can find both the phonon frequency ω q and the phonon linewidth ∆ω, which corresponds to the phonon inverse lifetime Γ = 2∆ω. Factor 2 takes into account that ∆ω corresponds to the decay of the amplitude of the vibration. The results for ω q are shown in Fig. 2.2 for three values of µ and q 100 . For sufficiently small values of wavevector q we see a nice linear dispersion curve ω q = vq, The Lorentz dispersion curves for different wave vectors q 100 direction and µ = 0.1. Closed diamonds correspond to the calculated values of S(q, ω) and lines are fitting curves according to Eq. (2.4). The number of atoms N = 50 3 (one realization). Insets: the Lorentz dispersion curves for q = 0.5 and q = 0.75. The frequency ω has units of Ω, wavevector q has units of 1/a 0 .

with the sound velocity v given by Eq. (1.36). It is independent of the q direction (i.e. the sound velocity is isotropic). With the 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 in Fig. 2.2 scales approximately as √ µ and the horizontal axis scales as µ 1/4 (sound velocity

v ∝ √ E ∝ µ 1/4
, and q max ≈ ω max /v ∝ µ 1/4 ).

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 the sharp maximum in VDOS g(ω) (see Fig. 1.19). Similar phenomenon exists in polariton physics (polaritons are quasiparticles resulting from strong coupling of electromagnetic waves with an electric or magnetic dipole-carrying excitations. They are an expression of the common quantum phenomenon known as level repulsion, also known as the avoided crossing principle. Polaritons describe the crossing of the dispersion of light with any interacting resonance. See for example [START_REF] Gurevich | Transport in phonon systems[END_REF]). The dip in the participation ratio P (ω) for µ = 0.001, µ = 0.01, and µ = 0.1 at ω ≈ ω max (see Fig. 1.20) 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. The dependence ω q on q for q 100 direction for various µ (1, 0.1, 0.01) in a cubic sample with N = 50 3 (one realization). Filled and open diamonds are the maximums of S(q, ω) as a function of ω for each discrete value of q n for frequencies below and above the Ioffe-Regel crossover correspondingly (see text below for details). Solid lines correspond to halves of the maximums. Dashed lines show ω = vq linear dependence with sound velocity v = E/ρ 0 . Horizontal dotted lines correspond to the maximum frequency ω max in g(ω) (taken from Fig. 1.19). Insets show the group velocity v g = dω/dq as a function of ω. The frequency ω has units of Ω units, wavevector q has units of 1/a 0 .

The negative dispersion of the group velocity v g is responsible also for the pronounced rise of the phonon VDOS above the ω 2 dependence, given by Eq. (1.35). It is clearly seen in Fig. 1.19. Indeed, taking the dispersion into account and disregarding weak anisotropy (taking place only for µ = 1) we can write instead of Eq. (1.35)

g ph (ω) = a 3 0 2π 2 q 2 (ω) v g (ω) . (2.5)
Here v g (ω) = dω/dq is the group velocity shown in Insets in Fig. 2.2. Taking for q(ω) and v g (ω) the data from Fig. 2.2 we obtain the points (filled and open diamonds) shown in Fig. 1.19. 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 linewidth ∆ω can be also found from fits similar to those shown in Fig. 2.1. The results are shown in Fig. 2.3. As follows from this figure, ∆ω ∝ ω 4 and does not depend on the direction of q. It is similar to Rayleigh scattering of phonons on a static disorder. However, 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 linewidth is due to strong resonant scattering of phonons by quasilocal vibrations responsible for the sharp peak in the VDOS, for q 100 , △ for q 110 , for q 111 . Filled and open symbols refer to excitations below and above the Ioffe-Regel crossover frequency ω ir correspondingly. similar to those introduced in [START_REF] Buchenau | Interaction of soft modes and sound waves in glasses[END_REF]]. The deep minimum in the diffusivity D(ω) around frequency ω max also supports this idea (see Fig. 2

.14).

With known value of ∆ω (and Γ), the phonon mean-free path l(ω) can be calculated as follows

l(ω) = v g Γ = v g 2∆ω . (2.6)
The phonons are well-defined excitations if their mean-free path l(ω) exceeds the phonon wavelength λ = 2π/q (the Ioffe-Regel criterion 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. 2.4 shows the ratio l(ω)/λ as a function of ω for several values of µ and different directions of the wavevector q.

The boundary between filled and open symbols (the full horizontal line) corresponds to frequency ω ir . Thus, filled and open symbols in Figs. 2.2, 2.3, 2.4, 2.6 belong to phonons with frequencies below and above the Ioffe-Regel crossover frequency correspondingly.

Usually, the Ioffe-Regel crossover frequency ω ir in glasses is correlated with the position of the boson peak ω b , see [START_REF] Gurevich | Theory of lowenergy Raman scattering in glasses[END_REF][START_REF] Parshin | Interaction of quasilocal harmonic modes and boson peak in glasses[END_REF][START_REF] Rufflé | Glassspecific behavior in the damping of acoustic-like vibrations[END_REF][START_REF] Rufflé | Boson peak and its relation to acoustic attenuation in glasses[END_REF][START_REF] Shintani | Universal link between the boson peak and transverse phonons in glass[END_REF]] and references therein. It is the frequency where the reduced VDOS g(ω)/ω 2 has a maximum. We also have a rather sharp boson peak in our disordered lattices (Fig. 2.5). As follows from Fig. 2.6 the left side of the boson peak is built from phonons having the negative dispersion of the group velocity dω q /dq. Similar conclusion was made recently for 2d and 3d Lennard-Jones glasses [START_REF] Tanguy | Vibrational modes as a predictor for plasticity in a model glass[END_REF][START_REF] Léonforte | Continuum limit of amorphous elastic bodies. III. Three-dimensional systems[END_REF][START_REF] Monaco | Anomalous properties of the acoustic excitations in glasses on the mesoscopic length scale[END_REF]. The right side of the boson peak consists of 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 2.1 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

N ph = ωir 0 g(ω)dω.
(2.7) These values are also given in the Table 2.1. We see that the relative number of phonons in the lattice is small for all investigated values of µ. It is in agreement with similar estimates for amorphous silicon (4%) [START_REF] Allen | Diffusons, locons and propagons: character of atomic vibrations in amorphous Si[END_REF]] and amorphous SiO 2 (0.2%) [START_REF] Taraskin | Ioffe-Regel crossover for plane-wave vibrational excitations in vitreous silica[END_REF]].

µ ω max /Ω ω b /Ω ω ir /Ω N ph 1 2.

Diffusons

In this section we are going to consider properties of diffusons (we introduced this notation in Section 1.4). 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.

Diffusion of momentum

First let us consider diffusion of momentum. Usually, the diffusion of momentum is related to viscosity in the system. If all atomic masses are equal (m i = m), the diffusion of momentum is equivalent to the diffusion of atomic 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 atomic displacements vanishes

i u i (t) = 0, (2.8)
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. One can complete the analogy by subtraction the mean displacement (the number of "particles" divided by the number of atoms) and taking the infinite number of "particles". In this case the displacements are not bounded, and the process is continuous.

By analogy with the diffusion of "particles" the information about diffusivity of displacements is absorbed in the displacement structure factor S(q, ω) (2.2). We remind that to calculate this structure factor we ascribed at the initial moment t = 0 the random displacements to all the atoms with Gaussian distribution (with zero mean and unit variance) and velocities equal to zero. So the condition (2.8) at t = 0 was satisfied. Therefore, let us analyze now this structure factor in the diffuson frequency range.

Firstly, let us consider the case of µ = 0 when phonons are absent, and only diffusons are present in the lattice. Fig. 2.7 shows the structure factor S(q, ω) as a function of wavevector q for three different directions in q space (symbols) and for q 100 , △ for q 110 , for q 111 . Full lines correspond to the structure factor S rw (q, ω) of the random walk on the lattice given by Eq. (2.9) with D rw = 0.7. Dashed line corresponds to the limit q ≪ 1 (see Eq. ( 2.12)). The frequency has units of Ω, the wavevector has units of 1/a 0 .

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 [START_REF] Kehr | Diffusion of particles on lattices[END_REF] for the case of the random walk on a lattice, S rw (q, ω) is given by expression

S rw (q, ω) = 2Γ(q) ω 2 + Γ 2 (q) . (2.9)
It is a Lorentz function, with a width Γ(q) given by

Γ(q) = D rw Q 2 (q) (2.10)
where D rw is a diffusion constant of the random walk. In a simple cubic lattice (with lattice constant a 0 ) the function Q(q) reads

Q(q) = 2 a 0 sin 2 q x a 0 2 + sin 2 q y a 0 2 + sin 2 q z a 0 2 . (2.11)
For small values of q ≪ 1/a 0 , Q(q) = q and in the continuum limit we have the well-known result for the diffusion structure factor [START_REF] Landau | Statistical Physics[END_REF], §89]

S rw (q, ω) = 2D rw q 2 D 2 rw q 4 + ω 2 .
(2.12)

Let us note that the structure factor (2.9) has a maximum at q values obeying the condition ω = Γ(q) = D rw Q 2 (q). (2.13)

We can specify it as a dispersion law for diffusons. The width of the maximum is Γ(q). For q ≪ 1/a 0 , Γ(q) = D rw q 2 .

A comparison of the displacement structure factor S(q, ω), (2.2), and the structure factor of the random walk S rw (q, ω), (2.9), is shown in Fig. 2.7. One fitting parameter was the diffusion coefficient D rw in Eq. (2.10). From the comparison of this data, we obtain D rw ≈ 0.7Ωa 2 0 . It means that the diffusion coefficient of atomic displacements D u ≈ 0.7Ωa 2 0 (see Section 2.6). Another fitting parameter was a height h(ω) of the random walk structure factor in the maximum. According to Eq. (2.9), in the maximum Γ(q) = ω and h(ω) = 1/ω, but to fit the data points in Fig. 2.7 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. ( 2.3)). As we can see from the figure, for the investigated frequencies the fit is perfect. With increasing frequency above ω ≈ 2Ω, the fitting becomes more and more poor since we approach the localization threshold at ω loc ≈ 5.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 atomic displacements at some frequency ω, expressed through eigenvectors u(R, ω) of the dynamical matrix M , reads

C(r, ω) = R u(R + r, ω)u(R, ω).
(2.14)

It is a Fourier transform of the displacement structure factor (2.1)

C(r, ω) = 1 8π 4 S(q, ω)e iqr dq.

(2.15)

Let us compare this correlation function with correlation function of the random walk. For distances bigger than the period of the lattice a 0 we can make use of the limit of small q ≪ 1/a 0 and integrate Eq. ( 2.12) for the random walk structure factor taken in the continuous media approximation. As a result, we derive

C rw (r, ω) = 1 2π 2 rD rw exp -r ω 2D rw cos r ω 2D rw .
(2.16) Fig. 2.8 shows a good agreement of our correlation function (2.15) with the correlation function of the random walk (2.16). For all investigated frequencies the numerical data collapse together and become indistinguishable from the theoretical prediction (2.16). We can see also on this figure the anticorrelation phenomenon (the (2.17)

It is also in a good agreement with our numerical results. Therefore, we can call a corresponding value of r found from Eq. ( 2.17) as a radius of diffuson. It is a typical size of the regions vibrating with frequency ω and having the same sign of all atomic displacements. According to (2.17), the radius of diffuson is given by 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 S n (q, ω) = S(q, ω)/ max q ′ S(q ′ , ω).

r d (ω) = π D rw 2ω ∝ ω -1/2 . ( 2
(2.19)

Figure 2.9. The normalized structure factor S n (q, ω) as a function of q (in units of 1/a 0 ) for some direction in q space and for each frequency ω (in units of Ω) for various values of µ (0, 0.01, 0.1, 1). The sample size is N = 50 3 . The averaging is performed over 100 realizations. Left sides of all plots are for q 111 , right sides are for q 100 . White horizontal dashes show the Ioffe-Regel crossover frequency ω ir . For µ = 1 the frequency ω ir is slightly different for different q directions. Black full line corresponds to Eq. ( 2.13) for the random walk on a simple cubic lattice with diffusion constant D rw = 0.7Ωa 2 0 .

The results are shown in Fig. 2.9 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 in Fig. 2.9a, which corresponds to diffusons. The maximum of the normalized structure factor S n (q, ω) (white regions) agrees well with Eq. (2.13) (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 atomic displacements is independent of µ. Deviations from S rw (q, ω) take place at high frequencies near the localization threshold.

For µ = 0 the radius of diffuson (2.18) takes a maximum value at ω ≈ ω ir . At smaller frequencies we have well-defined phonons. Since ω ir ≃ Ω √ µ and D rw ≈ a 2 0 Ω we can write for 0 < µ 1

r d (ω ir ) ≡ r c ≃ D rw /ω ir ≃ a 0 µ -1/4 .
(2.20)

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 The same normalized structure factor S n (q, ω) as in Fig. 2.9 but in q space in plane q x q y (q z = 0) for ω = 0.5Ω. The left picture corresponds to µ = 0 (a) and the right to µ = 0.1 (b).

with the Ioffe-Regel wavelength λ ir = 2π/q ir corresponding to frequency ω ir (see Section 2.5). Samples with a size smaller than r c have no phonon-like modes at all.

To compare phonon and diffuson structure factors, we show in Fig. 2.10 a cross section of the structure factor S n (q, ω) in q space for q z = 0 and frequency ω = 0.5Ω 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.

Diffusion of energy

We said above that there are two types of diffusion because we have two conservation laws. In this section we will consider the diffusion of the energy. In the harmonic approximation, all normal modes are independent so the energy conserves in the each small frequency interval. Therefore, the diffusivity of the energy D(ω) is a function of the frequency ω.

There are two common methods to determine the diffusivity D(ω). The first approach to calculate the diffusivity of energy D(ω) for vibrations with frequency ω is a direct numerical solution of Newton's equations for a cubic sample with N = L × L × L atoms 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 atoms, let us apply external forces with frequency ω and random phases ϕ i to all the atoms in the central layer x = 0 of our sample

f ext i (t) = sin(ωt + ϕ i ) exp - t 2 2τ 2 exc (2.

21)

Figure 2.11. The dependence of R 2 (t) in the case of µ = 0 for one sample with N = 100 × 100 × 100 atoms and 14 different frequencies (ω/Ω = 0.5, 1, 1.5, . . . , 7, from top to bottom). The numbers indicate integer frequencies. The slope of each line corresponds to each black dot in Fig. 2.12. Two points at ω = 2Ω and ω = 6Ω correspond to two distributions of energy E(x, t) over the sample for delocalized and localized modes correspondingly. They are shown in Fig. 2.13 (see below). The time has units of 1/Ω, the excitation width R has units of a 0 .

where ωτ exc ≫ 1. For a big sample with N = 100 × 100 × 100 = 10 6 atoms it is sufficient to excite only one atomic layer x = 0 with 100 × 100 = 10 4 atoms. The addition of two or more neighbor layers does not change the results. Increasing the width of the excited layer one should increase the length of the sample as well. The right and the left sides of the sample have coordinates x r,l = ±a 0 L/2. In such a way we excite vibrations with frequencies near frequency ω distributed in a small frequency interval (ω -1/τ exc , ω + 1/τ exc ). In calculations we used τ exc = 5/Ω for all frequencies ω. We started our calculations at time t 0 = -5τ exc 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

R 2 (t) = 1 E tot N i=1 x 2 i E i (t). (2.22)
Here x i is the x coordinate of the i-th atom, E i (t) is the vibrational energy of ith atom, and sum is taken over all atoms 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 > 5τ exc ).

Figure 2.12. The dependence of diffusivity D(ω) (in units of Ωa 2 0 ) on ω (in units of Ω) for µ = 0. Black dots are calculated by the direct solution of Newton's equations from Eqs. (2.22, 2.24) and Fig. 2.11 for N = 100 3 atoms (one realization). Full lines for N = 10 3 , 14 3 , 20 3 are calculated using formula of Edwards and Thouless (2.29) with c = 1 (see below). Averaging for lines is performed over frequencies in the small interval (ωδω, ω + δω) with δω = 0.25Ω and over several thousands realizations.

The energy of i-th atom E i (t) we define as a sum of the kinetic energy and a half of the potential energy of connected bonds

E i (t) = v 2 i (t) 2 - 1 4 j M ij u i (t) -u j (t) 2 .
(2.23)

Here v i (t) = ui (t) is an atomic velocity with the same notations as Eq. 1.2. Summation over all atoms in Eq. ( 2.22) we can divide into two steps. First we sum over all atoms 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 we have sample size L ≫ 1 in our case, we can change summation over different layers to integration over coordinate x for times where R(t) ≫ a 0 .

We will apply this method to the case of µ = 0 (i.e. for the lattice without phonons). The results are shown in Fig. 2.11. 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 ω. To calculate the slope, we take the time interval ∆t where, on the one hand, t > 5τ exc , and on the other hand, R ≪ a 0 L/2.

From the slope of R 2 (t), we can calculate the diffusivity of modes D(ω) using one-dimensional formula R 2 (t) = 2D(ω)t. (2.24) This diffusivity is shown by black dots in Fig. 2.12. At small frequencies it is approximately constant, and then it decreases with the frequency approaching zero at the localization threshold, ω loc ≈ 5.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 in Fig. 2.11. Black diamond corresponds to delocalized mode with frequency ω = 2Ω and has coordinates t = 234/Ω and R 2 = 166a 2 0 . The distribution of energy E(x, t) over the sample calculated numerically at this moment is shown by black diamonds in Fig. 2.13. The data are perfectly fitted by solid line drawn according to the solution of diffusion equation in 1d case

E(x, t) = E tot √ 2πR 2 exp - x 2 2R 2 , (2.25)
with the value of R 2 = 166a 2 0 . Black triangle in Fig. 2.11 corresponds to the localized mode with frequency ω = 6Ω and has coordinates t = 900/Ω and R 2 = 22a 2 0 . The distribution of energy E(x, t) over the sample calculated numerically at this moment is shown by black triangles in Fig. 2.13. This distribution is drastically different from the previous case. For localized modes we expect the usual exponential decay

E(x, t) = E tot √ 2R exp - √ 2|x| R . ( 2

.26)

The fit of the numerical data with this function and R 2 = 22a 2 0 is shown in Fig. 2.13. The fit is perfect except for the central point at x = 0 which lies noticeably above prediction of Eq. (2.26). The coefficients in Eqs. (2.25, 2.26) were taken to satisfy the obvious rules

∞ -∞ E(x, t)dx = E tot , 1 E tot ∞ -∞ x 2 E(x, t)dx = R 2 .
(2.27)

To find the diffusivity D(ω) for µ = 0, the method of numerical solution of Newton's equations is not accurate because we have phonons with long mean-free paths in this case. 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 [START_REF] Edwards | Numerical studies of localization in disordered systems[END_REF]]

D(ω i ) ≃ (a 0 L) 2 |∆ω i | (2.28)
where a 0 L is the length of the sample and ∆ω i is the sensitivity of the eigenfrequency ω i to a twist of boundary conditions. Physically, the frequency shift ∆ω i means the inverse time of the signal propagation from one boundary to another.

Allen et al. proposed a more precise formula with twisting by a small angle ϕ [START_REF] Allen | Diffusons, locons and propagons: character of atomic vibrations in amorphous Si[END_REF]:

D(ω) = c lim ϕ→0 (a 0 L) 2 ϕ 2 |∆ω(ω)| (2.29)
where c is some constant of the order of unity. The symmetric real matrix M was defined as usual (1.50) with periodic boundary conditions. The twisting of the matrix M by angle ϕ gives a new Hermitian matrix M ′ obtained as follows. For bonds between the left (l) and the right (r) boundaries of our cubic sample

M ′ lr = M lr exp(iϕ), M ′ rl = M rl exp(-iϕ). (2.30)
For all other bonds M ′ jk = M jk . So ∆ω i is the difference between i-th eigenfrequencies of matrices M and M ′

∆ω i = ω i -ω ′ i . (2.31)
Twisting of boundary conditions was performed for x direction only. For other two directions, the periodic boundary conditions were used.

The averaging in Eq. 2.29 is performed over frequencies ω in the small interval (ωδω, ω + δω) with δω = 0.25Ω and over several thousand realizations. For µ = 0 the results for D(ω) are shown in Fig. 2.12 for three different cubic samples (full ) as a function of ω (in units of Ω) for various µ (0, 0.01, 0.1, 1) for sample with N = 14 3 (crosses). The diffusivity was calculated using formula of Edwards and Thouless (2.29) with c = 1 and averaged over two thousand realizations. The arrows indicate frequencies ω max in the VDOS g(ω) for corresponding values of µ. Open symbols correspond to phonon diffusivity (2.32) below the Ioffe-Regel crossover frequency ω ir . lines). We compared these results with the numerical solution of Newton equations for µ = 0 (black dots) and get the constant c ≈ 1. Then we used this c value for µ = 0. The results are shown in Fig. 2.14. For µ = 0 we see clearly two different frequency regions in the function D(ω).

At low frequencies, the diffusivity increases with decreasing of ω. This range corresponds to phonons. Indeed, the diffusivity of phonons D(ω) can be calculated as follows

D(ω) = 1 3 l(ω)v g (ω).
(2.32)

Open symbols in Fig. 2.14 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 [START_REF] Xu | Energy transport in jammed sphere packings[END_REF][START_REF] Vitelli | Heat transport in model jammed solids[END_REF]]. The deep minimum in the diffusivity at ω ≈ ω max corresponds to strong scattering of phonons by the quasilocal vibrations near the sharp peaks in the VDOS g(ω) (see Fig. 2.6).

The thermal conductivity

One can find the thermal conductivity from the VDOS g(ω) and the diffusivity D(ω) [START_REF] Sheng | Heat conductivity of amorphous solids: simulation results on model structures[END_REF][START_REF] Feldman | Thermal conductivity and localization in glasses: Numerical study of a model of amorphous silicon[END_REF]]

κ(T ) = 1 a 3 0 ∞ 0 g(ω)D(ω)C(ω, T )dω (2.33)
where C(ω, T ) is the specific heat of the harmonic oscillator

C(ω, T ) = ω T 2 e ω/T
(e ω/T -1) 2 .

(2.34)

Diffusons have an approximately constant VDOS g(ω) ∼ 1/Ω and diffusivity D(ω) ∼ a 2 0 Ω. In this case from Eq. ( 2.33) we get a linear temperature dependence of the thermal conductivity

κ(T ) ∼ k 2 B T a 0 , k B T Ω. (2.35)
At high frequencies, all vibrational modes are excited, and the thermal conductivity saturates

κ(T ) ∼ k B Ω a 0 , k B T Ω.
(2.36)

Eqs. (2.35) and (2.36) with a 0 = 3 Å and Ω = 30 meV have a quantitative agreement with experimental data of the thermal conductivity of amorphous SiO 2 in the temperature range above 20 K (red line in Fig. 2). At the high-temperature limit, we get the experimental value κ ≈ 2 W/(m•K). Note, that the maximum frequency in SiO 2 is several times bigger than Ω (Fig. 4) as well as the maximum frequency in the random matrix model (Fig. 2.3).

Scaling relations

Finally, the concept of diffusons allows us to establish useful scaling relations between observable values and important parameters of our model. We have the dimensionless parameter µ. The second important parameter of the model is the variance of non-diagonal elements A ij of the random matrix

A A 2 ij = Ω 2 .
(2.37)

The parameter Ω has units of frequency and assigns the scale of typical frequencies in the system. In particular, the normalized density of states g(ω) for µ = 0 shown in Fig 1 .5 has the following scaling relation

g(ω) ≃ 1/Ω. (2.38)
Other dimensional parameters are the lattice constant a 0 and the atomic mass m. The diffusivity D does not depend on µ, therefore from the dimensional analysis we get D ≃ a 2 0 Ω.

(2.39)

From Eq.( 1.52) the Young modulus is

E ≃ mΩ 2 a 0 √ µ.
(2.40)

The Ioffe-Regel frequency is also proportional to √ µ ω ir ≃ Ω √ µ.

(2.41)

The corresponding length scale (2.20) is

λ ir ≃ l(ω ir ) ≃ q -1 ir ≃ D/ω ir ≃ D/v ≃ a 0 µ -1/4 . (2.42)
It is interesting to note that these scaling relations are identical to those found in jamming transition [START_REF] Vitelli | Heat transport in model jammed solids[END_REF]. Authors [START_REF] Vitelli | Heat transport in model jammed solids[END_REF]] study a model of amorphous packing of frictionless spheres interacting via the repulsive pair potential (see Chapter 4 for more details) (2.43) where the distance between the centers of atoms 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 atoms are just touching each other, and there is no overlap [START_REF] O'hern | Jamming at zero temperature and zero applied stress: The epitome of disorder[END_REF]]. At densities lower than φ c atoms 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

U (r ij ) ∝ (1 -r ij /σ ij ) α if r ij < σ ij , U (r ij ) = 0 if r ij > σ ij ,
∆z ≡ z -z c ∼ ∆φ 1/2 (2.44)
where z c = 2d, with d being the space dimension.

Chapter 2. Diffusion of vibrations 54 It was found that different quantities exhibit scaling behavior near the jamming point. According to [START_REF] Vitelli | Heat transport in model jammed solids[END_REF]], the Ioffe-Regel crossover frequency ω * and the shear modulus G behave as (below we use the notation of the paper [START_REF] Vitelli | Heat transport in model jammed solids[END_REF])

ω * ∼ ∆φ (α-1)/2 , G ∼ ∆φ (2α-3)/2 .

(2.45)

The transverse sound velocity v t and the diffusivity in the plateau region D 0 scale

v t ∼ ∆φ (2α-3)/4 , D 0 ∼ ∆φ (α-2)/2 .
(2.46)

The applied pressure p and the plateau in the density of states g 0 depend on the packing fraction as follows [START_REF] O'hern | Jamming at zero temperature and zero applied stress: The epitome of disorder[END_REF]]

p ∼ ∆φ α-1 , g 0 ∼ ∆φ (2-α)/2 .
(2.47)

Thus, if we put µ ∼ ∆φ, Ω ∼ ∆φ (α-2)/2 , (2.48) 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 g 0 in the jamming transition model respectively. In particular, the parameters µ and Ω in our model are equivalent to the parameter ∆φ and the inverse density of states 1/g 0 in the jamming transition model correspondingly.

In the Chapter we mainly considered a case of strong disorder, µ ≪ 1. Taking into account Eq. ( 2.48) 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 (2.42) and (2.48) is given by

l(ω ir ) ∼ ∆φ -1/4 , (2.49)
what also coincides with [START_REF] Vitelli | Heat transport in model jammed solids[END_REF]].

Conclusion

In this Chapter we have found that the delocalized vibrational excitations in the disordered lattice model 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 wavevector q.

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. [START_REF] Allen | Diffusons, locons and propagons: character of atomic vibrations in amorphous Si[END_REF]]. The diffusons again can be characterized by frequency ω, but have no well-defined wavevector q. Above ω ≈ ω IR the structure factor of atomic displacements S(q, ω) becomes very similar to the structure factor S rw (q, ω) of a random walk on the lattice. The corresponding vibrational linewidth Γ(q) ∝ q 2 was found in many glasses in the experiments on inelastic x-ray scattering, see for example [START_REF] Sette | Dynamics of glasses and glass-forming liquids studied by inelastic X-ray scattering[END_REF][START_REF] Ruocco | High-frequency vibrational dynamics in glasses[END_REF]] and references therein.

We have found that the diffusivity of the vibrational energy of diffusons is approximately constant (as well as the VDOS). The corresponding temperature dependence of the thermal conductivity have the same behavior as the experimental data: it is proportional to the temperature T above the plateau and then saturates.

Chapter 3 Vibrational properties of amorphous silicon-like materials

In previous two Chapters we consider a scalar model where vibrational modes have no polarization. In this Chapter we show that the polarization plays a crucial role in such materials as amorphous silicon. We begin with a detailed presentation of the numerical model of amorphous silicon. Then we compute the vibrational density of states (Section 3.2) and its transverse and longitudinal component. Then we analyze the spatial structure of the corresponding eigenmodes (Section 3.3). In Section 3.4 we study the dynamical structure factor, the dispersion laws, the corresponding phonon lifetimes and mean-free paths. Finally, in Section 3.5, we compare the results to the propagation of quasi-monochromatic wave packets with different frequencies allowing to measure the diffusivity of vibration energy in the materials, as a function of the frequency and of the bending rigidity. This allows to identify coherently well-defined vibrational domains, as it is summarized in the conclusion.

Numerical model

In this chapter we study the vibrational properties of a model amorphous silicon (a-Si) system consisting of N = 32768 atoms contained in a cubic box of lengths L x = L y = L z of approximately 87 Å (smaller systems of N = 8000 have also been studied to compare our results). The technical details of the preparation of the a-Si model are presented in Ref. [START_REF] Fusco | Role of local order in the small-scale plasticity of model amorphous materials[END_REF]. The Si-Si interaction in the system studied here is described by the Stillinger-Weber potential [START_REF] Stillinger | Computer simulation of local order in condensed phases of silicon[END_REF], where we have tuned the prefactor Λ of the three-body term as already done in [START_REF] Fusco | Role of local order in the small-scale plasticity of model amorphous materials[END_REF][START_REF] Fusco | Rheological properties vs. local dynamics in model disordered materials at low temperature[END_REF] the vibrational properties. The Stillinger-Weber potential is an empirical potential including two-body and three-body interactions, such that the total potential energy of the system is written as

U = i<j f (r ij ) + Λ i<j<k g(r ij , r ik , θ jik ) + g(r ji , r jk , θ ijk ) + g(r ki , r kj , θ ikj ) (3.1) with f (r ij ) = A B/r 4 ij -1 exp σ (r ij -a) -1 (3.2)
and

g(r ij , r ik , θ jik ) = (cos θ jik + 1/3) 2 exp α(r ij -a) -1 + α(r ik -a) -1 (3.3)
with the parameters proposed in [START_REF] Stillinger | Computer simulation of local order in condensed phases of silicon[END_REF]] A = 7.05, B = 11.60, α = 2.51 Å, σ = 2.06 Å, and a = 3.77 Å. The parameter Λ gives a measure of the bond's directionality: high values of Λ favor local tetragonal order (Λ = 21 is the original value proposed by [START_REF] Stillinger | Computer simulation of local order in condensed phases of silicon[END_REF] as an empirical model for a-Si). The atomic configurations corresponding to a-Si structures for different values of Λ have been obtained from the liquid state, using the open source LAMMPS package [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF] for classical Molecular Dynamics simulations, and following the procedure already described in [START_REF] Fusco | Role of local order in the small-scale plasticity of model amorphous materials[END_REF]]. Different configurations have been obtained, either by quenching in the NVT ensemble at a fixed density ρ = 2.303 g/cm 3 , giving rise to different final pressures as detailed in Table 3.1, either after pressure relaxation up to P ≈ 0 GPa.

In order to study the role of the local order on the vibrational properties of a-Si, we have calculated the dynamical matrices for different values of Λ. The dynamical matrix (1.3) derivative of the potential energy around the equilibrium atomic positions R i . Each realization of the amorhous silicon and its dynamical matrix M is different. In this sense the dynamical matrix M is random, however the matrix elements have complicated correlations with each other.

The elastic constants (shear and bulk modulus) are obtained as in [START_REF] Fusco | Role of local order in the small-scale plasticity of model amorphous materials[END_REF]] by measuring the quasi-static response of the system to a small deformation of the box. The corresponding values of the transverse c T and longitudinal c L sound velocities are summarized in Table 3.1 for the different values of Λ. Thick lines correspond to constant density configurations as described in Table 3.1. Dashed thin lines are relaxed configurations with P ≈ 0 GPa. Solid horizontal thin lines on the left show the low-frequency Debye predictions calculated from the static shear and bulk modules. (b) Boson peak as a function of the reduced wave vector q * = ω/c T . Vertical gray bands mark the position of q * 1 and q * 2 . Arrows show the position of the transverse Ioffe-Regel criteria (see Section 3.4).

Density of states

The dynamical matrix M has N = 3N eigenvalues that are squares of the corresponding eigenfrequencies ω j . The normalized vibrational density of states (VDOS) as a function of ω = 2πν reads

g(ω) = 1 N N j=1 δ(ω -ω j ). (3.4)
The full set of eigenvalues for a small system (with N < 10 4 ) can be obtained by standard numerical routines. We used the FEAST Eigenvalue Solver [START_REF] Polizzi | Density-matrix-based algorithm for solving eigenvalue problems[END_REF] for N = 8000 to get the full set of eigenfrequencies together with the eigenvectors of the dynamical matrix. However, this direct method requires too much time as well as random access memory for large enough systems (N > 10 4 ). For these purposes, it is necessary to use more powerful methods for the larger systems studied. In Appendix B we discuss the Kernel Polynomial Method (KPM) and velocity autocorrelation method.

The numerical results for VDOS obtained using KPM are presented in Fig. 3.1. They show the usual shape of VDOS obtained for amorphous silicon [START_REF] Kamitakahara | Vibrational spectrum of amorphous silicon: Experiment and computer simulation[END_REF], with a first peak related to transverse acoustic modes and a second well-defined peak at high frequencies that is reminiscent of optic modes in the crystal (Fig. 3.1a). The rescaling of the frequencies by the transverse sound velocity (Fig. 3.1b) allows drawing the density of states as a function of a reduced wave vector q * ≡ ω/c T . In this case, the low q-part of the spectra superimpose whatever the value of Λ, confirming the dominant transverse acoustic character of the lowfrequency vibrations, and suggesting the existence of a characteristic length at a wave vector q * ≈ 10 nm -1 independent of Λ, above which the rescaled VDOS split. The transverse sound velocities c T shown in the inset of Fig. 3.1 have a Λ 1/2 dependence at constant density, showing that the shear modulus is proportional to the threebody contribution to the total energy of the system. This effect is not shown in the longitudinal sound velocities c L because bending rigidity is not dominant for the propagation of compressive waves, contrary to shear waves. The parameter Λ thus allows tuning the transverse wave velocities independently of the longitudinal one.

The boson peak is visible after dividing the VDOS g(ν) by ν 2 (the Debye prediction) as shown in Fig. 3.2a. The shape of the boson peak (Fig. 3.2a) shows clearly a dependence on the bonds directionality quantified by the parameter Λ.

The boson peak appears to be magnified when the three-body interactions are low, and it decreases when the three-body interactions get more and more important as compared to the central interatomic forces. For Λ = 21 corresponding to a-Si, the initial very low-frequency peak is no more marked, but the boson peak is still visible with an excess of low-frequency vibrations as compared to the Debye prediction. As the value of Λ increases, the position of the peak is shifted to higher frequencies. This effect is clearly dominated by Λ. We have checked that pressure differences between the samples induce only a small change in the boson peak (thin lines in Fig. 3.2a) as compared to the role of Λ. In order to quantify the observed shift to higher frequencies, we again rescale the frequencies by the transverse sound velocity, as suggested in [START_REF] Léonforte | Inhomogeneous elastic response of silica glass[END_REF]]. The resulting reduced density of states is shown in Fig. 3.2b. The position of the boson peak as a function of the reduced wave vector q * appears now independent of Λ, suggesting a universal process dominated by transverse waves, that will be discussed later. Note, however, that the fine structure of the peak depends on the bonds directionality Λ: at a very low-frequency, a peak is visible for low values of Λ, located at q * 1 ≈ 2.7 nm -1 (corresponding to a wavelength ξ * 1 ≈ 23 Å). This very low-frequency peak disappears progressively, and a secondary peak appears at q * 2 ≈ 7.0 nm -1 (ξ * 2 ≈ 9 Å) when Λ > 21. The significance of these peaks will be discussed later.

Above we have shown that the low-frequency part of the VDOS has presumably a dominant transverse character. It would be very interesting to find a regular way to separate the VDOS into the longitudinal and the transverse components for the whole vibrational spectrum. In particular, it gives us a possibility to show more clearly that the boson peak has a transverse nature [START_REF] Schober | Vibrations and relaxations in a soft sphere glass: boson peak and structure factors[END_REF]]. In Appendix C we describe a generalized decomposition method without the notion of the wave vector, which is an ill-defined quantity in strongly disordered systems. This method is based on the volume variations of the Voronoi cells during the atomic motion. The atomic displacement of each atom can be decomposed into two components, one of them preserving the volume of the Voronoi cells. The displacements preserving the volume of each Voronoi cell are identified as transverse displacements, and the other as longitudinal displacements.

The separate contribution of longitudinal and transverse displacements to the total VDOS is shown in Fig. 3.3. In the low-frequency region (below 7 THz for Λ = 21), the transverse modes dominate the VDOS, thus confirming the transverse character of the vibrations in the region of the boson peak for amorphous silicon-like samples (Fig. 3.3b). However, at 7 THz there is a sharp transition from mostly transverse modes to mostly longitudinal ones. It corresponds to the maximum frequency (7.5 THz) of TA modes in crystalline silicon [START_REF] Tubino | Lattice dynamics and spectroscopic properties by a valence force potential of diamondlike crystals: C, Si, Ge, and Sn[END_REF]]. The maximum frequency of LA modes is much large (11 THz) due to the large difference between bulk and shear moduli. The elementary cell in the diamond-like crystalline silicon contains 2 non-equivalent atoms whose out of phase motion results in the three branches of high-frequency optic modes (one LO and two TO modes) with very low group velocity [START_REF] Kittel | Introduction to solid state physics[END_REF]]. In amorphous silicon out of phase motion of nearest atoms is similar to crystalline optic modes. These modes form the second peak of the VDOS, which is clearly seen in neutron scattering measurements [START_REF] Kamitakahara | Vibrational spectrum of amorphous silicon: Experiment and computer simulation[END_REF]]. Thus, the predominance of the longitudinal modes between 7 THz and 15 THz in amorphous silicon (Fig. 3.3a) corresponds indeed to the gap between the upper frequency of TA modes (7.5 THz) and the lower frequency of TO modes (13.9 THz) in crystalline Si. This frequency region in crystalline Si is totally occupied by LA and LO modes (without a gap). In amorphous Si in the same frequency region, the vibrations have a small transverse component (15-20%), in agreement with the results obtained in [START_REF] Marinov | Model investigation of the Raman spectra of amorphous silicon[END_REF] In order to complete this description, in the next Section we will calculate participation ratio and the correlation function to describe the geometrical structure of the eigenmodes that are obtained as eigenvectors of the dynamical matrix.

Participation ratio and spatial correlations

The exact diagonalization of the dynamical matrix is performed using FEAST Eigenvalue Solver [START_REF] Polizzi | Density-matrix-based algorithm for solving eigenvalue problems[END_REF]] on a system made with N = 8000 atoms. A series of Arrows are proportional to the displacements of the particles (×100). The 2D representation corresponds to a cut along the x-y plane (δz = 5 Å) that contains the particle supporting the largest displacement. N = 3N eigenmodes u i (ω j ) is then obtained with the corresponding eigenvalues ω 2 j . These eigenmodes are the normal modes of the amorphous material but they are not simple plane waves with a well-defined wave vector q. Examples of such eigenmodes are shown in Fig. 3.4. The low-frequency eigenmodes are a superposition of plane waves with softer regions supporting highly strained isolated vibrations (Fig. 3.4a). The modes supporting additional isolated vibrations are precursors of local plastic rearrangements when looking at the anharmonic mechanical response [START_REF] Tanguy | Vibrational modes as a predictor for plasticity in a model glass[END_REF]]. We identify them as soft modes because they occur only in the low-frequency part of the spectrum (as will be proved later) with soft spots due to very low local elastic stiffness. Other authors called the low-frequency modes quasi-localized modes [START_REF] Schober | Low-frequency vibrations in a model glass[END_REF][START_REF] Schober | Vibrations and relaxations in a soft sphere glass: boson peak and structure factors[END_REF]] in order to distinguish them from plane waves. At higher frequencies, the shape of the eigenmodes becomes more complex.

The amount of particles moving together in the vibrational eigenmodes is usually quantified by the participation ratio defined for each eigenmode j as Participation ratio as a function of the reduced wave vector q * = ω/c T for different values of the parameter Λ. Vertical gray bands show three characteristics reduced wave vectors q * 1 , q * 2 , and q * 3 . Inset: zoom on the low-frequency range.

P (ω j ) ≡ 1 N ( i u 2 i (ω j )) 2 i u 4 i (ω j ) . ( 3
For an isolated particle P ∼ 1/N , and for translational motion P ∼ 1 (see Section 1.4). The participation ratio is shown in Fig. 3.5 as a function of the reduced wave vector q * . Similarly to the VDOS, the low-frequency part of P superimposes for all Λ when plotted as a function of the reduced wave vector q * , suggesting the existence of a common geometrical origin involving mainly transverse vibrations. It can be schematized as follows: first an initial decay due to the wavelengths decrease of acoustic modes, together with very low values characteristic of soft modes. Then an increase up to a value close to P * = 0.5 (an example of such mode is shown in Fig. 3.4b). The value of P * is close to 0.6 for uncorrelated Gaussian random noise (see Section 1.4). After a secondary minimum (mode shown in Fig. 3.4c) the participation ratio decreases to zero at the mobility edge [START_REF] Allen | Diffusons, locons and propagons: character of atomic vibrations in amorphous Si[END_REF]] that follows the position of the high-frequency peak in the VDOS (see Fig. 3.1). A typical mode in this frequency range is shown in Fig. 3.4d. Quite remarkably, the position of the first minimum in the participation ratio corresponds for all Λ to the first maximum in the rescaled VDOS divided by ν 2 located at q * 1 , and the position of the common maximum P * = 0.5 is located at q * 2 corresponding to the second peak in the low-frequency rescaled VDOS divided by ν 2 . The departure from the plane waves participation ratio in this range, means that in all our systems, the boson peak is located in a frequency range where plane waves are no more the dominant contribution to the eigenmodes. This frequency range is limited by two characteristic distances ξ * 1 and ξ * 2 , that are independent on Λ. There is also a third characteristic reduced wave vector q * 3 = 11.7 nm -1 which corresponds to a secondary local minimum of the participation ratio for all values of Λ. It coincide with the sharp change of the nature of vibrations from almost transverse to almost longitudinal ones (Fig. 3.3). for the different values of Λ. The dashed line is r * = 0.449/q * , which corresponds to the first minimum of Eq. (3.10) for transverse modes. Vertical gray bands mark the positions of q * 1 , q * 2 , and q * 3 .

q 1 * q 2 * q 3 * C n (r, ) 0.1 0 0.1 r, nm
In order to detail the shape of the eigenmodes and compare with these characteristic lengths, we have computed their spatial correlation function C(r, ω). It is defined as in Ref. [START_REF] Tanguy | Continuum limit of amorphous elastic bodies: A finite-size study of low-frequency harmonic vibrations[END_REF] as

C(r, ω) = 1 N N j=1 u(r + r ′ , ω j ) • u(r ′ , ω j ) r ′ δ(ω -ω j ), (3.6)
where u(r, ω j ) is coarse-grained displacement field of jth eigenmode

u(r, ω j ) = N i=1 W(r -R i )u i (ω j ) (3.7)
with W a coarse-graining function normalized by W 2 (r)dr = 1. We used Gaussian coarse-graining function of width w CG = 0.5 Å. This length is less than the typical distance between atoms, so we can neglect the overlapping of different grains. In this case, normalization of the eigenmodes implies

u(r, ω j ) • u(r, ω j ) r = 1, (3.8)
which results in the property C(0, ω) = g(ω). For convenience we used the normalized correlation function

C n (r, ω) = C(r, ω) g(ω) . (3.9)
In order to calculate C(r, ω), we used the KPM (Appendix B). The amplitude of C n (r, ω) averaged over different directions of r for all the frequencies is shown in Fig. 3.6a for Λ = 21. Starting from C n (0, ω) = 1, it shows oscillations between positive and negative values characterizing a spatial flipping of the displacement field. For a three-dimensional plane wave with a given polarization (L or T) and the wavevector q the normalized correlation function has a form

C n (r, ω) = sin(qr) qr . (3.10)
The low-frequency behavior of the correlation function (Fig. 3.6) is indeed dominated by the wavelength of the plane wave: when plotted as a function of the reduced wave vector q * = ω/c T , it shows the characteristic behavior of transverse plane waves. Indeed, Fig. 3.6c shows the position of the first minimum r * of C n as a function of the reduced wave vector. In the low-frequency regime, it decays like r * = 0.449/q * in exact correspondence with the wavelength of the transverse plane wave. This means that C n is dominated by the collective dynamics of plane waves even in the presence of soft modes. However, the values of r * = ξ * 1 at q * 1 and r * = ξ * 2 at q * 2 , confirm the signature of a characteristic wavelength in the vibration modes. The origin of these lengths is not obvious. It was obvious from P (ω j ) that the eigenmodes are not simple plane waves at q * 1 , but have a very small participation ratio indicating the presence of isolated centers of enhanced vibrations. These centers are sufficiently few to not affect the long-range spatial correlations due to transverse plane waves in C n , but can affect the vibrational density of states through small frequency shifts. This description supports the fact that transverse plane waves and isolated centers of enhanced vibrations still coexist at q * 1 . ξ * 1 could be the distance between the isolated centers, that would correspond as well to the wavelength at q * 1 . At q * 2 a clear-cut change of behaviour appears in C n for all the systems studied indicating departure from transverse plane waves. This effect will be discussed again later. Finally, a common change appears at a larger reduced wave vector q * 3 = 11.7 nm -1 in C n , which was introduced in Fig. 3.5. It corresponds to the transition from transverse modes to longitudinal ones with bigger correlation radius.

We have shown in this Section, that the eigenmodes have characteristic features depending on the corresponding frequency range. In the low-frequency part of the spectrum, eigenmodes share common features independent of the bending rigidity of the modes: for example, the boson peak is bounded by two characteristic lengthscales ξ * 1 and ξ * 2 with a very low participation ratio in the first case, and a local maximum in the participation ratio in the second case. These two length-scales have a signature in the spatial correlation analysis of the modes. It is shown as well, that in the very low-frequency range, transverse plane waves coexist with local enhanced vibrations while the plane-wave character of the vibrations is questioned in the higher limit of the boson peak. We will now study the dynamical structure factor, in order to relate these observations to the study of density-density correlation functions, as can be tracked in neutron diffraction experiments for example [START_REF] Giordano | Fingerprints of order and disorder on the high-frequency dynamics of liquids[END_REF]Baldi et al. 2011b]. The analysis of the dynamical structure factor allows also to discuss the Ioffe-Regel criterion for waves scattering.

Dynamical structure factor

In this Section we analyze the dynamical structure factor in order to determine the dispersion law and the mean-free path for longitudinal and transverse phonons. We prove that it is an accurate method in the frequency range below the Ioffe-Regel criterion, where the mean-free path is still bigger than the half wavelength, and the notion of phonon dispersion is well-defined.

The dynamical structure factor is the self-correlation function of the mass currents [START_REF] Shintani | Universal link between the boson peak and transverse phonons in glass[END_REF] in the system at thermal equilibrium with some temperature T . The structure factor can be calculated by normal mode analysis. and various values of the wavenumber q (from left to right: 1.44, 2.89, 4.33, 5.77, 7.21, 8.66, 10.10 and 11.54 nm -1 ). The inset shows a full curve for q = 1.44 nm -1 .

Using a small displacement expansion of the density correlation function, combined with the projection of the displacements on the normal modes, and the classical approximation k B T ≫ ω for the equipartition of energy in the normal modes amplitude, it reads

S η (q, ω) = k B T m q 2 ω 2 F η (q, ω), (3.11)
where η denotes longitudinal (L) or transverse (T) component. In the above equation, F η (q, ω) is the longitudinal or transverse component of the Fourier transform of the eigenmodes

F L (q, ω) = N j=1 N i=1 q • u i (ω j )e iqR i 2 δ(ω -ω j ), (3.12
)

F T (q, ω) = N j=1 N i=1 q × u i (ω j )e iqR i 2 δ(ω -ω j ). (3.13)
Here q = q/|q| is a unit vector along q and u i (ω j ) is the displacement of the ith atom for jth eigenmode.

In order to calculate F η (q, ω), we also used the KPM (Appendix B). Fig. 3.7 shows the structure of eigenmodes F η (q, ω) in the reciprocal space averaged over all possible directions of q. For a better visual effect we divide F η (q, ω) by the magnitude of its maximum for each fixed value of ω. All color maps in Fig. 3.7 have two evident regions: low-frequency region with thin phonon branch, and highfrequency region without a certain relationship between the wavenumber q and the frequency ω.

In order to extract information about phonons in the low-frequency region we fit the structure factor S η (q, ω) using the DHO model (Fig. 3.8)

S η (q, ω) = A (ω 2 -ω 2 η (q)) 2 + ω 2 Γ 2 , η = L, T. (3.14)
We extract phonon dispersion ω η (q), the phonon inverse lifetime Γ(q) and a coefficient A from this fit.

The numerical results obtained by this method for phonon dispersion ω η (q) as well as the group velocity v η g = ∂ω η /∂q are presented in Fig. 3.9. With a known value of Γ and v g , the phonon mean-free path l(ω) can now be calculated as l(ω) = v g /Γ. The phonons are well-defined excitations if their mean-free path l(ω) exceeds the phonon half wavelength π/q (Ioffe-Regel criterion for phonons). Fig. 3.10 shows the value of Γ for all the samples with different bending rigidities Λ and Fig. 3.11 shows the position of the Ioffe-Regel frequency for longitudinal and transverse phonons for the parameter Λ = 21. From the similar figures for other values of the parameter Λ we find the remaining Ioffe-Regel frequencies (Table 3.2). Λ ν T IR , THz ν L IR , THz q T IR , nm -1 q L IR , nm Different comments are raised by these measurements. First, the sound wave velocities (Fig. 3.9b) are well-defined below the Ioffe-Regel criterion. It is not constant, but it varies with q. Transverse sound velocities show a small decay with q, starting at the low-frequency limit of the boson peak, as already measured in experiments [Baldi et al. 2011b]. Longitudinal sound velocities decay faster with a sudden increase at q ≈ 10 nm -1 (ν = 12.5 THz for Λ = 21), after transverse waves became strongly scattered in the sample.

The density of states of longitudinal and transverse phonons can be obtained from the dispersion laws q L (ω) and q T (ω) respectively

g L (ω) = L 3 6N π 2 q L (ω) 2 v L g (ω) , (3.15) g T (ω) = L 3 3N π 2 q T (ω) 2 v T g (ω)
. (3.16) and compared to the more general decomposition that was already discussed in Section 3.2. In the low-frequency limit, the ratio between them is

g L (ω)/g T (ω) = c 3 T /2c 3 L ≪ 1.
(3.17)

For Λ = 21 this ratio is equal to 0.057, which coincides well with the low-frequency part of this ratio is shown in the inset in Fig. 3.3a. The very low-frequency modes are naturally mainly transverse due to their lower sound velocity. The two estimations of longitudinal and transverse contribution to VDOS are compared in Fig. 3.3. The estimation from the dispersion law is close to the general estimation, but slightly lower. The main difference is on the boson peak Fig. 3.3b: the first low-frequency peak (attributed to soft modes in Section 3.2) is indeed completely absent in the dispersion law estimation of the VDOS that stays close to the Debye one in this frequency range. It means that this first peak results from a departure to the phonon-like behavior.

The inverse lifetime Γ (Fig. 3.10) is different for transverse and for longitudinal waves. The inverse lifetime of transverse waves varies approximately ∝ q 2 as discussed extensively in the literature [START_REF] Rufflé | Glassspecific behavior in the damping of acoustic-like vibrations[END_REF]], with a collapse for all Λ values at q * 2 . Longitudinal inverse lifetime is more sparse. It shows a sudden increase at q ≈ 10 nm -1 , that is after transverse waves are strongly scattered in the system and do not interfere anymore with longitudinal waves. The measurement of Γ allows computing a mean-free path l (Fig. 3.11) but only in the region where the sound velocity v g is well-defined. The computed longitudinal mean-free path is always larger than the transverse mean-free path. The measured values of l can overcome the size of the system, because it is computed from an estimation of the inverse lifetime Γ that does not result from a propagating process [START_REF] Damart | Nanocrystalline inclusions as a low-pass filter for thermal transport in a-Si[END_REF]] but only from a general fit of a geometrical function S(q, ω). The Ioffe-Regel crossover occurs at a well-defined wave vector q T IR = 2πν T IR /c T ≈ 7.5 ± 1.0 nm -1 for transverse waves, and q L IR = 2πν L IR /c L ≈ 9.7 ± 0.4 nm -1 for longitudinal waves (see Table 3.2), slightly larger than the upper limit q * 2 of the boson peak for the transverse one. This relation is only slightly sensitive to Λ, suggesting a universal mechanism for strong scattering in amorphous materials, independent of the specific interatomic interactions. We will now compare this estimation of the Ioffe-Regel criterion to the description of quasi-monochromatic wave packet propagation.

Diffusivity

In this Section we consider the diffusion of the vibrational energy. For this purpose, we excite a quasi-monochromatic wave packet in the middle thin layer of the sample around x = 0 in a small time interval around t = 0. In this Section we use a vector analog of the method described in Section 2.3.2. To excite vibrations in the sample we use the excitation force (3.18) where the phase ϕ iα is random for each atom i and each Cartesian projection α. The width of the excited layer is determined by the value of w = 3 nm and the duration of the excitation is given by τ exc = 1 ps. The latter determines the frequency resolution ∆ν = 1 THz. We start our calculations at time t 0 = -5τ exc when the external force is still negligible. In order to have a sufficiently large system size, the central sample with periodic boundary conditions and size L = 87 Å is duplicated into 4 images along x direction and 2 along y and z directions. As a result, in volume size we obtain a 16 times bigger sample. This allows a determination of large mean-free paths for phonons and diffusion of energy on longer distances. Indeed, the energy diffusion front reaches the boundaries of the original sample at t ≈ 1 ps when the excitation force is still active (see dashed lines in Fig. 3.12). The diffusivity of the vibrational energy in this extended sample is the same as in one big sample except a small region near the mobility edge. Localized modes with the localization length ξ > L look like delocalized in the repeated sample.

f ext iα (t) = sin(ωt + ϕ iα ) exp - x 2 i 2w 2 - t 2 2τ 2 exc
After applying the external force, the vibrational energy spreads in both directions from the central layer x = 0. The average radius squared of the energy diffusion front is defined as in Chapter 2 as:

R 2 (t) = 1 E tot i x 2 i E i (t). (3.19)
Here, x i is the x coordinate of the ith atom, E i (t) is the total energy of the ith atom, and the sum is taken over all atoms in the sample. E tot = i E i (t) is the total vibrational energy of the system. It is independent on time after the external force f ext iα (t) became negligibly small (i.e., for t > 5τ exc ). The energy of the ith atom E i (t) is the sum of the kinetic energy and a half of the potential energy of connected bonds with ith atom:

E i (t) = v i (t) 2 2 + 1 2 jαβ M iα,jβ u iα (t)u jβ (t). (3.20)
Here, v i (t) = ui (t) is the ith atom velocity with the same notations as Eq. 1.2.

Eq. (3.20) is the vector analog of Eq. (2.23).

The spatial vibrational energy distribution along the x direction is shown in Fig. 3.12 at different times t. Initial random phases ϕ iα allow keeping the center of mass of the energy in the central layer while the energy is progressively spread inside the sample.

To integrate the system with a given external force and zero initial conditions we used the Verlet method with a small enough time step δt = 0.6 fs and get the dependence R 2 (t) for different frequencies from ν = 2 THz up to ν = 20 THz. The results are shown in Fig. 3.14. We clearly see a linear temporal dependence in each curve. Their slope gives us the diffusivity by the equation for one-dimensional diffusion R 2 (t) = 2D(ω)t.

(3.21)

The resulting diffusivity is shown in Fig. 3.15a for different values of the parameter Λ. All curves have the same structure: (1) low-frequency modes with large diffusivity; (2) a flat region with relatively small diffusivity; (3) a prominent peak of the diffusivity; (4) a gradual decreasing of the diffusivity; ( 5) zero diffusivity for localized modes. The first two regions coincide for all values of Λ if we plot the rescaled diffusivity D/c T as a function of the reduced wave vector q * = ω/c T (Fig. 3.15b). After an initial decay, the diffusivity saturates at a minimum value. Whatever Λ, the flat region in the diffusivity occurs precisely between q * 2 and q * 3 , that is in the region between the upper bound of the boson peak (close to the Ioffe-Regel criterion for transverse waves) and transition from mostly transverse modes to mostly longitudinal ones. The lower boundary of the flat region is in perfect agreement with those obtained in Chapter 2 for a completely different random system. The relation between boson peak and Ioffe-Regel criterion was also suggested by experimental measurements [START_REF] Rufflé | Glassspecific behavior in the damping of acoustic-like vibrations[END_REF]] and molecular dynamics simulations in Lennard-Jones glasses [START_REF] Tanguy | Plastic response of a 2D Lennard-Jones amorphous solid: Detailed analysis of the local rearrangements at very slow strain rate[END_REF]]. It is shown here that the strong scattering gives rise to a very low diffusivity, and that it is possible to measure diffusivity in a purely harmonic model as soon as interactions are random. The shape of the instantaneous velocity field in the flat region is shown in Fig. 3.13 during the propagation of a wave packet. Rotational structures are clearly visible and responsible for the strong dephasing close to the Ioffe-Regel crossover. The flat region in the diffusivity is followed by a peak already discussed in Ref. [START_REF] Allen | Diffusons, locons and propagons: character of atomic vibrations in amorphous Si[END_REF].

The peak of the diffusivity is large for almost all values of the parameter Λ (Fig. 3.15). We can alternatively find the diffusivity of longitudinal and transverse phonons up to the Ioffe-Regel criteria using the approximate relation [START_REF] Kittel | Introduction to solid state physics[END_REF]]

D η (ω) = 1 3 l η (ω)v η (ω), η = L, T (3.22)
where ω < ω L IR for longitudinal phonons and ω < ω T IR for transverse phonons. We can see that the diffusivity of longitudinal and transverse phonons has monotonically decreasing behavior except the negligible peak for longitudinal phonons (Fig. 3.16). Eq. 3.22 cannot give the diffusivity beyond the Ioffe-Regel criteria, but we expect a small diffusivity decreasing down to 0 at the mobility edge. Therefore, the peak in the diffusivity cannot be explained by the diffusivity of longitudinal and transverse phonons separately. However, the total diffusivity depends on the ratio of density of states of longitudinal and transverse vibrations

D(ω) = g L (ω) g(ω) D L (ω) + g T (ω) g(ω) D T (ω). (3.23)
The resulting phononic diffusivity is shown in Fig. 3.16. It shows that the main peak located at q * ≈ 13.5 nm -1 is due to the large density of longitudinal modes g L (ω), enhancing the small diffusivity increase due to the absence of transverse modes in that frequency range. The rise of the diffusivity at 7 THz in amorphous silicon thus corresponds to the sharp change in the nature of vibrations from almost transverse to almost longitudinal ones having high sound velocity.

VDOS

Conclusion

In this Chapter we have proposed a coherent picture of the vibrational properties of harmonic amorphous solids with local tetrahedral order, by combining four independent approaches: the detailed study of the normal modes (resonant vibrational modes) and of the vibrational density of states, dynamic structure factor calculation and an analysis of propagation of a quasi-monochromatic wave packet. The bending rigidity of local interatomic bonds was used as a control parameter to tune the sound velocity. This allowed to get a coherent picture of the vibrational response of our model systems. Different regimes were highlighted. The results are summarized in Fig. 3.17.

The low-frequency vibrational response is dominated by transverse modes. In this region, the boson peak is visible and bonded by two characteristic wave vectors: the first is related to soft modes, and the second to the Ioffe-Regel limit for transverse waves. Remarkably, these two wave vectors are independent on the details of the interactions in the different systems studied here, and they define two characteristic mesoscopic length scales ξ * 1 and ξ * 2 having a signature in the spatial correlations of the normal modes. In silicon-like samples, the large difference in the transverse and longitudinal sound velocities yields a large gap between the Ioffe-Regel limit for transverse waves, and the Ioffe-Regel limit for longitudinal waves. In this gap, the vibrations sharply change the transverse character to longitudinal one near q * ≈ q * 3 , resulting in a deep increase of the diffusivity. As shown already by P.B. Allen at al. [START_REF] Allen | Diffusons, locons and propagons: character of atomic vibrations in amorphous Si[END_REF], the mobility edge and the transition to localized modes occur at higher frequencies. This transition was obtained as well in other disordered model materials, like lattice models [START_REF] Ludlam | Disorder-induced vibrational localization[END_REF]] and models of amorphous silica [START_REF] Ludlam | Universal features of localized eigenstates in disordered systems[END_REF], thus supporting its universal feature. The modes in the boson peak range preceding the Ioffe-Regel crossover for transverse waves have a characteristic random rotational structure yielding a dephasing of the wave front. In this frequency range, the participation ratio is at a maximum value, but the diffusivity is at a minimum value, and the diffusive inverse lifetime is proportional to q 2 . Similar results have already been obtained with Brillouin scattering measurements [START_REF] Rufflé | Glassspecific behavior in the damping of acoustic-like vibrations[END_REF]]. However, the precise sensitivity of the inverse lifetimes or of the vibrational density of states to the rescaled wave vector (or equivalently to the frequency) is system dependent [START_REF] Larkin | Thermal conductivity accumulation in amorphous silica and amorphous silicon[END_REF].

Chapter 4

Random matrix theory approach to the jamming transition

In this Chapter we show how one can apply the random matrix theory to estimate the vibrational density of states (VDOS) in jammed granular systems.

Model

Granular media, various emulsions (microdrops of one liquid in another immiscible liquid), and colloid suspensions (solid particles in the liquid) are widespread in nature, industry, and daily life. Such media demonstrate a wide variety of phenomena still poorly studied. In particular, they can flow as a liquid or have elastic properties as a solid, depending on external factors.

A transition between a solid phase, where all granules (microdrops, etc.) touch each other, and a phase of free particles is called the jamming transition (after traffic jam). Such a transition is described by a simple model, where N elastic granules enclosed in a certain volume are considered [START_REF] Liu | Nonlinear dynamics: Jamming is not just cool any more[END_REF]]. The most important parameter of this model is the ratio φ of the volume occupied by granules to the entire available volume. If φ is large and exceeds a certain critical value φ c , all granules touch each other and constitute something like a solid whose structure can withstand finite external loads (Fig. 4.1). In order to avoid the crystallization of the system, a mixture of granules with slightly different dimensions is taken. If φ < φ c , granules no longer touch each other and the system behaves as a gas. At φ = φ c , all granules touch each other, but the interaction between them is absent. In this work, we demonstrate how the random matrix theory allows the description of the VDOS in the solid phase, when φ is slightly larger than the critical value φ c . a b c The case where granules are spherical and friction between them is absent is of most interest. Despite seeming simplicity, this model is sufficient for qualitative description of the jamming transition. Such a model is mathematically described by a repulsive potential between each pair of granules touching each other [START_REF] Liu | Nonlinear dynamics: Jamming is not just cool any more[END_REF]]:

U (r ij ) ∝ (1 -r ij /σ ij ) γ , r ij < σ ij , U (r ij ) = 0, r ij > σ ij . (4.1)
Here, r ij is the distance between the centers of the ith and jth granules and σ ij is the sum of the radii of these granules. The exponent γ depends on the type of the interaction between granules. The harmonic potential (γ = 2) and Hertz potential (γ = 5/2), which corresponds to the interaction between three-dimensional elastic balls, are often used. Since we consider spherical granules in the absence of friction, rotational degrees of freedom of individual granules are of no significance.

We introduce the average coordination number z, i.e., the average number of contacts of each granule with neighboring granules. The larger the concentration of granules φ, the larger the number of the neighbors interacting with each granule. At the point of jamming transition φ = φ c , the average number of contacts is determined by the universal formula z c = 2d, where d is the dimension of the space, which is related to the Maxwell rule of counting bonds [START_REF] O'hern | Jamming at zero temperature and zero applied stress: The epitome of disorder[END_REF][START_REF] Maxwell | On the calculation of the equilibrium and stiffness of frames[END_REF]].

As the average number of contacts z decreases to the critical value 2d, various characteristics of the system satisfy a power law. In particular, the bulk modulus G and shear modulus B behave as follows [START_REF] O'hern | Jamming at zero temperature and zero applied stress: The epitome of disorder[END_REF]]:

G ∼ (z -z c ) 2γ-3 , B ∼ (z -z c ) 2γ-4 . (4.2)
We are interested in the behavior of the VDOS g(ω) as a function of the frequency ω. Numerical simulation shows that two characteristic frequency ranges can be identified: the range ω -< ω < ω + , where the VDOS is approximately constant, and the range 0 < ω < ω -with a relatively small number of oscillations (or a gap in the spectrum). The closer z to the critical value 2d, the closer is ω -to 0. At the critical value z = 2d, the VDOS g(ω) is approximately constant beginning with zero frequency. There is no simple theoretical explanation of such a behavior of the VDOS near the stability threshold. We will show that the random matrix theory can provide an adequate estimate of the VDOS near the jamming transition.

We consider the system close to a critical one when z > z c (φ > φ c ). In this case, granules touch its neighbors, but are insignificantly indented into each other. Then, the total potential energy is expanded near the stable equilibrium position as [START_REF] Wyart | Effects of compression on the vibrational modes of marginally jammed solids[END_REF]]

U (u 1 , . . . , u N ) = (ij) k ij 2 (u i -u j ) • n ij 2 . (4.3)
Here, u i is the displacement of the ith granule from the equilibrium position r

i , (ij) under the sum sign means summation only over the touching pairs of the ith and jth granules, and n ij is the unit vector along the direction connecting the centers of these granules r ir j . Owing to the repulsive potential under consideration, all k ij values are positive, k ij > 0. For certainty, we accept that the total number of pairs of touching granules is K = zN/2. The above formulas are valid for any dimensionality of the space d.

The dynamical matrix M is defined in terms of the second derivatives with respect to the potential energy of the system at the equilibrium position: .4) where α and β mean the projections of the displacements of granules on the Cartesian coordinates. The dynamical matrix has a dimension of N f × N f , where N f = N d is the number of degrees of freedom. In this case, the eigenvalues of the dynamical matrix are squares of the eigenfrequencies of the mechanical system under consideration. The eigenfrequencies include trivial zero frequencies corresponding to the translational and rotational motions of the system as a whole.

M iα,jβ = 1 √ m i m j ∂ 2 U ∂u iα ∂u jβ . ( 4

Decomposition of the dynamical matrix

Before the consideration of the general case of many interacting particles, we consider illustrative cases of two and three particles interacting with each other without the participation of other particles (Fig. 4.2). The potential energy of two particles (Fig. 4.2a) appears as a single term in Eq. ( 4.3): This energy is just the interaction energy of two particles connected by a spring with the longitudinal rigidity k. The corresponding dynamical matrix M (2d × 2d) can be written in the block form

U (u 1 , u 2 ) = k 2 (u 1 -u 2 ) • n 2 . ( 4
M = kn/m 1 -kn/ √ m 1 m 2 -kn/ √ m 1 m 2 kn/m 2 , (4.6)
where n is a d × d matrix with the elements nαβ = n α n β . We note that such a dynamical matrix can be represented in the form M = AA T , where A is the following 2d × 1 matrix:

A = k/m 1 n -k/m 2 n . ( 4.7)
The elements of the matrix A have the dimension of frequency and correspond to oscillations of the masses m 1 and m 2 connected by a spring with the rigidity k.

For three particles, the first of which touches the second and the second touches the third (Fig. 4.2b), the 3d × 3d dynamical matrix has the form

M =         k 12 n12 m 1 -k 12 n12 √ m 1 m 2 0 -k 12 n12 √ m 1 m 2 k 12 n12 + k 23 n23 m 2 -k 23 n23 √ m 2 m 3 0 -k 23 n23 √ m 2 m 3 k 23 n23 m 3         . ( 4.8)
This dynamical matrix is also represented in the form M = AA T , where the matrix A now has two columns according to two contacts and 3d rows according to 3d degrees of freedom existing in the system:

A =   k 12 /m 1 n 12 0 -k 12 /m 2 n 12 k 23 /m 2 n 23 0 -k 23 /m 3 n 23   . (4.9)
The positions of minus signs are arbitrary: each column of the matrix A can be multiplied by -1 without a change in the dynamical matrix M .

In the general case of a large number of interacting granules, the elements of the dynamical matrix M have the form

M iα,jβ = - k ij n α ij n β ij √ m i m j , i = j, (4.10) M iα,iβ = - j =i m j m i M iα,jβ . (4.11)
By analogy with the above examples, it can be represented in the form M = AA T , where A is a rectangular N f × K matrix. Here, as above, N f = N d is the number of degrees of freedom and K = zN/2 is the total number of pairs of interacting granules. The elements of the matrix A have the form

A iα,p = k p m i n pα (δ p 1 i -δ p 2 i ), (4.12)
where the subscript p enumerates pairs of touching granules and p 1 and p 2 are the numbers of granules contained in the pth pair. As a result, each row of the matrix A corresponds to a certain degree of freedom and each column corresponds to a certain pair of interacting granules. Furthermore, the representation AA T guarantees the stability of the mechanical system, because the matrix M = AA T is always positive definite for any rectangular real-valued matrix A [START_REF] Bhatia | Positive definite matrices[END_REF]].

We note that the eigenvalues of the dynamical matrix M do not change when the matrix A is multiplied from the left and right by arbitrary orthogonal matrices U and V , respectively. In other words, the matrix

U ij = V ij = 0, (4.13) U ij 1 U ij 2 = U j 1 i U j 2 i = δ j 1 j 2 1 N f , (4.14) V ij 1 V ij 2 = V j 1 i V j 2 i = δ j 1 j 2 1 K , (4.15)
because an individual column or row of a random orthogonal matrix is a randomly oriented unit vector. Consequently, the elements of the matrix

∼ A have the simple properties ∼ A ij = 0, ∼ A 2 ij = 1 N f K kl A 2 kl , (4.16)
Thus, all elements of the matrix for various filling densities φ. The solid lines are predictions of the random matrix theory and the dashed lines present the numerical data from [START_REF] Wyart | Effects of compression on the vibrational modes of marginally jammed solids[END_REF].

between the nearest granules. In this case, by definition, the matrices M = AA T and

∼ M = ∼ A ∼
A T have the same set of eigenvalues. We note that the elements of the matrix ∼ A have certain correlations. However, for simplicity and universality of the estimate of the distribution of eigenvalues, we neglect these correlations. Thus, we accept below that the matrix ∼ A is a N f ×K rectangular random matrix with the independent elements having the properties

∼ A ij = 0, ∼ A 2 ij = ω 2 0 N f . (4.17)
Here, ω 0 is the characteristic frequency of oscillations of touching particles. It is determined from Eqs. (4.12) and (4.16):

ω 2 0 = 1 K (ij) k ij 1 m i + 1 m j . ( 4
g(ω) = 1 πω 2 0 ω (ω 2 + -ω 2 )(ω 2 -ω 2 -), ω -≤ ω ≤ ω + (4.19)
where

ω ± = ω 0 K N f ± 1 , K N f = z 2d . (4.20)
It is seen that the value z = 2d is specific: in this case, the number of degrees of freedom N f is equal to the number of touching pairs K, the matrix A becomes square, and the VDOS has the form of a quarter of a circle. For values z > 2d, a gap in the VDOS appears in the frequency range 0 < ω < ω -. In this case, near the critical value 2d,

ω -≈ ω 0 |z -2d| 4d , ω + ≈ 2ω 0 . (4.21)
However, it is noteworthy that the case z < 2d in the model under consideration with elastic balls is not implemented, because the system at z < 2d completely breaks up into noninteracting granules and the average coordination number is z = 0.

The vibrational density of states

The VDOS obtained numerically in [START_REF] Wyart | Effects of compression on the vibrational modes of marginally jammed solids[END_REF]] for elastic spheres with the repulsive potential with γ = 2 is shown in Fig. 4.3 in comparison with the estimate by Eq. ( 4.19). To recalculate the density of filling φ to the average coordination number z, we used the relation z -2d = 7.5φφ c , which is quite accurately satisfied in the studied range of φφ c [START_REF] O'hern | Jamming at zero temperature and zero applied stress: The epitome of disorder[END_REF]]. We note that the characteristic frequency ω 0 can also depend on the difference z -2d. Indeed, the characteristic mutual penetration of particles δ is proportional to φφ c . For this reason, the characteristic rigidity of bonds, which is defined in terms of the second derivative of potential (4.1), is k ∝ δ γ-2 ∝ (φφ c ) γ-2 ∝ (z -2d) 2γ-4 . Correspondingly, the characteristic frequency is ω 2 0 = k/m ∝ (z -2d) 2γ-4 . However, for the case γ = 2, the frequency ω 0 is constant and the value ω 0 = 1.71 (in the same units as in the numerical calculation) was used for comparison with the numerical experiment.

The approach under consideration based on the random matrix theory appropriately predicts both an almost constant VDOS from ω -∼ z -2d up to the maximum frequency of the system and the almost complete absence of vibrational modes at ω < ω -. In a real system at ω < ω -, a number of vibrational states exist instead of the strict gap. The reason is that acoustic phonons are disregarded in the above approach. However, as is seen in the numerical experiment, the number of acoustic phonons rapidly decreases when z approaches the critical value 2d. Acoustic phonons in the gap were studied in Chapter 2. We note that the value ω 0 = 1.71 was the only fitting parameter (common for all curves). It simultaneously provides the correct frequencies ω -and an approximately constant density of states g(ω) ≈ 2/πω 0 = 0.37 in the range ω -< ω < ω + .

To conclude, we note that, strictly speaking, Eq. (4.3) is not exact. A more accurate formula has the form (4.22) where (u iu j ) ⊥ is the projection of the difference of the displacements u iu j on the plane perpendicular to the vector n ij . For the repulsive potential under consideration, k ij > 0 and e ij < 0. Thus, k ij and e ij make the stabilizing and destabilizing contributions to the potential, respectively. The destabilizing component is noticeable in the determination of the equilibrium position, but the ratio |e ij /k ij | at the equilibrium position is proportional to the mutual penetration of particles δ. It is equal to 0 at z = z c and, at z > z c , makes an insignificant change in the density of vibrational states, leading to a smoothing of the density of states at ω ≈ ω - [START_REF] Wyart | Effects of compression on the vibrational modes of marginally jammed solids[END_REF]].

U (u 1 , . . . , u N ) = (ij) k ij 2 (u i -u j ) • n ij 2 + e ij 2 (u i -u j ) ⊥ 2 ,

Conclusion

In this Chapter we have shown that the interaction between elastic granules in a granular system is described by the dynamical matrix M = AA T . Each row of the matrix A corresponds to a certain degree of freedom and each column corresponds to the elastic interaction between a certain pair of neighboring granules. The representation of the dynamical matrix in the form M = AA T , together with the random orthogonal transformation, allows to describe qualitatively the density of vibrational states making use the Wishart ensemble. In this case, the only significant parameter is the ratio of the total number of contacts K to the total number of degrees of freedom N f . The characteristic frequency ω 0 specifies only the scale of all frequencies and trivially appears in all formulas. If the total number of contacts K differs from the total number of degrees of freedom N f , the density of vibrational states has a gap whose width is proportional to K -N f . In reality, this gap in the vibrational spectrum is not perfect and contains a small number of acoustical phonons (see Chapter 2, where such a soft gap was called phonon gap). At K = N f , the gap is closed and the density of vibrational states is approximately constant starting from zero frequency. The results of this work are in agreement with the theoretical calculations performed by the effective medium method [START_REF] Wyart | Scaling of phononic transport with connectivity in amorphous solids[END_REF]].

General conclusion

In this work we have studied different vibrational properties of different amorphous systems, including a stable random matrix model (Chapters 1 and 2), a numerical model of amorphous silicon (Chapter 3), and a granular jammed solid (Chapter 4).

In Chapters 1 and 2 we have developed a stable random matrix approach to describe vibrations in strongly disordered systems. 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 [START_REF] Schirmacher | Harmonic vibrational excitations in disordered solids and the "boson peak[END_REF][START_REF] Taraskin | Origin of the boson peak in systems with lattice disorder[END_REF][START_REF] Grigera | Vibrations in glasses and Euclidean random matrix theory[END_REF] 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 . Since matrices AA T and M 0 are positive definite, such form of the dynamical matrix guarantees the mechanical stability of the system for any positive value of µ. The first term AA T is responsible for the disorder in the system, and the second term µM 0 describes the ordered part of the dynamical matrix. 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 µ 1 and fluctuating part of the dynamical matrix is bigger than the ordered part. In this case the Young modulus of the lattice E ∝ √ µ. As we have shown, this and other scaling relations map directly onto the scaling observed in jammed packings near the isostatic point. For example, in scaling relations the parameter µ plays the same role as the parameter ∆φ 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 acoustic phonons (plane waves) which can be characterized by frequency ω and wave vector q. However, with increasing of ω, due to the disorderinduced scattering, the phonon linewidth ∆ω increases rapidly as ∆ω ∝ ω 4 and at some frequency ω ≈ ω ir the phonon mean-free path l becomes of the order of the wavelength λ. Though this crossover is not sharp and has no critical behavior at ω = ω ir , the structure of the eigenmodes at higher frequencies quite soon becomes very different from the plane waves. As a result, at higher frequencies the original notion of acoustic phonons is lost and delocalized vibrational modes have a diffusive nature. They are similar to diffusons introduced by Allen and Feldman et al. [START_REF] Allen | Diffusons, locons and propagons: character of atomic vibrations in amorphous Si[END_REF]]. The diffusons again can be characterized by frequency ω, but have no well-defined wavevector 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 latter has a broad maximum as a function of q at q = ω/D u , where D u ≃ Ωa 2 0 is a diffusion coefficient of particle displacements. The vibrational line width is Γ(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 [START_REF] Sette | Dynamics of glasses and glass-forming liquids studied by inelastic X-ray scattering[END_REF][START_REF] Ruocco | High-frequency vibrational dynamics in glasses[END_REF]] and references therein. It was also found in molecular dynamics simulations of amorphous silicon [START_REF] Christie | Vibrational behavior of a realistic amorphous-silicon model[END_REF]].

The crossover between acoustic 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 [START_REF] Rufflé | Glassspecific behavior in the damping of acoustic-like vibrations[END_REF]], densified vitreous silica [START_REF] Rufflé | Observation of the onset of strong scattering on high frequency acoustic phonons in densified silica glass[END_REF]], vitreous silica [START_REF] Baldi | Sound attenuation at terahertz frequencies and the boson peak of vitreous silica[END_REF]Baldi et al. 2011b;Baldi et al. 2011a], glassy sorbitol [START_REF] Ruta | Communication: High-frequency acoustic excitations and boson peak in glasses: A study of their temperature dependence[END_REF]] and glycerol glass [START_REF] Monaco | Breakdown of the Debye approximation for the acoustic modes with nanometric wavelengths in glasses[END_REF]. The crossover frequency was found to be close to the boson peak position.

The experimental data and molecular dynamics simulations show that the diffusion coefficient of particle displacements D u and the diffusivity of energy D(ω) are of the same order D u ∼ D(ω) ∼ 1 mm 2 /sec in vitreous silica [Baldi et al. 2011b;[START_REF] Feldman | Realistic model calculations based on the Kubo theory for the thermal conductivity of amorphous insulators[END_REF][START_REF] Yu | Thermal conductivity computed for vitreous silica and methyl-doped silica above the plateau[END_REF], amorphous silicon [START_REF] Christie | Vibrational behavior of a realistic amorphous-silicon model[END_REF][START_REF] Allen | Diffusons, locons and propagons: character of atomic vibrations in amorphous Si[END_REF] and glycerol [START_REF] Ruocco | Nondynamic origin of the high-frequency acoustic attenuation in glasses[END_REF]. Our results for amorphous silicon are discussed in Chapter 3.

Since ω ir ∝ √ µ, we can vary the Ioffe-Regel crossover frequency and, therefore, the relative number of acoustic phonons N ph in the system, changing the parameter µ. It is zero when µ = 0 and there are no acoustic phonons in the lattice. In this case all delocalized vibrations are diffusons. If 0 < µ ≪ 1 we have acoustic 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 the silica glass we can estimate the relative number of acoustic phonons from the data [START_REF] Taraskin | Ioffe-Regel crossover for plane-wave vibrational excitations in vitreous silica[END_REF]. The Ioffe-Regel crossover frequency was estimated to be ν IR = 1 THz, and integrating density of states 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 acoustic phonons. As follows from Table 2.1 it corresponds to very small values of µ < 0.01. It means that small amount of acoustic phonons in disordered systems is a signature of a strong disorder.

We show that in the random matrix model functions g(ω) and D(ω) are approximately constant in some frequency interval (ω ir ω ω loc ), then we find that approximately κ(T ) ∝ T in the corresponding temperature range [START_REF] Xu | Energy transport in jammed sphere packings[END_REF]]. It explains a quasi-linear temperature dependence of the thermal conductivity above the plateau observed in glasses [START_REF] Cahill | Thermal conductivity of amorphous solids above the plateau[END_REF]. 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 a clear explanation for the temperature dependence of the thermal conductivity of glasses and other disordered systems above the plateau.

In Chapter 3 we have considered a numerical model of amorphous silicon. We have shown that properties of transverse and longitudinal vibrations in siliconlike amorphous materials are sufficiently different. It leads to a more complicated structure of the vibrational spectrum (Fig. 3.17). The low-frequency part of the VDOS has a dominant transverse character up to 7 THz. In this frequency region, there are the wide boson peak [START_REF] Beltukov | Density of states in random lattices with translational invariance[END_REF][START_REF] Beltukov | Diffusion of vibrations in disordered systems[END_REF][START_REF] Beltukov | Ioffe-Regel criterion and diffusion of vibrations in random lattices[END_REF][START_REF] Beltukov | Vibrations in amorphous solids beyond the Ioffe-Regel criterion[END_REF] and the transverse Ioffe-Regel frequency (4.5 THz). Modes in the frequency range 4.5-7 THz have a constant diffusivity and the linewidth Γ ∼ q 2 . Therefore, transverse diffusons dominate in this frequency range. At 7 THz there is a sharp transition from mostly transverse modes to mostly longitudinal ones. However, the longitudinal Ioffe-Regel frequency is much larger (13 THz). It explains the prominent peak in the diffusivity at ≈8 THz.

In Chapter 4 we have shown that dynamical matrix of the granular jamming system can be presented in the form M = AA T where A is a rectangular N f × K matrix. Here N f is the number of degrees of freedom and K is the number of elastic contacts between particles. We have shown that one can use the Wishart ensemble to qualitatively describe the density of vibrational states. In this case, the only significant parameter is the ratio K/N f . If the total number of contacts K differs from the total number of degrees of freedom N f , the density of vibrational states has a gap whose width is proportional to K/N f -1. In numerical experiments, this gap in the vibrational spectrum is not perfect and contains a small number of phonons. If K = N f , the gap is closed and the density of vibrational states is approximately constant starting from zero frequency. Thus, granular jamming systems are similar Here v i (t) is the instantaneous velocity of the ith atom and m is the atomic mass (all atomic masses are supposed to be the same). This method is much faster than a numerical diagonalization of the matrix M . However, it requires accurate integration of the equations of motion up to a large enough time t max with small enough time step δt ≪ 1/ω max and low temperature k B T < ∆E. Here ω max is the maximum frequency in the system, and ∆E the smallest energy barrier surrounding the referred equilibrium position. The resulting frequency resolution of the density of states g cvv (ω) in this method is δω ∼ 1/t max .

The kernel polynomial method [START_REF] Weiße | The kernel polynomial method[END_REF]] (KPM) is an alternative way to compute the VDOS for large systems. It is a more accurate and much faster method in comparison with the previous one, as will be discussed below. It makes it possible to find the VDOS using Eq. (3.4) with δ-function replaced by a series of polynomials. KPM was introduced in [START_REF] Silver | Calculation of densities of states and spectral functions by Chebyshev recursion and maximum entropy[END_REF] and detailed reviewed in [START_REF] Weiße | The kernel polynomial method[END_REF]]. Originally it was used for finding electronic DOS in disordered systems. It allows with controlled accuracy getting directly the distribution of the eigenvalues of some large matrix M , not calculating the eigenvalues itself. We will show how KPM can be adopted to find the VDOS, i.e. the distribution of the square roots of the eigenvalues of the dynamical matrix M without its exact diagonalization. In this method we use only moments of this matrix up to sufficiently high order which is controlled by the accuracy of the calculations. Below we shortly describe the KPM for our problem.

All eigenvalues ω 2 j of the matrix M are non-negative due to mechanical stability of the system and lie in some interval [0, ω 2 max ]. Usually, the precise value of the maximum frequency is unknown so ω max is an estimation of the maximum frequency which guarantees that ω j < ω max for all ω j . Let us introduce new dimensionless variable ε = 1 -2ω 2 /ω 2 max in order to rescale all eigenfrequencies squared ω 2 j to the interval [-1, 1] for variable ε j . Thus, we can transform Eq. (3.4) as

g(ω) = 4ω N ω 2 max N j=1 δ(ε -ε j ) (B.2)
where ε i = 1 -2ω 2 i /ω 2 max are eigenvalues of the matrix ∼ M = I -2M/ω 2 max where I is the unit matrix.

For -1 < ε < 1 and -1 < ε j < 1 we can expand the δ-function in Eq. (B.2) in second kind Chebyshev polynomial series .3)

δ(ε -ε j ) = 2 π √ 1 -ε 2 ∞ k=0 U k (ε)U k (ε j ). (B

Figure 1 .

 1 Figure 1. A two-dimensional illustration of atomic order in crystalline (a) and amorphous (b) SiO 2 . Lines show the direction of chemical bonds, dots and circles are Si and O atoms respectively [Kittel 1949].

Figure 2 .

 2 Figure2. The thermal conductivity of crystalline and amorphous SiO 2 . Filled symbols are the first results of Eucken[START_REF] Eucken | Über die Temperaturabhängigkeit der Wärmeleitfähigkeit fester Nichtmetalle[END_REF]]. Open symbols are more recent Pohl's results[START_REF] Zeller | Thermal conductivity and specific heat of noncrystalline solids[END_REF][START_REF] Raychaudhuri | Thermal conductivity of neutronirradiated silica[END_REF][START_REF] Vandersande | Simple apparatus for the measurement of thermal diffusivity between 80-500 K using the modified Ångström method[END_REF][START_REF] Cahill | Thermal conductivity of amorphous solids above the plateau[END_REF][START_REF] Cahill | Lattice vibrations and heat transport in crystals and glasses[END_REF]. The black bottom line is the existing theory of the thermal conductivity of amorphous solids[START_REF] Buchenau | Interaction of soft modes and sound waves in glasses[END_REF]]. The straight red line shows the linear dependence (it will be obtained in Section 2.4).

Figure 3 .

 3 Figure 3.Comparison of the thermal conductivity of different amorphous dielectrics: GeO 2 , Se, SiO 2 , solid varnish GE 7031, polymer PMMA[START_REF] Zeller | Thermal conductivity and specific heat of noncrystalline solids[END_REF].

Figure 1 . 1 .

 11 Figure 1.1. VDOS of sparse 1000 × 1000 dynamical matrices for different values of n. The line n = ∞ is the theoretical prediction (1.13). The frequency has units of Ω.

Figure 1 . 2 .

 12 Figure 1.2. Tree topology.

Figure 1

 1 Figure 1.3.A simple cubic lattice.

Figure 1 . 4 .

 14 Figure 1.4. Schematic diagram illustrating the interaction of atoms in a cubic lattice. Shown are atoms interacting with the central (black) atom with random stiffness. Different colors mark random bonds with different stiffness distribution.In total, the central atom interacts with 24 surrounding atoms (nearest and next nearest neighbors).

Figure 1 . 5 .

 15 Figure 1.5. VDOS for 1d, 2d, 3d, and 6d simple cubic lattices. The system size is N = 1000 1 , 32 2 , 10 3 , 3 6 respectively. The frequency has units of Ω.

Figure 1 . 6 .

 16 Figure 1.6. The spatial eigenmode structure of random matrix M = AA T for the lowest frequency ω min in two-dimensional square lattice 400 × 400.

Figure 1 . 7 .

 17 Figure 1.7. The numerical calculation of the participation ratio for cubic lattices in different dimensions d. a) N ≈ 1000 atoms (1000 1 , 32 2 , 10 3 , 3 6 ) b) N ≈ 16000 atoms (16000 1 ,126 2 , 25 3 , 5 6). The black horizontal line corresponds to the theoretical value P (ω) = 1/3. The frequency has units of Ω.

Figure 1 .

 1 Figure 1.8.The numerical calculations of the participation ratio of the amorphous SiO 2 using molecular dynamics methods[START_REF] Jin | Dynamic structure factor and vibrational properties of SiO 2 glass[END_REF]].

Figure 1 .

 1 Figure 1.11. A schematic representation of the Young modulus measurement for a system with N = 5 3 atoms.

Figure 1 . 13 .

 113 Figure1.13. The distribution of Young modulus (in units of mΩ 2 /a 0 ) for 10 5 different realizations of the system. The results for different system sizes are almost indistinguishable from each others. The black line is the exponential distribution(1.41).

Figure 1 .

 1 Figure 1.14. a) Non-affine displacements of the surface atoms (in arbitrary units) for N = 20 3 sample. b) The distribution of displacements.

Figure 1 .

 1 Figure 1.15. A schematic diagram illustrating the interaction of atoms in a square lattice. It is a two-dimensional analog of Fig. 1.4.

Figure 1 .

 1 Figure 1.18. Young modulus E (in units of mΩ 2 /a 0 ) as a function of µ for dynamical matrix M = AA T + µM 0 built on a simple cubic lattice with N = 100 × 100 × 100 particles (one realization). Black dots are calculated values, and the line is the best least-square fit.

Figure 1 .

 1 Figure 1.19. The normalized VDOS g(ω) for dynamical matrix M = AA T + µM 0 and five different µ (0, 0.001, 0.01, 0.1, 1) calculated with precise numerical KPM solution for a simple cubic lattice with N = 200 3 (full lines). Straight lines correspond to the Debye prediction (1.35). Inset: dependence ω max (µ) ∝ √ µ. The frequency has units of Ω.

Figure 1 .

 1 Figure 1.20. Participation ratio for different µ as a function of ω (in units of Ω) for N = 27 3 (one realization). The arrows indicate positions of ω max in g(ω) for corresponding values of µ (see Fig. 1.19).

Figure 1 .

 1 Figure 1.21. The normalized VDOS g(ω) for dynamical matrix M = AA T + µM 0 and different µ (0, 0.001, 0.01, 0.1, 1) calculated with precise numerical KPM solution for a simple cubic lattice with N = 200 3 (full lines). The sum rule (1.14) is violated. Inset: dependence ω max (µ) ∝ √ µ. The frequency ω has units of Ω.

Figure 1 .

 1 Figure 1.22. The normalized VDOS g(ω) for dynamical matrix M = AA T + µM 0 with µ = 1 and different percentage 100%p of cut out springs calculated with precise numerical KPM solution for a simple cubic lattice with N = 200 3 (full lines). Straight lines are calculated according to Eq. (1.35) with sound velocity v = E/ρ 0 . The Young modulus E is calculated in the same way as in Section 1.8.

Figure 1 .

 1 Figure 1.23. The normalized VDOS g(ω) for dynamical matrix M = AA T + βBB T with different β calculated with precise numerical KPM solution for simple cubic lattice with N = 100 3 (full lines). Straight lines are calculated according to Eq. (1.35) with sound velocity v = E/ρ 0 . The Young modulus E is calculated in the same way as in Section 1.8.

Figure 2 .

 2 Figure 2.1.The Lorentz dispersion curves for different wave vectors q 100 direction and µ = 0.1. Closed diamonds correspond to the calculated values of S(q, ω) and lines are fitting curves according to Eq.(2.4). The number of atoms N = 50 3 (one realization). Insets: the Lorentz dispersion curves for q = 0.5 and q = 0.75. The frequency ω has units of Ω, wavevector q has units of 1/a 0 .

  2.14).

Figure 2 .

 2 Figure2.2. The dependence ω q on q for q 100 direction for various µ (1, 0.1, 0.01) in a cubic sample with N = 50 3 (one realization). Filled and open diamonds are the maximums of S(q, ω) as a function of ω for each discrete value of q n for frequencies below and above the Ioffe-Regel crossover correspondingly (see text below for details). Solid lines correspond to halves of the maximums. Dashed lines show ω = vq linear dependence with sound velocity v = E/ρ 0 . Horizontal dotted lines correspond to the maximum frequency ω max in g(ω) (taken from Fig.1.19). Insets show the group velocity v g = dω/dq as a function of ω. The frequency ω has units of Ω units, wavevector q has units of 1/a 0 .

Figure 2 . 3 .

 23 Figure 2.3. The phonon linewidth ∆ω as a function of ω (in units of Ω) for different µ in a cubic sample with N = 50 3 (one realization). Different symbols correspond to different q directions.for q 100 , △ for q 110 , for q 111 . Filled and open symbols refer to excitations below and above the Ioffe-Regel crossover frequency ω ir correspondingly.

Figure 2 . 4 .

 24 Figure 2.4. The ratio l(ω)/λ as a function of ω (in units of Ω) for different µ. Different symbols correspond to different q directions as explained in Fig. 2.3. The full horizontal line (separating filled and open symbols) corresponds to Ioffe-Regel crossover l(ω) = λ/2.

Figure 2 . 5 .

 25 Figure 2.5. The reduced VDOS g(ω)/ω 2 is a function of ω/ω b where ω b is the boson peak position.

Figure 2 . 6 .

 26 Figure 2.6. The normalized VDOS g(ω) for dynamical matrix M = AA T + µM 0 and five different µ (0, 0.001, 0.01, 0.1, 1) calculated with precise numerical KPM solution for a simple cubic lattice with N = 200 3 (full lines). Straight lines correspond to the Debye prediction(1.35). Filled and open diamonds correspond to phonon contribution to the VDOS below and above the Ioffe-Regel crossover frequency ω ir correspondingly. The frequency has units of Ω.

Figure 2 . 7 .

 27 Figure 2.7.The displacement structure factor S(q, ω), Eq. (2.1) (symbols) for µ = 0 and for three different frequencies. Different symbols correspond to different q directions. for q 100 , △ for q 110 , for q 111 . Full lines correspond to the structure factor S rw (q, ω) of the random walk on the lattice given by Eq. (2.9) with D rw = 0.7. Dashed line corresponds to the limit q ≪ 1 (see Eq. (2.12)). The frequency has units of Ω, the wavevector has units of 1/a 0 .

Figure 2 . 8 .

 28 Figure 2.8. The correlation function C(r, ω) for µ = 0 and six different frequencies (ω/Ω = 0.14, 0.31, 0.49, 0.66, 0.84, 1.01) for sample with N = 50 3 atoms. The full lines are our numerical results obtained from Eq. (2.2). Each line starts from r = r min which is about 2.5 interatomic distances (marked by arrows). The dashed line corresponds to Eq. (2.16) with D rw = 0.7Ωa 2 0 . The distance r has units of a 0

  .18) At ω = 0 the correlation function(2.16) decays slowly as 1/r. In disordered systems at critical point the correlation function decays as C(r) ∝ 1/r d-D 2 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.

Figure 2 .

 2 Figure2.10. The same normalized structure factor S n (q, ω) as in Fig.2.9 but in q space in plane q x q y (q z = 0) for ω = 0.5Ω. The left picture corresponds to µ = 0 (a) and the right to µ = 0.1 (b).

Figure 2 .

 2 Figure 2.13. Black points (diamonds and triangles) show the distribution of energy E(x, t) contained in the layer x as a function of x for two different frequencies ω = 2Ω and ω = 6Ω at times t = 234/Ω and t = 900/Ω, respectively, calculated numerically with Newton method. Full lines are theoretical predictions for delocalized (diffusive) and localized modes given by Eqs. (2.25, 2.26) with R 2 ≈ 166a 2 0 and R 2 ≈ 22a 2 0

Figure 2 .

 2 Figure 2.14. The diffusivity D(ω) (in units of Ωa 2 0) as a function of ω (in units of Ω) for various µ (0, 0.01, 0.1, 1) for sample with N = 14 3 (crosses). The diffusivity was calculated using formula of Edwards and Thouless (2.29) with c = 1 and averaged over two thousand realizations. The arrows indicate frequencies ω max in the VDOS g(ω) for corresponding values of µ. Open symbols correspond to phonon diffusivity (2.32) below the Ioffe-Regel crossover frequency ω ir .

3 Figure 3 . 1 .

 331 Figure 3.1. (a) The VDOS for different values of the parameter Λ. (b) The VDOS as a function of the reduced wave vector q * = ω/c T . Inset: Λ-dependence of the sound velocities c T (green symbols) and c L (red symbols). Green line shows fit with c T ∝ √ Λ dependence.

Figure 3 . 2 .

 32 Figure 3.2. (a) The VDOS divided by ν 2 , that shows the boson peak for different Λ.Thick lines correspond to constant density configurations as described in Table3.1. Dashed thin lines are relaxed configurations with P ≈ 0 GPa. Solid horizontal thin lines on the left show the low-frequency Debye predictions calculated from the static shear and bulk modules. (b) Boson peak as a function of the reduced wave vector q * = ω/c T . Vertical gray bands mark the position of q * 1 and q * 2 . Arrows show the position of the transverse Ioffe-Regel criteria (seeSection 3.4).

  .

Figure 3 . 3 .

 33 Figure 3.3. (a) The decomposition of the total vibrational density of states (Σ) to longitudinal and transverse components for Λ = 21. Vertical arrows show the transverse and longitudinal Ioffe-Regel frequencies. Inset: the relative number of the longitudinal modes g L (ω)/g(ω) (green line between hatching regions). The relative number of the transverse modes g T (ω)/g(ω) = 1g L (ω)/g(ω) is shown by red hatching between green line and the value 1. (b) The boson peak (Σ) and its longitudinal and transverse components for Λ = 21. Thin horizontal lines show the low-frequency Debye prediction calculated from the static shear and bulk moduli. Dashed lines are estimations from Section 3.4 of the phononic contribution, below the Ioffe-Regel limit. The vertical arrow shows the transverse Ioffe-Regel frequency.

Figure 3 . 4 .

 34 Figure 3.4. Vibration modes corresponding to different frequency range, for Λ = 21.Arrows are proportional to the displacements of the particles (×100). The 2D representation corresponds to a cut along the x-y plane (δz = 5 Å) that contains the particle supporting the largest displacement.

  Figure 3.5. Participation ratio as a function of the reduced wave vector q * = ω/c T for different values of the parameter Λ. Vertical gray bands show three characteristics reduced wave vectors q * 1 , q * 2 , and q * 3 . Inset: zoom on the low-frequency range.

Figure 3 . 6 .

 36 Figure 3.6. (a) The spatial correlation function of the atomic displacements C n (r, ω) as a function of the frequency ν for Λ = 21. The amplitude of the correlation function is indicated by color. The negative correlation is marked by blue color. All amplitudes above 0.1 are shown as 0.1. (b) The same correlation function for different frequencies ν = ω/2π. (c) Position in the first minimum of the correlation function as a function of the reduced wave vector q *for the different values of Λ. The dashed line is r * = 0.449/q * , which corresponds to the first minimum of Eq. (3.10) for transverse modes. Vertical gray bands mark the positions of q * 1 , q * 2 , and q * 3 .

Figure 3 . 7 .

 37 Figure 3.7. The longitudinal (L) and transverse (T) components of eigenmodes in the reciprocal space as a function of the wavenumber q and the frequency ν for the parameter Λ = 21 and Λ = 40.

Figure 3 . 8 .

 38 Figure 3.8. Fits of the dynamical structure factor S L (q, ω) to Eq.(3.14) for Λ = 21 and various values of the wavenumber q (from left to right:1.44, 2.89, 4.33, 5.77, 7.21, 8.66, 10.10 and 11.54 nm -1 ). The inset shows a full curve for q = 1.44 nm -1 .

Figure 3 . 9 . 2 TFigure 3 .

 3923 Figure 3.9. (a) Phonon dispersion curve obtained from the fitting (3.14) for the parameter Λ = 21 (full line) and Λ = 40 (dashed line). (b) Group velocity for the parameter Λ = 21 (full line) and Λ = 40 (dashed line). All curves are shown up to the corresponding Ioffe-Regel frequency. Longitudinal and transverse phonons are denoted by L and T respectively.

Figure 3 .

 3 Figure 3.11. The mean-free path l as a function of the frequency (points), compared to the half-wavelength π/q given by DHO fit (solid lines). The crossing points determine the Ioffe-Regel criterion (shown by arrows). Inset: the Ioffe-Regel frequencies for longitudinal (red symbols) and transverse (green symbols) phonons for different values of the parameter Λ. The solid line shows the trend ν T IR = c T /ξ 2 .

Figure 3 . 12 .

 312 Figure3.12. The spatial distribution of the vibrational energy for different time t after the maximum of the exciting pulse for ν = 8 THz. Solid lines were calculated for the repeatedly extended sample with 4 × 2 × 2 periodic blocks. Dashed lines show the spreading of the energy in one periodic block only. Vertical thin black lines show boundaries of one periodic block. The width of the extended sample is larger than the horizontal plot range.

Figure 3 .

 3 Figure 3.13. Snapshot of rotons with wave packet for Λ = 21, ν = 4 THz and t = 2 ps.

Figure 3 .

 3 Figure 3.14. Spreading of the vibrational energy in space R 2 (t) for different frequencies for Λ = 21. Numbers near curves represent frequencies in THz. Solid lines were calculated for the repeatedly extended sample with 4 × 2 × 2 periodic blocks. Dashed line shows the spreading of the energy in one periodic block only (for ν = 8 THz).

Figure 3 .

 3 Figure 3.15. (a) The diffusivity as a function of frequency for different values of the parameter Λ. Upward and downward arrows show the transverse and the longitudinal Ioffe-Regel criteria respectively. (b)The rescaled diffusivity D/c T as a function of the reduced wave vector q * = ω/c T for the same values of Λ. Vertical gray bands mark the positions of q * 2 and q * 3 .

Figure 3 . 16 .

 316 Figure 3.16. Diffusivity of longitudinal phonons (L) and transverse phonons (T). Dotted line (L+T) shows the estimation of phonon contribution to the total diffusivity by Eq. (3.23). Solid black line is the total measured diffusivity. The vertical arrows show the transverse and longitudinal Ioffe-Regel frequencies.

3 *Figure 3 .

 33 Figure 3.17. Schematic description of the different crossovers in the vibrational properties of harmonic amorphous solids. The VDOS and the diffusivity of vibrations of model amorphous silicon with Λ = 21 are shown in the background.

Figure 4 . 1 .

 41 Figure 4.1. States of the granular media. (a) Free particles, φ < φ c . (b) Critical state, φ = φ c . (c) Solid phase, φ > φ c [O'Hern et al. 2003].

Figure 4 . 2 .

 42 Figure 4.2. Illustration of the interaction between two and three particles.

A

  T has the same eigenvalues as the matrix M = AA T if ∼ A = U AV . The dimensions of the orthogonal matrices U and V are N f × N f and K × K, respectively. For arbitrary (i.e., random) orthogonal matrices U and V ,

∼AFigure 4 . 3 .

 43 Figure 4.3. The VDOS for 1024 spheres with the repulsive potential with γ = 2for various filling densities φ. The solid lines are predictions of the random matrix theory and the dashed lines present the numerical data from[START_REF] Wyart | Effects of compression on the vibrational modes of marginally jammed solids[END_REF]].

A

  is the so-called Wishart ensemble in the random matrix theory. The corresponding VDOS is described by the Marchenko-Pastur distribution[START_REF] Marĉhenko | Distribution of eigenvalues for some sets of random matrices[END_REF] (see Appendix A):

Table 2 .

 2 1. The frequency of maximum in VDOS ω 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.

		5	2.4	2.2 *	0.12
	0.1	0.78	0.74	0.62	0.027
	0.01	0.23	0.23	0.19 0.0066
	0.001 0.072	0.07

Table 3 .

 3 , to quantify here the effect of local order on 1. Transverse and longitudinal sound velocities obtained from the elastic moduli for different values of the parameter Λ with N = 32768.

	Λ	ρ, g/cm 3 P , GPa c T , m/s c L , m/s
	17 19 21	2.303 " "	-1.82 -0.096 0.638	3334 3570 3854	7833 7750 7965
	23.5	"	1.38	4133	8226
	26.25	"	2.1	4386	8484
	40	"	5.07	5305	9490
	17 21	2.339 2.295	-0.011 0.013	3312 3714	8436 8367
	40	2.248	-0.114	5118	9350

Table 3 . 2 .

 32 Transverse and longitudinal Ioffe-Regel criteria for different values of the parameter Λ.

The notion of "phonons" has different meanings in different communities. We use the notation phonons for vibrations, which can be described by a wavevector q and a finite mean free path l (as a result, there is an uncertainty in the wavevector ∆q ∼ 1/l) if ql 1.

to the random matrix model considered in Chapter 1. However, they have different control parameters.

Acknowledgments

I would like to express my sincere gratitude to my advisors Prof. Dyakonov, Prof. Kozub, and Prof. Lorman for the continuous support of my Ph.D. study and related research, for their patience, motivation, and immense knowledge. Their guidance helped me in all the time of research and writing of this thesis. I sincerely thank my co-author Prof. Parshin, who helped me a lot in preparation of publications and this thesis. I am very grateful to my external reviewers Prof. Barrat and Dr. Skipetrov for their insightful comments in this thesis. My sincere thanks also go to Prof. Tanguy, who provided me an opportunity to join their team. I also want to thank the French Government for the financial support. Last but not the least, I would like to thank my family: my parents, my wife Dina and her parents for their love and supporting me spiritually throughout writing this thesis and my life in general.

The random matrix theory has many interesting results. It gives properties of eigenvalues of random matrices of certain symmetry classes (so-called random matrix ensembles). Unfortunately, the derivation of many results is mathematically rather difficult for non-specialists. Therefore, below we present the derivation of these results in a some simplified form. We obtain the density of states in the Wishart ensemble. Then we consider the Brownian motion of eigenvalues and demonstrate the level-repulsion effect.

Density of states

Let us consider the Wishart ensemble

where A is a random N × K matrix with independent and identically distributed elements having the properties

In the general case, the matrix A can be a real, complex, or quaternionic matrix.

In the random matrix literature, these three ensembles are called the Laguerre orthogonal ensemble (LOE), Laguerre unitary ensemble (LUE), and Laguerre symplectic ensemble (LSE) respectively [START_REF] Forrester | Log-Gases and Random Matrices[END_REF]].

Let us consider a perturbation δA

The transformation

93

does not change the variance of matrix elements. Therefore, the matrix M ′ = A ′ A ′ † belongs to the same ensemble as the matrix M . The transformation of the matrix M can be written as .5) where δM = (A + δA)(A + δA) † -M = δA A † + A δA † + δA δA † .

(A.6)

The matrix M = AA † is Hermitian (and positive definite) matrix. Therefore, there is an orthogonal (unitary, symplectic) matrix U such that

is a diagonal matrix. According to the perturbation theory, the change of an eigenvalue ε i of the matrix M (due to the addition of the matrix δM ) can be written in the form

where δ

Let us average Eq. (A.8) over the perturbation δA. The matrix A we keep fixed. Therefore, the matrices U and M , and eigenvalues ε i are also fixed. According to the perturbation theory, the change of an eigenvalue of the matrix M (due to the addition of the matrix δM ) can be written in the form

The matrix A can be presented as a singular value decomposition (A.11) where U and V are orthogonal (unitary, symplectic) matrices, and D is diagonal

where δ

The matrix δ ∼ A has the same mean and variance as the matrix δA

Indeed, the matrices U and V are orthogonal (unitary, symplectic) so

From Eqs. (A.12) and (A.14) we get the properties of the matrix δ

where β is the Dyson's index, β = 1, 2, 4 for orthogonal, unitary and symplectic ensembles respectively. Therefore, the change of an eigenvalue of the matrix M (A.10) is

The matrix M ′ belongs to the same ensemble as the matrix M , therefore, δε i = 0. The matrix size is large N, K ≫ 1, therefore, we can write the sum as an integral of a continuous eigenvalue distribution ρ(ε)

Thus, we get the integral equation

It has a solution

where we introduce the notation (z) + = max(z, 0) and

The corresponding VDOS is

where

Brownian motion and the level repulsion

Let us consider a repeated transformation (every time with a new realization of the perturbation δA)

According to Eq. (A.17)) each transformation slightly changes eigenvalues ε i of the matrix M = AA † . We can ascribe a small time interval δt = 4w 2 to each transformation. In this case, eigenvalues ε i (t) depend on the time t. This time evolution looks like Brownian motion (Fig. A.1). The idea of Brownian motion of eigenvalues was introduced by Dyson for Wigner ensemble [Dyson 1962a]. He also showed that the eigenvalue statistics is the same as the statistics of one-dimensional Coulomb gas [Dyson 1962e;[START_REF] Dyson | Statistical theory of the energy levels of complex systems[END_REF][START_REF] Dyson | Statistical theory of the energy levels of complex systems[END_REF][START_REF] Dyson | Statistical theory of the energy levels of complex systems[END_REF]Dyson 1962a]. Below we will show that eigenvalues "repulse" from each other in the Wishart ensemble.

After each transformation (A.24) the distance between two consecutive eigenvalues S = ε i+1ε i changes by δS = δε i+1δε i . One can show that for 96 small values of S δS = 4w 2 ε S , (A.25)

where angle brackets denote averaging over the perturbation δA. In the continuum limit δt → 0 the distance between two consecutive eigenvalues S(t) performs a Brownian motion with the diffusivity D = ε/β, mobility µ p = ε and temperature T = 1/β (according to the Einstein relation). The Eq. (A.25) means the repulsion force F = 1/S for S → 0.

After a large time t, there is a certain distribution Z(S) of the distance between two consecutive eigenvalues S. This distribution obeys the differential equation .28) This distribution goes to 0 if S → 0. Therefore the probability of close eigenvalues is small (smaller than one for independent eigenvalues, see Section 1.5). It is the so-called level repulsion effect.

Appendix B

The kernel polynomial method

The eigenvalue and eigenvector analysis of large sparse matrices is a common problem in physics. The time of the full diagonalization of N × N matrix is O(N 3 ). In three dimensions the number of atoms is N = L 3 and the computation time grows extremely fast as O(L 9 ). Thus, the full diagonalization is limited by the size L ≈ 30 which is not enough for accurate analysis of long-wave modes.

However, the precise values of all eigenvalues and eigenvectors sometimes are unnecessary information. Many measurable physical quantities (the VDOS, dynamical structure factor) are averaged over many eigenvalues and eigenvectors. The kernel polynomial method (KPM) can calculate such quantities for L 300 with a high accuracy. The KPM was introduced in [START_REF] Silver | Calculation of densities of states and spectral functions by Chebyshev recursion and maximum entropy[END_REF] and detailed reviewed in [START_REF] Weiße | The kernel polynomial method[END_REF]. In this Appendix we present a modification of the KPM for vibrational analysis and compare it with a classical correlation method.

Vibrational density of states

There are several methods for calculation of the VDOS without exact diagonalization of the dynamical matrix. Following standard statistical mechanics methods [START_REF] Dove | Introduction to lattice dynamics[END_REF]], one can make a molecular dynamics simulation of a large system excited at a small temperature T . Then assuming the equipartition of the energy, we can compute VDOS as a Fourier transform of the velocity autocorrelation function [START_REF] Sampoli | Numerical study of the low-frequency atomic dynamics in a Lennard-Jones glass[END_REF]]

99

Chebyshev polynomials of the second kind are defined by recurrence relations

They have an equivalent trigonometric definition

From Eqs. (B. 2) and (B.3) the density of states can be expressed in terms of the sine Fourier transform

where ϕ depends on ω as ϕ = 2 arcsin(ω/ω max ) and µ k is the k-th Chebyshev moment

It is not possible to calculate the infinite number of the Chebyshev moments µ k , so we can cut off the series (B.3) and (B.8) at some K-th degree which is controlled by the desired accuracy of the calculations. For the δ-function it gives the following approximation

The damping factors γ k were introduced to avoid Gibbs oscillations. With increasing k these factors decrease gradually from 1 to 0 (for k = K + 1). One of the best choice for γ k are Jackson damping factors [START_REF] Weiße | The kernel polynomial method[END_REF]]. The finite number of moments leads to the finite-width approximation of the δ-function [START_REF] Weiße | The kernel polynomial method[END_REF]]

The width δε = π √ 1ε 2 /K corresponds to the frequency resolution δω = π ω 2 maxω 2 /2K. The greater the degree K of the polynomial, the closer it to the delta function. As a result, we can approximately calculate the VDOS as

(B.12) This sum can be calculated now by the Fast Fourier Transform (FFT) which is implemented in many mathematical libraries.

We turned the calculation of the VDOS g(ω) into the calculation of Chebyshev moments µ k . Their definition (B.9) can be written as

where we use ket notation |j for the jth eigenvector of the matrix ∼ M (the eigenvectors of the matrices ∼ M and M are the same). For a sufficiently large matrix ∼ M the sum in Eq. (B.13) can be replaced by the averaging over several realizations of a Gaussian random vector |u 0 with unit norm

Indeed, let us expand the random unit vector |u 0 over eigenvectors |j of the matrix

The random vector |u 0 is normalized, so |β j | 2 = 1/N . As a result, we have

Chebyshev moments µ k for k = 0, . . . , K can be easily found by recurrence matrix-vector multiplications like (B.4) -(B. 6)

After averaging these projections over several number of realizations R we obtain Chebyshev moments µ k = m k . Then the resulting VDOS is calculated making use of Eq. (B.12). ). We can see that K = 30 is not enough because both peaks of the resulting VDOS are sufficiently broadened. On the other hand, the value of K = 3000 is unnecessarily big, and we can see peaks from distinct eigenfrequencies (Fig. B.1b). We have found that the optimal value of K is around K = 300. For this value the KPM takes about one minute on a modern computer for calculation the VDOS.

Chebyshev polynomials of the first kind have similar recurrence relations (B.4)-(B.6). So KPM can be implemented with the first kind polynomials as well as with the second kind polynomials. However, since usually g(0) = 0, the second kind polynomials give better approximation in the low-frequency region.

We would like to emphasize here that the KPM has interesting physical meaning. The recurrence relations (B.19) indeed reveals a connection between KPM and CVV method. By definition, with δt = 2/ω max and |u k is the vector at step k in the KPM. The Eq. (B.21) has the same form as the first step Verlet integration of equations for atomic displacements [START_REF] Verlet | Computer "experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules[END_REF]]

where the acceleration üi (t) of the ith atom is defined by the Newton's law (1.2) and u i (t) is the atomic displacement. Therefore, we can consider the integer variable k as a discrete time t = k δt. Usually, the time step δt should be much less than 1/ω max for reasonably small error in the integration procedure. The KPM relaxes this requirement to δt = 2/ω max . The Chebyshev moments µ k = u 0 |u k for k = 0, . . . , K can be considered as auto-correlation functions of atomic displacements µ(t) = u(0)u(t) for 0 ≤ t ≤ t max = K δt. The resulting VDOS (B.12) is the Fourier transform of the Chebyshev moments. The finite frequency resolution δω corresponds to the finite-time limit 1/t max . Thus the KPM is similar to the CVV method (B.1), however, the remarkable properties of the Chebyshev polynomials allows to take a big time step δt = 2/ω max instead of a much smaller time step δt ≪ 1/ω max .

Eigenvector analysis

The correlation function (Sec. 3.3) and the dynamical structure factor (Sec. 3.4) are two of the main eigenvector characteristics. The former shows the spatial correlations of vibrations and the latter shows the structure in the reciprocal space. The direct eigenvector analysis with full diagonalization of the dynamical matrix takes too much computational time. In this Section we show how to properly modify the KPM for eigenvector analysis.

103

The definition of the correlation function (3.6) is similar to the definition of VDOS (3.4), but it contains in addition (as multiplier) the bilinear form of the eigenmodes, which depends on the external parameter r. In this case we can use the whole evaluation procedure of the KPM with modified projection (B.20) of the state |u k to the initial state |u 0

Thus, for each fixed r we can efficiently calculate the correlation function as a sum of Chebyshev polynomials.

The same idea is applicable to the dynamical structure factor. It depends on the spatial Fourier transform of the eigenmodes F η (q, ω), which also contains a bilinear form of the eigenmodes. Therefore, we can use modified projections (B.20), which are slightly different for longitudinal and transverse components

A more general and detailed information about the calculation of eigenvector characteristics and Green's function can be found in the review [START_REF] Weiße | The kernel polynomial method[END_REF]].

Appendix C General decomposition in transverse and longitudinal vibrations

In Chapter 3 we have shown that the polarization of vibrational modes plays a crucial role in amorphous silicon-like materials. However, the classical notion of transverse and longitudinal modes are not applicable for amorphous systems because it relies on the well-defined wavevector q . In this Appendix we generalize the notions of transverse and longitudinal vibrations for disordered systems like amorphous silicon.

Longitudinal and transverse components

Low-frequency vibrations below the Ioffe-Regel criterion are well-defined plane waves (see Section 3.4 for details). In continuous medium approximation, the displacement fields u(r) for longitudinal (L) and transverse (T) waves have a form

However, above the Ioffe-Regel criterion, the wave vector q is ill-defined so we cannot use the definition (C.2) in a general case.

The transverse displacement field u T (r) has zero divergence. Therefore, it conserves the local volume. A natural analog of the local volumes in amorphous media are Voronoi cells constructed around each atom. By definition, the Voronoi cell V i , associated with atom i is the set of all points in the surrounding space whose distance to the atom i is not greater than their distances to the other atoms j [START_REF] Aurenhammer | Voronoi diagrams a survey of a fundamental geometric data structure[END_REF]]. This type of cells is also known as Wigner-Seitz cells in Displacements of atoms u i in amorphous media may (or may not) change volumes of Voronoi cells. We will call the displacement of atoms u i to be transverse if it does not change the volumes of all Voronoi cells. For that let us introduce a matrix A (not be confused with the matrix A from the main text) which is responsible for the relative change of the ith Voronoi cell volume V i under jth atom displacement in the direction α

The explicit formula for the matrix A will be derived in the next Section. Using this matrix the displacement of jth atom in the direction α, u jα results in the following relative change of the Voronoi cell volumes A i,jα u jα . Summing over all j and α gives the relative change of the ith Voronoi cell volume

In the bra-ket notation, this equation reads |ε = A|u where A is a rectangular N ×3N matrix (with N being the number of atoms), and |u is a displacement vector with 3N elements. The matrix A is a discrete analog of the divergence operator.

By definition the transverse component |u T of an arbitrary |u satisfies to equation A|u T = 0, i.e. |u T is the projection of the displacement |u to the null space of the matrix A. The longitudinal component |u L is a remaining orthogonal component of the displacement field, and it is the projection of |u to the row space of the matrix A. These projections have the following forms [Meyer 2000, Eqs. (5.13.3) and (5.13.6)]

where

One can easy check that A|u T = 0 and u L |u T = 0.

Thus, P L |j and P T |j are projections of the eigenmode |j to longitudinal and transverse components respectively. Therefore, the total VDOS g(ω) can be decomposed into the longitudinal and transverse components in general case independently on frequency ω

where eigenfrequency ω j corresponds to the eigenvector |j . In three dimensions we have g L (ω)dω = 1/3 and g T (ω)dω = 2/3.

The definition of longitudinal and transverse components of the VDOS (C.9) contains the bilinear form of the eigenmode (as well as the correlation function (3.6) and the Fourier transform (3.12), (3.13)). Therefore, one can apply the KPM with the modified projection (B.20)

The results of this method are discussed in Section 3.2.

Derivation of the matrix A

Let us show how matrix A can be derived from geometry only. Let we shift only the atom j by the vector u (Fig. C.2). The vector u is small, so the Voronoi cells after shifting have almost the same structure, but their facets are slightly shifted and rotated in space. The shifting u of the atom j can change the volume V i of a nearest neighbor cell i. By definition the cells i and j are nearest neighbors if they have a common facet S ij .

By definition of the Voronoi cell, the facet S ij lies in the plane P ij , which has the normal n ij = r ij /r ij , where r ij = r jr i is the vector connecting atoms i and j and go through the point b ij = (r i + r j )/2. Therefore, the equation for this plane reads After shifting of the atom j, the new facet S ′ ij lies in the plane P ′ ij which satisfies the equation

where

The signed distance from an arbitrary point r to the plane

This distance has the sign "+" if r lies on the same side to the plane P ′ ij as the atom i and the sign "-" if r lies on the opposite side to the plane. If r lies in the plane P ′ ij , then the distance d ′ ij (r) is equal to 0, and we get the Eq. (C.12). In the linear approximation on u, the change of the volume V i is the integral of d ′ ij (r) over the surface of the initial facet S ij

where S ij and c ij are correspondingly the area and the centroid of the facet S ij . In the linear approximation on u, the Eq. (C.14) reads

where we have taken into account that

If all neighboring atoms shift, the change of the volume V i can be written in the matrix form δV i = jα W i,jα u jα .

(C.17)

The nondiagonal elements of the matrix W we have already found above

The diagonal element W i,iα means the change of the ith Voronoi cell volume under shifting of the ith atom itself. Shifting of the all atoms by the same vector does not change the volumes of the Voronoi cells. Therefore

After dividing by the volume, we finally get the matrix

It is notable that finite-elements methods also often use Voronoi cells and have a similar definition for the finite differences for the divergence operator [START_REF] Mishev | Finite volume methods on Voronoi meshes[END_REF]].

Abstract It is well known that various amorphous solids have many universal properties. One of them is the temperature dependence of the thermal conductivity. However, the microscopic mechanism of the heat transfer above 20 K is still poorly understood. Recent numerical simulations of amorphous silicon and silica show that vibrational modes in the corresponding frequency range (called "diffusons") are delocalized, however they are completely different from low-frequency acoustic phonons.

In this work we present a stable random matrix model of an amorphous solid. In this model one can vary the strength of disorder going from a perfect crystal to extremely disordered soft medium without macroscopic rigidity. We show that real amorphous solids are close to the second limiting case, and that diffusons occupy the dominant part of the vibrational spectrum. The crossover frequency between acoustic phonons and diffusons is determined by the Ioffe-Regel criterion. Interestingly, this crossover frequency practically coincides with the boson peak position. We also show that, as a function of frequency, the diffusivity and the vibrational density of states of diffusons are practically constant. As a result, the thermal conductivity is a linear function of temperature up to rather high temperatures and then saturates. This conclusion is in agreement with numerous experimental data.

Further, we consider a numerical model of amorphous silicon-like materials and investigate the role of disorder for longitudinal and transverse vibrations. We also show that the random matrix theory can be successfully applied to estimate the vibrational density of states of granular jammed systems.

Keywords amorphous solids, vibrations, random matrices.

Résumé Il est bien connu que divers solides amorphes ont de nombreuses propriétés universelles. L'une d'entre elles est la variation de la conductivité thermique en fonction de la température. Cependant, le mécanisme microscopique du transfert de chaleur dans le domaine de température supérieure à 20 K est encore mal compris. Simulations numériques récentes du silicium et de la silice amorphes montrent que les modes de vibration dans la gamme de fréquences correspondante (dits diffusons ) sont délocalisés. En même temps ils sont complètement différents des phonons acoustiques de basse fréquence.

Dans ce travail, nous présentons un modèle stable de matrice aléatoire d'un solide amorphe. Dans ce modèle, on peut faire varier le degré de désordre allant du cristal parfait jusqu'au milieu mou extrêmement désordonné sans rigidité macroscopique. Nous montrons que les solides amorphes réels sont proches du deuxième cas limite, et que les diffusons occupent la partie dominante du spectre de vibration. La fréquence de transition entre les phonons acoustiques et diffusons est déterminée par le critère Ioffe-Regel. Fait intéressant, cette fréquence de transition coïncide pratiquement avec la position du pic boson. Nous montrons également que la diffusivité et la densité d'états de vibration de diffusons sont pratiquement constantes en fonction de la fréquence. Par conséquent, la conductivité thermique est une fonction linéaire de la température dans le domaine allant à des températures relativement élevées, puis elle sature. Cette conclusion est en accord avec de nombreuses données expérimentales.

En outre, nous considérons un modèle numérique de matériaux de type de silicium amorphe et étudions le rôle du désordre pour les vibrations longitudinales et transverses. Nous montrons aussi que la théorie des matrices aléatoires peut être appliquée avec succès pour estimer la densité d'états vibrationnels des systèmes granulaires bloqués.

Mots clés solides amorphes, vibrations, matrices aléatoires.