System Size and Quench Temperature Dependencies in the Stiffness of Amorphous Solids

(1)

System Size and Quench Temperature

Dependencies in the Stiffness of Amorphous

Solids

Nigel van Herwijnen

July 7th, 2017

Institute for Theoretical Physics

University of Amsterdam

Report on the Bachelor Project Physics and Astronomy (15 EC) Conducted between 03-04-2017 and 07-07-2017

Written by Nigel van Herwijnen (student number 10330879) Supervised by dr. Edan Lerner

(2)

(3)

Abstract

Amorphous solids are known for their disordered structure. It is this disorder that gives rise to local excitations when the system is supercooled from a liquid state. Unfortunately, these excitations are difficult to study because of the presence of plane waves. In this paper, dipole forces have been applied to systems of soft spheres to simulate these local excitations. The stiffness as a result of this dipole force was calculated and was found to scale asκ ∝ L−2/5 in 2D andκ ∝ L−3/5in 3D, where L is the length of the system. Furthermore, a length of 11 particle diameters was found to be the size of the simulated local excitations in both 2D and 3D. Subsequently, the stiffness was measured for excitations generated in solids that were prepared with different quench temperatures and was found to be significantly lower for high quench temperatures.

(4)

(5)

Popular summary (Dutch)

Chaotische ballenbak

Hoewel sommige vaste stoffen een prachtige kristalstructuur hebben, zijn er ook vaste stoffen die een ongelooflijk chaotische structuur hebben. Over vaste stoffen met een kristalstructuur is al een hoop bekend, maar over het gedrag van de chaotische vaste stoffen zoals glas en gel, die amorfe stoffen worden genoemd, is nog veel onbekend. Als je kijkt naar een groot oppervlakte van een dergelijke stof, dan zie je er trillingen doorheen bewegen. Beeld je dit in als een golf door het water. Dit gedrag zien we ook in een vaste stof met kristalstructuur dus hier weten wij gelukkig al een hoop over! Maar als je de deeltjes van dichtbij bekijkt in een amorfe stof, dan zie je kleine ongeordende uitspattingen. Deze uitspattingen worden veroorzaakt door de chaotische structuur van de stof, maar wanneer ze precies ontstaan en hoe de rest van de stof erop reageert is een groot raadsel waarvan wij een klein stukje hebben kunnen oplossen.

Een voorbeeld van een gesimuleerde amorfe stof.

Om deze uitspattingen beter te begrijpen, hebben we een simulatie gemaakt van een bak met balletjes van twee verschil-lende maten. Dit simpele model gedraagt zich natuurlijk niet precies hetzelfde als glas, maar het lijkt er dusdanig veel op dat we er wel van kunnen leren! Door twee balletjes die tegen elkaar aan liggen uit elkaar te duwen, simuleren we een kleine uitspat-ting die we vervolgens kunnen bestuderen.

Wij vroegen ons af hoe de rest van de stof zou reageren op deze uitspatting, dus hebben wij gekeken naar hoe hard de andere balletjes terugduwden nadat wij een zetje hadden

gegeven. Lieten ze zich makkelijk verplaatsen en is het systeem heel soepel of was er heel veel energie voor nodig en is het systeem heel stijf? Hoe stijf een systeem is, is een belangrijke eigenschap omdat het ons bijvoorbeeld vertelt hoe energie zich door de stof verspreid. Door de grootte van de bak met balletjes te veranderen en de stijfheid te meten, hebben we geleerd hoe een amorfe stof reageert op een uitspatting. Maar misschien nog wel belangrijker, door de bak steeds kleiner te maken hebben we gevonden hoe groot zo een uitspatting eigenlijk is! Dat was namelijk nog helemaal niet bekend. Dat laat maar weer eens zien dat er over amorfe stoffen een hoop open vragen zijn die we één voor één proberen op te lossen.

(6)

(7)

CONTENTS

1 Introduction

Crystalline solids are well-known for the ordered structure their particles are in. Amorphous solids, on the other hand, have a very disordered structure [1]. It is this disorder that gives rise to local soft excitations when a shear is exerted, which, in turn, affects the properties of the solid [2]. It is already known that amorphous solids, such as glasses or gels, behave like con-tinuum elastic solids at large length scales [3]. The effect of local excitations at much smaller length scales, however, is still open for research. Recent research has, for example, shown that these soft modes are highly correlated with the relaxation processes in deeply supercooled liquids [4]. It is clear that studying these local excitations as a result of the disordered sys-tem in amorphous solids is of great importance in understanding the origin of the syssys-tem’s properties.

The problem in studying these very localized excitations is the presence of plane waves, or Goldstone modes [5]. The density per unit volume for these modes follow Debye’s theory as D(ω) ∼ ωd −1, whereω is the frequency of the mode and d is the dimension of the system [6]. This means that there is an abundance of plane waves covering the frequency range where local excitations are otherwise found. In addition, recent studies on structural glasses have found a density of low frequency vibrations that grows as D(ω) ∼ ω4, independent of spatial dimensions [7]. As visualized in figure1, the presence of these phonons make it hard to study the local excitation. Instead, this paper focuses on studying a dipole force, which creates an event similar to local excitations [8].

Studying the response of systems of soft spheres to a dipole force, enables us to study the behaviour of local excitations without the presence of plane waves. This paper shows the findings of a system size dependency in the stiffness of an amorphous solid as a result of a local excitation. Furthermore, a length scale has been found for the excitation itself. These have been found by measuring the stiffness of systems of different sizes. Subsequently, we discuss the influence of the parent temperature of the simulated systems on the elasticity of the systems.

(10)

INTRODUCTION

Figure 1: A visualized displacement field in a system of soft spheres caused by a shear on the system. The plane waves are clearly visible, but hide the presence of a local excitation which finds its origin in the disordered structure of the system. Lerner et al. [9]

(11)

THEORETICAL BACKGROUND

2 Theoretical background

In order to find a deeper understanding of amorphous solids at small scales, the elasticity of such solids has been studied numerically on systems of various sizes containing soft spheres. By means of a dipole force, a displacement field is generated in the solid from which the stiff-ness can be determined. The reason for calculating the stiffstiff-ness of the system is that it shows us how much energy is required for a displacement caused by a dipole force. By comparing the stiffness of systems of different sizes, a size dependency might be found in the energy shift of the system.

2.1 Energy and displacement

As a start, the total energy G of the system is formulated. In the following formula U denotes the potential energy, ~f is an applied force and~x −~x0 is the displacement of the particles.

The force can then be written in terms of a unit vector for the direction and a scalar for the magnitude.

G = U − ~f · (~x −~x0) (1)

~f = f ˆf (2)

The total net force on all particles in the system can be derived by taking the derivative of the energy with respect to the displacement.

~F = −∂G

∂~x (3)

= −∂U

∂~x + ~f (4)

The system is in equilibrium when the force is applied, so the total net force has to remain zero under variation of the magnitude of the applied force ~f .

(12)

THEORETICAL BACKGROUND ∂~F ∂f = 0 (5) ∂ ∂f (− ∂G ∂~x) = 0 (6) ∂2_G ∂f ∂~x+ ∂2_G ∂~x∂~x· ∂~x ∂f = 0 (7)

In the last step of the last equation the chain rule has been applied. The result is an extra term containing_∂f∂~x. This derivative of the displacement with respect to the magnitude of the applied force, shows us how the particles will move if the applied force is changed. Filling in equation4into equation7gives us

0 = ∂ 2_G ∂f ∂~x+ ∂2_G ∂~x∂~x· ∂~x ∂f (8) 0 = ˆf − ∂ 2_U ∂~x∂~x· ∂~x ∂f (9) ∂~x ∂f = ³ ∂2U ∂~x∂~x ´−1 · ˆf (10) ≡ M−1fˆ (11)

whereM is defined as the dynamical matrix. To find the actual displacement, which is called ~z from now on, the derivative of the space vector is multiplied by the magnitude of the force.

~z =_∂f∂~x· f (12)

= M−1f ˆf (13)

This equation is of use because the applied force is known and can therefore be solved by using the conjugate gradient method.

(13)

THEORETICAL BACKGROUND

2.2 Stiffness

To find the intrinsic property of the system that tells us how the system reacts to a dipole force, we determine the stiffness of the system. Firstly, the displacement vector~z and space vector ~x are redefined as follows,

~z = z ˆz (14)

~x =~x0+~z = ~x0+ z ˆz (15)

where~x0denotes the original place of the particles and~z the place the particles will move to

as a result of the dipole force. It should be noted that the derivative of the space vector with respect to the magnitude of the displacements only leaves the direction of the displacements.

∂~x

∂z = ˆz (16)

As explained before, we are actually interested in how the energy changes as a result of the dipole force. By taking the Taylor expansion of the potential energy, an expression of the energy in terms of the displacement is found.

U (z) ' U (0) +X i ∂U ∂~xi ·∂~xi ∂z z + 1 2 X i , j ∂~xi ∂z · ∂2_U ∂~xi∂~xj ·∂~xj ∂z z 2_{+ · · ·} ₍₁₇₎

The higher order terms in the expansion are left for now. The first order term vanishes, be-cause the derivative is evaluated at mechanical equilibrium. This results in only the zero’th order term and the second order term. The zero’th order term holds information on the en-ergy of the system at equilibrium, but since we are interested in the enen-ergy shift as a result of the dipole force it can be moved to the other side of the equation.

(14)

THEORETICAL BACKGROUND U (z) −U (0) ' 1 2 X i , j ∂~xi ∂z · ∂2_U ∂~xi∂~xj · ∂~xj ∂z z 2 ₍₁₈₎ =1 2( ˆz · M · ˆz)z 2 ₍₁₉₎ ≡1 2κz 2 ₍₂₀₎

The result is an equation for the energy shift as a function of the magnitude of the displace-ment. This is a simple second degree polynomial as shown in figure2where the curvature is determined by the stiffnessκ. If a system has a very high stiffness, the curvature will be large and a small displacement will cost a great amount of energy whilst a system with a very low stiffness will be easily displaced.

Figure 2: An example of the shift in potential energy as a function of the displacement of the particles in the system and the corresponding Taylor approximation. The curve has the shape of a parabola in which the curvature around zero is determined by the stiffness of the system.

Since we are not interested in the actual energy in the system, but how the energy shift changes when the system size changes, the stiffness has been determined for different system sizes. By taking the hessian of the normalized displacement field, we find the stiffness of the system defined asκ ≡ ˆz · M · ˆz.

(15)

ANALYSIS THROUGH SIMULATIONS

3 Analysis through simulations

In order to be able to find the size dependency of the stiffness of amorphous solids, the re-action of a system of soft spheres to a local dipole force is simulated. The computer model allows for accurate measurements and full control over the system that would otherwise be impossible.

3.1 The system

The systems used are computer-generated structural glasses. These binary systems consist of spheres of two sizes and were brought to zero temperature by a rapid quench after a short equilibration run in the liquid phase [2]. Instead of simulating solid boundaries at the edges of the box, periodic boundaries have been used. This prevents the boundaries from interfering with the system but instead is expected to behave like an infinite medium for large systems.

To find the size dependency, systems of different sizes have been used. In 2D these range from 25 particles per system to 1.6 million particles per system; in 3D the size ranges from 100 particles per system to 1 million particles per system. This wide variety of system sizes enables us to measure asymptotic behaviour if it is present. Because of the random nature of the packings, data has been collected from multiple systems per system size. Averaging over numerous systems, allows us to find accurate and correct data.

3.2 Measuring the stiffness

In this paper, we describe the research of the stiffness as a result of a dipole force. Two particles are chosen and on both the same, but opposite, force is exerted. This will result in the particles moving in opposite direction, simulating a perturbation in the system. As described in the theoretical background, the stiffness of the system can be obtained from the displacement field that results from this dipole force.

Only particles that are adjacent should be used as dipoles. Unlike particles that are neigh-bours, particles that are not adjacent would not naturally be excited and would therefore not be suitable to simulate a reaction with. The first step in finding all neighbours is finding the radial distribution functions. This tells us from which distance from the centre of a particle we can assume there to be no direct neighbours anymore. After finding all eligible pairs, the pairs used in finding the displacement fields are chosen at random.

(16)

Figure 3: An example of a system used for the numerical analysis. This particular 2 dimensional binary system contains a total of 400 particles. Lerner et al. [9]

One pair at a time, a dipole is made and the resulting displacement field is computed using the conjugate gradient method. After the displacement field is normalized, the stiffness is calculated using the dynamical matrix, as shown in equations19and20. This is done for a minimum of 1,000 dipoles per system, or all dipoles if there are less than 1,000 dipoles suitable for a dipole force. This process is repeated for multiple systems per system size, after which all stiffness values are averaged over.

3.3 Preparation protocol

While using rapidly quenched glassy configurations, one might question the influence of the temperature at which the system was quenched on the elasticity of the system. The liquids are supercooled before an instantaneous quench to zero temperature [7]. While the system is in a supercooled state, the particles in the system will be able to move around and interact with each other. If one looks at a macroscopic level, one will see the system naturally reorgan-ising, trying to find a lower energy state. Compared to higher temperatures, systems at lower temperatures will be moving around much slower and will not be able to overcome high en-ergy barriers. This will cause the system to slowly find an equilibrium before being quenched

(17)

Figure 4: An example of a displacement field (right) as a result of a dipole force (left) in a 2D binary system. Lerner et al. [9]

to zero temperature. Contrary to that, systems with higher temperatures will have particles moving around at higher speeds resulting in a very chaotic structure.

The stiffness has been measured of samples quenched at different temperatures to study the effects of their parent temperature T0on the elasticity of the system. All samples are 3

dimensional binary systems of 2,000 particles. To study the change in macroscopic elasticity, the bulk modulus B end the shear modulus G have been measured as well. These can be used to compute the Debye frequency using:

cL=p(B + 4G/3)/ρ (21) cT =pG/ρ (22) ωD = Ã 18π2(N /V ) 1/c3_L+ 2/c_T3 !1/3 . (23)

In the systems used, the mass densityρ and number density N_V are equal to 0.82. These quan-tities will be compared to study how the parent temperature affects the system at zero tem-perature.

(18)

RESULTS

4 Results

4.1 Radial distribution function

In order to find the maximum distance of a neighbour, radial distribution functions have been made. The results for 2D systems and 3D systems can be found in figure5. Although the nor-malization seems to be wrong in the 3D distribution functions, the shape of the distribution in all systems is clear and as expected. The cut-off values used in 2D and 3D for finding the neighbours of particles are shown in table1.

Table 1: Cut-off values in 2D and 3D as found in figure5.

2D 3D gss 1.3 1.3

gsl 1.5 1.5

gl l 1.7 1.7

Figure 5: Radial distribution functions of particles in a 2D binary system (left), and in a 3D binary sys-tem (right). gssdescribes the function for the distance between small particles, gsl between different

sized particles and gl l between large particles. Data has been collected from systems of 1600 particles

and average values have been calculated from 10 samples.

(19)

RESULTS

4.2 System size dependency

To find the system size dependency on the elasticity of the system, the stiffness has been calculated for systems of different sizes. The results of this can be seen in figure6. A decay in stiffness is clearly visible as the system size grows.

Figure 6: Stiffnessκ as measured in binary systems in 2D (left) and 3D (right) for different system sizes

L.

For very large systems the stiffness seems to converge to an asymptote. In order to find a qualitative description for the stiffness, a power function has been fitted of the form

κ = κ∞+C1L−β (24)

where C1andβ are constants and κ∞is the value the stiffness converges to for an infinitely

large system. The result is plotted in figure7.

The graphs show that in 2D the stiffness followsκ ∝ L−2/5for L > 11 and that in 3D the stiffness followsκ ∝ L−3/5for L > 11. In both cases, the stiffness deviates from the fit for small system sizes smaller than a length L = 11 particle diameters.

4.3 Parent temperature dependency

In order to find the influence of the parent temperature on the elasticity of the system, the stiffness has been measured for systems with different parent temperatures. The result can

(20)

RESULTS

Figure 7: Stiffnessκ as measured in binary systems in 2D (left) and 3D (right) for different system sizes

L. In 2D the asymptote has been found to beκ∞= 16.3; in 3D κ∞= 31.6. A line has been plotted of the

form of equation24.

be seen in figure8. The stiffness that has been measured for systems with a higher parent temperature is significantly lower than that of systems with a lower parent temperature.

In addition, the shear modulus and bulk modulus have been measured to calculate the Debye frequency ωD following equation23. The ratioκ/ω2_D for systems at different parent

temperatures is plotted in figure9.

(21)

RESULTS

Figure 8: The stiffnessκ as measured in 3D binary systems which contain 2000 particles for different parent temperatures T0.

Figure 9: The ratio between the stiffness and the Debye frequency squared as measured in 3D binary systems which contain 2000 particles for different parent temperatures T0.

(22)

DISCUSSION

5 Discussion

5.1 System size dependency

There is a significant difference in stiffness measured between small and large systems in both 2D and 3D. The lowest and highest points differ in 2D with a factor of 4 and in 3D with a factor of 1.3, showing a clear elastic decay further from the dipole force. The asymptotic behaviour we see in figure6is what we would intuitively expect from an elastic solid. From a certain distance from the dipole force the particles would no longer undergo any interaction as a result of this force, since it would already have decayed.

In figure 7this asymptote has been filtered out to find the slope at which the stiffness decays. Although the 2D systems undergo a much larger decay, the 3D systems decay faster. Where the 2D systems decay with L−2/5, the 3D systems decay with L−3/5. The reason for this steeper decay can originate in the amount of places a particle can move to. There are simply more empty spaces around particles in 3D systems than in 2D systems, which will result in more degrees of freedom and a system that will find a favourable state more easily.

For smaller systems, the stiffness deviates from this function and explodes. Remember that for a small region around the dipole force, the continuum description breaks down. As can be seen in figure 4, a disordered core is created which does not behave elastically like the rest of the system. When the system is made smaller, the boundaries of the system will clash with the core of the dipole force at some point. Below this system length L the elastic behaviour of the system is no longer calculated, but instead measurements are made to the disordered core itself. The size of the system at this length could therefore be interpreted as the size of the core and is found to be around 11 particles in both 2D and 3D.

To find the stiffness of a system, the average value of multiple dipoles per system has been calculated. For each data point in the figures in chapter 4: results, the average stiffness has been calculated from multiple systems. This means that for the smaller systems millions of dipoles have been used to calculate the average stiffness. Unfortunately, this amount of data points was not possible for the larger systems because of the heavy computational power needed to perform the calculations. This can explain the small deviation of some data points for large system sizes in figure7and can be solved by simply using more dipoles to calculate the mean stiffness. Evaluating even larger systems could help as well, especially in 3D, to confirm the slope in the decay of the stiffness. To pinpoint the exact size of the core of a local excitation, more system sizes need to be evaluated around the size measured at 11 particle

(23)

DISCUSSION

diameters. These are small systems and are therefore easy to do calculations on. The exact value is now hard to determine, but more data points would help enormously.

5.2 Parent temperature dependency

From figure8can be read that the stiffness of the system drops with a factor of 2 at high parent temperatures with respect to lower parent temperatures. This significant change is due to the configuration of the system after the cooling sequence and before the instantaneous quench. Particles in systems with high temperatures are moving at higher velocities creating a more chaotic solid. Because of the high energy in the system, the particles do not need to find a low energy state before the quench. When a dipole force is put onto such a system, the system is still in a high energy state and will therefore easily find a relatively low energy state. A system quenched at a lower temperature, on the other hand, will have had the time to find a low energy state during the cooling process. It will therefore be harder for these systems to find a low energy state when undergoing a dipole force.

Although we have been studying the effects of local events on the energy of the system, it is interesting to know how the parent temperature influences the energy at a macroscopic level as well. This has been done by plotting the ratio of the stiffness with respect to the Debye fre-quency squared as shown in figure9. The Debye frequency is the highest allowed frequency in a solid, so a change in Debye frequency means a change in the allowed phonon energies in the solid [6]. What can be noticed in the figure is that the ratio difference at different par-ent temperatures is significant. From this it can be concluded that the local excitations are influenced more by the parent temperature than the phonons.

5.3 Continuum calculation

In order to confirm the size dependency that has been found, a continuum calculation should be performed. This has been tried, but unfortunately has not resulted in a conclusive expres-sion for the stiffness of the system.

(24)

CONCLUSION

6 Conclusion

By numerical analysis on simulations of systems of soft spheres, a system size dependency in the stiffness of the system has been found when a dipole force is applied. By applying a dipole force, a local excitation has been simulated. The system follows a stiffness that goes asκ ∝ L−2/5in 2D andκ ∝ L−3/5in 3D, where L is the length of the system. Furthermore, a length scale for the core of the excitation has been found to be L = 11 particle diameters.

It is advised to calculate more data points for larger systems to confirm the slope of the decay of the stiffness. Although these are computational very heavy, it would be greatly ben-eficial to the fit if more data was available. More data at smaller system sizes would help as well to find a more accurate value for the length scale of the core of the excitation. Naturally, it would be best to confirm the functions found for the stiffness with a continuum calculation.

Furthermore, it has been found that the parent temperature of the system has a great in-fluence on the stiffness of the system. This is because the system will find a low energy state when it is cooled to a low temperature before being quenched, whilst quenching at high tem-perature will leave the system in a high energy state. This will result in a system that is either easily moved, or hard to move by a dipole force. Finally, it has been found that the system is influenced more at small wavelengths than at large wavelengths by the parent temperature by looking at the ratio of the stiffness and the Debye frequency.

(25)

CONCLUSION

Acknowledgment

I would like to thank my supervisor dr. Edan Lerner for his guidance during this project and the never ending stream of new ideas and interpretations he provided with much enthusiasm. This made the project both interesting and enjoyable to do.

(26)

REFERENCES

References

[1] E. Lerner, E. DeGiuli, G. Düring, and M. Wyart. Breakdown of continuum elasticity in amorphous solids. arXiv:1312.2146v1, 2013.

[2] E. Lerner, G. Düring, and E. Bouchbinder. Statistics and properties of low-frequency vi-brational modes in structural glasses. Physical Review Letters, 2016.

[3] F. Leonforte, R. Boissière, A. Tanguy, J. P. Wittmer, and J.-L. Barrat. Continuum limit of amorphous elastic bodies. iii. three-dimensional systems. Phys. Rev. B, 72:224206, Dec 2005. doi: 10.1103/PhysRevB.72.224206.

[4] A. Widmer-Cooper, H. Perry, P. Harrowell, and D. R. Reichman. Irreversible reorganization in a supercooled liquid originates from localized soft modes. Nature Physics, 4(9):711–715, 2008.

[5] J. Goldstone, A. Salam, and S. Weinberg. Broken symmetries. Physical Review, 127(3):965, 1962.

[6] P. Debye. Zur Theorie der spezifischen Wärmen. Annalen der Physik, 344:789–839, 1912. doi: 10.1002/andp.19123441404.

[7] E. Lerner. Effect of instantaneous quenches on the density of vibrational modes in model glasses. arXiv:1705.01037v1, 2017.

[8] F. Leonforte, A. Tanguy, J. P. Wittmer, and J.-L. Barrat. Continuum limit of amorphous elastic bodies ii: Linear response to a point source force. Physical Review B, 70(1):014203, 2004.

[9] E. Lerner, G. Düring, E. DeGiuli, and M. Wyart. Bond-space perspective on the elasticity of disordered solids. Lecture Roskilde University, March 2014.

(27)

System Size and Quench Temperature Dependencies in the Stiffness of Amorphous Solids