Colloidal Simulations of Equilibrium Gels and Analysis of Their Structure and Properties

Colloidal Simulations of Equilibrium Gels and

Analysis of Their Structure and Properties

Jurriaan van den Berg

August 14, 2020

Abstract

In this thesis the temperature dependence of the structural and the elastic properties of equilibrium gels are studied using theory and simulations. Gels, which are semi-solids consisting of a three dimensional network of clustered particles inside a medium liquid, are studied using N-body simulations of colloidal systems, where the underly-ing particle-particle interactions are modelled by a potential consistunderly-ing of short range attraction and longer range screened electrostatic repulsion, as introduced by [1]. This work confirms a strong temperature dependence of the formation of clusters of par-ticles, where lower temperatures allow larger clusters to form. At the lowest studied temperature T = 0.09 (system units) almost all the particles combine into one net-work spanning the entire system, while at the highest studied temperature T = 0.25 only small clusters consisting of a few particles are present. Furthermore, the elastic properties, here quantified via the shear modulus and the bulk modulus, of the solid gels were measured. Qualitatively the shear modulus shows the same temperature dependence as the maximum cluster size, rising for lower production temperatures as the system gets more clustered. The bulk modulus rises for both low and high temperatures, showing a minimum at the temperature T = 0.12, where the system transitions from a state with many small clusters to one containing a single large clus-ter. Finally, stress correlation functions of the gels have been studied, indicating short time temperature dependant oscillations.

This thesis was completed as a bachelor project for the bachelor Physics and As-tronomy at the University of Amsterdam (size 15 EC), and was conducted between March 30, 2020 and August 14, 2020.

Student number 11880082

Supervisor Dr. Edan Lerner

Second examiner Dr. Sara Jabbari Farouji Daily supervisor Corrado Rainone

Period March 30, 2020 to August 14, 2020. Study Bachelor Physics & Astronomy Institute Institute for Theoretical Physics

Contents

Preface 4

Acknowledgements . . . 4

Populair-wetenschappelijk samenvatting (in Dutch) . . . 4

1 Introduction 5 2 Preliminaries 6 2.1 Overview . . . 6

2.2 Liquids and solids . . . 6

2.2.1 Shear modulus . . . 7

2.2.2 Bulk modulus . . . 8

2.2.3 Viscosity . . . 9

2.3 Temperature . . . 9

2.4 Going supercooled . . . 10

2.4.1 Equilibrium and stability . . . 10

2.4.2 Glasses . . . 11 2.4.3 This research . . . 11 2.5 Colloidal suspensions . . . 12 2.5.1 Colloids . . . 12 2.5.2 Particle interactions . . . 12 2.5.3 Gels . . . 13 3 Simulating Gels 14 3.1 Particles . . . 14 3.2 Potentials . . . 14

3.2.2 The Yukawa potential . . . 15 3.2.3 This research . . . 16 3.2.4 Forces . . . 17 3.3 Integration . . . 18 3.3.1 Calculating forces . . . 18 3.3.2 Verlet’s method . . . 18 3.4 Regulating temperature . . . 19 3.5 Experiments . . . 20 3.6 Making solids . . . 20 3.6.1 Energy minimisation . . . 21 4 Results 22 4.1 Energy Evolution . . . 22 4.2 Structure . . . 23 4.2.1 Clustering . . . 23

4.2.2 Radial distribution function . . . 24

4.2.3 Structure factor . . . 24

4.3 Correlations . . . 27

4.4 Elastic properties . . . 28

5 Discussion 29

Preface

Acknowledgements

I would like to thank my daily supervisor Corrado Rainone, without his help and guidance this project would not have been possible, and I would like to thank my supervisor dr. Edan Lerner, giving me the opportunity to work with code from his group and granting me access to the SURFsara LISA computing cluster. Throughout this project my friends and family have supported me immensely, for which I would like to express my deepest appreciation. Moreover, I would like to Annika, Steven and Maaike in particular for going out of their way to proofread this thesis, you have been a tremendous help!

Populair-wetenschappelijk samenvatting (in Dutch)

Dit onderzoek houdt zich bezig met gels, puddingachtige substansies die vandaag de dag overal om ons heen te vinden zijn. Natuurlijk heb je cosmetica zoals haargel, maar ook veel voedselproducten, zoals pudding, en zelfs ook kaas, zijn gels. Daarnaast worden gels ook veel voor medicijnen in de farmacie gebruikt en hebben ze veel toepassingen in de elec-tronica, zo worden ze bijvoorbeeld gebruikt voor het maken van efficiente batterijen. Een gel is vaste stof die is opgebouwd uit een driedimensionaal netwerk van samengeclusterde deeltjes en een vloeistof waarin dit netwerk van deeltjes gesuspendeerd is. Dit netwerk geeft gels de elastische eigenschappen van een vaste stof. De term ’gel’ wordt ook vaak ge-bruikt voor vloeistofachtige substansies, maar wetenschappelijk gezien zijn gels alleen vaste stoffen, die niet wegstromen na een lange tijd. Gels hebben dus veel toepassingen, daarom is het erg belangrijk dat we precies weten wat gels zijn en hoe ze worden gemaakt, om zo ook verbeteringen te kunnen maken in zowel hun eigenschappen als de productieprocessen waarmee ze gemaakt worden. In deze scriptie wordt precies hiernaar gekeken. Zo wordt er uitgelegd wat gels zijn en hoe hun eigenschappen, zoals de stevigheid, kunnen worden beschreven en gemeten. Dit wordt gedaan met behulp van computersimulaties, waarmee de gels worden nagebootst. Door gesuspendeerde deeltjes te simuleren met een specifieke onderlinge interactie zullen deze deeltjes samenclusteren en zo een gel vormen. In welke mate dit gebeurt, hangt af van de temperatuur waarop deze deeltjes worden gesimuleerd. Er is gevonden dat voor lage temperaturen, waarbij de deeltjes dus minder bewegen, er meer clusters ontstaan dan bij hoge temperaturen, waarbij de deeltjes relatief snel bewe-gen. Precies voor de lage temperaturen zijn de gels ook veel steviger, wat dus aangeeft dat de stevigheid van gels komt door het samenclusteren van de deeltjes.

1.

Introduction

Gels have seen a lot of uses in the past decades, from pharmaceutics to cosmetics to food products, and as such they have been widely studied in research. Studying their properties allows for better gels to be made, with better and better tunable properties. Gels are solids that are essentially liquids with a rigid network of microscopic particles inside, giving the gels their solidity. Because of this, they can be sturdy and flexible, all the while still having relatively light liquid mass. However, making regular gels is an unstable process, which involves phase separating a suspension into the gel and some leftover waste. Moreover, it is hard to tune the properties of the gel using this process, and the gels that are created during such a process are generally not stable. This means that such gels settle over time, losing their gel properties in the process. This is where equilibrium gels come in. Equilibrium gels are a very useful kind of gel which does not phase separate over time, and as such remains a stable gel for forever. Moreover, the production process of equilibrium gels allows for better control of its properties, as the gels are created from a liquid suspension that changes continuously to a gel state, instead of relying on phase separation mechanisms. The uses of equilibrium gels include applications in super light batteries, allowing to store higher amounts of power with less weight than other batteries, and in medicine coating, which can be used to efficiently deliver drugs into the body. [2]

This thesis will study these equilibrium gels, using computational simulations to both simulate equilibrium gels and analyse them. To achieve this goal this thesis consists of two parts. The first part serves as a general introduction to the field of soft matter physics and gels, in which some of the basic concepts used in this field will be introduced. These concepts will then be used in the second part, where a study of equilibrium gels will be presented. Here the simulation and experimental procedures will be explained, the results of which will be given and discussed. The experiments will consist of three dimensional N-body simulations of supercooled colloidal systems at different temperatures, with an underlying particle-particle potential consisting of short range attraction and longer range screened electrostatic repulsion, following the simulations done in [1]. The effects of using this potential will be described and analysed, and the gels created in these experiments will be studied to find their structural and elastic properties related to the temperature they have been prepared at.

This, however, raises the question: why are these systems studied using simulations, and not using experiments on real equilibrium gels? These equilibrium gels exist in real life of course, and as such can be studied there. However, simulations allow for much more precise control over the system, and make it possible to do things that would be impossible to do in a laboratory. For example, regulating temperature precisely, homogenously and perhaps instantaneously over a whole system is often very difficult in real experiments, but can easily be done using a computer. Moreover, using simulations, every aspect of the system is fully known, and measurements are easily done. For these reasons, simulations are an important tool in physics, and as such are a logical device to use to study equilibrium gels.

2.

Preliminaries

2.1 Overview

Before it is possible to properly look into the effects of this potential a short introduction to the basic concepts used to describe liquids, glasses and gels is in order. In this sec-tion relevant terms and definisec-tions are given with short explanasec-tions. Furthermore, more detailed explanations are presented for the more essential concepts used in this thesis. Gels are semi-solids that can be described as liquids with a cross-linked network within, giving the gel its ’hard’, solid structure. A more precise definition of gels will be given later, but for now it is important to understand that gels and liquids are closely related. In fact, the simulations used in this thesis start out in a liquid state, and it is only due to progressive clustering of particles during the time evolution of the system that gels are formed. Moreover, many of the properties used to describe liquids can also be used to describe the properties of gels, which is why those will be introduced here.

2.2 Liquids and solids

First it is important to ask: what is the difference between solids and liquids? For starters, solid do not flow, while liquids do. However, as is almost always the case, in reality this line is more blurred. What does it mean for a liquid to flow, and how slow must a substance ’flow’ to be considered not flowing at all? Of course, there is matter which can be only be described as a solid, like solids with a specific periodic lattice. Leave these undisturbed for a long time, and they will still remain in the same shape. However there are also solids which do not exhibit a strict lattice, but are considered solids nonetheless. This has everything do to with the relaxation time τR of the substance, an important concept in soft matter physics. The relaxation time identifies the time it takes for a perturbed system to return, or relax, back to its equilibrium state. This concept is linearly related to the viscosity η of a liquid, which measures the resistance to deformation. A liquid with a high viscosity flows slowly, resulting in a high relaxation time after any disturbance. Even more so, taking η → ∞ would result in τR → ∞; the liquid does not flow, and as such could better be described as a solid. Caution has to be taken however, as here another concept introduces itself: whether a substance is a liquid or a solid depends on how long the system is observed for. Namely, if the system was to be observed over an experimental time of t → ∞, then at one point even an ’infinite viscosity solid’ would still flow. This may seem like a contradictory definition then, why not restrict solids to substances which do not flow on any timescale? But it is not! The experimental time-scale of the experiment is actually essential. How much a liquids flows depends on how long it is being observed. On time scales shorter than the relaxation time of liquid, the liquid does not flow, and it behaves like a solid. As such, describing it as liquid would make no practical sense. These concepts can be defined more rigorously by introducing the shear modulus, for reasons which will become apparent later.

2.2.1 Shear modulus

A F

y θ

∆x

Figure 1: A solid cube of height y to which a shear force F is applied to an area A, causing a shear displacement ∆x.

Understanding the qualities and implications of the concepts given above it is time to give more rigorous definitions of the viscosity and the re-laxation time. For that we first have to go back to solids, and work our way to liquids from there. Consider a solid cube, to which a shear force F is applied, resulting in a shear by an angle θ, see Figure 1. This shear gives a shear displacement ∆x at a height y of the cube, where these two values are related as

∆x = γy, (1)

where

γ = uxy (2)

is the shear strain in this case, which quantifies how much the cube is sheared. In general the shear strain is given by

uxy ≡ ∂∆x

∂y , (3)

so for this linear shear

uxy = γ = tan(θ) ≈ θ. (4)

The solid will resists this push, as it wants to stay in the same shape. This causes a responding force equal and opposite to the applied force F , which results in a shear stress σxy = F A = Guxy, (5) where G ≡ σxy uxy = F y A∆x (6)

is the shear modulus of the solid [3]. The shear modulus defines the elastic response of a material to an applied shear, and as such G can be seen as analogous to the spring constant in Hooke’s Law. If the shear is applied at a time t0 _{and the response is measured} at a time t, then it is possible to generalise this formula to the time-dependent shear modulus G(t − t0), given by

G(t − t0) ≡ σxy(t − t 0₎ uxy

. (7)

For a finite G(t) that does not decay in time the material tries to stay in place when a small shear force is applied, pushing itself back into its original shape. This is of course precisely a solid: when you push on a table the table does not move out of the way. Then, what happens for materials where G(t) = 0? It turns out that this actually describes liquids.

When such a substance is pushed there is little to no elastic response; the material just flows out of the way, and does not push back with an opposing force. In other words, for liquids G(t) → 0 as the time t → ∞. This then leads to the fundamental definition of the shear relaxation time τR, as defined in the Maxwell model for viscoelastic materials, which states that

G(t) = G∞e−t/τR, (8)

with G∞ the zero time shear modulus [4]. τR defines how long it takes for the shear modulus to decay with a factor of 1

e, and as such represents the time-scale on which a substance is considered a liquid or solid. It is therefore reasonable to say that for times t < τR the substance responds with finite shear stress to a shear, and as such behaves like a solid. On the other hand, for times t > τR the stress response of the substance has relaxed, and the substance flows as a liquid. In general real materials do not have singular relaxation times, and the shear modulus may be a summation of multiple elements from Equation 8, such that

G(t) =X i

Gi_∞e−t/τRi. (9)

This does not change the previous discussion however, as the largest relaxation time τi R still defines the relaxation time-scale of the substance, although more nuanced analysis may be done when more or all the relaxation times are considered [5].

2.2.2 Bulk modulus

Figure 2: A substance that is being compressed to a smaller size, illustrating the affect of the bulk modulus of the substance. Image taken from [6]. The shear modulus is one of several quantities

used to describe the stiffness or elasticity of ma-terials. The other quantities include the bulk modulus, Young’s modulus, and Poisson’s ratio [7]. In this section the first one of these will be shortly discussed.

The bulk modulus describes how much a sub-stance responds to compression. When a pres-sure P is applied to a normal substance, the substance will naturally shrink, as can be seen in Figure 2. How much the substance gets com-pressed compared to this pressure P is described by the bulk modulus K. The bulk modulus is more formally defined as the ratio between an

infinitesimal pressure increase and the relative decrease in volume the pressure increase causes. This means that the bulk modulus is given by

K ≡ −V dP

dV , (10)

where V is the volume of the solid, and P the pressure [8]. The bulk modulus is an important quantity in describing the behaviour of liquids and solids alike. However, it does not relate to the relaxation time fundamentals in the same way the shear modulus does.

2.2.3 Viscosity

It is useful to expand the previous definitions of shears to include the viscosity of a liquid. Consider a time-dependent shear strain uxy(t). Here each infinitesimal small strain duxy(t0) at a time t0 _{induces a infinitesimal small shear stress}

dσxy(t) = G(t − t0) duxy(t0) = G(t − t0) ˙uxy(t0) dt0, (11) where ˙uxy is the shear rate: the time derivative of the shear strain. Each of these small stresses add up over time, building up the total stress. Say that, without loss of generality, uxy(t0) = 0 for some time t0 < 0, then it is possible to write the total shear stress

σxy(t) = Z dσxy(t) = Z t 0 G(t − t0) ˙uxy(t0)dt0. (12) For a constant shear rate, or equivalently a linear shear strain, this gives

σxy(t) = ˙uxy Z t

0

G(t − t0)dt0, (13)

where, if the integral on the right converges, then in the limit t → ∞ we get σxy = lim t→∞u˙xy Z t 0 G(t − t0)dt0 = ˙uxy Z ∞ 0 G(s)ds = ˙uxyη, (14) with η ≡ Z ∞ 0 G(s)ds. (15)

Here η is the viscosity of the liquid: it is the ratio between the shear stress and the shear rate, in some respects analogous to the shear modulus. For normal liquids, such as the ones described by Equation 8, the time-dependent shear relaxation modulus decays faster than 1/s, and therefore the viscosity converges to a fixed value specific to the liquid. In this case η = Z ∞ 0 G∞e−s/τRds = G∞ τR , (16) so τR= G∞η, (17)

from which we recover the previously mentioned relation

τR∝ η. (18)

2.3 Temperature

The temperature of a substance is directly related to the movement of the particles making up the substance. For example, for a monoatomic ideal gas the temperature T and the kinetic energy Ek are related by

3

2kbT = Ek = 1 2mv

where kbis the Boltzmann constant. This is a specific form of the equipartition theorem, where in this case there are three degrees of freedom, one for the translation movement in each dimension. Liquids are often more complicated however. Other than that there are usually more degrees of freedom present other than the three translation ones, it is here the case that the particle-particle interactions cause the equipartition theorem to break down. Nevertheless it is still generally possible to use Equation 19 to relate the liquids temperature and its average translational kinetic energy, and, to an approximation, the relation between particle movement and temperature it describes remains true. [9]

In that spirit, from Equation 19 it can be seen that for lower temperatures it must be the case that also the average velocities of the particles in the liquids must be lower in magnitude. In other words, for lower temperatures the dynamics of the liquid slows down, and as such it takes longer for the liquid the return to equilibrium after any applied strain. From this then follows an important relation; as the temperature of a normal substance goes down, its relaxation time increases. In addition, the relaxation time the temperature also influences the relative effects of the particle-particle interactions. The higher the temperature of the system, the more kinetic energy the particles have, and how relatively less impact the underlying particle interactions influence the system. The particles, having higher energies, are more easily able to overcome the effects of potential energy barriers. This in turn, ignoring the thermodynamic details, gives rise to phase transitions. Below the melting temperature Tm a liquid crystallises into a solid, while above the boiling point Tb the liquid vaporises into a gas.

2.4 Going supercooled

2.4.1 Equilibrium and stability

Even if a liquid is cooled below its melting point Tm, that does not necessarily mean that the liquid will crystallise. This has everything to do with nucleation, which is the phenomenon that the solid phase does not come into existence homogenously, but instead grows from a finite amount of already crystallised points. The full details will not be described here, for more information see [10], but for now it is important to know that if this nucleation does not, or is not allowed to, happen, then the liquid will not crystallise. If this is the case then what is left is what is known as a supercooled liquid: a liquid below its melting temperature. In principle this liquid phase is not stable, as the solid crystalline phase is the thermodynamic minimum energy configuration, or the equilibrium configuration, of the system. However, it can still be said that the supercooled phase is in some kind of equilibrium, as it is possible to experimentally ’equilibrate’ a supercooled liquid such that time-translation invariance holds, meaning that the system is relaxed and that its observable quantities only fluctuate around a static point over time. If this is the case, then there is no way to know that the system is not in true thermodynamic equilibrium just by observing it. Because of this we expand the definition of equilibrium for our purposes to include this kind of equilibrium, where the solid is the stable and the supercooled liquid phase is the metastable equilibrium phase.

2.4.2 Glasses

Now, given that a supercooled liquid can be in equilibrium, what does it mean for one to not reach equilibrium? Say that the relaxation time of the substance is too large for the liquid to equilibrate in any experiment. In this case we are of course actually dealing with a solid. This is what is known as a glass, an off-equilibrium solid. Important to note is that because of this, glasses and liquids are structurally indistinguishable; it is the shear relaxation time, and as such also the presence of a shear modulus, that distinguishes the two states. In both glasses and liquids there is no long range order, unlike of course the equilibrium stable crystalline phase. The previous section described that lowering the temperature causes the relaxation time of the liquid to smoothly increase. While this is true, glassy states are often formed by a much sharper rise in the relaxation time, even increasing up to 14 order of magnitude when the temperature is only lowered a relatively tiny amount [11]. Such ’transitions’ are called glass transitions, and can be defined as

τR(T < Tg) > τexp, (20) where τRis the (shear) relaxation time of the liquid, τexpthe experimental time frame, and Tg the glass transition temperature. Glass transitions seem to be very general, appearing in liquids with a wide range of compositions and underlying particle interactions. However, it is important to note that glass transitions are not phase transitions, even though it may look like they are. Unlike regular phase transitions, glass transitions are purely dynamical in nature, and form a glass state that is not thermodynamically stable.

While supercooled liquids and glasses are closely connected - the deeper the degree of supercooling the easier it is to form a glass - Tg is not required to be below the melting temperature Tm, and as such glasses can form even before supercooling is reached [12]. The phenomenon of glasses, and also importantly the glass transition, is a whole research field on its own, which this thesis will not cover in further detail. For more information on this interesting topic, see [12], [13].

2.4.3 This research

For the purposes of this research it is sufficient to know that both the crystalline phase and the glass phase are to be avoided. Of course the crystalline phase is not wanted, as in this solid phase the characteristic clustering of gels cannot take place. Moreover, the glass phase has to be avoided, as the sharp increase in relaxation time would prevent clusters from being formed during our experimental simulation time, and the interesting physics of gels would also be lost. From this, it follows that too low temperatures cannot be reached during the experiments. However, at too high temperatures the effects of the potential will be drowned out by the kinetic energy of the particles, which means it is also not preferable to go to too high temperatures. It is therefore possible to hypothesise that the clustering takes place only in a specific temperature range.

2.5 Colloidal suspensions

2.5.1 Colloids

Figure 3: An illustration of a 2d colloidal suspension, consist-ing of particles suspended in a solvent.

Now that we know in which circumstances the clustering could happen, we have to ask ourselves how the cluster-ing actually takes place. To understand this we first need to look at what gels are, which will be done in this sec-tion. Gels are a specific form of colloidal suspensions, or colloids for short. Colloidal suspensions are mixtures consisting of particles suspended in a medium fluid, as can be seen illustrated in Figure 3. These phase sepa-rated particles, also sometimes called colloids by them-selves, can be solid or liquid in nature, and can be ap-proximately between 1 and 1000 nanometers in size [14]. Note that only the term colloidal suspensions unambigu-ously refers to the mixture of the particles and the sol-vent together, and colloid can refer to the particles or the whole mixture. This thesis will furthermore use colloid to only refer to the particles themselves.

It’s insightful to make the distinction between a solution,

a regular suspension, and a colloidal suspensions. In a solution the solute and the solvent combine together to form a singular phase, while in a (colloidal) suspension the solute are distinct dispersed particles, and are in a separate phase as the medium solvent they reside in. On the other end of the spectrum lie regular suspensions, where the solute contains particles of any size larger than 1 micrometer. Because of this the particles in suspensions generally settle after some time. This is not the case for colloidal suspensions however, as here the particles are small enough such that random collisions with the medium fluid and between the particles themselves prevent settling, or cause the particles to only settle after a long time [14].

2.5.2 Particle interactions

The movement of the particles suspended in the colloidal suspension is governed by the interactions they exchange with each other. As mentioned in the previous section, colloidal suspensions do not settle, and as such external forces such as gravity can be ignored. Due to their small size the colloids interact much like atoms would in a liquid. In fact, the kind of simulations used for this thesis are also used to model liquids, and to study the glass transition. As such, the same kind of forces play a role in both cases. These forces can be split up into the following four parts: excluded volume repulsion, electrostatic interactions, van der Waals forces, and entropic forces [15]. Volume repulsion simply states that it is not possible for the particles to overlap, and as such there are forces which prevent this from happening. Electrostatic interactions describe the forces caused by the electrical charge of the particles and of the solvent the particles reside in. This can cause the particles to attract or repel each other, where the charge of the solvent may influence (for example, screen) such interactions. The van der Waals forces are caused by the interactions arising from the polarization of the electrons of the particles. While

mostly intermolecular in nature, these effects are also present in colloidal suspensions, even if the colloids themselves do not carry a permanent dipole [16]. Finally there are the entropic forces, which are forces resulting from entropy maximization, following the second law of thermodynamics. These forces are present when shaped colloids or polymer suspensions are considered, but arise even if the colloids are simply spherical in shape, see [17] for more information.

2.5.3 Gels

Figure 4: An illustration of a 2d gel, consisting of clustered particles suspended in a solvent. Then we are left with the question which started this

sec-tion: what is a gel? As mentioned before, a gel is a col-loidal suspension, but with a twist. In a gel the interac-tions described in the previous section combine to cause the particles to cluster together, forming a cross-linked network structure inside the solvent. A 2d illustration of this is given in Figure 4. This three dimensional network causes gels to be stiff and behave like solids, while still having mostly the weight of the liquid medium. Impor-tant to note is the difference between glasses: the solid like behaviour of gels is caused by the network inside it, while the solid like behaviour of glasses has no apparent structural origin. More formally a gel is defined as a nonfluid colloidal network or polymer network that is ex-panded throughout its whole volume by a fluid. [18] A subset of gels are equilibrium gels, which are gels

which do not settle, or phase separate, over time. This means that in equilibrium gels the structured network is consistent, and the gel is considered stable.

For the purposes of this thesis only colloidal systems will be considered. Here simple spherical, monoatomic particles are used, which means that it is possible to use Equation 19 to find the temperature of the system. Furthermore, the interactions of the colloids in real gels will be modelled using a single particle-particle potential, which allows the clustering to take place.

3.

Simulating Gels

3.1 Particles

Figure 5: A simplified illustra-tion showing a colloidal system containing 16 particles.

The colloidal gels are modelled using a 3D N-body sim-ulation, where each dispensed particle can move around freely. These particles interact which each other via a direct particle-particle interaction, using a radial poten-tial which models the interactions present in real gels. In the simulations studied in this thesis the indirect effects of the solvent are assumed to be negligible and exter-nal forces such as gravity are also ignored. Moreover, to model a boundless system, and to ignore surface in-teractions, periodic boundary conditions are used. The system then is a cube shaped box of volume V = L3_, in which N particles are contained. These two quanti-ties are related by the number density ρ of the system, given by

ρ = V

N. (21)

The particles are modelled as simple spherical particles of diameter σ, which implies that the volume density φ is given by

φ = ρπσ 3

6 , (22)

a value commonly used instead of the number density. Each of the particles has a position r = (x, y, z) and a momentum p = (px, py, pz), where the periodic boundary conditions state that

(x, y, z) = (x + L, y, z) = (x, y + L, z) = (x, y, z + L). (23) For simplicity we assume m = 1, such that the velocity v = p and that the particle acceleration is equal to the force on the particles. The whole system is then described by two 3N vectors: one describing the positions of the particles and one their momenta.

3.2 Potentials

To empirically model all the interactions present in gels, a radial potential consisting of two parts is used. One part is the short range Lennard-Jones potential, and the other part is a repulsive, long range Yukawa potential. These potentials will be introduced in this section.

3.2.1 The Lennard-Jones potential

The Lennard-Jones (LJ) potential consist of two parts: a short range strong repulsive part, which represents volume exclusion, combined with a more long range attractive part, which

0.8 1.0 1.2 1.4 1.6 1.8 2.0 r/σ -1.0 0.0 1.0 2.0 ULJ ( ² ) α = 6 α = 12 α = 18

Figure 6: The Lennard-Jones potential, given by Equation 24, for different values of the interaction strength α.

0.0 1.0 2.0 3.0 4.0 r/ξ 0.0 1.0 2.0 3.0 UY ( A ) k = 0 k = 1 k = 2

Figure 7: A repulsive version of the Yukawa potential, given by Equation 25, for different values of k.

models electrostatic attraction. In general the potential is given by ULJ = 4  σ r 2α −σ r α , (24)

where  is the base energy strength of the potential and σ is the distance at which the LJ potential is zero, which corresponds to 2 times the radius of the interacting particles. α describes the range, or strength, of the interaction, for which commonly α = 6 is used. A plot of the LJ potential for different values of α is given in Figure 6. Of note is the strong repulsion for r < σ, caused by the first term in Equation 24. For large r the attractive longer range second term becomes dominant instead. In between the two terms combine to form a potential energy well, which causes particles to stay, or ’bond’, together.

The Lennard-Jones potential is commonly used in molecular simulations, and as such has been researched intensively. Its advantages include its computational simplicity, while still being able to be closely fitted to many real world systems. Besides this it is rather hard for systems using the Lennard-Jones potential to crystallise, which is very useful when studying the glass transition where the crystal phase is to be avoided.

3.2.2 The Yukawa potential

The Yukawa potential, which is also known as the screened Coulomb potential, describes exactly that: screened electrostatic interactions. Screening is the effect that electric fields are damped due to the presence of mobile charge carriers, which arrange themselves in such a way for this to take place. The Yukawa potential is given by

UY = A e−kr/ξ

r/ξ , (25)

where A is the interaction strength of the potential and k and ξ are scaling factors. The range of the potential is then given by roughly r = ξ

0.0

1.0

2.0

3.0

4.0

r/σ

-0.5

0.0

0.5

U

(

²

)

Figure 8: The particle-particle potential used in this paper, consisting of a short range Lennard-Jones potential and a long range repulsive Yukawa potential.

just the regular Coulomb potential multiplied by an extra exponential damping term. A plot of this potential is given in Figure 7 for different k. Here A is positive, which gives a repulsive force in this case. Apart from condensed matter physics, the Yukawa potential is often used in particle physics, although often in a slightly different, but equivalent, form [19].

3.2.3 This research

As mentioned before, the potential describing the particle-particle interaction used in the simulations consists of two parts: a short-range Lennard-Jones potential ULJ (Equation 24) and a repulsive, long range, Yukawa-potential UY (Equation 25), where

ULJ= 4 _ σ r 36 −σ r 18 , (26) UY = 4 e−2r/σ r/σ . (27)

Here again r is the distance between the two particles,  is the base energy strength of the potential and σ is the distance at which the LJ potential is zero, which corresponds to the diameter of the particles. This thesis focuses on the case where α = 18 for the Lennard-Jones potential, and A = 8, ξ = σ/2 for the Yukawa part. For simplicity energy units of  will be used, so  will be set to 1 in the following equations. Moreover, the actual simulations will use monodisperse particles, and as such σ = 1, but the system could be expanded to use differently sized particle as well. This gives the total particle-particle potential U (r) = 4 " σ r 36 −σ r 18 + 4e −2r/σ r/σ # , (28)

for which its form is shown Figure 8. The notable difference between the regular LJ potential is the bump after the LJ potential well, which will cause particles to repel each other when their distance r > 1.28σ. Theoretically this change would cause the particles to cluster together and form bonds consisting of a few particles, as there is only space for

Figure 9: A one dimensional Bernal spiral structure of particles. This configuration is the ground state energy of the potential given in Equation 28. Image taken from [1]. a few particles in the potential well. These clusters then repel particles not in the cluster, which makes sure the system does not crystallise. Due to these effects the particles form one dimensional spirals at low enough temperatures, shown in Figure 9. These Bernal spiralsform the groundstate configuration of the potential [1].

3.2.4 Forces

To calculate the force on a particle j due to a particle i we have to calculate the gradient of the potential from Equation 28 with respect to the distance r = rij from particle i to particle j. We define rij ≡ rj− ri =   xj− xi yj− yi zj− zi  , (29) rij ≡ |rij| = q (xj− xi)2+ (yj − yi)2+ (zj − zi)2 = r. (30) The force on a particle j due to a particle i is then given by

Fij = −∇ijU (r) = −∇ijULJ− ∇ijUY, (31) where −∇_ijUL= −4 −36σ36 r_ij37 ˆrij+ 4 −18σ18 r_ij19 ˆrij = " 144σ36 r_ij38 − 72σ18 r_ij20 # rij and, using the quotient rule,

−∇ijUY = −4 −_σ2rij σ e −2r_ij/σ₋ 1 σe −2r_ij/σ (rij/σ)2 ˆ rij = " 8e−2rij/σ r_ij2 + 4σe−2rij/σ r_ij3 # rij. This gives the total force on each particle i for an N particle system

Fi = N X j6=i (−Fij) = − N X j6=i " 144σ36 r_ij38 − 72σ18 r_ij20 + 8e−2rij/σ r_ij2 + 4σe−2rij/σ r_ij3 # rij. (32)

3.3 Integration

3.3.1 Calculating forces

Using the analytical expressions for the force it is possible to simulate the Newtonian dy-namics of the colloids. Starting from some initial configuration, the system gets simulated, or integrated, for some time trun, using time steps of length ∆t. Here each time step the forces on all the particles are calculated using Equation 32, and integrated to find the particle velocities. The velocities are then used to update the positions of the particles. However, it would be quite computationally expensive to calculate the forces between all combinations of two particles, especially since the potential becomes negligible for large r. For this reason the potential is cut off a distance rc, reducing the need to calculate the value of the forces when they would be close to zero. Moreover, using this, each particle only has to check the neighbours that fall in the range of rc, resulting in a lot less itera-tions that have to be done each time step. However, it is not possible to simply say that U (r) = 0 when r > rc, as this would break the continuity of the potential, meaning that the force and the higher order derivatives would be indefinite at r = rc. This problem is solved by adding a correction potential Ucorr to Equation 28, which forces the potential and its derivatives to be 0 at r = rc. In this case

Ucor(r) = C0+ C2r2+ C4r4+ C6r6, (33) which gives Usim(r) = U (r) + Ucorr(r) = 4 " σ r 36 −σ r 18 + 4e −2r/σ r/σ # + C0+ C2r2+ C4r4+ C6r6. (34) Here C0,2,4,6 are obtained by solving a system of linear equations consisting of the new Usim(r) and its first three derivatives, forcing these four equations to be zero at r = rc. Adding more correction terms would allow for a greater number of derivatives to be put to zero at rcutoff, but for the purposes of this thesis the first three are enough. For completeness, this gives the final force on a particle i:

Fi = − N X j " 144σ36 r 38 ij −72σ 18 r 20 ij +8e −2rij/σ r 2 ij +4σe −2r_ij/σ r 3 ij − 2C2− 4C4rij2− 6C6rij4 # rij. (35) 3.3.2 Verlet’s method

The integration technique used in this thesis is called the velocity Verlet algorithm, a second-order integration method designed to integrate Netwonian equations of motion [20]. Velocity Verlet integration is often used in molecular simulations as it provides important physical properties, like time reversal symmetry and preservation of energy, all the while it has relatively low computing costs, only having to compute the forces on the particles once per time step. This means that it is approximately as fast as the simple Euler method, while also providing better numerical stability. The velocity Verlet algorithms differs from

regular Verlet integration in that is calculates the velocities of the particles explicitly, which is not done in the Verlet algorithm [21]. The velocity Verlet algorithm is used so that it is possible to easily calculate the kinetic energy of the particles. Moreover, Verlet integration is used in favour of other, higher order integration techniques for its efficiency, while still providing accurate, physical integration.

Velocity Verlet integration is able to integrate second-order differential equations of the form

¨

x(t) = ˙p(t) = a(t) = F(x(t)) (36) where x, p, a, and F are the positions, the momenta, the accelerations and the forces on the particles respectively. Note that here for simplicity m = 1. Starting at some starting conditions x(0) and p(0), with a(0) = F(x(0)), time is advanced one step by performing the following algorithm:

1. Calculate p(t +1 2∆t) = p(t) + 1 2a(t)∆t. 2. Calculate x(t + ∆t) = x(t) + p(t +1 2∆t)∆t.

3. Calculate a(t + ∆t) = F(x(t + ∆t)) using Equation 35. 4. Calculate p(t + ∆t) = p(t +1

2∆t) + 1

2a(t + ∆t)∆t.

This algorithm is then repeated for as many time steps as required.

3.4 Regulating temperature

During the studying of molecular dynamics it is important that the temperature may be kept to a desired, adjustable value. In real experiments this is commonly done with a heat bath kept at a constant temperature, with which the system that is being studied can freely exchange energy, such that the temperature of the system becomes equal to that of the bath. Such systems that modulate the temperature of the system are called thermostats, and are also important in molecular simulations. In principle it is possible to simulate a whole gel and bath system, but this is unwieldy and costly. For this reason various algorithms have been devised to efficiently mimic such thermostats, each with their ups and downs [22]. The algorithm that will be used in this thesis is the Berendsen thermostat, which has the advantage that it allows for global temperature fluctuations, which more closely represents real world systems [23]. Here the (instantaneous) temperature Ti, given by Equation 19, is regulated by scaling the velocities of the particles at each time step with with a factor ξT = s 1 + T Ti − 1 ∆t τT , (37)

given the desired temperature T and some coupling parameter τT which regulates the strength, or time scale, of the thermostat. This then ensures that

dTi dt =

T − Ti

τ , (38)

which can be solved to give

creating exponential decay from the temperature Ti to the desired temperature T . τT can be freely chosen to suit the needs of the experiments.

3.5 Experiments

Using the techniques described in the previous sections it is possible to simulate colloidal gels at various temperatures. Simulation will be done using N = 4096 colloidal particles in a box with number density ρ = 0.3, using a time step of ∆t = 0.005. Six temperatures will be studied, namely T = 0.09, 0.1, 0.12, 0.15, 0.2, 0.25, which will be regulated using a Berendsen thermostat with τT = 1.0. These temperatures are chosen to be in the range which allows clustering to take place, while going as low as is possible to simulate in our time frame. The experiments will proceed in three steps: generation, equilibration, and production, which are, for statistical reasons, repeated for a total of 96 runs.

In the generation step the particles are first initialised in a cubic grid formation and are given random velocities. The system is then run for some time t at high temperature, resulting in a configuration with randomised, but physical, particle positions. From this a ’snapshot’, containing the particle positions, is saved and used in the next step: equili-bration. In this step the system is cooled down to one of the configuration temperatures T, and is run until equilibrium is reached. This relaxation is completed if the total energy of the system does not decay any more over time, and only fluctuates around an average value. From this the relaxation time τR can be estimated per temperature. The configu-ration after equilibconfigu-ration is then stored and used as the starting point for the production step. In this step the system is run at the same corresponding temperature for t = 100τR, where after each relaxation time a snapshot is saved, or ’produced’. Although the system is relaxed during the production run, the particles still move around, causing clusters to break and form. Simulations are run for a maximum of 5 continuous days. During simu-lation certain quantities are saved every few time steps, which include the stress and the pressure on the system, and the average potential and kinetic energies per particle. See Figure 10 for a schematic overview of these steps.

3.6 Making solids

The systems simulated in the experiments are still liquid in nature, which is indicated by their relaxation times being smaller than the experimental time scale. To create actual solid gels the snapshots created in the production run have to be quenched. During a quench all the thermal energy is instantaneously removed, and potential energy of the system is minimised. Physically this means that the movement of the colloidal particles becomes negligible, and that the system is put into its (local) ground state energy. As there is no particle movement and the forces on the particles all balance each other in the potential energy minimum, the relaxation time becomes infinite, describing a solid. Given these solid gels it is possible to find the shear and bulk moduli of these configurations, giving their dependence on the temperature T the system was equilibrated in.

3.6.1 Energy minimisation

The process of finding the local ground state energy of a system of particles, such that all the forces on the particles are close to zero, is called energy minimisation. In essence, energy minimisation is a mathematical optimization problem, and as such many algorithms devised for these purposes can be used. In this thesis a combination of two algorithms are used, namely a variable step size gradient descent algorithm, and the conjugate gradient method. In gradient descent the system just moves particles in the negative direction of the potential energy gradient, so following the force on the particles. The conjugate gradient method is a more sophisticated algorithm which achieves the same results more efficiently, although this thesis will not go into further detail of how it works. For more information about this algorithm, see [24].

Note that the process of energy minimisation itself does not have a physical meaning; there are no Newtonian dynamics being simulated. However the resulting configurations do describe a physical solid gel.

Generation:

Equilibration:

Production:

For serial = 0, 1, 2, . . . , 95 do:

Give particles initial r, random p and set T = 1.0 Run dynamics for trun = 100

For Tconf= 0.09, 0.1, 0.12, 0.15, 0.2, 0.25 do: Set Tconf and t = 0.0

Run until equilibrated, e.g. trun≈ τR Run for trun = 100τR, storing 100 snapshots

4.

Results

4.1 Energy Evolution

In Figure 11 the total energy evolution of the equilibration and the production steps for the six studied temperatures are displayed, here averaged over all the 96 runs. As can be seen, the systems starts out at the high energy from the generation step, after which the system relaxes to its equilibrium energy. Note that the higher temperatures start at a higher energy, even though the same starting snapshots are used. This is due to the velocities not being stored in the snapshots, but being initialized randomly according to the Maxwell-Boltzmann distribution for that specific temperature [9]. The transition between the equilibrate run and the production run is marked for each temperature with a dotted vertical line. Notably, for T = 0.12 the system moves out of equilibrium during production, indicating that it was not fully equilibrated. Moreover, for the lowest temperature T = 0.09 equilibrium is not fully reached at all. For the lowest two temperatures not all 100 relaxation times could be completed in the 5 available days of computation, and as such less snapshots have been saved for these temperatures.

Figure 11: The evolution of total energy, consisting of the potential energy U and the kinetic energy Ek, of the equilibrate and production steps for different temperatures T , averaged over 96 runs. The vertical lines indicate the transition from equilibration to production.

(a) T = 0.09 (b) T = 0.1 (c) T = 0.12

(d) T = 0.15 (e) T = 0.2 (f ) T = 0.25

Figure 12: Some examples of the largest clusters for each of the studied temperatures present in the production snapshots. For the lowest two temperatures almost the entire system is spanned by a single cluster.

4.2 Structure

4.2.1 Clustering

During and after the system is equilibrated the particles will cluster together. How much clusters form, and how large the clusters are, depends on temperature of the run. A few examples of the largest clusters for different temperatures are shown in Figure 12. For the lowest two temperatures it is possible to see strands of particles form a network through the entire system. For the larger temperatures the cluster size decreases with temperature, resulting in smaller and smaller clusters. Two particles are considered to be in the same cluster if a chain of bonded, or nearest-neighbours, particles connect them. Two particles are bonded if the distance between the two is smaller than rnn = 1.05. This is chosen in favour of rnn = 1.28, which corresponds to the maximum of the potential energy barrier (see Figure 8), as this more closely captures if particles are bonded, even at higher temperatures.

A more quantitative overview of the amount of clusters n and the size of the largest cluster smax with respect to the configuration temperature are given in Figure 13 and Figure 14 respectively. Of course as more clusters combine the maximum cluster size increases while the cluster count decreases. For the lower temperatures nearly all of the particles end up in the one cluster network which penetrates the whole system.

0.10

0.15

0.20

0.25

T

0.00

0.05

0.10

0.15

0.20

n

/N

Figure 13: The average amount of clus-ters in the produced gels over temperature, normalised by the system size. Particles are considered bonded when their distance r < 1.05.

0.10

0.15

0.20

0.25

T

0.0

0.2

0.4

0.6

0.8

1.0

s

/N

Figure 14: The average largest cluster size in the produced gels over temperature, normalised by the system size. Particles are considered bonded when their distance r < 1.05.

4.2.2 Radial distribution function

The structure of the gels can more precisely be described using the radial distribution function g(r). This is a useful quantity in describing the structure of liquids and gels alike, and is given by

g(r) = 1 N 1 4πρr2 N X i N X j6=i δ(r − rij). (40)

The radial distribution describes how much particles there are at a certain distance r from a center particle, normalised for volume. The radial distribution function per temperature can be seen in Figure 15a, averaged over 96 production snapshots. Of course there are no particles for r < 1.0σ due to volume exclusion, and as such g(r) is zero there. The peak at r ≈ 1.04 is the nearest neighbour peak, of which a zoomed in version can be seen in Figure 15b. This peak indicates bonded particles, and rises as the temperature goes down. This confirms that at the lower temperatures more particles are bonded or clustered together. This is also strengthened by looking at the rest of the distribution, showing large peaks for lower temperatures, while the higher temperatures are described by a more flat distribution, indicating less to no long range structure. Moreover, the absence of high r peaks indicate that the system has not crystallised for any of the temperatures. Notably, small but different peaks emerge for T = 0.12.

4.2.3 Structure factor

Another quantity which is often used to study molecular systems is the structure factor S(q), the Fourier transform of the radial distribution. It has as advantages that is is more easily calculated, and can often be directly measured in experiments, as it represents how

1.0 1.5 2.0 2.5 3.0 3.5 4.0

r

(

σ

)

0

1

2

3

4

5

6

7

g

(

r

)

T

= 0

.

09

T

= 0

.

1

T

= 0

.

12

T

= 0

.

15

T

= 0

.

2

T

= 0

.

25

(a) The full radial distribution function.

0.95

1.00

1.05

1.10

1.15

r

(

σ

)

0

5

10

15

20

25

30

35

40

45

g

(

r

)

T

= 0

.

09

T

= 0

.

1

T

= 0

.

12

T

= 0

.

15

T

= 0

.

2

T

= 0

.

25

(b) The distribution zoomed in at the nearest neighbour peak.

Figure 15: The radial distribution g(r) function per temperature, averaged over 96 pro-duction snapshots. Notice that the lower the propro-duction temperature T the more peaks form in the distribution, which indicate that more clustering takes place for these temper-atures.

0

2

4

6

8

10 12 14

q

(

1

/σ

)

0

1

2

3

4

5

6

S

(

q

)

T

= 0

.

09

T

= 0

.

1

T

= 0

.

12

T

= 0

.

15

T

= 0

.

2

T

= 0

.

25

Figure 16: The structure factor S(q) of the simulated gels for different temperatures, av-eraged over 96 production snapshots.

0

2

4

6

8

10 12 14

q

(

1

/σ

)

0

1

2

3

4

5

6

S

(

q

)

T

= 0

.

09

T

= 0

.

1

T

= 0

.

12

T

= 0

.

15

T

= 0

.

2

T

= 0

.

25

Figure 17: The structure factor S(q) of the created solid gels for different temperatures, averaged over 96 minimised snapshots.

0

2

4

6

8

10

12

14

q

(

1

/σ

)

0

1

2

3

4

5

6

S

(

q

)

t

= 5

.

0

10

t

= 1

.

0

10

t

= 3

.

0

10

t

= 5

.

0

10

t

= 1

.

0

10

t

= 3

.

0

10

t

= 5

.

0

10

t

= 1

.

0

10

t

= 2

.

0

10

t

= 3

.

0

10

t

= 4

.

0

10

t

= 5

.

0

10

Figure 18: The evolution of the structure factor over time for T = 0.08 and ρ = 0.24 during equilibration. The increasing spikes indicate the growth of clusters.

a material scatters incident radiation [25]. The structure factor is given explicitly by S(q) = 1 N N X i N X j e−iq·(ri−rj)_, (41)

which can for isotropic materials be rewriten as S(q) = 1 N N X i N X j sin(qrij) qrij . (42)

The structure factors for the gels generated in the production step are given in Figure 16, here averaged again over 96 snapshots. As can be seen the structure factor peaks grow with lower temperature, as the troughs become deeper. This is especially clear for the low qpeak. The structure factors provide evidence for a transition from a flat liquid to a more structured, clustered phase.

To show how the structure changes during the solid gel generation in the minimisation process, the structure factor for the minimised snapshot is given in Figure 17. The lower temperatures do not change much during minimisation, as they where already low in their potential energy minimum. However the higher temperatures do change significantly, becoming more structured after minimisation. Notably the height of the low q peaks for T = 0.2 and T = 0.25 becomes nearly the same, but both peaks remain at their original and different wavelengths q.

In Figure 18 the time evolution of the structure factor during equilibration is illustrated. Here a low temperature system with T = 0.08 and ρ = 0.24 is run for its relaxation time, showing that clustering increases with time during equilibration.

(a) The full stress correlation function evolution for each temperature.

(b) The first part of the correlations, zoomed in on the dip.

(c) The second part of the correlations, zoomed in on the decay after the dip.

Figure 19: The correlation functions of the stress σ over the time difference t for different temperatures, normalised by the system volume over the temperature.

4.3 Correlations

In the equilibration step the relaxation time was estimated using the energy evolution of the system. However, this was just an approximation. Ideally the relaxation time can be measured more precisely. Given a correlation function C(t) of an observable quantity ψ(t) and using time translation invariance it is possible to write that

C(t) = hψ(t)ψ(0)i = 1 tmax− t

Z tmax−t

t0₌₀

ψ(t0+ t)ψ(t0)dt0 (43) evolves, for a liquid, as

C(t) = C0e−t/τ, (44)

where τ is the relaxation time of the system [12]. In other words, the correlations within the system decay exponentially in time, and become negligible for time differences higher

0.10

0.15

0.20

0.25

T

0.00

0.05

0.10

0.15

0.20

0.25

0.30

0.35

G

(

²

σ

)

Figure 20: The shear modulus G of the gels over the production temperature T .

0.10

0.15

0.20

0.25

T

4

5

6

7

8

9

10

11

N

σ

/G

Figure 21: The standard deviation of the shear modulus normalised by its mean, plotted over the production temperature T . than the relaxation time. ψ(t) can in principle be any observerable, as it can be expected that there is only one fundamental time scale of the system (such as the shear relaxation time τR) to which all other time scales are just rescalings. Using this it is possible to find the relaxation time τ = τR per temperature by looking where these correlations of the observable fall off with a factor of 1/e.

In this thesis the correlations of the stress σ on the system is chosen. The correlations of the stress generated using the production runs for each temperature are displayed in Figure 19a. It can be seen that no exponential decay is taking place. For small t the correlations dip below zero, indicating short oscillations, see Figure 19b. The width of the dip is approximately the same for each temperature, but it gets deeper with lower temperature, until below T = 0.12, after which it becomes less deep again. Note that the resolution for the lower temperatures becomes low, but this is due to less values being saved for this relatively shorter time fraction, and has nothing to do with the physics of the system. After the dip the correlations rise and fall again, after which some kind of logarithmic decay seems to occur, with again a strong temperature dependence. For lower Tthere are stronger correlations, although they fall off faster than the higher temperatures.

4.4 Elastic properties

After the snapshots generated during the production run have been minimised using the algorithms described in Section 3.6.1, the shear and bulk moduli of the created solid gels can be measured. By calculating the eigenvalues of the Hessian matrix1 _{of the potential of the} system, the shear and bulk moduli may be found. The exact details are beyond the scope of this thesis, so for more information see [26]. Using this algorithm the elastic properties with respect to the production temperature T were calculated, which are displayed in Figure 20 and 22. As can be seen in Figure 20 the shear modulus increases the lower the production temperature T , meaning the ’harder’ the generated gel becomes. This is to be expected, as the system becomes more clustered for these lower temperatures, in turn strengthening

1

The Hessian matrix is a square matrix of all the second-order partial derivatives of a scalar valued function, which is the potential in this case.

0.10

0.15

0.20

0.25

T

0.20

0.25

0.30

0.35

0.40

0.45

0.50

0.55

K

(

²

σ

)

Figure 22: The bulk modulus K of the gels over the production temperature T .

0.10

0.15

0.20

0.25

T

0

1

2

3

4

5

6

7

N

σ

/K

Figure 23: The standard deviation of the bulk modulus normalised by the its mean, plotted over the production temperature T . the solid network inside the gel. The bulk modulus K is different however, as it has a minimum at T = 0.12, where K is nearly the same for T = 0.09 and T = 0.2. The bulk modulus seems to increase both for higher temperatures and for lower temperatures than T ≈ 0.12. Moreover, the standard deviations σG,K over the mean of the shear and bulk moduli, normalised by √N, are displayed in Figure 21 and 23. This dimensionless and N-independent value represents the mechanical disorder in the system. The fluctuations fall off the higher the temperature for both the shear and the bulk moduli. As can be seen, for the shear modulus at T = 0.12 the fluctuations rise even above those of the lowest two temperatures.

5.

Discussion

Figures 12, 13 and 14 confirm that the particle-particle potential given by Equation 28 does cause particles to cluster together, forming gel-like networks. The amount of clustering that takes place grows as the temperature of the system goes down, and disappears almost completely for temperatures above T = 0.25, where the kinetic energy drowns out the effects of the potential. However, the lower the temperature the more the dynamics of the system slows down, as can be seen in Figure 11. This shows that the relaxation time of the system increases drastically for even a small decrease in temperature.

Of note is the energy evolution for T = 0.12, shown in Figure 11. As can be seen, after seemingly having relaxed, it goes out equilibrium again, linearly lowering its total energy in time. This comes paired with the interesting minimum in the bulk modulus for T = 0.12, as can be seen in Figure 22, which is not present for the shear modulus in Figure 20. Figure 14 seems to indicate that T = 0.12 lies at the transition point where the system evolves from containing a lot of small clusters to a system where these small clusters combine to form larger clusters, eventually spanning the whole system. To better study what is happening here more temperatures around T = 0.12 could be studied, perhaps revealing

more of how this effects evolves with temperature.

What is also clear from Figure 11 is that the lowest two temperatures T = 0.1, 0.09 have not completely equilibrated after the production has been started, with T = 0.09 not reaching full equilibrium at all. For the production run to create accurate snapshots, and for the stress correlations to be properly calculable the system needs to be in equilibrium, such that time translation invariance holds. However for these low temperatures this was not possible in the five days of computation time, as the dynamics have slowed down too much. This has also caused that the system could not be run for 100 production snapshots, only reaching about 28 relaxation times for T = 0.09, and 43 for T = 0.1. Further studies concerning these low temperatures should be able to do computations for longer times, allowing the system to fully relax and gather more production data. Going to even lower temperatures should then also be possible, although this would maybe not give significantly different results, as for T = 0.09 almost the full system has already clustered to a singular network. This indicates a limit in the formation of the gels. However, it would be interesting seeing how the properties of the system evolve with temperature after this limit has been reached.

That the system was not fully relaxed can not cause the dips in the stress correlations however, as the system has been fully relaxed for the higher temperatures, where the dip is mostly present. Again a turning point can be seen at T = 0.12, where the dip is at its minimum. These negative correlations indicate oscillations within the stress, such that at some time t apart the stress is negative what is was before. This maybe indicates numerical instability, but tests using a smaller (and a larger) time step did not change the shape of the dip at all, revealing that this was not a problem. Moreover, it can not possibly be an artifact of bad statistics caused by a lack of data, as the dip exists for every of the 96 production runs per temperature. Moreover, averaging less correlation functions does not change the existence or depth of the dip. This makes it seem like this dip is physical in nature, and that the oscillations may be caused by the formation of smaller clusters. That is, until the system has fully clustered, and the dip disappears. It may be some artifact of using the stress correlations, and using correlation functions of other quantaties, such as the dynamical scattering function, may give better results.

In Figure 20 the shear modulus increases with lowering the production temperature. This was to be expected, as the lower the temperature the more the particles in the system cluster together, which gives the gels their solid properties. For higher temperatures the shear modulus seems to be little temperature dependent, as here the clustering disappears altogether, while it exponentially increases for lower T . Further research could look into the question if this exponential low T growth keeps up for temperatures even lower than T = 0.09, as at this temperature the aforementioned clustering saturates the system. For the bulk modulus, given in Figure 22, the same rise for low T as for the shear modulus can be seen. However, for temperatures higher than T = 0.12 the bulk modulus actually also increases. This perhaps indicates two limits. For T = 0.12 there are a lot of small clusters in the system. Each of these clusters contains closely packed particles, leaving more room in the system for these clusters to move around. As such this system may be compressed more easily than if the system consists of just particles on their own, not in clusters. If the particles are not closely packed, which is the case for higher temperatures, and their mutual distances are larger than r = 1.28σ, then the potential is repulsive, pushing the particles away from each other. This then resists compression, and as such

the bulk modulus goes up. In the other low temperature limit the colloids form a complete network through the entire system, which supports itself. This network is stiff and acts as a solid, pushing back any compression. This then again causes the bulk modulus to increase.

Furthermore, in this thesis only ρ = 0.3 was studied. It would be interesting to find how changing the density would affect the cluster formation and in turn the elastic properties of gels. Understanding this better would allow manufacturers to better tune the properties of the gels than just using the production temperature.

6.

Conclusion

This thesis has given a basic introduction to the theory of soft matter physics and how this theory relates to equilibrium gels. An overview was presented for colloidal simula-tions, which was applied to simulate equilibrium gels with a particle-particle potential consisting of short range attraction and longer range screened electrostatic repulsion. The results presented in this work showed a strong temperature dependence on the formation of clustered particle networks within the gels, from little to no clustering at temperatures T = 0.25 and higher to an almost fully clustered network of particles for temperatures T = 0.1and below. This was confirmed by both directly looking at the cluster count and sizes inside the gels, and by studying the radial distribution function and the structure factor of the systems. Furthermore, the elastic properties of the gels after they have been quenched show the same relation; the shear modulus shows significant increases the lower the temperature, varying between G = 0.0393(28)  σ−3 _{for T = 0.25 to G = 0.27(4)  σ}−3 for T = 0.09. These results combined indicate that the stiffness of gels is indeed caused by the formation of a network of clustered particles through the system. The bulk modulus shows similar results but only for low temperatures below T = 0.12. At this transitional temperature point the system transitions from having a lot of small clusters to the particles forming a singular network spanning the entire gel. This is indicated by a minimum in the bulk modulus at this temperature, where K = 0.240(9)  σ−3 _{for T = 0.12, after which it} rises again for higher temperatures, reaching K = 0.4851(28)  σ−3 _{for T = 0.25. Here the} repulsive part of the potential causes this to be higher than for the lowest temperature T = 0.09, where K = 0.45(4)  σ−3. Finally, the correlation functions of the systems shear stress have been studied to find the temperature dependence of the relaxation time of the gels. While the relaxation times could not be determined using this method, the stress correlations showed an interesting negative dip for small time differences, indicating some form of oscillations in the gels.