Elasticity of disordered systems with rigid backbones

Bachelor Physics and Astronomy

Bachelor Thesis

Elasticity of disordered systems with rigid

backbones

Marten Folkertsma

July 8, 2017

Supervisor:

Edan Lerner

Assessor:

Philippe Corboz

Abstract

Fibre networks can be described as networks with very strong springs and weak angle in-teractions. The average number of connections per node in the fibre networks is below to 2d isostatic point. To investigate the physics behind these networks we looked into the limiting case where stretching is not allowed, calling it the rigid backbone network. First we describe how we can find the vibrational modes. From numerical analysis we determine that the energy of the modes is not strongly dependent on the coordination and the delocalisation increases with increased coordination. We also find that the rigid backbone network gives a good description of the bending modes in the fibre networks. An analytic description of the shear modulus of such a network is given. We derive that if the coordination comes close to the isostatic point the shear modulus G increases with G ∼ 1

δz where δz is given as δz = 2d − z, twice the dimension minus

the average coordination. Numerical data are shown which agree with this analysis. We conclude by describing the equation of motion and signalling a breaking of the constraints at second order approximation.

Contents

1 Dutch summary 2

2 Introduction 3

3 The model 5

3.1 Generating the networks . . . 5

4 Finding the vibrational modes 6

4.1 Numerical analysis of the vibrational modes . . . 8

5 Deriving the shear modulus 12

5.1 Numerical Analysis . . . 15

6 Deriving equations of motion 16

7 Discussion 17

7.1 Network size . . . 17 7.2 Breaking constraints at second order . . . 17 7.3 Scaling in the energies of the system . . . 18

1

Dutch summary

De physica van goed geordende stoffen, de kristalstructuren, is al reeds bekend. Denk hierbij aan materialen als metalen, kristallen en zouten. Vanwege hun welgeordende molecuulroosters is het relatief makkelijk om de physische eigenschappen van deze materialen te beschrijven. We kunnen precies uitrekenen wat de buigzaamheid van bijvoorbeeld metaal is.

Buiten deze welgeordende stoffen vinden we ook stoffen in de natuur die geen geordende mole-cuulstructuur hebben, maar waarbij de moleculen bijna random verspreid liggen in het materiaal. Dit noemen wij amorfe stoffen. Hierbij is het een stuk moeilijker om de fysica van de materialen te beschrijven. In dit onderzoek zijn we dieper in deze stoffen gedoken om te kijken hoe we de buigzaamheid van zulke stoffen kunnen beschrijven.

Een voorbeeld van die soort stoffen zijn netwerken van vezels. Vezels vinden we overal in het menselijk lichaam, bijvoorbeeld in spieren, en ze geven stabiliteit aan het menselijk lichaam. Daarom is het belangrijk om meer over de buigzaamheid van deze vezelnetwerken te weten. Met netwerken bedoelen we hier systemen die beschreven worden met punten en verbindingen. Een voorbeeld van een netwerk zijn metaalroosters, dan zijn de moleculen de punten en de molecuulverbindingen de verbindingen tussen de punten.

Bij vezelnetwerken gaat het om random georganiseerde netwerken. Er is geen duidelijke structuur in waar de punten liggen en in welke punten met elkaar verbonden zijn. Een belangrijke eigenschap die we kunnen bekijken is de gemiddelde connectiviteit van een netwerk. De gemiddelde connectiviteit wordt gegeven door het gemiddeld aantal verbindingen dat aan een punt vast zit. Dit noemen we ook wel de co¨ordinatie van het netwerk.

Een vezel netwerk is een voorbeeld van zo’n chaotisch netwerk, waarvan de verbindingen vaak beschreven worden als stijve springveren. In zo’n netwerk zijn dus alle punten met springveren aan elkaar verbonden. Deze springveren zijn echter zo stijf dat het heel moeilijk is om ze in te drukken of uit te rekken, dit kost veel energie. Naast het uitrekken en indrukken van de veren is er echter een andere manier hoe de punten zich in het netwerk kunnen verplaatsen. Ze kunnen om elkaar heen draaien. Hierdoor verandert de lengte van de verbindingen niet maar de hoek tussen de verbindingen wel. In een vezelnetwerk kost zo’n beweging veel minder energie dan een beweging waarbij veren worden uitgerekt. Hier noemen we dit soort bewegingen draaibewegingen van het netwerk.

Deze draaibewegingen hebben een grootte impact op de buigzaam van een materiaal. Als er meer mogelijke draaibewegingen zijn in het netwerk kan het materiaal makkelijker buigen. Om de impact van deze buigbewegingen te onderzoeken bekijken we hier het limietgeval waarbij de verbindingen ¨

uberhaupt niet kunnen worden uitgerekt of ingedrukt. We vervangen als het ware alle springveren door solide staven.

Voor dit soort materialen geldt dat als de co¨ordinatie boven een bepaald punt komt, namelijk twee keer de dimensie, er geen draaibewegingen meer mogelijk in het netwerk. Dit punt noemen we het isostatische punt. Omdat er geen draaibewegingen mogelijk zijn kan het hele netwerk niet meer buigen.

In dit onderzoek geven we een analytische beschrijving van hoe de buigzaamheid afhangt van de co¨ordinatie van zo’n netwerk. Daarnaast hebben we zo’n systeem gesimuleerd. Uit beide analyses vinden we dat de buigzaamheid van zo’n materiaal heel hard afneemt als de co¨ordinatie dicht bij het isostatische punt komt.

2

Introduction

The mechanics of well ordered systems, such as crystalline solids, are well understood. In contrast, the mechanics of amorphous solids are a lot harder to grasp because of lack of order in the system. Examples of highly irregular, or even randomly organized systems are found in materials like glasses and proteins. To better understand these systems they are usually subdivided into two categories: packings of spheres and networks. In this research we will only be looking at systems which resemble networks. A characteristic that all many body networks share is the average connectivity z, the average amount of bonds per node. We call the average connectivity the coordination of network. For randomly organized systems this is a very important characteristic, for there are no symmetries in the system that can be used to describe the physics of the system. If a physical network has a coordination below a certain threshold zc = 2d, twice the dimension, the system becomes unstable,

i.e. has some degrees of freedom. Maxwell [1] describes how frames with higher coordinations lose all degrees of freedom and therefore become rigid. Networks act with respect to coordination in the same way as Maxwells frames. We find the same rigidity if the average connectivity is higher than this threshold z > zc, assuming homogeneous distribution of connections.

If we consider a network with springs as bonds and with a coordination lower than z < zcwe find

extra degrees of freedom in the system. These extra degrees of freedom give rise to vibrational modes with zero energy cost, the network is able to deform without stretching any bond. We call these modes floppy modes. These floppy modes utilize the degrees of freedom, created by the lack of extra bonds, by rotating particles in such a way that none of the bonds stretch. This causes instabilities in the network, because it can freely deform without any energy cost. This clearly has an enormous impact on the stiffness of materials.

To stabilize these kind of spring networks with low coordinations we can add bending potentials to the network. This bending potential is added between every pair of bonds connected to the same node. This adds an energy cost to rotation around the nodes, causing the floppy modes to disappear. This stabilizes the system, because there are no energy free movements left. The energy costs of bending are usually much lower than the cost for stretching.

Figure 1: On the left showing stretching movement, where the bonds need to be stretched, and on the right bending movement, where the bonds rotate around the node.

Biological materials like fibre networks and protein structures can be described as networks of springs with coordinations below the threshold zc [2]. When you look at how the stiffness of the

material increases with the coordination we find a very non-linear pattern. These systems can best be described as networks consisting of strong springs. We typically find that bending in such networks has a much lower energy cost than stretching. These systems will therefore first try to utilize bending as much as possible before stretching. By increasing the coordination of the network less bending is possible and more stretching is needed. This causes a non-linear behaviour in the stretching of the material [3]. To better understand this behaviour we will be looking at a limit case in which the network does not consist of springs but is a network of rigid bones. Using this we want to find the dependence of the shear modulus on the coordination.

The rigid backbone networks are similar to the spring networks only we exchange the springs for rigid bones, with constraint that the rigid bones cannot be stretched at all. These networks disallow stretching of the bonds completely. This is a hard constraint on the movement of the particles. As in the spring network a bending potential is added to stabilize the system. The only deformations allowed will be those involving bending. We call this network the rigid backbone network.

Figure 2: Image on the left showing a spring network (a) and on the right the rigid backbone network (b). Network (a) contains bending and stretching energies, network (b) only bending.

This paper will continue as follows: first of all in chapter three there will be a description of the model and how it works. Secondly in chapter four we will look at the density of states of such a system trying to find the vibrational modes. This involves applying the hard constraint of the system on the movements of the nodes. Thus the vibrational modes will only allow movements which change the angles between the bonds and not allow movements which stretch or compress the bonds. To find these modes we will first do a mathematical description followed by a numerical analysis of the states. This will be compared to a system of very stiff springs to see if this is a proper way of describing the transition from below to above the isostatic point(z = zc). Next in chapter five we will analyze

how the system behaves under a shear stress. Looking at how the shear modulus increases versus the relative coordination δz, the difference from the isostatic point. Here we say: δz = z − 2d thus δz is the difference from the isostatic point. Again we will approach this with a mathematical analysis followed by numerical calculation. Finally in chapter six we look at how such a system would evolve over time if we would put it in a thermal bath, seeing how the different vibrational modes evolve over time.

3

The model

The model we use has Nn nodes and Nb bonds. These nodes and bonds make up a 2d network. The

network consists of rigid backbone, i.e. the length of a bond cannot change: ∂rij

∂Rk = 0. The potential

energy of the system will only be dependent on the change of the angles between the bonds. We define it as follows: U = Na X i 1 2(θi− θ 0 i) 2 ₍₁₎

Here we sum over all the angles between every two bonds connected to the same node and next to each other which are labeled by i. Where θ0

i is an initial angle between two connections, θi is the

angle after a deformation and Na is the total amount of angles. Important to note is that there is no

initial stress in the system. From here we will start deriving the formulas for the vibrational modes and the shear modulus.

3.1

Generating the networks

To do this research it is necessary to generate homogeneous random networks with a particular mean connectivity. In this research we will be using a method developed by Edan Lerner and Jacques Zylberg. To get our desired networks we first generate random networks with high coordinations. This is done by generating a well-compressed disordered packing of soft spheres. For a 2D network N soft discs are generated. N/2 of these are small discs and N/2 are large discs with the ratio between the radii of 1.4. When the discs are packed into a space of L2_{,the size of the 2d network, the centers}

of the touching discs are connected. After this, the discs are removed and all that is left is a highly coordinated network with coordination zl.

To get to the desired coordination zf, which is smaller than zl, bonds are iteratively removed. A

bond between the ith and jth node is denoted by eij and the number of bonds connected to the ith node by zi_{. At every iteration, the set of Γ of bonds e}ij _{with the highest sum of ranks z}i_{+ z}j _is

identified. Then, a subset of Γ is built by finding the bonds eij _{∈ Γ that have the smallest absolute}

difference of rank |zi_{− z}j_{|. A random member of this subset is removed from the network. This}

procedure is repeated until the target mean connectivity zf is reached.

For the numerical analysis we used a combination of c, c++, matlab and python using the formulas derived below. The networks all had 1024 particles at different coordinations z.

Figure 3: Constructing packing-derived initial random networks. (a) start with an athermal packing, system at zero kelvin, of soft discs. (b) connect the centers of interacting discs by bonds and (c) then get rid of the particles all together.

4

Finding the vibrational modes

Vibrational modes are displacements of the system that cause a force in the direct opposite of the displacement. This causes an oscillation in the direction of the displacement. Every displacement of the system can be described as a linear combination of vibrational modes. If we look at networks with springs that have energies in stretching as well as bending of the bonds we get the following potential: U = Na X i 1 2l0 κi(θi− θi0) 2₊ Nb X α 1 2l2 0 µα(rα− rα0) 2 ₍₂₎

Here all the information of interaction between the particles is in the energy of the system. To find the vibrational modes we have to look at the following expansion:

|F i = −∂U ∂Ri − ∂ 2_U ∂Ri∂Rj |Rii (3)

Because there is no initial stress in the system and it is in mechanical equilibrium we know _∂R∂U

i = 0.

What we are left with is what we call the Hessian operator M [4]:

M = ∂

2_U

∂ ~Ri∂ ~Rj

(4)

Because of the calculus involved we can see that this matrix is symmetric and so diagonalizable. The eigenstates that we find from the Hessian are displacements that give rise tot restorative forces. If we would diagonalize the Hessian we would find the following possible solutions:

M |χi = ω2|χi = −|F i (5)

The eigenmodes of the Hessian give us a differential equation just like that of an ordinary spring. Just as a spring this differential equation gives the set of modes in which the system can vibrate. Only these vibrations are movements of all the nodes insteat of the movement of one spring.

Using this method we get an easy differential equation and clearly an oscillating system with frequency ω emerges.

In this research however we are dealing with a network consisting of a rigid backbone and a potential containing only bending energies. As described in the chapter about the model the potential energy is given as follows:

U = Na X i 1 2(θi− θ 0 i)2

The explicit formula for M for bending becomes:

Mnm= ∂2_U ∂ ~Rn∂ ~Rm = κX ij ∂θijk ∂ ~Rn ∂θijk ∂ ~Rm + (∆θijk) ∂∂θijk ∂ ~Rn∂ ~Rm = κX ij  (δin− δjn) rij  _−ny ij nx ij  +(δkn− δjn) rkj  _−ny kj nx kj  ∗ (δim− δjm) rij  −ny_ij nx ij T +(δkm− δjm) rkj  _−ny kj nx kj T! +(∆θijk)2(δin− δjn)(δim− δjm) " 1 r2 ij !  _nx ijn y ij −nxijnxij+ 1/2 ny_ijny_ij− 1/2 −ny_ijnx ij # +(∆θijk)2(δkn− δjn)(δkm− δjm) " 1 r2 kj !  _nx kjn y kj −n x kjn x kj+ 1/2 ny_kjny_kj− 1/2 −ny_kjnx_kj # (6)

Here the energetics of a system do not contain enough information to calculate the vibrational modes in the same way. The hard constraints of having rigid bones are not found in the potential energy. If we would use the same method, this would give rise to non-allowed stretching modes. Just expanding the system would cost no energy according to the bending energies and would generate a zero mode, but is not allowed in the system with rigid bones. This constraint must still be applied on the system to ensure that modes involving stretching disappear. To do so, we need to measure the change of

the lengths of the bonds dependent on the variation in the position of the nodes. So, in case of a displacement the lengths of the bonds must stay the same. Using Taylor expansion we find the following:

rα≈ rα+

∂rα

∂ ~Ri

( ~Rj− ~Rj0) (7)

Here rαis the length of a bond and α is the index for all the bonds. Using that rα can’t change

we get the following:

∂rα

∂ ~Ri

( ~Rj− ~Rj0) = 0 (8)

We call this operator S:

Siα=

∂rα

∂ ~Ri

(9)

The zero modes of this matrix are the movements of nodes which leave the lengths between the bonds constant and thus respect the constraints of the rigid bones. Its dimensions are given by: d ∗ Nb

and d ∗ Nn with Nb the amount of bonds, Nn the amount nodes and d the dimension of the system.

S is not square, but ST_{S is, it has dimension d ∗ N by d ∗ N ; the same dimensions as the Hessian.}

This gives us the following set of zero modes of S :{ψ : ST_{S|ψi = 0}. If we now find eigenmodes of}

the Hessian consisting of linear combinations of these modes, we find the vibrational modes:

χ =X

l

alψl (10)

We now want to find only the allowed frequencies. For this we use the Hessian with the energetics of the system and the allowed movements from the S operator:

∂U ∂al = ∂U ∂ ~Ri ∂ ~Ri ∂al = ∂U ∂ ~Ri ψil (11) ∂ ∂ak  ∂U ∂al  = ψj_kMijψil+ ∂U ∂ ~Ri ∂ψi_l ∂ak (12) ∂U ∂ ~Ri ψi l ak = 0 (13) Akl:= ψjkMijψil (14)

Here the ψs are dimensionless per construction thus Akl has the same units as the Hessian. The

eigenvalues of Akl must therefore also be frequencies:

Akl|νqi = ω2q|ν

q_{i (15)}

We can see Akl as a reduced version of the Hessian in which only allowed movements are left.

Akl contains the energetics of the system and the constraints on the movements of the nodes. The

frequencies given from Aklbelong to displacements of the system which leave the lengths of the bonds

constant and so are allowed eigenfrequencies in our system. Now we want to find the vibrational modes. For this we need an extra operator to remove the modes from M that are not allowed, but not reduce M in size:

P := I − ST SST
−1S (16)

This operator is zero when multiplied by a non-zero mode of S. It only contains the zero modes of S. This can be derived in the following way:

SST|µi = λ2_|µi

Where the different λ are the non-zero eigenvalues of ST_S.

STS|ψi = λ2|ψi S|µi = λ|ψi

SST
−1= 1 λ2|µihµ|

I =X

λ

|ψλihψλ| With these equations it is shown that:

P = I − ST SST
−1S =X λ |ψλihψλ| −X λ>0 |ψλihψλ| =X λ=0 |ψλihψλ| (17) Taking the eigenstates of the reduced matrix Akl and the zero modes of S we can construct the

following vector:

φq=X

l

ν_lqψl (18)

Now it is possible to show the following:

(I − ST SST
−1S)M (I − ST SST
−1S)|φqi = ω2 q|φ

q_{i (19)}

By now combining equations (14), (15), (17), (18) and (19) we can show the following:

P M P |φqi = (I − ST SST
−1S)M (I − ST SST
−1S)|φqi =X l |ψlihψl|MX k |ψkihψk|X k0 |ψk0i|ν_kq0i =X l |ψl_iA kl|ν_kqi = X l ω_q2|ν_lqi|ψl_{i = ω}2 q|φ q_i

These |φq_{i are the vibrational modes that respect the constraints of the rigid backbone.}

4.1

Numerical analysis of the vibrational modes

While exploring the vibrational modes we did find a projection operator P which projects out all the states that involve stretching and leaves us with only allowed states of the Hessian. Here we used dimensional analysis to find that the frequencies from Akl are the same as the frequencies from

the Hessian. However the vibrational modes found from the P M P operator are not the same as the eigenstates of the Hessian. Therefore it is important to compare the limiting case, with no movements in the bones, with a case of very strong springs. For this we use the model of equation (2) with κi

l0µi = 10

−6_{. If we look at the frequency distribution of a system like this we get the following:}

Figure 4: Display of the eigenstate distribution of the Hessian of a network with springs and bending interaction with z = 3.2

This is a system where stretching the bonds is allowed but it costs a lot of energy. We clearly see that this causes the modes to split into two categories: modes with stretching and modes without stretching. The modes with high frequencies are the modes that cost the most energy because of the stretching of the bonds. The lower part frequencies are modes of low energy cost. These modes utilize the degrees of freedom associated with rotating the bonds instead of stretching them. If we now compare this with the system of rigid bars we see that they are quite similar:

Figure 5: Comparison between a network containing springs and a network containing rigid bones. The same network with 1024 nodes and average connectivity z = 3.2 was used. In one case we allowed stretching and added a stretching interaction in the Hessian. In the other case stretching was not allowed. From the stretching case we only see the low energies.

This graph clearly shows that the frequency distribution is quite similar in the case of rigid bones and the stiff spring network. We clearly see a frequency gap between the translational zero modes and the first bending modes in both cases. The peak of the distribution is also almost the same and the upper bound is as well. We can conclude that the vibrational modes of the network with rigid bones are quite similar to the the vibrational modes of the stiff string network.

We also find a high correlation between the delocalisation and the coordination δz, as also shown by D¨uring et al. [5]. At the highest coordinations we find complete delocalisation in the material.

Figure 6: Showing delocalisation increase with increased coordination. In all cases the networks had 1024 nodes but different coordinations. Shown are the vibrational modes corresponding with the highest frequencies found in the networks.

It is clear that the states get less localized when the coordination increases. Displacements spread themselves over a larger part of the system, involving more nodes. With an increase of δz the degrees of freedom decrease. With less degrees of freedom deformations must spread more over the entire network causing the system to be less localized. At coordination 3.9 we see that the delocalization of the vibrational modes becomes comparable to the size system .

Figure 7: Frequency distribution of networks at different coordinations. When the coordination in-creases less vibrational modes are possible in a network. Therefore, to get balanced statistics we used more networks at higher coordinations, all with 1024 nodes. We used 100, 200, 600, 800 networks for coordinations 3.2, 3.6, 3.8, 3.9 respectively.

This displays the density of frequencies at different coordinations. The zero modes have been removed for clearer statistics. Per network there are always two zero modes, one of every node moving to the right and one of every node moving upwards. From this graph we clearly see that the peak of the distribution moves to the left if δz increases. The frequencies stay relatively the same. There is a small movement of the location of the peak and a small decrease of the width. In other words, the peak becomes more and more narrow and there are no singularities forming:

Figure 8: Fit of highest and lowest frequencies found in the frequency distribution shown in figure 4.1.

It is clear that the lowest frequency increases and the highest frequency decreases with increasing δz. The energies in the system do not increase much even though there is a delocalization in the system.

5

Deriving the shear modulus

Materials above the isostatic threshold are marginally stable with small deformation, because all the degrees of freedom are balanced by the constraints on the system. For lower coordinations this is not the case. There are still degrees of freedom left which can be utilized for some small deformations. A measurement for the rigidity of these materials is the shear modulus. The shear modulus is a measurement of the stiffness in a material, i.e. the resistance of a material to a small shear. Biological materials like fiber networks and protein structures usually have lower coordinations than the isostatic threshold. Typical coordinations for these kind of materials are between 3 and 4, well below the 2D and 3D threshold [6].

Broedersz et al. show a very non-linear correlation between the coordination of a spring network and the shear modulus [7]. Broedersz et al. used spring networks with energy in stretching and bending of the bonds. To get more insight in this non-linear behaviour we use networks with only bending and no stretching at all. We first derive an analytical formula for the shear modulus and than use this description to test it numerically.

For the derivation we consider a bulk material in mechanical equilibrium and apply a small shear on it. The small shear can be seen as a transformation of coordinates ~Rk → H ∗ ~Rk with H being:

H =1 γ

0 1



(20)

Full derivatives with respect to shear are given as follows:

d dγ = ∂ ∂γ+ ˙ ~ R ∂ ∂ ~R (21)

First we want to find the shear modulus dependent on the energetics and the kinematics of the system. Here we start wit the shear stress which is given by:

σ = 1 V dU dγ = 1 V  ∂U ∂γ + ∂U ∂R| ˙Ri  (22)

In which V is the volume of the system. The shear modulus is defined as the derivative of the shear stress to a shear:

G ≡ dσ

dγ (23)

This gives the following shearmodulus:

G = 1 V  ∂2_U ∂γ2 + 2  _∂2_U ∂γ∂R     ˙ R  + +h ˙R| ∂ 2_U ∂R∂R| ˙Ri +  ∂U ∂R     ¨ R 

Here ∂U_∂R is zero because the system is initially without any stress:

G = 1 V  ∂2_U ∂γ2 + 2  _∂2_U ∂γ∂R     ˙ R  + h ˙R| ∂ 2_U ∂R∂R| ˙Ri  (24)

We see that the shear modulus depends on the energetics of the system, derivatives of the potential energy, and the relocation of the nodes. From the analysis of vibrational modes we know that there is no singularity in the energies. Everything involving the energies, the derivatives of U, behaves normally. The energy does not depend in any special way on δz. The asymptotic behaviour of the shear modulus must be in the movement of the nodes, which does involve the hard constraints. The leading terms will be the terms involving | ˙Ri. We need to find an expression for | ˙Ri.

Mechanical equilibrium causes the net force to be zero, so the forces in the rods are exactly in balance: ~ F = −∂U ∂ ~R + ~F (rods)_{= −}∂U ∂ ~R+ X α ∂rα ∂ ~Rfα (25)

Here α denotes a bond between the ith and jth particle and fα is the force in a specific bond.

Furthermore we have that ∂U

∂ ~R = 0, because the system is created without any initial stress, and so

P

α ∂rα

∂ ~Rfα= 0 as well.

d ~F dγ = 0 = − ∂2U ∂ ~R∂γ − ˙ ~ R ∂ 2_U ∂ ~R∂ ~R+ X α ∂2rα ∂ ~R∂γfα+ X α fα ∂2rα ∂ ~R∂γ ˙ ~ R +X α ∂rα ∂ ~R ˙ fα (26)

Here we know that fαis zero and we can write it in matrix notation:

− ∂ 2_U ∂ ~R∂γ − ˙ ~ R ∂ 2_U ∂ ~R∂ ~R+ X α ∂rα ∂ ~R ˙ fα= 0 −     ∂2_U ∂ ~R∂γ  − M | ˙Ri + ST| ˙f i = 0 (27) −SM−1     ∂2_U ∂ ~R∂γ  + SM−1ST| ˙f i = S| ˙Ri (28)

We also know that the length of the bones should not change when we apply the shear on the system: drα dγ = 0 = ∂rα ∂γ + ˙ ~ R∂rα ∂ ~R     ∂r ∂γ  + S| ˙Ri = 0 S| ˙Ri =     ∂r ∂γ  (29)

Now we can combine equations (28) and (29) to find a solution for | ˙f i:

SM−1ST| ˙f i = SM−1     ∂2U ∂ ~R∂γ  +     ∂r ∂γ  | ˙f i = SM−1ST
−1 SM−1     ∂2_U ∂ ~R∂γ  + SM−1ST
−1     ∂r ∂γ  (30)

Using equations (28) and (30) we can solve for | ˙Ri:

M | ˙Ri = −     ∂2_U ∂ ~R∂γ  + ST SM−1ST
−1 SM−1     ∂2_U ∂ ~R∂γ  + ST SM−1ST
−1     ∂r ∂γ  | ˙Ri = −M−1     ∂2U ∂ ~R∂γ  + M−1ST SM−1ST
−1SM−1     ∂2U ∂ ~R∂γ  + M−1ST SM−1ST
−1     ∂r ∂γ  | ˙Ri = −M−1  ST SM−1ST
−1  SM−1     ∂2_U ∂ ~R∂γ  −     ∂r ∂γ  +     ∂2_U ∂ ~R∂γ  (31)

To find the scaling of | ˙Ri we need to compute h ˙R| ˙Ri. For this we must look at what happens with the individual matrices and vectors when δz goes to zero.

From M we know that it has modes with very small, or even zero, eigenfrequencies. Just com-pressing or expanding the system will not change the angles between the bonds and thus will not cost any energy. These modes heavily utilize the fact that there is no energy in stretching, In other words, these modes are the ones that break the constraints most heavily. They give a force field that only stretches and compresses the bonds and thus will not be present in our system. Because these modes are disallowed we can disregard them for the scaling of | ˙Ri. The modes of M that respect constraints are the modes that do change the angles of the system and thus cost a finite amount of energy. The eigenfrequencies of these modes do not display any significant scaling properties dependent on δz as shown in chapter 3. Therefore we know that M and M−1 will not have any interesting scaling properties dependent on δz. This is expected because M only contains information about the energies and not about the constraints of the system. It does know about z, for increasing z is adding extra angles to the system. However, it does not know anything about δz, because it doesn’t contain any information about the constraints of the system.

SST _{however contains the information of the constraints on the system. Therefore we expect that}

Figure 9: Spectrum of SST at increasing coordinations. 40 networks were used at all shown coordi-nations.

Here we clearly see that lower frequencies are allowed when we decrease δz, i.e. increase the coordination closer to the isostatic point, as shown by D¨uring et al. [5]. This is due to the effect that ST_{|ξi is a factor of balancing the forces on the bonds. Only when z = z}

c we find complete balancing

on the bonds and we find exactly one solution for ST_{|ξi = 0. If z < z}

c we have too few constraints

and can’t find a finite solution for ST_{|ξi = 0. The closer we get to z = z}

c, the more constraints are

added to the system and we can optimize ST_{|ξi = 0 better and better. We find that the optimization}

of this solution scales with δz. For the vector SM−1

   ∂2_U ∂ ~R∂γ E −  ∂r ∂γ E

we now know that M−1 physically does nothing interesting with increasing coordination. S has some modes which go to zero as shown above, but because ∂2U

∂ ~R∂γ

is only dependent on the energies of the system, we would not expect it to couple only to those modes. So we can assume that SM−1

 ∂2_U ∂ ~R∂γ E −  ∂r ∂γ E

behaves like a random vector coupling on all the modes. Let us for simplicity call it |Y i:

SM−1     ∂2_U ∂ ~R∂γ  −     ∂r ∂γ  = |Y i

Now we get the following equation for | ˙Ri:

| ˙Ri = −M−1  ST SM−1ST
−1|Y i +     ∂2U ∂ ~R∂γ  (32)

And for h ˙R| ˙Ri:

h ˙R| ˙Ri =  _∂2_U ∂ ~R∂γ     M−1+ hY | SM−1ST
−1 SM−1   M−1ST SM−1ST
−1 |Y i +     ∂2_U ∂ ~R∂γ 

If we now look at the separate terms and take into account that (SST₎−1 _{has modes that go to}

infinity and M−1 does not, we get:

 _∂2_U ∂ ~R∂γ     M−2     ∂2_U ∂ ~R∂γ   _∂2_U ∂ ~R∂γ     (ST SST
−1     Y 

 Y     SST
−1 S     ∂2_U ∂ ~R∂γ  hY | SST
−1 SST SST
−1 |Y i

Here the last term contains most (SST₎−1 _{and will therefore go fastest to infinity, taking into}

account that |Y i is a random vector and thus couples as much to the low energy states as to the high energy states. We get the following scaling for h ˙R| ˙Ri:

h ˙R| ˙Ri ∼ hY | SST
−1

|Y i (33)

Decomposing in the system of SST _{we get the following:}

hY | SST
−1 |Y i =X l hψl|Y i 2 ω2 ∼ Z inf δz D(ω) ω2 dω ∼ 1 δz (34)

5.1

Numerical Analysis

Figure 10: Shear modulus of networks at increasing connectivity. All data points are averaged over five networks.

Here we clearly see the predicted pattern, that the shear modulus goes like 1

δz. As we expected,

the energies show a singular behaviour dependent on the distance from the critical point. If the coordination gets closer to the critical point the shear modulus increases very quickly. When we get to the critical point the nodes cannot rearrange anymore and therefore the shear modulus goes to infinity. The fast increase might also partly be a consequence of the delocalization of the modes. The delocalization of the modes increases with increased coordination as shown in chapter three. A more coordinated movement is necessary to re-organize the nodes causing an increased movement of the nodes and thus increasing the shear modulus.

6

Deriving equations of motion

We found the vibrational modes of the system, now it is interesting to look at how the system evolves over time if we would put it in a thermal bath. In a thermal bath the temperature of the system will rise and this gives rise to vibrational movements. To look at these vibrations we need the equation of motion. To derive an equation of motion we first need to look at the constraints in the dynamics of the system. We know that the particles can only move in directions which do not break the constraints. These constraints are in the S operator as shown before. The constraint equation that we find is as follows:

S| ˙Ri = 0 (35)

With | ˙Ri defined as follows:

|Ri = Ri_{− R}i 0

| ˙Ri = d dt|Ri

If we derive our constraint equation with respect to time we get the following equation:

∂S

∂R| ˙Ri + S| ¨Ri = 0 | ¨Ri = −S−1h ˙R|∂S

∂R| ˙Ri (36)

Now using Newton’s second law we can derive a differential equation from which it is possible to derive the equation of motion. For simplicity we assume all masses equal to one:

Fext+ STf = | ¨Ri (37)

Fext≈ −

∂U

∂R− M |Ri (38)

Using equations (36), (37) and (38) we can derive the following equation for f:

f = − SST
−1 h ˙R|∂S

∂R| ˙Ri + SS

T
−1

SM |Ri (39)

Using equation (39) and (37) we get the following equation of motion:

−P M P |Ri − ST _SST
−1 h ˙R|∂S

∂R| ˙Ri = | ¨Ri (40)

Interesting is the rise of term dependent on the velocity ST _SST
−1 h ˙R|∂S

∂R| ˙Ri in the equation of

motion. This term came into existence by applying the constraints of the system on the velocity. It will cause vibrations to spread over multiple vibrational modes. A displacment of the system in the direction of a vibrational mode will not make the system stay in one vibration. The term dependent on velocity will cause energy to spread itself over different vibrational modes and eventually will cause the system to vibrate in all vibrational modes.

7

Discussion

While doing this research we found a few problems which we will discuss in this chapter.

7.1

Network size

In this research of the behaviour of backbone networks we used networks of size N = 1024. We found that the uncertainty of the shear modulus at the higher coordinations became really high. This is because of the small network size we used. In a coordination of z = 3.99 there are just 10 degrees of freedom left. Here we clearly lose the strength of statistics. At this coordination about five bonds are being removed from the isostatic point. With so few bonds being removed we can’t assume that there is no average rotation of the bonds. There is net rotation of the bonds. Hence there will be a preferred direction of movement and this will influence the calculation of the shear modulus. For better analysis of the shear modulus on coordinations this close to the isosatic point we need to look at larger systems.

7.2

Breaking constraints at second order

When we calculated the equation of motion, describing the acceleration, we only used the first order approximation for our constraints. However, in the equation describing the acceleration we find a term dependent on the velocity and the second order approximation of our constraints.

−P M P |Ri − ST _SST
−1 h ˙R|∂S

∂R| ˙Ri = | ¨Ri (41)

The second order term here is: ST SST
−1h ˙R|∂S_∂R| ˙Ri. This equation tells us in which direction the particles will accelerate. Of course the particles will only be able to move in the constrained displacement, in other words, in the direction of the displacement in the zero space of S. Therefore the acceleration must also be in the zero space of S, to make sure the particles will only move in allowed directions. However, if we multiply the equation of motion by S we get this:

−SP M P |Ri − SST SST
−1 h ˙R|∂S

∂R| ˙Ri = S| ¨Ri −h ˙R|∂S

∂R| ˙Ri = S| ¨Ri (42)

h ˙R|∂S_∂t| ˙Ri clearly does not have to be zero, so the acceleration is not in the zero space of S and changes the lengths of the bonds. This is a violation of the hard constraints. At second order in our approximation we find a breaking of our hard constraints and we can’t evaluate the movements of the particles. For a complete solution there must be no breaking of the hard constraints in any order of the approximation.

A possible solution may be the use of an alternative coordinate system. This must be a system in which the lengths of the bonds are not dependent on the variables in the coordinate system. One such coordinate system could be one in which we describe the position and orientation of the bonds instead of the positions of the particles. To describe the system like this we would need two variables: ~

λα and θα. Here ~λα describes the position of the centre of a bond and θα the rotation of the bond

compared to the x-axis. In the vector θα is written as follows: |θi =

cos(θα)

sin(θα)



Now we write the coordinates in a vector as follows:

|λ, θi =|λi |θi 

An illustration of this coordinate system would be like this:

Figure 11: Illustration of coordinate system using λ and θ.

In figure 11, i and j are two particles that are connected by one bond. The vector for the location of the bond is λ and the orientation is given by θ. To get from the coordinate system of positions of particles to the one above, the following matrix is used:

Vnµ= 1 nν (δni+ δnj) + 1 2nu (δnj− δni)ru (44)

Here ν goes from 1 to half max(µ) and u goes form half max(µ) to max(µ). µ loops twice through all the bonds and is defined as: µ = ij where µ is de bond between particle i and j. Now the following is true:

|Ri = Vnµ|λ, θi (45)

In this coordinate system the lengths of the bonds do not depend on the variables. Hence the lengths of the bonds are conserved with changing of the variables. In this way the hard constraints are never broken.

7.3

Scaling in the energies of the system

In chapter four we made a physical argument why M , the Hessian, should not influence the singular behaviour of the shear modulus. It contains only the information about the energy in the system. However, we also calculated M numerically and it turned out that M actually does show a singularity with respect to δz. By increasing the coordination we find that modes of M start to approach zero:

Figure 12: Spectrum of the Hessian with increasing coordination. 40 networks were generated and analyzed for every coordination.

These modes that are going to zero should not be a problem, because these are the modes that break the hard constraints. However, if we combine these modes with

   ∂2U ∂ ~R∂γ E , as with SM−1    ∂2U ∂ ~R∂γ E ,

we see a singularity with respect to δz. We call this combined vector F = SM−1 

∂2_U ∂ ~R∂γ

E :

Figure 13: Fit describing the dependence of vector F on δz. At every coordination 50 networks were analyzed and the averaged result is shown.

The result we found, of the shear modulus depending on δz with G = _δz1, should still hold. We also found it numerically. A physical argument that this extra singularity in M does not change the dependence of the shear modulus on δz is as follows: If we would replace the bending interaction

with weak springs we would find the same behaviour of the network. However in that case M does not depend on δz anymore. It would not contain any information about the critical point.

Eventually the singularities in M we find here should cancel out, but this is something that still should be further investigated.

8

Acknowledgements

I want to express my utmost appreciation of Edan Lerner, my supervisor, for assisting me during this project. Every meeting we had he spoke with much enthusiasm about the project. Every meeting concluded in new interesting ideas of how to approach the problems at hand. He gave me a great insight into his research field.

I also want to thank Robbie Rens, who is doing a PhD with Edan, for his immense support during my project. He helped me a great deal in understanding parts of the physics behind the project. He also helped me a great deal by sharing his program for analysis and helping me whenever I got stuck. It was great working alongside him.

References

[1] J. Clerk Maxwell F.R.S. L. on the calculation of the equilibrium and stiffness of frames. Philo-sophical Magazine Series 4, 27(182):294–299, 1864.

[2] M Sheinman, CP Broedersz, and FC MacKintosh. Nonlinear effective-medium theory of disordered spring networks. Physical Review E, 85(2):021801, 2012.

[3] Abhinav Sharma, AJ Licup, R Rens, M Vahabi, KA Jansen, GH Koenderink, and FC MacKintosh. Strain-driven criticality underlies nonlinear mechanics of fibrous networks. Physical Review E, 94(4):042407, 2016.

[4] Hui Li and Jan H Jensen. Partial hessian vibrational analysis: the localization of the molecular vibrational energy and entropy. Theoretical Chemistry Accounts: Theory, Computation, and Modeling (Theoretica Chimica Acta), 107(4):211–219, 2002.

[5] Gustavo D¨uring, Edan Lerner, and Matthieu Wyart. Phonon gap and localization lengths in floppy materials. Soft Matter, 9(1):146–154, 2013.

[6] Stefan B Lindstr¨om, David A Vader, Artem Kulachenko, and David A Weitz. Biopolymer network geometries: characterization, regeneration, and elastic properties. Physical Review E, 82(5):051905, 2010.

[7] Chase P Broedersz, Xiaoming Mao, Tom C Lubensky, and Frederick C MacKintosh. Criticality and isostaticity in fibre networks. Nature Physics, 7(12):983–988, 2011.