Acoustic waves in glasses

MSc Physics

theoretical physics

Master’s Thesis

Acoustic waves in glasses

by

Thomas W. Hiemstra

6220894

July 2017 60 ECTS 1-8-206 to 1-8-2017 Faculty of Science (FNWI)

Supervisor:

dhr. dr. E.Lerner

Examiner:

dr. C.J.M.Coulais

Contents

1 Introduction 1 2 Numerical methods 4 2.1 Dynamics . . . 5 2.2 Autocorrelation functions . . . 7 2.3 Programming details . . . 9 2.4 Results . . . 10

3 Comparison with experimental data 12 4 Coarse-graining and auto-correlations 13 4.1 Programming details . . . 13 4.2 results . . . 14 5 Nonlinear dynamics 16 5.1 Results . . . 17 6 Density of states 18 6.1 Programming details . . . 18 6.2 Results . . . 19 7 Soft spheres 20 7.1 Results . . . 21 8 Hybridation coefficients 23 8.1 Results . . . 24 9 Conclusion 26 A Coarse graining 28 References 30

1

Introduction

Glasses are known as amorphous solids which mean that, unlike a crystal, they lack long range order. These type of solids are of growing interest to physicists. This is due to the technical applications in for example amorphous semiconductors. But the experimental and theoretical aspect of recovering solids state physics without the periodic lattice is also of great interest

Decreasing the temperature of a liquid causes it undergo a first order phase transition to a crystal. However it is possible to avoid the crystallization and keep the system in a meta stable phase. The system is now in what is known as the supercooled phase these supercooled liquids are still poorly understood. For instance, how deeply can the liquid be supercooled before a metastable limit is hit and how is the crystallization prevented?

The most interesting feature of supercooled liquids is the dynamic glass transi-tion. Below a certain temperature the relaxation time of the system increases by such an amount that a significant dynamical arrest intervenes and it becomes impossible to equilibrate the system within a reasonable experimental time. The nature of the glass transition and the physical processes governing it are still an open question[1] [2].

Figure 1: A very slow cooling rate leads to a discontinuous change in volume to a crystal state (purple curve). A rapid quench leads to a continuous change in volume to a glass (blue curve).

described by relatively simple physics. Instead the microscopic arrangement of the system has to be taken into account which presents many different prob-lems on several length scales. For instance the disorder of the system leads to so called soft spots. Areas of the solid which can move more or less freely requiring little energy to move.

Many of the interesting properties of Glasses occur at low energies where the be-haviour of glasses differs strongly from that of crystals. For example the density of states of a glass shows some remarkable features at low frequencies. In recent years there has been major progress in understanding the low energy properties of glasses. This is due to improved scattering experiments [3][4][5], numerical simulations [6][7] and effective medium theories[8][9].

This work focuses on the results and predictions of a recent paper called: ”Anomalous phonon scattering and elastic correlations in amorphous solids”. [10] Written by S. Gelin, H. Tanaka and A. Lemaˆıtre. In this paper the acoustic properties of glassy solids are studied.

In crystals the vibrational excitations are generally well understood. If the crystal is perfect the vibrational motion of the atoms can be decomposed into independent normal modes. Normal modes are plane-wave excitations, called phonons, and are characterized by a frequency, wave-and-polarization vector. In a crystal anharmonic effects cause coupling between the modes which in turn leads to a finite life time of the phonons. When a crystal is not perfect the defects in the material also give rise to scattering also limiting the lifetime of phonons. In other words: the phonons are damped and will not continue forever in the material.

In a glass the situation is much more complex however. Since in a glass it is no longer possible to make the assumption of periodicity as in a crystalline solid. The normal modes are no longer plane waves but are much more complex. A consequence is that phonons are damped in a glass just like in an anharmonic crystal. This effect becomes more and more important as the wavelength of the phonons reaches the inter atomic spacing. A remarkable feature of glasses is that numerically a harmonic approximation still leads to scattering in a glass[11], this would not happen in a crystal. Even though topological disorder exists in glasses, phonons are still present in the system at very high frequencies (in the THz range). The physics underlying acoustic waves in a glass are still poorly understood. An attempt is made in this work to uncover some of the mechanism governing this phenomenon.

Classically the scattering of phonons in a solid is attributed to scattering due to random impurities in the solid. This in turn then leads to Rayleigh scattering which means that the sound attenuation (or damping) should scale as kd+1_[12]. In the paper a deviation from the Rayleigh scattering is found for acoustic waves in a glass using a numerical approach. They do this by exciting a standing wave in a numerical model of a glass and looking at the damping of the wave. They then show this by plotting the sound attenuation vs kd+1_:

Figure 2: Log-lin plots of Tp(k)/kd+1 versus k in 2D for a, longitudinal and b, transverse waves.[10]

They then argue that all the data is consistent with Tp(k)/kd+1 increasing linearly with ln k for decreasing k. Thus instead the relation Tp(k) ∼ kd+1ln k is proposed. However it is also suggested that there might be a crossover between Rayleigh scattering at low frequencies and a different kind of scattering at high frequencies:

[13].

Figure 3: acoustic excitations of a 3 dimensional system showing 2 different regimes.

The question is then what causes this deviation from the predicted Rayleigh scattering? In the paper they argue that in the derivation of the kd+1 be-haviour the long range correlations in the material are ignored. They present

evidence that this is incorrect and that the correlations are the reason for the log correction.

We have found however that altering the model by removing the internal forces from the system causes the log correction to disappear. The correlations in the material are also sown to remain the same after removing the forces thus proving they are unrelated to this correction term. A look is also taken at a different kind of system also displaying this log correction without having any kind of long range correlations further proving the correlations are unrelated to the correction term. The density of states is also examined for a two dimen-sional system. Finally a new and promising approach is proposed to find the theoretical explanation for the log correction. Namely the behaviour of the so called ”hybridation coefficients”[10] in the paper.

2

Numerical methods

The system under consideration is described by a potential energy U ( ~R). Here the location of the kth particle is denoted by ~Rk and the vector of all particles as ~R and the equilibrium position of all the particles as ~R0. Before continuing the definition of a solid is needed. For the current work a solid is simply a state that corresponds to a local minimum of the potential energy:

∂U ∂ ~Rk    _{R= ~}~ _R(0) = 0 (1)

known as mechanical, equilibrium. This means that the net force per particle is zero. Next a computationally convention form of the potential energy is needed. As stated in the introduction the harmonic approximation of the potential is used which reads:

U ( ~R) ' U    _{R= ~}~ _R(0) +X k ∂U ∂ ~Rk    _{R= ~}~ _R(0) · ( ~Rk− ~R (0) k ) (2) +X k,` 1 2( ~Rk− ~R (0) k ) · ∂2U ∂ ~Rk∂ ~R`    _{R= ~}~ _R(0) · ( ~R`− ~R (0) ` ) . (3)

Denoting the displacement from the equilibrium position as ~∆l= ~Rl− ~R0l and the energy at the mechanical equilibrium state as U0= U |_{R= ~}~ _R(0). Finally the matrix of second derivatives which is known as the hessian or dynamical matrix is written as: ↔ Mk`≡ ∂2<sub>U</sub> ∂ ~Rk∂ ~R`    _{R= ~~ R}(0) . (4)

Using equation 1 and the definitions given above the potential energy becomes: U ( ~∆) ' U0+1₂∆~k·

↔

where repeated indices are summed over. Since the energy of the system is always defined up to a constant the first term on the right hand side can be omitted:

U ( ~∆) ' 1₂∆~k· ↔

Mkl·~∆l, (6)

A quick note on the condition of solidity: There are two conditions that must be satisfied before a system described by the harmonic approximation can be considered a stable solid. First there is the mechanical equilibrium given by equation 1.

Before the second condition is discussed first the eigenvalues of the dynamical matrix need to be discussed. The dynamical matrix is real and symmetric thus diagonalizable. The eigenvalues of the dynamical matrix are denoted as λmand the eigenvectors as ~Ψ(m) _{which then leads to:}

~

Ψ(m)·M ·~↔ Ψ(m)= λm. (7)

where the eigenvectors are orthonormal: ~

Ψ(m)· ~Ψ(n)= δnm. (8)

The second condition then states that all the eigenvalues must satisfy λm> 0. Or for any displacement the energy is positive:

~

∆·M ·~↔ ∆ > 0 (9)

The dynamics in the harmonic approximation are only valid below 300K [11] meaning we can only look at the solid phase under this approximation.

2.1

Dynamics

The force in a system is simply minus the gradient of the potential: ~

Fk= − ∂U ∂ ~Rk

(10) The displacement about the mechanical equilibrium of the system is defined as ~

∆ = ~Rk− ~R (0)

k which means that we can freely change the derivative according to:

∂ ∂ ~R =

∂

∂ ~∆ (11)

Which then leads to:

~ Fk = −

∂U ∂ ~∆k

And so in the harmonic approximation the force can be written as: ~ Fk≈ − X ` ↔ Mk`·~∆` (13)

From this it is clear that an expression is needed for the dynamical matrix act-ing on a displacement.

So far everything has been completely general. In order to work out an expres-sion for the force a specific potential has to be chosen. The model used here assumes that the energy is a sum of pairwise contributions:

U =X

α

ϕα (14)

Just for completeness the specific form of the the potential is:

ϕij =    0 _ σij rij n +Pq l=0c2l _r ij σij 2l _r ij σij < xc 0 rij σij > xc (15)

rij being the distance between the ithand the jthparticle,  is and energy scale and xc is the dimensionless distance for which the potential vanishes continu-ously up to q derivatives. Here xc= 1.48, n = 10 and q = 3. c2lis given by:

c2cl= (−1)k+1 (6 − 2k)!!(2k)!! (β + 6)!! (β − 2)!!(β + 2k)· r −(β+2k) c (16)

So ϕαis the pair potential which only depends on the distance from the ithand jth particles. And α denotes the specific pair of interacting particles. So the sum over alpha can be written asP

α≡

P

i<j i.e. the sum over all interacting particles. The force on the kth_{particle is now:}

~ Fk = − ∂U ∂ ~Rk = −X α ∂ϕα ∂ ~Rk = −X α ∂ϕα ∂rα ∂rij ∂ ~Rk =X i<j fij ∂rij ∂ ~Rk , (17) Where fij≡ − ∂φij

∂rij are the pairwise forces. The dynamical matrix can now be written as: ↔ M`k= ∂2U ∂ ~R`∂ ~Rk = −X i<j ∂fij ∂ ~R` ∂rij ∂ ~Rk −X i<j fij ∂2rij ∂ ~R`∂ ~Rk = −X i<j ∂fij ∂rij ∂rij ∂ ~R` ∂rij ∂ ~Rk −X i<j fij ∂2<sub>r</sub> ij ∂ ~R`∂ ~Rk = X i<j κij ∂rij ∂ ~R` ∂rij ∂ ~Rk −X i<j fij ∂2<sub>r</sub> ij ∂ ~R`∂ ~Rk , (18)

Where κij= − ∂fij ∂rij= ∂2_φ ij ∂r2 ij

is called the stiffness. What we need now is: ∂rij ∂ ~Rk = (δjk− δik) ~ Rij rij = (δjk− δik)ˆnij, (19)

Where the unit vector is defined as ˆnij ≡ ~ Rij

rij, the differences ~Rij ≡ ~Rj− ~Ri and the distances rij =

q ~

Rij· ~Rij. What we need is the contraction with the displacement vector ~∆ so we need:

X

k ∂rij ∂ ~Rk

· ~∆k = ~∆ij· ˆnij. (20)

And similar to the previous exercise we obtain: ∂2_r ij ∂ ~R`R~k = (δjk− δik)(δj`− δi`) "↔ I rij −nˆijnˆij rij # , (21)

Which then leads to the following relation: X `,k ~ ∆`· ∂2rij ∂ ~R`R~k · ~∆k= ~ ∆ij· ~∆ij rij −( ~∆ij· ˆnij) 2 rij . (22)

Using all that we finally arrive at an expression necessary to evaluate the dy-namics of the system:

X k ↔ M`k ·~∆k = X i<j  κij( ~∆ij· ˆnij)ˆnij− fij rij h ~_{∆ij− (~∆ij· ˆnij)ˆnij}i = X i<j  κij+ fij rij  ( ~∆ij· ˆnij)ˆnij− fij rij ~ ∆ij  . (23)

In the code the 2 terms in the above equation are calculated for every pair of interacting particles: −fij rij and  κij+ fij rij 

. Removing the forces means the first term is removed and the second term simply reduces to κij. This effectively turns the system into a network of masses connected by springs.

2.2

Autocorrelation functions

As stated in the introduction we are looking at the damping of standing waves in a solid. In order to study the dynamics of the system the discrete wave equation has to be solved:

¨ ~ ui(t) =

↔

The system is at rest for t < 0 i.e. ~u = ˙~u = 0 here u is the displacement vector which in the code is called z. Here k is the probed wave vector, k is then re-stricted according to k =2πn_l where l is the length of the system. This problem is the same as the elastic response without a source term with initial conditions ~ui= 0 and ˙~ui= a sin(~k · ~ri).

Now all that is left to do is solve equation 24. This can be done by either a full diagonalization of the dynamical matrix or by leapfrog integration. Diag-onalizing the dynamical matrix becomes almost impossible as the system size increases since this matrix is of size N × N where N is the amount of particles. The alternative then is leapfrog integration. Simply put leapfrog integration is equivalent to updating the position and velocity of all the particles at interleaved time points. In terms of formulas this would be:

xi= xi+ vi−0.5∆t (25)

ai= F (xi) (26)

vi+0.5= vi−0.5+ ai∆t (27)

Where the acceleration due to the force is done by equation 23. Of course in order to calculate x1 we need v0.5. An educated guess for v0.5 would be v0.5 = ai∗ 1/2∆t which should be accurate enough if the time step is small enough.

After a sufficient amount of the time initial wave has completely lost it’s struc-ture due to scattering in the system. Visually this looks like this:

Figure 4: the velocity of all the particles. On the left the initial velocity dis-placement, on the right the system after a certain amount of elapsed time The observable of the system is then the velocity auto-correlation function:

C(t) = ˙ ~ u(t). ˙u0 p ˙~u0_{. ˙~u}0 (28) i.e. the overlap of the time evolved velocity at time t with the initial velocity. Since the system used in the model has N particles the result will never be

perfect since N should be infinity. So in order to get result which are robust we need N as large as possible. This in turn will increase the time it takes for the calculation to complete. One way around this is to take N reasonably small and average the result of many different systems of the same size. The average of the velocity auto-correlation function of 100 systems is then fitted to exp(Γ(k)t/2) cos(Ω(k)t) for a range of k vectors. The lowest k is determined by the system size and the largest k is chosen to be 1. The result for a specific N and k then looks like this:

0 0.5 1 1.5 2 2.5 3 t -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 C(t)

Figure 5: the velocity autocorrelation fitted to a damped harmonic oscillation Where the red line is the fitted damped harmonic oscillation. The damping coefficient Γ(k) can then be extracted from the fit.

2.3

Programming details

In order to perform the leapfrog integration one has to specify the sys-tem size, whether the initial displacement is either transverse or longitudinal, whether or not to include the internal forces, the minimum and maximum mode (the n in k = 2πn_l ) and the amount of solids to include.

In the case of N = 409600 and larger we only have 10 different solids rather than 100 to work with. In order to still average over 100 different solids the initial displacement is simply displaced: ˙~ui= a sin(~k · ~ri+ φ). where φ is the displacement. The initial wave is displaced 10 times where φ is kept smaller than half the wavelength. This provides us with 10 different results to average over per solid.

2.4

Results

First the results are plotted for the case with the internal forces for both the transverse and longitudinal wave. Also plotted are Γ(k)/k3 _{to show the} deviation from the predicted behavior. These result should be the same as in [10] which they appear to be:

10-2 10-1 100

k

10-4 10-3 10-2 10-1 100 101 102

Γ

T

(k

)

k

3 10-2 10-1 100

k

0 20 40 60 80 100 120

Γ

T

/k

.

3 L=690 L=345 L=172 10-2 10-1 100

k

10-4 10-3 10-2 10-1 100 101 102

Γ

L

(k

)

k

3 10-2 10-1 100

k

0 20 40 60 80 100 120

Γ

L

/k

.

3 L=690 L=345 L=172

a)

b)

c)

d)

Figure 6: damping coefficient Γ(k) with the internal forces included a. transverse wave. b. longitudinal wave. Γ/k3 _{for several solid sizes to show that the data} does not collapse to k3 _{c. transverse waves d. longitudinal waves}

Next the result for the case without internal forces is given which clearly shows that it follows the classical prediction of k3_:

10-2 10-1 100

k

10-5 10-4 10-3 10-2 10-1 100

Γ

T

(k

)

1

3

k

3 10-2 10-1 100

k

0 5 10 15 20

Γ

T

/k

.

3 L=690 L=345 L=172 10-2 10-1 100

k

10-4 10-3 10-2 10-1 100 101

Γ

L

(k

)

1

3

k

3 10-2 10-1 100

k

0 5 10 15 20 25

Γ

L

/k

.

3 L=690 L=345 L=172

a)

b)

c)

d)

Figure 7: damping coefficient Γ(k) without the internal forces a. transverse wave. b. longitudinal wave. Γ/k3 for several solid sizes to show that the data does collapse to k3c. transverse waves d. longitudinal waves

This clearly shows that removing the internal forces causes the system to behave according to Rayleigh scattering thus proving that the forces are responsible for the log correction.

3

Comparison with experimental data

Here a short overview follows of the comparison between the damping coef-ficients obtained in 3D systems for both experimental and numerical data. The experimental data covers: sorbitol at 80 K [14], glycerol at 150 K [15], densified vitreous silica at 570 K [16], sodium silicate at 100 K [4], vitreous silica at 1620 K[17], 300 K [3] and 1 K[18].

The data of the first 4 systems are plotted as Γl/k4 in figure a. Only in the last system was the temperature varied to asses the conditions of emergence of anharmonic effects. And consequently showing their absence at 100 K for this system. This is plotted in panel b.Finally the numerical data is plotted in the last 2 panels clearly showing that the model fits the real world data.

Figure 8: Damping of sound measured in various glasses. a and b represent experimental data where for a: glycerol(circles)[15], sorbitol(crosses)[14], densi-fied vitreous silica (squares)[16] and sodium silicate (triangles)[4]. b: virtreous silica at 1620 K (circles)[17], 300 K (triangles)[3] and 1 K (crosses)[18]. c and d show numerical with circles representing longitudinal initial waves and triangles transverse ones.

4

Coarse-graining and auto-correlations

The authors of the paper [10] claim that the log correction is due to the cor-relations of elastic stiffnesses in the material. If the corcor-relations in the material would change when removing the internal forces then this explanation could still hold. So to investigate this we repeat the procedure from the paper with the forces both intact (to see if our results match up) and without the forces. If the results stay unchanged when removing the forces then their explanation should be incorrect.

In order to look at at elasticity at the local scale a coarse-graining approach is used. It can be shown that the problem of elasticity at the local scale can be reduced to a fluctuating elastic problem based on effective elastic constants SαβκΞ, which each pair of interacting atoms (i,j) contributes:

S_ijαβκχ= C_ijαβκχ+ δακσ_ijβχ (29) Greek subscripts refer to Cartesian coordinates, σβχ_ij is the pair contributing to stress and C_ijαβκχ are the microscopic elastic constants:

C_ijαβκχ= hijnαijn β ijn κ ijn χ ij (30) where hij = ∂2φ(rij) ∂r2 ij r2 ij − ∂φ(rij)

∂rij rij with φij the pairwise interaction. Writ-ing ~nij= (cos(θij), sin(θij)) the elastic constants then take the form Cijαβκχ= hijcosn(θij) sin4−n(θij). They specifically look at the coarse grained field of C3

ij= hijsin(2θij) and it’s spatial auto-correlation field.

Their claim relies on the fact that the spatial auto-correlation of C3 shows a clear cos(4θ) symmetry. They also claim that it decays as 1/r2_{in space (they do} this by taking a cut of the correlation function along π/4 see [10]). They then go on the state that in the work of John and Stephen [19] a 1/rd_{correlation of} disorder leads to Γ(k) ≡ −kd+1_{ln(k). They note however that their calculation} was based on a treatment of mass rather than elastic disorder. However they claim that in scattering theories both types of disorder are essentially equivalent [20].

4.1

Programming details

A derivation of the coarse-grained Cαβκχ _{is given in appendix A. A list is} generated of all the pairwise interactions hij as well as a version without forces (hij=

∂2_φ(r ij) ∂r2

ij

r_ij2). Now equation 39 can be evaluated.

With the coarse-grained map of C3 in hand we can now proceed wit the second map provided in the paper; the spatial autocorrelation map. In order to get this correlation map we start with a point i0, j0 on the coarse-grained map which I’ll call map[i 0][j 0]. Now image the correlation map with it’s centre at i0, j0 like so:

A point on the correlation map (i,j) is now:

map[i 0][j 0]· map[i 0+∆i ][j 0+∆j]. Of course the point map[i 0+∆i ][j 0+∆j] may be outside of the course-grained map but that is fixed using periodic bound-ary conditions. Now we repeat this process for every possible point on the coarse-grained map and sum the results. This will yield the spatial autocorre-lation map.

4.2

results

The results for the correlation map are not expected to change since removing the forces reduces hij=

∂2_φ(r ij) ∂r2 ij r2 ij− ∂φ(rij) ∂rij rij to hij = ∂2_φ(r ij) ∂r2 ij r2 ij. Since the leading order in the potential is of the form r−10 the second derivative is an order of magnitude bigger than the first derivative. So removing the forces, which results in removing the first derivative of the potential here, should not make a noticeable difference.

a)

b)

c)

d)

Figure 9: a. Local coarse grained field C3 _{with forces. b. Local coarse grained} field C3 _{without forces. c.Spatial auto-correlation of the field C}3 _{with forces.} d. Spatial auto-correlation of the field C3 _{without forces.}

It is clear to see that there are no differences between the maps with and without forces, as expected. Thus removing the internal forces removes the log correction but the correlations are left in tact. This is a first proof that the correlations and the log corrections are in fact unrelated.

5

Nonlinear dynamics

An attempt has been made at a numerical approach beyond the harmonic approximation, expanding the Potential further up to a 4th order term. The idea is that the harmonic approximation should be sufficient since it agrees with experimental data [10][11]. However investigating the higher order terms might lead to interesting results. To give an idea of the contraction of the higher order term with the displacement the 3rd order is given here:

∂3U ∂ ~Ri∂ ~Rj∂ ~Rk : ~VjU~k= X α Γαi _ϕ000 α r3 α − 3ϕ00_α r4 α + 3ϕ0_α r5 α  ( ~Rα· ~Vα)( ~Rα· ~Uα) ~Rα +X α Γαi _ϕ00 α r2 α − ϕ0_α r3 α  h ( ~Rα· ~Vα) ~Uα+ ( ~Rα· ~Uα)~Vα+ (~Vα· ~Uα) ~Rα i , (31)

Where Γαk ≡ δjk− δik. Where we contracted with the displacement twice since in equation 13 we have 1 free index on the left side.

It should be noted here that the initial displacement does not have to be normal-ized anymore. Instead the initial velocity is first normalnormal-ized and then multiplied by s√N .

And this is where we run into the first problem. Picking s > 1 causes the dy-namics to be unstable and produces nonsensical results.

Having the dynamics blow up at s > 1 even for a much smaller time step than normal means this part is unfinished and might need some work. Preliminary results however indicate that the higher order terms in the expansion do not contribute significantly.

5.1

Results

For s < 1 however we get a complete agreement with the harmonic approx-imation, the only difference can be seen at the finite sized effects end of the spectrum. But since we are not interested in that part of the numerics (or we fix it by simply increasing N) the results are rather boring. (here again Γ(k) is plotted)

6

Density of states

As mentioned before we are interested in the low energy properties of glasses. Of course when talking about vibrations in a system one has to mention the density of states (DoS). The DoS describes the probability distribution of vibra-tional frequencies in a solid. The DoS of glasses has attracted a lot of attention during the last few decades[21][22]. The density of states in a solid is described by the Debye model. This can be derived for a periodic system but it also holds in the continuum case. A glass however is disordered and displays different behaviour. For instance in 3 dimensions one can plot the DoS of a glass vs ω2 and an excess of vibrational modes (called the boson peak) will be found at low frequencies:

The origin of this strange phenomenon is still an open question.

To study the behaviour of the DoS of a glass the low-frequency excitations of the system have to be exposed. In simulations we can achieve this by decreasing the system size L. This causes the vibrational frequencies of waves to increase but it does not effect the frequencies of the vibrational modes that stem from the disordered structure.

6.1

Programming details

In order to get proper numerics we will need a lot of systems of a small size. From these system we will then need the eigenvalues of the dynamical matrix (eigenfrequencies). Finally we’ll need to get a plot of the probability distribution

of these eigenfrequencies. For the DoS we’ve used 1.105 systems of N = 256, 5.104 systems of N = 484, 2.104 systems of N = 1024 and 2.101 systems of size N = 2116.

6.2

Results

After the binning is complete the results for the case with forces are plotted on a loglog scale. We clearly see the bump moves backwards as the system size increases so in an even bigger system this information would not be visible:

0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 ω 100 101 102 D( ω ) 256 484 1024 2116

Figure 11: The DoS of a 2D system with the internal forces The results here should overlap but they are not (yet) normalized.

In the case without the internal forces a loglog scale is no longer necessary and a normal binning will suffice. Here it is clear to see that doubling the system size causes the second phonon of the smaller system to overlap with the first phonon of the larger system:

0 1 2 3 4 5 6 ω 0 1 2 3 4 5 6 D( ω ) ×104 256 484 1024 2116

Figure 12: The DoS of a 2D system without the internal forces

7

Soft spheres

Another system that can be used to illustrate the point that the log correc-tion and correlacorrec-tions are unrelated is a system composed of soft spheres. This is the simplest model of amorphous solids [23] which can be seen as randomly jammed particles at zero temperature with the following pairwise potential:

ϕij = ( 1 − rij σij α /α rij < σij 0 rij > σij (32)

Where  is the characteristic energy scale of the interaction, rij the separation between the centers of the particles i and j and σij is the sum of the radii of the particles i and j. Here α = 2, repulsive harmonic springs. Particles only interact in this system if they overlap.

7.1

Results

First we look at the damping coefficients of the system with the forces in-cluded since the authors of the paper [10] claim that this system should have a log correction as well. They also claim that this system should display cor-relations giving further evidence for their claim. The log correction is indeed present in this system:

10-2 10-1 100

k

10-4 10-3 10-2 10-1 100

Γ

T

(

k)

k

3

Figure 13: sound damping in a system composed of soft spheres Which shows the same behaviour as the system with the pairwise potential did.

When looking at the correlation function (with the forces included) what we get is something that, at first glance, does not show any correlations at all. Not only that but the auto-correlation function does not show the cos(4θ) symmetry which is crucial for their claim:

a)

b)

Figure 14: a. The local coarse grained field C3 for the soft spheres b. the auto-correlation function for the field C3

We can also look at a comparison between the cut along the correlation function along π/4 which in the glassy system shows a clear 1/r2 _{relation (in 2D) and} (according to their claim) should be present in the soft spheres system as well: The glassy system does display the 1/r2_{behavior while the soft spheres system} just shows statistical noise thus further proving the spheres system lacks any kind of long range correlations.

100 101 102 r 10-4 10-2 100 102 104 106 <C 3(Q)C 3(r)> spheres forces no forces

Figure 15: A comparison between the cut along the correlation function for the glassy system and the soft spheres system

8

Hybridation coefficients

So now we have a pretty solids understanding of what is not the reason for the log correction. A first step towards a possible explanation is detailed here. We start by looking at the formal solutions in terms of eigen modes and eigen values of the system. Any vibration of the solid can be expanded in the eigenen-frequencies of the system. Doing this for the velocity of the displacement we get:

˙ ~ ui = X m ˙ ~ u0ψmcos(ωmt)ψm (33)

Where the sum runs over the relevant modes (i.e. after eliminating the trans-lations associated with periodicity). From this we can write the velocity auto-correlation function as:

C(t) = ˙ ~ u(t). ˙~u0 p ˙~u0_{. ˙~u}0 =X m ξ2cos(ωmt). (34)

with the coefficients:

ξ = ψm. ˙~u 0 p ˙~u0_{. ˙~u}0

Which are called the ”hybridation coefficients” [10].

So the imposed wave is decomposed on the eigenmodes ψmof the system. The distribution of the hybridation coefficients follows a normal distribution around a natural oscillation. Now that the wave is composed of many waves each with a different frequency. This leads to so called phase de-coherence a phenomenon known as oscillator phase noise [24]. This in turn leads to a damping of the wave.

If we regard equation 34 as a continuous problem we can recognize a Fourier transform: C(t) =X m ξ2cos(ωmt) ≈ Z G(ω) cos(ωmt)dω (36)

and since we already know that the velocity auto-correlation should be a damped harmonic what we have here is simply the Fourier transform of a damped har-monic oscillation: exp(Γ(k)t/2) cos(Ω(k)t) = Z G(ω) cos(ωmt)dω (37) Where [24]: G(ω) = Γ Γ2_{+ (ω − Ω)}2 (38)

So the behaviour of the distribution of the hybridation coefficients can tell us something about the camping coefficient.

8.1

Results

The results are, as always, plotted for the case with and without forces choosing an appropriate bin size. Firs the result from the paper is reproduced: Next the results without the internal forces are plotted:

0 1 2 3 4 5 6 ω 0 0.5 1 1.5 2 2.5 ξ 2 ×10-3

Figure 16: distribution of ξ2_{plotted against the eigenfrequencies of the system,} forces included 0 1 2 3 4 5 6 7 8 ω 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 ξ 2

Figure 17: distribution of ξ2plotted against the eigenfrequencies of the system, forces not included

These results suggest that in the case without forces the coefficients are far more sharply peaked around a certain frequency. Indeed plotting the results on the same graph reveals this difference:

This shows that the behavior of the hybridation coefficients differs for both cases. We therefore think it is this quantity that should be looked at when trying to construct a theoretical framework for the log correction.

0 1 2 3 4 5 6 7 8 ω 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 ξ 2 Forces No forces

Figure 18: distribution of ξ2 plotted against the eigenfrequencies of the system

9

Conclusion

The model used in this work perfectly reproduces the results found in [10] meaning that the log correction is present in glassy systems. Removing the in-ternal forces from the system however also removes the anomalous behaviour of the damping coefficients. This leads us to believe the internal forces are respon-sible for the log correction rather than the long range correlations in the system as proposed by the paper[10]. In order to falsify the claim made by the authors of the paper it is shown that the correlations remain unchanged when removing the internal forces from the system. Thus the correlations in the material are unrelated to the log correction found in sound dampening.

To further emphasize this point a system of soft spheres is used which also shows a log correction on the sound dampening but has no correlations at all. So a system with correlations can be manipulated to lose the log correction (remove the forces). And a system without correlation still has the log correction in it. Further proving that the correlations in the system and the anomalous be-haviour of the damping coefficients are unrelated.

Finally a theoretical direction is suggested by looking at the so called hybrida-tion coefficients. It is shown that the distribuhybrida-tion of these coefficients changes when removing the forces from the system. It is therefore suggested that the

hybridation coefficients might play a crucial role in understanding the origin of the anomaly in the damping coefficients.

Further research can be focused on improving the results from the nonlinear approach since that still does not work properly. And the theoretical base for this work is still missing for which so far only a direction has been suggested.

A

Coarse graining

We define the coarse-grained Cαβκχas:

Cαβκχ(~r) = X i,j>i C_ijαβκχ Z 1 0 ω |~r − ~ri− s~rij| Λ  ds (39)

where ~rij ≡ ~rj− ~ri and ω(x) is some coarse-graining function that peaks at x = 0 and decays to zero as x → ±∞ We choose:

ω(x) = (

A(1 − x2₎4 _{, |x| < 1}

0 , |x| > 1 (40)

Where the coefficient A is determined by requiring that Z

dV dsω |~r − ~ri− s~rij| Λ

 = 1

This specific choice of (2) enables one to calculate the integral in (1) exactly. We begin with |~r − ~ri− s~rij| 2 = r2_ij " s2− 2s(~r − ~ri).~rij r2 ij +|~r − ~ri| r2 ij # = r2ij  s2− 2s(~r − ~ri).~rij r2 ij + (~r − ~ri).~rij r2 ij !2 − (~r − ~ri).~rij r2 ij !2 +|~r − ~ri| r2 ij   = r2_ij   s − (~r − ~ri).~rij r2 ij !2 +|~r − ~ri| 2 ~ rij− ((~r − ~ri).~rij)2 r4 ij   (41) using the identity |~a|2_|~b|2_{− (~a ∗ ~b)}2_{= |~a × ~b|}2_{, this becomes}

|~r − ~ri− s~rij|2= rij2   s − (~r − ~ri).~rij r2 ij !2 +|(~r − ~ri) × ~rij| 2 r4 ij   (42) Notice that: |(~r − ~rj) × ~rij| = |(~r − ~ri) × (~rj− ~ri)| = |~r × ~rj− ~r × ~ri− ~ri× ~rj+ ~ri× ~ri| = |~r × ~rj− ~r × ~ri− ~ri× ~rj|

and

|(~r − ~rj) × ~rji| = |(~r − ~rj) × (~ri− ~rj)|

= |~r × ~ri− ~r × ~rj− ~rj× ~ri+ ~rj× ~rj| = |~r × ~rj− ~r × ~ri− ~ri× ~rj| = |(~r − ~ri) × ~rij| So the last term on the RHS of (4) is

|(~r − ~ri) × ~rij|2 r4 ij =|(~r − ~ri) × ~rij||(~r − ~rj× ~rij| r4 ij (43) which is symmetric with respect to interchanging i and j. For simplicity we denote Cij(~r) = (~r − ~ri).~rij r2 ij and Eij(~r) = Eji(~r) = |(~r − ~ri) × ~rij| r2 ij (44) such that |~r − ~ri− s~rij|2= rij2 (s − Cij)2+ EijEji  (45) To calculate the integral in (1), we must first find expression for the lower and upper limits s±; these are obtained by solving

|~r − ~ri− s~rij|2= Λ2 (46)

the solutions are

s±= Cij± s Λ2 r2 ij − EijEji (47)

we now calculate the integral in (1) Z 1 0 ω |~r − ~ri− s~rij| Λ  ds = Z s+ s− ds  1 − |~r − ~ri− s~rij| 2 Λ2 4 = Z s+ s− ds 1 −r 2 ij(s − Cij)2+ EijEji Λ2 !4 =r 8 ij Λ8 Z s+ s− ds Λ 2 r2 ij − EijEji− (s − Cij)2 !4 (48)

Notice that s± must only be considered in the range 0 ≤ s± ≤ 1. We define Gij ≡ Λ 2 r2 ij − EijEji, and solve Z 1 0 ω |~r − ~ri− s~rij| Λ  ds = r 8 ij Λ8  (C_ij2 − Gij)4s − 4Cij(Cij2 − Gij)3s2− 4 3(C 2 ij− Gij)2(Gij− 7Cij2)s 3 −2Cij(7Cij4 − 10Cij2Gij+ 3G2ij)s4+ 2 5(35C 4 ij− 30Cij2Gij+ 3G2ij)s5 +4 3Cij(−7C 2 ij+ 3Gij)s6− 4 7(−7C 2 ij+ Gij)s7− Cijs8+ 1 9s 9  s+ s− . (49)

References

[1] Gary L Hunter and Eric R Weeks. The physics of the colloidal glass tran-sition. Reports on Progress in Physics, 75(6):066501, 2012.

[2] Andrea Cavagna. Supercooled liquids for pedestrians. Physics Reports, 476(4):51–124, 2009.

[3] Giacomo Baldi, Valentina M Giordano, and Giulio Monaco. Elastic anoma-lies at terahertz frequencies and excess density of vibrational states in silica glass. Physical Review B, 83(17):174203, 2011.

[4] G Baldi, VM Giordano, B Ruta, R Dal Maschio, A Fontana, and G Monaco. Anharmonic damping of terahertz acoustic waves in a network glass and its effect on the density of vibrational states. Physical review letters, 112(12):125502, 2014.

[5] Giulio Monaco and Valentina M Giordano. Breakdown of the debye ap-proximation for the acoustic modes with nanometric wavelengths in glasses. Proceedings of the national Academy of Sciences, 106(10):3659–3663, 2009. [6] J¨urgen Horbach, Walter Kob, and Kurt Binder. High frequency sound and the boson peak in amorphous silica. The European Physical Journal B-Condensed Matter and Complex Systems, 19(4):531–543, 2001.

[7] Fabien Leonforte, Anne Tanguy, JP Wittmer, and J-L Barrat. Inhomoge-neous elastic response of silica glass. Physical review letters, 97(5):055501, 2006.

[8] SN Taraskin, YL Loh, G Natarajan, and SR Elliott. Origin of the boson peak in systems with lattice disorder. Physical review letters, 86(7):1255, 2001.

[9] W Schirmacher, B Schmid, C Tomaras, G Viliani, G Baldi, G Ruocco, and T Scopigno. Vibrational excitations in systems with correlated disorder. physica status solidi (c), 5(3):862–866, 2008.

[10] Simon Gelin, Hajime Tanaka, and Ana¨el Lemaˆıtre. Anomalous phonon scattering and elastic correlations in amorphous solids. Nature Materials, 15(11):1177–1181, 2016.

[11] J¨urgen Horbach, Walter Kob, and Kurt Binder. Specific heat of amor-phous silica within the harmonic approximation. The Journal of Physical Chemistry B, 103(20):4104–4108, 1999.

[12] PG Klemens. The thermal conductivity of dielectric solids at low tem-peratures (theoretical). In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, volume 208, pages 108– 133. The Royal Society, 1951.

[13] Giulio Monaco and Stefano Mossa. Anomalous properties of the acoustic excitations in glasses on the mesoscopic length scale. Proceedings of the National Academy of Sciences, 106(40):16907–16912, 2009.

[14] B Ruta, G Baldi, F Scarponi, D Fioretto, VM Giordano, and G Monaco. Acoustic excitations in glassy sorbitol and their relation with the fragility and the boson peak. The Journal of chemical physics, 137(21):214502, 2012.

[15] Giulio Monaco and Valentina M Giordano. Breakdown of the debye ap-proximation for the acoustic modes with nanometric wavelengths in glasses. Proceedings of the national Academy of Sciences, 106(10):3659–3663, 2009. [16] G Baldi, M Zanatta, E Gilioli, V Milman, K Refson, Bj¨orn Wehinger, B Winkler, A Fontana, and G Monaco. Emergence of crystal-like atomic dynamics in glasses at the nanometer scale. Physical review letters, 110(18):185503, 2013.

[17] G Baldi, VM Giordano, G Monaco, and B Ruta. Sound attenuation at terahertz frequencies and the boson peak of vitreous silica. Physical review letters, 104(19):195501, 2010.

[18] W Dietsche and H Kinder. Spectroscopy of phonon scattering in glass. Physical Review Letters, 43(19):1413, 1979.

[19] Sajeev John and Michael J Stephen. Wave propagation and localization in a long-range correlated random potential. Physical Review B, 28(11):6358, 1983.

[20] B Velicky. Sound in granular matter: a case of wave propagation in random media. 2006.

[21] Antina Ghosh, Vijayakumar K Chikkadi, Peter Schall, Jorge Kurchan, and Daniel Bonn. Density of states of colloidal glasses. Physical review letters, 104(24):248305, 2010.

[22] HR Schober and C Oligschleger. Low-frequency vibrations in a model glass. Physical Review B, 53(17):11469, 1996.

[23] Corey S O’hern, Leonardo E Silbert, Andrea J Liu, and Sidney R Nagel. Jamming at zero temperature and zero applied stress: The epitome of disorder. Physical Review E, 68(1):011306, 2003.

[24] Melvin Lax. Classical noise. v. noise in self-sustained oscillators. Physical Review, 160(2):290, 1967.