Theoretical Physics
Master Thesis
credited 60 ECTS, conducted in the period Sep. 2017 to July 2018
Collective mechanics of a dense suspension of active
particles
by
Ricardo Reyes Garza
11391863
Supervisor :
Assessor :
I would like to thank Dr Edan Lerner for all what he shared with me and his mentorship. I would also like to thank Dr Gustavo D¨uring, from the Pontifical Catholic University of Chile, for giving access to his computer cluster so I could run the numerical simulations that were needed to accomplish this project.
Suspensions are systems composed of particles immersed in a fluid. In athermal suspen-sions, the packing fraction determines the system’s mechanical responses to stresses applied on its boundary. It has been observed and widely documented that above a certain critical value, known as jamming point, the system responds elastically to small applied stresses while below this point, it behaves as a viscous fluid. By considering a suspension made of active particles, we switch the dynamics’ driving mechanism from an external force acting in the boundaries to an inner property of the system. The aim of this project is to describe two-dimensional dense suspensions of hard active disks close to their jamming point; par-ticularly, we focus on the scaling of the particles’ velocity and the contact forces between interacting neighbours with respect to the system’s connectivity. This project comprises a theoretical study, based on the so-called Affine Solvent Model, and numerical simulations of active suspensions and physical networks of nodes and rigid rods, the latter serving as a first approximation to an actual suspension. For the physical networks, the numerical results for the forces between connected nodes are in strong agreement with the theoretically expected scaling, while the nodes’ velocity seems to vanish, as we approach the critical connectivity, slightly faster than the model’s prediction. For the active suspensions, the simulated data of particular ratios of the quantities of interest seem to have an overall behaviour in accord with the theory while the individual scaling of these same observables presents a strong anomalous behaviour except in a small interval of the connectivity. This small interval spans the values of the connectivity that are accessible to most of the simulated suspensions regardless of their packing fraction. These results suggest that the self-organization in active suspensions can give rise to some configurations whose physical properties are consistent with the predictions but also allows other states whose characteristics are beyond the reach of the considered model.
Introduction 1
1 Jamming phenomenon in suspensions 5
1.1 How is mayonnaise made? . . . 5
1.2 The diluted suspension regime . . . 7
1.3 The dense suspension regime . . . 10
1.3.1 Jamming and its critical behaviour . . . 12
1.3.2 The contact network . . . 17
2 Active and dense suspensions- making sense out of the welter 21 2.1 The system’s overdamped mechanics . . . 22
2.2 Theoretical formalism . . . 25
2.2.1 Equations of motion . . . 25
2.2.2 The N and M operators . . . 29
2.2.3 A first approximation to active suspensions- physical networks . . . . 33
2.2.4 A more detailed approach . . . 40
3 Numerical approach to the active case 47
3.1 Preliminary aspects . . . 47 3.1.1 Methodology and simulation scheme . . . 49 3.2 Results and discussion . . . 53
4 Conclusions and suggestions for future research 65
4.1 Suggestions for future research . . . 67
Mixtures of substances that remain as distinguishable phases have an important role and presence in our daily life but most of the time they are hardly ever noticed. This kind of systems usually have the form of discrete elements of matter dispersed in a continuous medium, coining the term dispersions, and if the discrete elements happen to be very small, we will not be able to distinguish them from the continuous phase. Blood is an example of this: The red and white cells have a well-defined boundary with respect to the plasma in which they are transported but are so small that at simple sight the mixture seems as one-component substance. Sand in a jar, on the other hand, is a dispersion (sand grains and air) for which we can distinguish both phases. Clouds, mist, foam, milk, paint and the very air and water surrounding, filled with bacteria, pollen and pollutants, are typical examples of this kind of systems.
From a physical perspective, we may ask how different are these mixtures from each other? Do seemingly similar systems present the same properties as, for example, sand and dry granulated sugar? What about systems of shared components but different concentration? If we add a handful of sand into a bucket of water does the overall behaviour of the systems differs from that of a bucket with pure water? I am confident that we can agree that there is a diametrical difference between a bucket of water with a handful of sand and a bucket of sand with a tablespoon of a water but, how much sand do we need to add to our bucket of water before the system starts behaving just as sand? Is there a point in between these states?
These last questions become particularly interesting for suspensions, this is, for systems that consist of solid particles immersed in a viscous liquid. It has been observed [9,12-19] that the relative density of the solid particles can have a drastic impact on the suspension’s behaviour compared to that of the pure liquid phase. These differences can be as dramatic as a complete loss of the fluidity properties of the system as it undergoes the so-called jamming transition from a fluid-like state to an emergent solid-like state. Intuitively, we can
understand this transition by picturing a system so densely crowded that all its particles are constrained to very small regions of space, so it effectively behaves like a solid. This transition and its similarity to other critical phenomena is by itself a huge topic of debate and research nevertheless, there is still no accepted microscopic theory of this phenomenon. Suspensions are the general topic of this thesis, but we are not interested in an arbitrary kind of suspended particles. Instead, we are drawing our attention to self-propelled, or active, particles which are particles that posses some intrinsic mechanism that allows them to swim or walk through the available space. These particles are widely used in numerical simulations to study the motion of living organisms and even non-living systems as Jannus particles nevertheless, these systems show many collective properties that are still not well understood. The present project aims to provide a description of the contact forces and the velocity of an ensemble of active particles that form a dense suspension near its jamming point, i.e., the critical density at which the system jams.
Before diving into the actual research, we start the first chapter by presenting some previ-ous works regarding jamming and the physics of dense and dilute suspensions. Particularly, we emphasize the critical behaviour of the jamming transition and how its physics can be described by power-laws. After this, we introduce the concept of contact network, which is nothing more than the geometric arrangement described by the particles in contact. The im-portance of this concept cannot be overstated since its topology can be used to characterize the jamming point, so the physical behaviour of the system is linked to a purely topological feature of the particles configuration. This bibliographical research serves as a motivation for this project and also provides us with the conceptual notions necessary to understand the physics of our active suspension.
The second chapter is dedicated to the theoretical description of active suspensions from a microscopic perspective. For this, we first need to specify the characteristics of the systems we are interested in. We have already said that our suspension is made of active particles but there are some other properties that we are interested in: The particles must be perfectly hard, should not describe any sort of Brownian motion and their dynamics must be restricted to an overdamped limit. After broadly discussing these properties and their implications, we dive into the theory’s formalism: The Affine Solvent Model [33]. We also introduce an auxiliary system, a physical network of nodes connected by rods, that serve as the link between the active suspensions and their corresponding contact network. In this chapter, we also derive some theoretical prediction for the system’s pressure, contact forces and particles’ velocities. These predictions are shown as power-law functions that depend on the average number of contacts in the system and its distance from the jamming point.
Along with the theoretical work, we numerically simulated active suspensions and physical networks. In the third chapter, we present the methodology behind these simulations as well as the data obtained for the relations between several measurable quantities. The numerical results are then compared to our theoretical model’s predictions. Finally, the fourth chapter is dedicated to a short overview of the whole work, final conclusions, remarks and some suggestions for future research.
With this, we conclude the description of the topics that are presented in this written work. So, let us start from the beginning: the first chapter.
Jamming phenomenon in suspensions
1.1
How is mayonnaise made?
Any traditional mayonnaise recipe contains two main ingredients: egg yolks and vegetable oil. There are other ingredients, like whole eggs, mustard, salt, lemon juice, vinegar and so on, that can be used to season the sauce and help to stabilize it. What do we mean by stabilizing it ? Well, egg yolks are made of very complex organic molecules but their main compound (that can make up to 51% of the egg yolk) is water [1] and, as it is commonly known, mixing water and oil together is a rather difficult achievement. Usually, if the recipe is not followed properly or if the techniques are not adequate, our oil and water emulsion will separate, this is, our mayonnaise as a mixture was not stable.
The regular procedure for making mayonnaise sauce starts by whisking the egg yolks together with the seasoning ingredients if any. After obtaining a bright yellow mixture the oil is added in a very thin stream while still whisking. The slower the beating speed (a blender is faster than a hand with a fork) the thinner the stream should be. The mixture will become thicker and thicker and when the desired texture is obtained the mayonnaise is finished. It is frequently unnoticed that through these simple steps we have obtained a remarkable thing: We have combined two low-viscous fluids, which are typically immiscible, and obtained a high-viscous substance. In fact, mayonnaise is considered to be a shear-thinning fluid [2], this is, a fluid whose viscosity decreases under a shear stress and can have an elastic respond to small stresses [4].
In order to understand why the resulting substance behaves so different compared to 5
its original constituents, we need to take a look into the microscopic phenomenon. As we pointed out, an important part of the process is to add the mix while vigorously whisking. This will break the oil into small droplets so the resulting system consists of low-viscosity drops immersed in a strong flow of egg yolk. Then, as we keep stirring, two things are going to happen: 1) The amphiphilic lipids 1 get in contact with the oil droplets creating a hydrophilic external layer which prevent them to merge and become large. 2) The flow will deform the droplets and eventually burst them into smaller ones which, under de adequate stirring, can be of microscopic size [3, 5]. So, up to now, we have a stable, aqueous solution of very small oil droplets that, thanks to the lipids in the egg yolk, will not merge together or create a different phase from the water. As we add more and more oil, we create more and more drops while keeping the amount of the rest of the ingredients constant. This means that we are increasing the droplet volume fraction which eventually will cause the system to jam, i.e., the system as a whole will develop a shear modulus high enough so that it will behave as an elastic, disordered solid at least for small stresses.
Our daily life is filled with systems that go from freely moving to mechanically stable states or vice versa. Back to the cooking examples, the reader might be familiar with the experience of trying to pour some flour from a bag into a bowl by gently shaking it, failing in doing so, we usually proceed by shaking the bag a little bit harder and ending up with a mess all over the table. In this example, we see the contrary of what happened while making the mayonnaise: we started with a static configuration of flour and, after mechanically perturbating the system, it lost its stability and started to flow. Glass, foams, pastes, cars in a highway and the users of public transport are other examples of systems that can be encounter in their static or fluid phase depending on the applied stresses, their temperature or the density of the suspended particles [4, 12]. In all these cases, the solid phase lacks the long-range periodic order that characterizes crystalline solids. Instead, their constituents are restricted to a very small part of the phase space in a disordered configuration either because of the proximity of their neighbours or its low thermal energy. They have simply jammed.
For the mayonnaise, the jamming transition is characterized by the density of the number of oil drops or, more precisely, by its packing fraction, φ, which is defined as the ratio of the volume occupied by the solute (oil drops in this case) to the system’s total volume. If we do not add enough oil or if we do not create enough droplets, the packing fraction will be too
1Anphiphilic substances are chemical compounds that have hydrophilic and lipophilic properties at the
same time. Usually one of its ends has the hydrophilic property while the other is lipophilic. This kind of compounds are commonly found in organic matter and egg yolk is not the exception. As a friendly tip, adding salt will make the lipids and the proteins of the egg yolk to separate more easily so the emulsion stabilizes quicker.
low and the result will be a runny mixture rather than the thick sauce we are used to and, on the other side of the spectrum, by adding too much oil the final product will resemble a gelatine.
It is important to notice that this process is happening at a certain lengthscale. If we study mayonnaise at atomic or molecular scales we would observe incessant moving molecules, i.e., a liquid-like behaviour. The jamming transition is happening at lengthscales of the order of 1 µm [6], which is roughly the size of the oil droplets. At this scale and for a high packing fraction, the oil droplets movement is constrained by their neighbours to a very small region. The higher φ is, the smaller the configuration space they can explore until we observe a jamming transition. Despite our interest in high values of φ, very interesting things already occur in the low packing regime. Before adventuring into the physics of high densities, it might be useful to briefly discuss the low density case since it will provide us with an intuitive understanding of the phenomenon.
1.2
The diluted suspension regime
The 18th and 19th-century were marked, among many other important things, by a heated debate: The existence of atoms2. Rudolf Clausius, Michael Faraday, James Maxwell, Ludwig
Boltzmann and Max Planck are just a few names of the great minds that took a side either in favour or against the discrete nature of matter [7, 8]. During the first half of the 18th-century, thinkers, like Daniel Bernoulli in his book Hydrodynamica (1733), acknowledged that it was possible to derive the experimentally known macroscopic laws of gases by assuming they are made of small, hard and fast spheres. It was nearly two centuries later, in 19053,
when Albert Einstein used this assumption in his article On the movement of small particles suspended in a stationary liquid demanded by the molecular-kinetic theory of heat (published originally in German in Annalen der Physik, 17) to obtain a theoretical description of a phenomenon that could not be understood using classical thermodynamic theory: the motion
2It is worthy to clarify that during this time the notion of atom was far from our current understanding of
these quantum fellows. Naturally, their inner structure was indeed an important question of the debate that would become crucial for electrodynamics and quantum mechanics but, ultimately, the deepest consequences of accepting their existence transcends physic itself. Atomism, as a philosophical posture, encourage us to accept that there is more in nature than what our senses can perceive.
3The year 1905 has an almost mystic aura around it. During this year, Einstein published five
ground-breaking works: One related to the photoelectric effect (for which he was awarded the Physics Nobel prize in 1912), two more about statistical mechanics and Brownian motion, one dedicated to special relativity and lastly one where he stated his famous relation between matter and energy.
that microscopically-visible bodies carry out when suspended in a liquid, i.e., Brownian motion [9]. It was until 1913 when the French physicist and chemist Jean Baptiste Perrin put an end to the debate by providing experimental confirmation of Einstein’s predictions in his book The Atoms.
Even though Einstein would have to wait eight years before seeing Perrin’s results and conclusions, he certainly never doubted the truth within the atomistic view of matter. Back in 1905, Einstein presented his doctoral thesis A new determination for molecular dimension. As its title already suggests, Einstein was trying to find a way to determine the size of the molecules of the substances, particularly liquids, that he was sure made matter up and whose size was so small that any direct measurement was impossible. Einstein started by considering a mixture of a solute and a solvent that fulfilled two conditions: 1)The size of the solute’s particles, even though invisible for the naked eye, was big compared to the solvent’s molecule. Notice that by particle we are referring to objects of a characteristic length much bigger than the size of the molecules of the solvent. 2) The solvent’s total volume was small compared to that of the solvent. Einstein wanted these assumptions to hold because they allowed him to describe the system through the theory of hydrodynamic but there is something more in them, they already resemble the first stages of our mayonnaise example when we had just started to add a small amount of oil into our egg yolks, i.e., they are sufficient conditions for a diluted suspension4
We will not explicitly discuss the details of Einstein calculation here since they do not provide a significant insight towards our goal nevertheless, a very nicely written overview can be found in [10]. In general terms, Einstein’s approach consisted of two steps. First, he considered an incompressible fluid provided with a stationary velocity profile, v. Because of the fluid’s viscosity, there is a intrinsic stress that arises from the flow. Given some arbitrary region Ω, these stresses will perform some work (per unit time), W , in its boundary, ∂Ω, that will dissipate in the form of heat due to the fluid’s viscosity. This work is given as
W = Z
∂Ω
viσijnjdA, (1.2.1)
where σ is the system’s stress tensor and n is a unit vector normal to the infinitesimal area element dA and we are summing over repeated indices. In this environment, Einstein placed a tiny sphere and determined the perturbation that it would cause to v. From this, by taking Ω to be a sphere of volume V enclosing the tiny sphere, where V is much bigger than the
4In chemistry, there is a technical difference between suspensions, colloids and emulsions depending on
the properties of the suspended particles. For our purpose, there is no need to such distinction so we will grant ourselves the freedom to misuse the language and talk about suspensions when referring to a liquid medium with small particles suspended in it.
particle’s volume, W becomes W = Cη V + Φ 2 = CηV + CηΦ 2, (1.2.2)
where Φ is the volume of the tinny sphere, η is the fluid’s dynamical viscosity (this is, the parameter that quantifies the system’s internal friction produced by shear stresses) and C in some constant that depends on the particular velocity profile. We can immediately see that, in the absence of the tiny sphere, the work is simply W = CηV
Done this, Einstein’s second step consisted in determining what would happen if instead of having one tiny sphere we have n tinny spheres of the same size, homogeneously distributed in the fluid and separated by much bigger distances than their radii, i.e, Einstein explored the case of a diluted suspension composed of a fluid and suspended tiny spheres. In this case the work, W will be
W = CηV 1 + 5φ 2 , (1.2.3)
where φ is a number smaller than one and is equal to the ratio between the total volume occupied by the n tiny spheres and the volume V , i.e., the system’s packing fraction.
Let us now suppose that we are presented with a fluid to which some tiny spheres have been secretly added. In our innocence, we will suppose that this fluid is completely clean, so we are expecting that if it starts flowing, its viscosity will produce a work C ˜ηV for some ˜η. At the same time, the value of this product has to coincide with that of the expression 1.2.3 which can only happen if the viscosity that we measure, ˜η, is related to the original fluid’s viscosity as ˜ η = η 1 + 5φ 2 , (1.2.4)
this is, the apparent viscosity of the system as a whole has increased by an additive term proportional to the packing fraction of the system!
Einstein continued his thesis by describing how does the system’s diffusion coefficient changes with respect to the pure solvent case and finished it by calculating Avogadro’s number using some experimental data. With all these elements, Einstein completed his goal since, by knowing the mass of the solute, we can know of how many molecules it is composed (with the help of Avogadro’s number), then we can measure the ratio between the viscosities to obtain a numerical value of φ and, from it, we get the radius of the molecules assuming they perfect sphere which, in the worst case scenario, gives us at least the order of magnitude of the molecule’s size. As we have said, Einstein developed these ideas in a time
when atomism was still considered a heuristic curiosity among important scientists but do not worry, Einstein’s dissertation was accepted and approved.
The great contribution of Einstein to our cause is the expression 1.2.4. Just by adding small particles into our fluid, the macroscopic behaviour of the system as a whole changes! The suspended particles hinder the flow of the fluid so, effectively, what we observe is an increment in the viscosity. The more particles we add to our mixture, the bigger φ and, consequently, ˜η will be. Nevertheless, if we keep adding suspended particles, eventually we will move away from the diluted regime and the particles will start noticing and interacting with each other. Even in this case, we might be tempted to say that viscosity should keep increasing but, as we already know for our mayonnaise example, there is a drastic change in the macroscopic behaviour of the system, there is a jamming transition. Let us, then, move forward into the discussion of the physics of dense systems.
1.3
The dense suspension regime and its jamming
tran-sition
As we have already discussed through our mayonnaise example, a system of suspended particles undergoes a transition between a fluid and a solid phase if the suspension’s density is high enough. This solidity is an emergent property of the system that strongly depends on the interaction between particles, otherwise, the system would continue to behave as a fluid of increasing viscosity. The solids that arise from jamming do not form an ordered arrangement so, in this transition, there has not been any continuous symmetry breaking as opposed to other kinds of phase transitions that form crystalline solids. In this sense, we might ask how different is jamming to other critical phenomena and even if it can be considered such kind of process.
In 1998, Andrea Liu and Sydney Nagel proposed in their article Jamming is not just cool any more [15] the idea that jamming might be used to describe a wide variety of systems re-gardless of the type of interaction between their constituents and characteristic size. Powders and sand can be driven into jamming by an applied pressure or stress [18] and liquids might jam into a disordered solid through a glass transition5 controlled by the system’s
tempera-5Glasses, in contrast to crystal, are disordered materials that behave as solids. The most common way to
obtain a glass is to cool a liquid below its melting temperature while avoiding crystallization. The resulting system is a highly viscous substance know as supercooled liquids and as the temperature keeps decreasing the systems solidifies and becomes a glass. We can already see that the glass transition sounds pretty much like
ture. With this in mind, Liu and Nagel proposed that the jamming of these diverse systems can be tied together through a phase diagram as the one shown in the figure 1.1. As we can see in this figure, the space is divided into two regions: The jammed region surrounded by the unjammed region. For granular materials, bubbles, droplets and grains, the particles usually do not have a kinetic energy that we can interpret as the system’s temperature so the jamming transition can be described in the Σ − φ−1 plane, while the T − φ−1 plane would be the appropriate one for the glass transition of a supercooled liquid.
For the authors, this diagram also suggests that the three parameters T, Σ and φ are relevant for every kind of system. For example, by shaking a granular material, one should obtain a liquid-like behaviour since we are adding an effective temperature. It also insinuates that the physics of the transition might be controlled by its proximity to the jammed region, which would explain why different systems behave so similarly. Even more, for the particular case where both the temperature and the applied shear are zero, which is an easy system to work with, it is possible to get into the jammed region only by increasing its packing fraction! The value of φc at which the system jams, marked as J in the diagram, is known
as the jamming point.
Figure 1.1: Scheme of the especulative jamming phase diagram proposed by Liu and Nagel. Here, T indicates the system’s teperature, φ its packing fraction, Σ is the applied shear stress and J is the jamming point. Retrived from [16].
The existence of the jamming point was numerically observed years before Liu and Nagel proposed their phase diagram. Through numerical simulations, it was observed that for a two-dimensional system of a bidisperse mixture6 of perfectly hard, frictionless spheres (this
jamming. The relation between them is still a fascinating topic for research and discussion, nevertheless, it goes beyond the scope o the present work.
is, a mixture of impenetrable and undeformable spheres characterized by two different radii), the jamming point occurred at φ = 0.84 regardless of the ratio of the two different radii and the relative number of spheres of each size [21]. Later, in 2002 and 2003, Liu and Nagel, together with other authors, published several articles [16, 19] showing, again through numerical simulations but now of soft frictionless spherical particles, that the jamming point exists and has a value of φc ≈ 0.84 for two-dimensional systems and φc ≈ 0.64 for
three-dimensional systems independently of the number of particles in the system and the particular way they interact. This interaction is usually modelled as a repulsive potential that depends on the overlap between the spheres and has a form
V (rij) = (1 − rij/Rij)α/α if rij < Rij, 0 if rij ≥ Rij, (1.3.1)
where rij is the distance between the two centers of the ith and jth sphere, Rij is the sum of
their radii, is an energy scale and α is a number that determines the potential [4, 16]. So, indeed, there is a particular point in the φ−1 axis that separates the liquid and solid phase of a suspension. Starting with a liquid-like system of suspended particles characterized by some packing fraction φ < φc, we will observe an increment in its viscosity as we approach
the jamming point until it diverges and the system becomes a solid. This behaviour has been already experimentally observed for two-dimensional suspensions of hard spheres [17]. Starting from a solid configuration, i.e. φ > φc, the shear modulus of the system will decrease
as φ −→ φ+
c until it completely vanishes and the system becomes a fluid [13].
Now that we are sure that the jamming point exists, how much does the system’s physics actually depend on its proximity to it? With this we mean: Given a system of suspended particles characterized by a packing fraction φ, can we describe its macroscopic observables as mathematical functions that depend on (φ − φc) ? Let us proceed by answering this question.
1.3.1
Jamming and its critical behaviour
Phase transitions are part of our daily life yet the physics involved can be rather com-plicated. The behaviour of macroscopic systems near their critical point (their critical be-haviour ) arises as an emergent property coming from the intricate interaction between its many constituents. Remarkably, near the critical point of a transition, the details of these
an ordered state, i.e. the system crystallizes, and then we cannot talk about a jamming transition since we would not be describing an amorphous solid. [21, 22]
interactions are completely washed away so many different systems behave in the exact same way! For example, it has been theoretically predicted and experimentally observed [23] that for a liquid-gas transition, the system’s heat capacity scales as
C(T ) ∼ |T − Tc|−0.11..., (1.3.2)
where T is the system’s temperature and Tc is the critical temperature. Amazingly enough,
the same power law, up to its proper approximation, has been numerically observed for the 3D Ising Model [24]. Even though we might think of these two as unrelated systems, their critical behaviour is identical given that the scaling of their heat capacities shares the same critical exponent, ≈ −0.11. This feature of shared behaviour is called universality and the systems that share the same critical exponents are said to belong to the same universality class.
Whether jamming can be studied under the scope of traditional phase transitions theory is still a matter of discussion, although there is numerical and experimental evidence that suggests there are many similarities between the two phenomena [25, 26]. If we are going to consider the jamming as a critical phenomenon, we would expect that, close to φc, some of
the macroscopic observables scale as power laws characterized by their exponents.
In 2007, Peter Olsson and Stephen Teitel published an article where they numerically proved a statement that we have been mentioning over and over again: the divergence on the viscosity at φc [25]. In this article, they studied two-dimensional binary mixtures of soft
particles characterized by the repulsive potential 1.3.1 with α = 2. They found that, as they increase the system’s packing fraction, the viscosity, η, scales as
η ∼ |δφ|−1.65, (1.3.3)
where δφ ≡ φc− φ. So, indeed, the viscosity diverges at φc.
In the same article where Liu, Nagel, et al. obtained the value of the jamming point [16], they also reported quite interesting scaling relations. Starting from the solid-like phase, the system’s shear modulus, G, scales as
G ∼ |δφ|0.50 for α = 2 |δφ|1.0 for α = 5/2, (1.3.4)
where, once more, α is the characteristic number of the potential 1.3.1. This result does not depend on the dimension of the system or the polydispersity of the particles’ radii but it does depend on the details of the interactions, i.e., on α!
Another macroscopic observable discussed in [16] is the pressure of the system associated with the collectiveness of jamming. In order to understand what we are talking about, let us think for a moment that we are inside a bus that has stopped at a crowded bus stop. Before all the additional passengers get inside, there is an intrinsic pressure acting on us, namely, the atmospheric pressure. As more and more people get in, we all start getting closer and closer. Eventually, we will have to start pushing each other so that more people can get in. Intuitively, we can visualize that our bodies will be subjected to a higher pressure than before due to the contacts that have formed. If we manage to measure the contact pressure between each interacting passenger, we can define the collective pressure of the system as the sum of all these individual pressures. The same happens in our mayonnaise example where there is a collective pressure due to the contacts between oil droplets in addition to the intrinsic fluid’s pressure and the same happened in the simulations we are discussing. This collective pressure is right at the core of our project, so let us take a moment to understand it better. Let’s assume that our system of suspended particles has a packing fraction φ, a temper-ature T = 0 and shear stress σ = 0. In order to avoid dealing with the physics in the surface of our system, let us concentrate in a volume V in the solid away from the boundaries. We can define the linear momentum within V as
PV =
X
i∈V
pi,
where i denotes the ith particle inside V . If we have a homogeneous system, we would expect
the net flux through V to be zero and, assuming there is no external forces acting on V , a change in PV can only be due to the contact forces. So,
˙ PV =
Z
V
∂βσαβ,
where σαβ is a stress tensor associated to V such that the right-hand side in the above
expression becomes the average contact force in V . Then, it turns out that σαβ = − 1 V X i>j in V dΦ drij (~rij)α(~rij)β rij , (1.3.5)
where we are summing over pairs in contact, Φ is the interaction potential given in 1.3.1 (with a change in notation to avoid confusion between the symbol of volume and potential), rij is the distance between the ith and jth particles in contact, (~rij)α is the α component of
the vector ~rij between the centers of the particles in contact and the same for (~rij)β. The
In order to proceed, we need to observe that the factor −dΦ/drij is the value of the
force between the particles in contact, fij and that this force is parallel to ~rij provided that
~
fij = −∇Φ. Now, we will define the pressure tensor, Pαβ, as minus the stress tensor 1.3.5.
Pαβ = 1 V X i>j in V fij (~rij)α(~rij)β rij = 1 V X i>j in V fij ~rij rij · ˆeα(~rij)β = 1 V X i>j in V ( ~fij)α(~rij)β. Or simply P = 1 V X i>j in V ~ fij ⊗ ~rij. (1.3.6)
From the theory of elasticity we already know the the pressure of a system is p = −T r(σ)/d, where d is the system’s dimension. So, using the pressure tensor, we can write it as p = 1 dT r(P) = 1 dV X i>j in V ~ fij · ~rij = 1 dV X i>j in V fijrij. (1.3.7)
Going back to our main topic of discussion, Liu, Nagel, et al. found that, for the systems they studied, the pressure scales as
p ∼ |δφ|1.0 for α = 2, |δφ|1.5 for α = 5/2. (1.3.8)
Be aware that in this last expression α corresponds to the parameter of the interaction po-tential and not to a component of a tensor. Again, these results depend on the particular interaction but now we can understand why. As the authors cleverly observed, we can synthe-size the results in 1.3.8 as p ∼ |δφ|α−1. Now, if the system’s response under a compression is
perfectly affine, or at least approximately, the pressure must scale with the same exponent as the force, which we already know scales as something to the α − 1 (just derive Φ with respect to rij!). Another reported scaling relation is that of the system’s bulk modulus, defined as
B ≡ φdp/dφ. Just by its definition, we expect it to scale as B ∼ |δφ|α−2, which is exactly
What about the shear modulus? Well, the authors took a particular shear strain, γ, for which the shear modulus took the form G = d(−Pxy)/dγ. Since, as for B, there are two
derivatives of the potential involved, we would expect an α−2 exponent instead of the α−1.5 exponent of 1.3.4. Here is where the subtleties of jamming make their entrance. By adding a shearing to the system, we are driving it out of equilibrium so it will have a zero-time response, i.e. a change in its macroscopic quantities that occur before any particle changes its position, and an infinite-time response, this is, the value of its macroscopic quantities after the particles had reorganized and the system’s energy had reached a local minimum. The scaling in 1.3.4 corresponds to the infinite-time response while, for the zero-time modulus, they found G0 ∼ |δφ|α−2 as expected. About the difference between the exponents of these
two limits, the authors argued that it is originated in an effective non-affine deformation caused by the shearing since it changes the relation between the contact forces and the pressure.
There is one more relation from which we can get a deeper insight into jamming: the scaling of the connectivity. When we wrote the expression 1.3.5 we stated that only the particles in contact have a role in the collective physics of the system, which might seem obvious since the interaction potential 1.3.1 is non-zero only for particles in contact. We might then ask, can we describe jamming just by knowing the pairs in contact and their overlap? As we will see, for the case of frictionless spheres the answer is ‘yes’. Just as a side note, we will not discuss in the present work either the case of non-spherical or frictional particles nevertheless, the curious readers might want to take a look at [4] and [32].
1.3.2
The contact network
Given a system of suspended particles, we define the contact network of the system as the set of all the contacts, characterized by their direction with respect to the laboratory frame of reference. In figure 1.2 we can see an example of such set of contacts for a binary mixture of hard disks. By looking at figure 1.2b we can appreciate that describing them as a network it is not just a linguistic whim, they do form a network whose nodes and edges correspond to the centers of the disks and lines connecting the centers of particles in contact contacts respectively. Naturally, for a system of hard spheres, the length of the edges is given by the radii of the particles in contact whereas for soft particles it also depends on the overlap.
(a) (b)
Figure 1.2: Example of the contact network of a binary mixture of hard spheres at φ = 0.80.(a) Both, contacts and particles are shown, the smaller ones in blue and the bigger ones in purple.(b) Only the network is visible. In both images the width of the contact in proportional to the repulsive force of the contact.
How can we use the notion of a network to study the jamming transition? Well, let us think of four points connected by four rigid rods as in figure1.3a. Even though the rod’s length cannot be changed, we can change the relative distance between some of the points as shown in figure 1.3b. In this sense, the arrange of four points and four rigid rods does not work as an analogy of a rigid body. On the other hand, if we add a fifth rod as in figure 1.3c then the relative distance between the points is impossible to change unless we change the length of the rods. Just by adding an additional constraint, we created a frame that
resembles a rigid body. In the case we substitute the rigid rods by springs, we will reach the same conclusion: Five springs guarantee that the only way to modify the distances between the points is by compressing or stretching them. Ultimately, the rigid rod network serves as an analogy of a rigid body while the spring network models a solid that has an elastic response.
(a) (b) (c)
Figure 1.3: (a) & (b) Two different configurations of four points, labelled as A, B, C and D, connected by four rigid rods. (c) Four points connected by five rigid rods.
Let us assume we have set of N points in a d-dimensional space. Do we need to explicitly fix de dN degrees of freedom of the system in order to obtain a rigid body? No, in fact, just by giving dN − d(d + 1)/2 constraints in the form of fixed distances between the points we guarantee that the system will become rigid [29]. As in figure 1.3c, even though we had eight degrees of freedom, we just needed five rods to solidify it. It is easy to convince our self of this by looking at the two-dimensional case. Note that the d(d + 1)/2 term comes from the rigid body degrees of freedom (translations and rotations) that do not compromise its rigidity.
With this in mind, the idea of a network of contacts become more useful. If it has the adequate geometry, i.e. number of edges, the particles involved in its formation will collectively behave like a solid since they cannot change their relative position unless we modified the network itself by taking contacts out of it. At the same time, if there are not enough contacts in the network, the system will manifest a liquid-like behaviour provided that the particles, or at least some of them, can change their relative position without affecting the set of contacts and, consequently, the overall interaction potential. Under this scope, jamming becomes a topological property of the contact network! Just by distinguishing between a rigid and a floppy contact network we can determine if the system is in its solid or liquid phase respectively.
Let us assume we have a contact network that has exactly dN − d(d + 1)/2 edges. It is clear that we must be right at the jamming point since having one contact less will allow us to deform the network without changing the edges’ lengths. As we have stated, even as a solid, the system preserves its continuous symmetries. In particular, this means that the system must be homogeneous (despite its inevitable statistical fluctuations) which implies that the contact network should also be homogeneous. Thus, every particle might have, on average, the same number of contacts, zc. Then, the next expression must hold:
dN −d(d + 1)
2 =
zcN
2 , (1.3.9)
where we need the one-half factor on the right-hand side so we do not overcount the number of contacts. From this last expression, and assuming that N 1, we obtain that
zc' 2d. (1.3.10)
In general, we can define the network’s connectivity, z, as the average number of contacts per particle,
z ≡ 2Nc
N , (1.3.11)
where Nc in the number of contacts and, as we have learned a couple of lines above, the
jamming point corresponds to a connectivity of 2d. When talking about networks, instead of jammed network the usual term is isostatic network. So we have found that the jamming point corresponds to the isostatic point of its contact network. Even more, we can use the connectivity as a control parameter of the jamming transition.
In their article, Liu, Nagel, et al. found that in the rigid regime,
|δz| ∼ |δφ|0.50, (1.3.12)
where δz ≡ z − zc, for every single studied system, regardless of its dimensionality or
inter-action potential. So, indeed, we can express all the scaling relations we have seen in term of δz, confirming that the jamming transition and its criticality can be understood as emergent properties of the geometry of the contact network.
From all that we have already discuss about jamming in suspensions, we can distinguish that this transition is very similar to other critical phenomena in that its macroscopic physics is governed by the distance between the system’s state and a critical point, and it can be expressed in the forms of power laws. Nevertheless, the exponents that appear in this scaling relations depend on the particular form of the interaction between particles but not in the dimensionality of the system, as opposed to traditional critical phenomena. Furthermore, we
have learned that jamming is a manifestation of the geometry of the network formed by the particles in contact. Having all this in mind, let us then conclude this chapter and proceed to study the system we actually want to describe: a suspension active particles.
Active and dense suspensions- making
sense out of the welter
After reviewing some general aspects of the jamming transition, is time for us to tackle the objective of this work. The system we are interested in is a dense suspension of non-Brownian, active, frictionless hard disks. In this chapter, we will develop a microscopic theoretical description for this kind of systems as well as predictions for the scaling behaviours of its collective pressure. These relations were tested through numerical simulations but we will dedicate the third chapter of this work to that part of the research. Before presenting this theoretical description, it is worth the time to precisely describe the characteristics of the system we will be studying.
First of all, we will consider that our suspended particles are perfectly hard spheres. As we have already stated in the previous chapter, this means that particles are infinitely stiff so they cannot be deformed and there cannot be overlaps between particles in contact. A more formal way to state this is by setting the interaction potential as follows:
Φh.s.(rij) = lim→∞(1 − rij/Rij)2/2 if rij < Rij, 0 if rij ≥ Rij, (2.0.1)
where we are using the same notation as in the expression 1.3.1 although we have chosen a repulsive harmonic potential (α = 2) jut for convenience and the subscript h.s. refers to hard spheres.
Secondly, they are active particles. This means that they possess an intrinsic mechanism that allows them to self-propelled through the medium they are immersed in. The features,
virtues and vices of this mechanism are not important for us but we can always think of them as living organisms crawling in the bottom of a pond or bees in a beehive although active particles are not exclusively used to simulate animals since some Janus particles might show a self-propelled behaviour. In any case, this auto-propulsion is characterized by a particular orientation and magnitude. We will be interested in the case where:
• The orientation and magnitude do not change in time for any of the particles. • The self-propulsion between the suspended particles is uncorrelated.
When talking about movement, orientation and magnitude, we think, almost as a Pavlo-vian reflex, about equations of motion. Indeed, we will implement these two points into the particles’ equations of motion but before that, we need to know what kind of dynamics we will be dealing with.
2.1
The system’s overdamped mechanics
Most of our daily life experiences occur in what in physics is called a high Reynolds number regime. Even though we are perpetually immersed in a fluid, the air around us, we barely notice any change in our movement due to its viscosity. This is to say, that the friction between our body and the air is so low that it is not involved in the dissipation of our momentum. Following the same line of thought, if the viscosity of the air happened to be much higher, we would actually experience a noticeable loss in our inertia and the higher the viscosity, the shorter the time scale in which this would happen. Mathematically, for spherical objects these effects can be written into their equations of motion as a drag force that, when neglecting turbulence effects, is proportional to the speed of the object but acts in the opposite direction, i.e.
Fdrag = −kv, (2.1.1)
where k is some positive number that depends on the fluid’s viscosity and the characteristic length of the object. The drag force written in this fashion is known as Stoke’s law. Let us consider an object moving through a fluid and subjected to a constant force of magnitude f . Then, its equation of motion, in one dimension, has the form
˙v = f − kv, (2.1.2)
and its solution is
v(t) = f k − Ce
−kt
where C is determined by the initial conditions.
We can observe that 1/k sets a timescale for the dissipation of the initial momentum of the particle. For long times with respect to this scale, v(t) ≈ constant so the object has lost all its inertia and its movement is only driven by the force acting on it. The higher the viscosity, the shorter the time scale 1/k. So, in the limit of very high viscosity, we would only observe that an object moves when a force is acting on it and that remains at rest otherwise1.
Interestingly, these effects do not exclusively depend on the medium’s viscosity, η. It turns out that the characteristic length of the object, Lo, its density, ρ, and speed, Vo, determine
the relevance of inertia in its mechanics. All these quantities are combined in one single dimensionless parameter, the Reynolds number :
Re ≡
ρLoVo
η . (2.1.4)
For systems with high Reynolds number, the effects of viscosity are negligible while for low Reynolds number systems they become important [11]. For example, an ice skater does not need to worry about the air viscosity since Re ∼ 105 while for a fairyfly (an insect whose
length is of the order of 0.1 mm) Re ∼ 10−3, so the air seems to it like a thick syrup2.
We will assume that our suspension of active particles has a very low Reynolds number so that at every moment, the velocity of the particles is equal to the force acting on them. Effectively, this means that, when writing their equations of motion, we will not add any term that contains a time derivative of the velocity3. The implication of this is that, if we want
the particles to self-propelled through the medium, the active mechanism must manifest as a force acting of the particles, the active force. So, for each one of the suspended particles, their equation of motion reads
F − ξV = 0, (2.1.5)
where F consists of all the forces acting on the spheres. Notice also how we have already taken into account the non-Brownian nature of the particles since there are no stochastic terms in the equation. Assuming that there is no friction between particles in contact, no lubrication, capillarity effects, sedimentation or any other hydrodynamic effects, the F term contains only the active force plus the forces arising from each contact. In the same way there
1If we suddenly set f=0, the object’s equation of motion would be of the form ˙v = −kv with solution
v(t) = Ae−kt. If 1/k is shorter than any of our measurement times, we would observe that the object stopped moving exactly when f was turned off.
2A broader (and very nice) discussion regarding systems of low Reynolds number can be found in Edward
Purcell’s article Life at low Reynolds number, American Journal of Physics, 45 (1977).
is a normal force between a table and a book placed on top of it due to the book’s weight, two hard spheres pushing each other will exert a repulsive force in one another. Thus, we can write our equation of motion as
Fact+ Fcont− ξV = 0, (2.1.6)
where act and cont refer to the active force and contact forces respectively. Going back to the bullet points at the beginning of this chapter, the first point is fulfilled by setting for each particle Fact = constant. Regarding the second bullet point, we will consider that each of the components of the active force is a uniformly distributed random number between (−1, −1) so that there is no correlation between the active forces driving any of the particles.
There is a subtlety in all that we have just said. How are the contacts formed? In a crowded environment of active spheres, we can imagine them moving around and colliding with each other. What is more difficult to notice is that once two particles have collided, they do not bounce back because of the overdamped dynamics. They remain stick together while sliding around each other until they get a free path to continue their movement.
As shown in figure 2.1, when two particles collide we can decompose their active forces into two components: One parallel to the straight line the centers of the two particles, Fn
and a second one orthogonal to the first one, Ft. Notice that the sum of these two forces is
not equal to the active force since Fn contains already the contact force, which by geometric
arguments, also has to act parallel to the line between the centers of the two spheres. We can visualize the movement during the next instant of time as two steps: First, the particles will infinitesimal move away from each other in the Ft direction and then they will get back
together by moving in the Fndirection. Because there is no friction between the spheres, each
of the active forces’ direction remain unchanged nevertheless, the line between the centers did change its orientation with respect to the observer, it becomes slightly flatter. Thus, the new decomposition will be such that the Ft components are bigger than in the previous
step. Repeating this steps multiples times, we will reach to a configuration in which the particles will continue their motion away from each other. Of course, this particular scenario corresponds to that shown in figure 2.1. Surely we can think of other configurations that develop differently but in all cases, the idea of decomposing the force into two will tell us how they will move.
Figure 2.1: Representation of the collision between two hard disks. The active force, Fact, the forces
parallel to the line between the centers (depicted as a dotted line), Fn, and the forces orthogonal
to his same line, Ft, are represented by black arrows.
Thanks to the overdamped limit of our system’s mechanics, the collision between two particles become a contact that persists for some time. This will become crucial for us given that, as we have seen, studying the network of contacts can tell us quite a lot of the system’s collective properties. Now that we have clarified all these details, let us proceed by studying the interactions between particles through their equations of motion.
2.2
Theoretical formalism
2.2.1
Equations of motion
The theoretical formalism that we will develop in this section is based on a microscopic framework for non-Brownian suspensions called the Affine Solvent Model [33]. Let us consider then a system of N hard, frictionless active particles, suspended in a medium that provides them with overdamped dynamic. Then, for the ith particle, its equation of motion has the form of the expression 2.1.6, this is
Facti + Fconti − Vi = 0, (2.2.1)
where i = 1...N . Since our objective is to find scaling relations, we will not pay much attention to the units of our equations, that is why we have dropped out the ξ term that
appeared in the expression 2.1.6.
Let us take a closer look to the contact force term. Suppose that the ith particle is in
contact with several others. We will denote the magnitude of the repulsive force between the ithand the jthparticle as fij (note that, by Newton’s third law, fij = fji). Before proceeding,
let us introduce some notation: Let Ri be the coordinates of the ith particle’s center so the
directional distance between two particles is Rij = Rj− Ri and has a magnitude rij =
q R2ij. With these last two quantities in mind we can introduce the normalized direction between the two particles: nij = Rij/rij. Notice that nij = −nji. So, the total contact force exerted
on the ith particles can be written as
Fconti = X
j in contact with i
njifji.
This last sum can be expressed also as X
all pairs k, j in contact
(δji− δki)nkjfkj.
A clever observation made by Lerner et al. in [34] is that the deltas can be written as a derivative of the distance between the particles so we end up with the next expression:
Fconti = X all pairs k, j in contact ∂rkj ∂Ri fkj. (2.2.2)
Let us stop for a moment and look closely at this last equality. If we are able to know the forces between all the particles in contact and how does every distance between those pairs change with respect to a change in the position of the ith one, then we can know the total contact force the ith particle is experiencing. Following the formalism presented in [34,35] we
can build a vector |Fconti whose elements are the components of the vector Fcont
i for every
i ∈ N . In the same way, we can build a vector |f i whose elements are the fij magnitudes for
every j in contact with i without overcounting (if fab is listed, then fba is not). These new
elements allow us to rewrite equation 2.2.2 for every particle in the compact form
|Fconti = ST|f i . (2.2.3)
By their definition, we can see that |F i has dimension 2d (where d is the dimensionality of the system) and that the dimension of |f i is equal to the number of contacts, Nc. In
what follows, we will use capital letters for vectors of dimension dN and lowercase letters for vectors of dimension Nc. In this sense, ST is a linear operator from the space of contacts to
the space of particles so it has dimensions dN × Nc. Just to fix ideas, let’s take a look to this
small example:
Figure 2.2: System of five rigid disks with a particular set of contacts.
In figure 2.2 a particular arrange of five disks are shown. In this example the disk number 1 is contact with the disks 3 and 5 while disks 2 and 4 form a separated contact. So, the equations for the contact forces between the disks are
F1 = f31n31+ f51n15,
F3 = f13n13,
F5 = f15n15,
F2 = f42n42,
F4 = f24n24,
where, for the sake of smplicity, we have omitted the cont superscript. Since fij = fji and
F1,x F1,y F2,x F2,y F3,x F3,y F4,x F4,y F5,x F5,y = n31,x n51,x 0 n31,y n51,y 0 0 0 n42,x 0 0 n42,y n13,x 0 0 n13,y 0 0 0 0 n24,x 0 0 n24,x 0 n15,x 0 0 n15,y 0 f13 f15 f24 ,
This last expression is nothing more that the expanded form of equation 2.2.3. We can see that for this case the 10 × 3 matrix correspond to the linear transformation ST and our vector |f i has only three rows since our system has three contact. In contrast, the |F i vector has has ten rows since we have five disks that can move in two dimensions.
Returning to our main discussion, if we apply the same formalism to equation 2.2.1 and using equation 2.2.3 we get the next expression for the velocities of the disks:
|V i = |Facti + ST |f i . (2.2.4)
Intuitively, we can argue that the contact and active forces are not independent of each other since a particle will follow the trajectory imposed by its self-propulsion until it collides with some of its fellow companions. In order to obtain this relation, we need to understand first how do the pairwise distances rij of disks in contact are changing in time. In the most
general case, the distance between the ith and the jth particle will be a function of the
positions of the rest of the particles. Then ˙rij = X k ∂rij ∂Rk Vk= (Vj − Vi) · nij. (2.2.5)
If we rewrite equation 2.2.5 in our bra-ket notation, we will end up with a linear transfor-mation from the space of particles to the space of contacts through a similar operator than that one from equation 2.2.2 (the derivative in 2.2.2 is the same as that in 2.2.5! The only difference is the way the sum is taken). Thus,
| ˙ri = S |V i , (2.2.6)
where S has dimensions Nc × dN . Notice also that now we have an intuitive picture of
what does the operator S does: Given the velocities vector |V i, S returns the relative radial velocities of those particles in contact.
Now, given a group of particles in contact, we can build the contact network of the system. As we have discussed, because of the overdamped dynamics, there must exist an interval of time in which the number of contacts in this network remains unchanged either because no new contacts have been added or because the original ones have not broken yet. This implies that, in this interval of time, the pairwise distance between particles in contact is constant, i.e., | ˙ri = 0 and therefore S |V i = 0. Using this result on equation 2.2.4 we obtain that
|f i = −N−1S |Facti , (2.2.7)
where N ≡ SST.
So equation 2.2.4 becomes
|V i = Id − STN−1S |Facti . (2.2.8)
Once more, just to fix ideas, if we go back to our example system displayed in figure 2.2, its corresponding N matrix is
N = 2n2 13 n13· n15 0 n13· n15 2n215 0 0 0 2n2 24
We can observe that the diagonal terms of N tell us which disks are in contact: 1 and 3, 1 and 5, 2 and 4. At the same time, the off-diagonal terms contain the information of the relative orientation between contacts that share a node. This means that the entire contact network is completely described through the N matrix.
When writing the expression 2.2.7, we did not stop to think if something as N−1 actually exists. From our tiny example we might be inclined to think that the existence of such inverse is not a wild assumption but, as we will see, we need to be a little bit more careful with this matter.
2.2.2
The N and M operators
Throughout the first chapter of this work, we defined a solid as a system of particles whose relative distances do not change at least during the time scales of our observation. In fact, we used this idea to build a bridge between contact networks and amorphous solids. Nevertheless, there is another way to talk about solids (that we briefly mentioned at the end of the subsection 1.3.1 of chapter one, when describing the way Liu, Nagel, at al. measured
the elastic modulus of their systems): particles in contact whose configuration minimizes, at least locally, the interaction energy of the system. For soft particles systems, where the interactions are given by the overlap between particles, the interaction energy, U , has the form
U = X
all contacts
Φ, (2.2.9)
where Φ is the interaction potential. Notice how U = U (R1, · · · , RN). To simplify notation,
we will write U = U (R).
So, given a particular configuration of the suspended particles R0 = (R0,1, · · · , R0,N), the
system is a solid accordingly to our new definition if ∂U ∂R R0 = 0 (2.2.10)
and all the eigenvalues of
Mij ≡ ∂2U ∂Ri∂Rj R0 (2.2.11) are strictly positive except for the zeros that arise from the symmetries of U (remember that we are dealing with disordered arrangement that do not break any continuous symmetry when becoming a solid). The M operator is know as the dynamical matrix of the system. [36]
Let us suppose that our system has Nccontact and that we have enumerate them. Then
we can write the interaction energy 2.2.9 as
U = Nc X µ=1 Φ(rµ) = Nc X µ=1 1 2(rµ− Rµ) 2 , (2.2.12)
where µ is the index that labels the contact number, rµ is the distance between the centers
of the particles that form the µth contact, Rµ is the sum of their radii and we have choose
an energy scale such that = 1. Using this expression we can explicitly write the elements of the dynamical matrix as follows:
Mij = ∂2 ∂Ri∂Rj Nc X µ=1 1 2(rµ− Rµ) 2 = ∂ ∂Ri Nc X µ=1 (rµ− Rµ) ∂rµ ∂Rj = Nc X µ=1 ∂rµ ∂Ri ∂rµ ∂Rj + (rµ− Rµ) ∂2r µ ∂Ri∂Rj .
By setting rµ= Rµ, this is, in the isostatic case that applies for both soft and hard spheres,
the second term in the right-hand side of the last expression vanishes.
Using our new labelling for contacts and looking back to the expressions 2.2.2, 2.2.3, 2.2.5 and 2.2.6, we can see that we can express the S and ST operators as
Sµi= SiµT =
∂rµ
∂Ri
, (2.2.13)
µ ∈ {1 · · · Nc}, i ∈ {1 · · · dN }. So, we can write the dynamical matrix as
M = STS. (2.2.14)
This last expression suggests that there might be special relations between the matrices M and N since they are formed by the product of the same operators but in different order. In order to reveal this properties, let |ψωi be an eigenvector of M with eigenvalue ω2(≥ 0
according to our definition of solid) then
M |ψωi = ω2|ψωi
⇒ SM |ψωi = ω2S |ψωi
⇒ N S |ψωi = ω2S |ψωi ,
i.e., S |ψωi is an eigenvector of N with eigenvalue ω2 for arbitrary ω2. Thus, every eigenvalue
of M is an eigenvalue of N . We have to be a little bit careful with this assertion because if somehow S |ψωi = 0 then it cannot be an eigenvector and |ψωi would have zero as its
eigenvalue. Analogously, starting from N |φi = λ |φi we will conclude that every eigenvalue of N has to be an eigenvalue of M, therefore they are positive or zero, and that ST |φi is an eigenvector of M.
On the other hand, if the system is such that Nc< dN , M has dN − N c more eigenvector
that N . We know that every S |ψωi has to be an eigenvector of N (or a zero vector). Let’s
take one of those extra eigenvectors, say |ψνi, and suppose it maps to some eigenvector of N
under S. Since we are taking one of the extra eigenvectors there must exist a |ψωi outside
the set of extras that maps to the same eigenvector of N . This is, S |ψωi = S |ψνi
⇒ M |ψωi = M |ψνi
⇒ ω2|ψ
ωi − ν2|ψνi = 0
Therefore, all the dN − N c extra eigenvector of M have a zero eigenvalue. We will call then the zero-modes of M or floppy modes. Using the same reasoning, if Nc > dN then N has
Nc− dN zero-modes.
Now, given that the eigenvector of N form an orthonormal basis of the contact space and both operators share the same spectrum, if S |ψωi is an eigenvector of N with eigenvalue ω2
then it has to be equal to α |φωi, for some coefficient α and where N |φωi = ω2|φωi. Then
ω2|ψωi = αST |φωi
⇒ hψω| ω2|ψωi = α hψω| ST |φωi
⇒ α = ω. Therefore, we can conclude that
S |ψωi = ω |φωi (2.2.15)
and following a similar reasoning
ST |φωi = ω |ψωi . (2.2.16)
Returning to the initial question from which all this operators’ discussion arose, when can we talk about N−1? Well, we need that all its eigenvalues are non-zero and if this happens to occur, then we can guarantee that Nc ≤ dN . So, if Nc > dN , N has zero-modes and we
cannot invert it. Now, notice that we can rephrase this in term of the connectivity of the network: If N has no zero-modes then z ≤ 2d ≡ zc. And similar for the dynamical matrix:
If M has no zero-modes then z ≥ zc. This shows the consistency of our definitions of a solid!
If we have a suspension of particles collectively behaves like a solid, then M does not have zero modes, then z ≥ zc which implies in return that our suspended particles behave like a
solid. Therefore
Rigid phase ⇔ z > z
c⇔
M has no zero-modesN has zero-modes⇔ Rigid phase,
Floppy phase ⇔ z < z
c⇔
N has no zero-modesM has zero-modes⇔ Floppy phase
and the isostatic configuration, z = zc, correspond to a marginal solid where both N and M
Incidentally, what is the physical meaning of these eigenvectors? Before picking up the thread we left in the previous subsection, let’s answer this question. Suppose we have dense suspension of particles with positions |R0i such that the system is in its solid phase. So, for
small deviations of the equilibrium configuration, |Ri, we can write the system’s energy as U (|Ri) ' U (|R0i) +
1
2hR| M |Ri . Then, the forces on each particle, |F i = −∂RU , is
|F i = −M |Ri . (2.2.17)
Since the eigenvector of M form a basis for the vectors in the space of particles, we can choose our displacement to be proportional to a particular eigenvector, say |Ri = c |ψωi. So,
|F i = −ω2|Ri .
Since the eigenvalues of M are positive, the eigenvectors represent particular directions that, for displacements along them, will cause a set of response forces pointing back to the equilibrium configuration4 The zero-modes, on the other hand, are directions for which the
correspondent displacements can be done without any collective opposition.
2.2.3
A first approximation to active suspensions- physical
net-works
Before carrying the discussion of the N and M operators out, we had obtained an expres-sion for the contact forces, |f i, and the particles velocity, |V i, as functions of the particles’ active force and the geometry of the contact network, namely expressions 2.2.7 and 2.2.8. These two quantities, the contact forces and the particles’ velocities, will be right at the center of our theoretical predictions for this active system since their values strongly depend on the collective behaviour of the system and are physical observables.
4For systems where inertia plays a role in their mechanical behaviour we can write
| ¨Ri = −M |Ri .
After writing the displacements as a linear combination of the eigenvector, |Ri =P
ωcω|ψωi, we would find
that
cω(t) = Aωei(ωt+Bω),
this is, the small perturbation in the configuration causes a superposition of oscillation whose frequencies are given by the eigenvalues of the dynamical matrix. Since we are interested in overdamped dynamics, these oscillations are suppressed and will not appear in our description.
