Direct imaging of contacts and forces in colloidal gels

A longstanding aim for studies of soft solids is to understand the mechanisms by which they fail or yield, either due to internal stresses or imposed shear or other external fields.^1–6 Theoretical approaches can be limited since these materials are often far-from-equilibrium and their properties depend on the details of the preparation protocol and mechanical history. This is problematic because yielding processes are often heterogeneous, and in tackling this challenge, it is useful to think of localized irreversible (or plastic) rearrangement events driven by stress at the microscopic scale.^4,5,7 On the micro-scale, stress is a fluctuating quantity that is intrinsically linked to packing of particles, i.e., the local structure. However, on passing to the macro-scale, fluctuations in stress are no longer apparent, and the system obeys a constitutive relation that relates applied forces (stress) and material response (strain). To connect the microscopic and macroscopic pictures, thus, emerges as a major challenge, and in addressing this, measurement of the microscopic stress is crucial.

An important class of soft solids is colloidal gels,⁸ which are encountered in numerous foods,⁹ cosmetics, coatings, crop protection suspensions, and pharmaceutical formulations. In addition to colloidal systems, a wide range of materials also exhibit gelation, including proteins,^10,11 phase-demixing oxides,¹² and metallic glassformers.¹³ In gel-forming systems, failure under load can take the form of crack propagation.¹⁴ In colloidal gels, gravitational stresses can become important, leading the system to collapse under its own weight.^15,16 This last effect is an important determinant of the shelf-life of industrial products, such as agrichemicals. Among the most challenging aspects of gel collapse is that, prior to collapse, the elastic modulus of the gel increases, so it becomes harder before it fails.¹⁷ 

A promising way to address phenomena such as gel failure is to use particle-resolved studies where the coordinates of individual particles are tracked.^8,18 In soft amorphous solids, such as colloidal glasses, this technique has been used to image local re-arrangements of colloidal particles that may be precursors to large-scale failure.^19,20 In colloidal gels, particle resolved studies have revealed the rich nature of their local structure.^21–25 So far, while investigations of gel failure have related yielding to local crystallization²⁶ and ingenious combinations of rheological methods and simulation and scattering have revealed the role of local plasticity²⁷ and bulk two-point structure,²⁸ direct imaging of particle rearrangements has largely focused on colloidal glasses^8,19,20 rather than gels.²⁹ 

However, rather than imaging of particle coordinates, an alternative route to understanding gel failure is to consider the local stress, as one expects that regions of high stress are where failure may occur. Now, the local stress is manifested in the forces between the particles. While using particle-resolved studies to obtain the coordinates of the particles is useful,^8,19,20 it is clear that a major development would be some means to determine the force that each particle is under. This is, in principle, possible from the measurement of the coordinates and knowledge of the interaction potential between the particles. While the latter can be estimated to a good approximation,^8,30 the inevitable errors in the determination of particle positions and polydispersity of the particles mean that it is very hard to convert the distance between two particles into a potential energy or force. This means that this kind of measurement has hitherto only been possible in very special circumstances where the force varies slowly as a function of particle separation and the particles are far apart such that their positions, relative to the lengthscale over which the force varies, can be very accurately determined.³¹ Thus, from coordinate data only, stress correlations are typically inferred indirectly.³² Computer simulations, of course, also provide access to coordinate data and the forces between particles.^33–35 Often, data are obtained similar to those of particle-resolved experiments,³⁶ particularly if hydrodynamic interactions are included.^23,37,38 However, while some aging effects are accessible³⁴ most phenomena pertinent to failure in colloidal gels,²⁸ particularly delayed collapse,¹⁷ the timescales and system sizes lie beyond those accessible to particulate simulation.

In granular systems, forces have been characterized for particles with diameters of at least 10 μm^39–43 (and often cm^44,45), and potential has been demonstrated for a scaled up version of a popular colloidal model system.⁴⁶ Unlike athermal granular systems, the thermal motion exhibited by colloids leads to a multitude of new phenomena, such as the emergence of long-lived networks in the colloidal gels that we are interested in here. Identifying contacts and forces in colloidal systems is challenging due to sub-resolution length scales relevant to obtain forces between colloids. Here, we take a first step to address this challenge by investigating interparticle contacts and forces in colloidal gels via high-resolution optical microscopy. Here, we focus on the methodology to identify local stress-bearing regions and consider the static properties of a system prior to failure. In particular, we use an emulsion system with a solvatochromic dye, which is sensitive to the compressive forces between droplets.⁴⁰ In this way, we obtain force contacts between droplets and measure the local stress. We correlate these via structural quantities and compare with computer simulations. We find that droplets in local structures associated with rigidity are more likely to be under higher pressure. Our work, thus, presents a method that can be used to identify potential failure points at the particle level in amorphous colloidal solids.

This paper is organized as follows: In Sec. II, we explain our experimental protocol to identify contacts between particles and proceed to describe how we may determine a local measure of compressive stress and how these are connected to form force chains. We detail the computer simulation methodology that we use to compare with our experimental results. In Sec. III, we compare the experimental results for the number of contacts per droplet and their coordination with simulation. We go on to consider the distribution of compressive forces. We, then, investigate the length of force chains and compare these with computer simulations. Finally, we consider correlations between some of the quantities we have investigated. We discuss our findings in Sec. IV.

Colloidal polydimethylsiloxane (PDMS) emulsion droplets were synthesized following Elbers et al.⁴⁷ Sodium dodecylbenzenesulfonate (SDBS) surfactant (2 mM) and potassium chloride salt solution (20 mM) were added in order to stabilize PDMS emulsions and screen charges on droplet surfaces, respectively. The solvatochromic dye Nile Red was employed to fluorescently label PDMS emulsions. Glycerol was then added to obtain a refractive index matched emulsion with a weight ratio of water to glycerol around 51%:49%. The droplets have a mean diameter of d = 3.2 μm, which is determined from the first peak of the radial distribution function obtained from particle tracking. The Brownian time to diffuse a diameter is

where η is the solvent viscosity and k_BT is the thermal energy.

The non-absorbing polymer utilized to induce depletion attraction is hydroxyethyl cellulose (Natrosol HEC 250 G Ashland–Aqualon) with molecular weight 3 × 10⁵ g mol⁻¹. Colloid–polymer mixtures are prepared by adding stock solutions of HEC polymers (10 gl⁻¹) to concentrated emulsions with volume fractions around random close packing that we take to be ϕ_rcp ≈ 0.64. All colloid–polymer mixture samples are observed under a confocal microscope at about 16 τ_B after loading the sample cell. Our system is not density matched between droplets and the solvent. In particular, colloidal systems, including depletion gels, are known to undergo batch settling (or, here, creaming) such that the local volume fraction in the bulk of the system is largely unaffected at short times.^16,48–50 We ensure that we analyze data from the bulk of the sample where little change in the volume fraction due to sedimentation is observed. At much longer timescales than those we consider here, a dense sediment is found.^16,51 Further sample details are listed in the  Appendix.

We used a Leica SP8 confocal microscope with a continuously tunable white light laser. We use two-channel imaging with excitations 514 and 580 nm that correspond approximately to the absorption peaks of Nile Red in non-aqueous and aqueous environments, respectively. Nile Red emits at differing frequencies in non-aqueous and aqueous environments with peaks of 545 and 645 nm, respectively. We exclude particles whose centers are closer than one diameter to the edge of the image to mitigate the boundary effects.

To obtain information on the interdroplet contacts and forces, we developed a method to mitigate the challenges resulting from the limited spatial resolution of the microscope. Previous work^39,40 considered much larger droplets, but, here, we must contend with the challenge to optical microscopy posed by rather smaller colloidal droplets. The size of the contacts, particularly their separation from one another, is comparable to the resolution of the microscope. We proceed by tracking the droplet coordinates⁵² in the droplet images [green channel in Fig. 1(b)]. Since our system is reasonably monodisperse (polydispersity $≈8$%), we know that the contacts should be approximately equidistant between the centers of two neighboring droplets. The set of midpoints between neighboring droplets, thus, gives a trial set of candidate force-bearing contacts. Each of these is populated with a sphere, which we term a blob. From this, we determine the magnitude of the force in the image by comparing with the number of pixels within this spherical volume and their intensity in the contact image [red pixels in Fig. 1(c)].

To obtain a measure of the force, we apply a threshold to the contact image. The contacts are identified on the basis of the number of pixels in the contact image that correspond to the “blobs” that are potential contacts. Our analysis gives a measure of the relative magnitude of the compressive force at contact points on each droplet. We compare our results to those which are approximately matched to the experimental system. Further details of our analysis are given in the  Appendix.

The interaction between emulsion droplets is complex and depends on the local geometry.^53,54 Here, we seek an estimate of the energy scales involved. Now, the surface tension γ = 9.2 m Nm⁻¹,³⁹ which amounts to an energetic cost comparable to the thermal energy for a microscopic change in the surface area of the droplet. Therefore, in the case of our mesoscopic droplets, we expect deformations to be very small. For such small deformations, we assume that two interacting droplets are deformed such that the surface in contact between them is a circle and determine the change in the surface area with respect to two undeformed droplets of the same total volume. To leading order, the interaction energy

for r ≤ d. Here, β = 1/k_BT. For our parameters, we expect that very small deformations around 0.1% are sufficient to result in an interaction energy of many times that of the thermal energy (Fig. 2). Our droplets, therefore, approximate closely hard spheres.⁵⁵ In fact, in comparison to measurements of the softness of other experimental systems that approximate hard spheres, such as sterically stabilized poly-methyl methacrylate (PMMA),⁵⁶ Eq. (2) indicates that these emulsion droplets are even closer to hard spheres than the well-known PMMA system. Note that some other emulsion systems exhibit rather lower surface tension and, therefore, more deformation is found.^41,57

The polymer size is much smaller than that of the droplets, such that our system is toward the “sticky sphere” limit of short-ranged attraction strength. We presume that the effective attractions between the droplets are of the Asakura–Oosawa form,

Here, q = 2R_g/d is the polymer–colloid size ratio and $zpr$ is the polymer fugacity in a reservoir in thermodynamic equilibrium with the colloid–polymer mixture, which we assume to be equal to the polymer number density in the reservoir, as would be the case for ideal polymers. Here, R_g is the radius of gyration of the polymer. We neglect the contributions from electrostatics due to the Debye screening length that we estimate as 2 nm, which is much smaller than the range of depletion attraction. Furthermore, using Derjaguin, Landau, Verwey and Overbeek (DLVO) theory, we arrive at a contact potential due to electrostatic interaction between two droplets less than k_BT.

To estimate the interactions between the droplets, we assume that the attractive interaction remains for small compressions of the droplets r < d,

The interaction potential is plotted in Fig. 2, where it is seen that the AO attraction is swiftly overwhelmed by the strong repulsion u_drop(r).

To proceed, we require the polymer radius of gyration and we estimate it from the gelation boundary. The phase diagram of our system is given in Fig. 9 in the polymer reservoir representation.⁵⁸ We map our experimental values of the polymer concentration to the reservoir representation using Widom particle insertion.⁵⁸ The polymer reservoir concentration corresponding to gelation $cpr, gel$ is then 0.71 ± 0.1 gl⁻¹. We express the polymer concentration as a ratio of this value. The polymer radius of gyration is then fixed by requiring that the reduced second virial coefficient at the gelation boundary $B2*=−3/2$⁠.⁵⁹ While this holds for criticality, in fact, for gels undergoing arrested spinodal decomposition (as is the case for colloid–polymer mixtures⁸), such short-ranged attractions lead to a very flat liquid–gas spinodal^36,60 such that the polymer concentration for gelation varies very little across a wide range of colloid volume fractions. In this way, we arrive at a polymer with a radius of gyration R_g = 21.2 nm and a polymer–colloid size ratio q = 0.013. This is close to the value quoted for HEC 250 G in Ref. 61. The resulting effective droplet–droplet interaction potential is shown in Fig. 2. We are interested in the compressive forces between the droplets. We, thus, interpret these as −d[βu_exp(r)]/dr for r ≤ d.

We now outline our method to obtain a measure of the local stress. Consider a reference particle p₀, for example, with three neighbors p₁, p₂, and p₃ that touch through contacts c¹, c², and c³, as shown in Fig. 3(a). The compressive force from the neighboring particle p₁ is f¹, with magnitude f¹, which is determined from the size of contact c¹. Our single-particle stress measurement, σ, is calculated by summing stress contributions from each neighbor on all element axes, as indicated in Eq. (4). This yields a 3 × 3 matrix σ whose elements are denoted by σ_ij, where i, j label Cartesian components. Similarly, $fic$ is the ith Cartesian component of the force on particle c,

where n^c is the number of contacts of p₀. Dividing this quantity by a suitable volume gives the Cauchy stress tensor, but here assigning the volume presents a challenge. As Fig. 1 shows, gels are heterogeneous materials. Thus, partitioning space according to a Voronoi decomposition leads to unphysically large separations. On the other hand, using the droplet volume does not fill space, as the volume fraction ϕ_c < 1. Here, we consider normalized quantities in reduced units where the mean particle diameter is set to unity. We shall, therefore, refer to σ as a reduced stress tensor, noting that we apply it at the single-particle level. For each particle, we obtain σ by analogy to the stress tensor, and diagonalization generates three eigenvalues and eigenvectors, which represent principal stresses and principal directions, respectively. After diagonalization, the sum of all principal stresses

Note that the quantity −tr(σ) is analogous to the local pressure.

Here, to identify force chains, we consider a quasilinear assembly of at least three particles where stress is concentrated.⁶² Based on this definition, a method was developed by Peters et al.,⁶² which we illustrate schematically in Fig. 3(b). If the minor principal stress of a particle is larger than the average compressive stress in the system, these particles are candidates for force chains. After selecting these particles with large stresses, we require that the particles with concentrated stresses should be quasilinear, allowing only small amounts of rotation. Given a reference particle 1, from its center, we define a region that deviates an angle of ±θ from the direction of stress σ_A. Here, we set θ = π/4, and we also require that the direction of symmetry on the second particle to be within θ.

As noted above, to verify and calibrate our experimental data, we perform Brownian dynamics computer simulations. We use point particles interacting via a spherically symmetric potential with Hertzian repulsive forces and a short-ranged attractive term, which we shift and truncate at a range r_cut. The repulsive Hertzian contribution to the potential is

while the attractive term is

so that the full potential is

The resulting u_sim(r) is plotted in Fig. 2.

To accurately model the short ranged, highly repulsive interaction between the droplets of the experiment is highly challenging for conventional computer simulations. While novel methods have been developed for Monte Carlo simulations,⁶³ here we are interested in dynamical behavior. We, therefore, set βA = 1000, δ = 0.02, and r_cut = 1.3 d, which results in a very short range attraction with a soft core (see Fig. 2). The interaction strength ɛ and the number density of the system ρ characterize the state points. Using the Barker–Henderson effective hard sphere diameter,

we map number densities to effective volume fractions $ϕeff=πdeff3ρ/6$⁠. We presume that the effective volume fraction corresponds to the absolute droplet volume fraction in the experiments.

Simulations are performed in the isothermal–isochoric ensemble (NVT) solving the Langevin dynamics

for particles of equal mass m in the presence of a zero-mean unit-variance random force ξ. For this purpose, we employ a suitably modified version of the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) molecular dynamics package.

We use a velocity-Verlet integrator with timestep dt = 0.001τ₀ with $τ0=mβd$⁠. Fluctuations of the temperature are allowed to damp on a relatively short timescale of $τd=100dt=0.1mβd$⁠, a setup for which results are similar to the overdamped limit.⁵⁰ This timescale also sets the Brownian time $τB=γd2/24kBT=τ02/24τd$⁠, which allows us to compare the numerical results with the experiments via Eq. (1).

As non-equilibrium systems, gels coarsen over time.^8,64 Now, there is a significant difference in this coarsening process between experiments and Brownian dynamics simulations due to hydrodynamic interactions in the former that are not fully accounted for in the latter.^23,37,38,65 Therefore, even if a precise matching of timescales were to be carried out, in fact, one would still expect considerable differences between experiment and simulation, as has been found previously.²³ For our purposes, then, we select a point in the time-evolution of the simulations of 12τ_B, in which we find that a number of time-dependent properties are comparable to those in the experiments (see Sec. III). However, we emphasize that Brownian dynamics simulations provide a robust description of the phenomenology of colloidal gelation even if the comparison between experiment and simulation is not fully quantitative.^23,37,66 In particular, gels simulated in this manner are composed of percolating networks that exhibit a yield stress when probed on timescales much less than those associated with aging,³⁴ consistent with experiments.¹⁷ 

To map the state point between simulation and experiment, we set the volume fraction ϕ_eff = 0.2, and to compare the interaction strength with the experiments, we scale by the value corresponding to gelation ɛ^gel. Following the experiments, we set ɛ^gel as that at which the reduced second virial coefficient $B2*=−3/2$⁠.^36,59 In addition to the approximations noted above, given that the interaction potential in the simulations is also rather softer than that of the experiments, we regard the comparison between our two approaches to be semi-quantitative, rather than the simulations being an accurate reproduction of the experiments.

We organize this section as follows. We begin by discussing the structural properties, the number of neighbors, and the number of contacts. We then move on to consider the forces between droplets inferred from the contacts, leading to quantities such as the reduced stress tensor. We then consider force chains. Throughout, we compare our experimental results with those from computer simulations.

We consider, schematically, the imaging methodology and method for interparticle force extraction in Fig. 1(a). Representative data of each fluorescent channel is shown in Figs. 1(b) and 1(c), and their combination is shown in Fig. 1(d). We render the droplets actual size and, following the identification of contact analysis outlined in Sec. II C and described in more detail in the  Appendix, the contacts are shown as pink sticks in Fig. 1(e). This constitutes our basic data. Having demonstrated the principles of our method, we consider the quantities of interest.

We proceed to show renderings of the properties of particular interest in Fig. 4 for a polymer concentration of $cp/cpr, gel=1.5$⁠. Other gel state points appear similar. We show the number of contacts for each particle n_c, which appears to be rather heterogeneous throughout the system. Before considering the other quantities, we move to a quantitative discussion of the coordination and the number of contacts around Fig. 5.

In Fig. 5(a), we show the distribution of the number of neighbors, Z, which are defined as being closer than 1.2d, which is close to the first minimum of the radial distribution function g(r). The number of neighbors requires knowledge of only the droplet coordinates, and thus, a comparison can be made to other work with particle-resolved studies, and indeed, similar behavior can be found, e.g., in Fig. 2 of Ref. 67. The number of contacts n_c for the same data is shown in Fig. 5(b). This has a smaller value than the number of neighbors. The distributions of neighbors and contacts in our simulations show very similar behavior, as shown in Figs. 5(c) and 5(d). In the simulations, the ɛ/ɛ^gel = 1.3 state point has fewer neighbors and contacts than the others we have shown. However, it is worth noting that this is rather closer to the gelation boundary than the others (the next closest being the experimental state point at $cp/cpr, gel=1.5$⁠), which could account for the difference.

We can estimate the relative compressive forces on each droplet. At the level of our analysis, we determine the number of pixels in the contact image within a “blob” (see Fig. 11). We take the sum of these contact pixels as a measure of the contact volume v_c, which we plot in Fig. 6. For granular systems, the force has been identified with the contact area,^39,43 which should scale as $vc2/3$⁠. We, therefore, plot the distribution of $vc2/3$ in Fig. 6(b). This is much sharper than the distribution of volumes. In our simulations, we have direct access to the compressive forces, and we plot these in Fig. 6(c). The distribution from the simulations is rather broader. Indeed, except for the smallest forces, the experimental data are roughly compatible with a Gaussian distribution. However, there is some evidence in the simulations for an exponential decay [black dashed line in Fig. 6(c)]. As noted above, although the time-evolution in the experiments and simulations differs, the structural quantities in Fig. 5 are rather similar across the two systems. Therefore, it is possible that the difference in the interaction potential between the two (Fig. 2) may underlie the difference between the force distributions that we obtain. Since the measured contact volumes also depend on the particle dynamics and the imaging process, it is likely that they do not respond to very fast force fluctuations. For example, to image a “contact” whose size is 0.1σ corresponds to 2.2s ≈ 0.11τ_B, during which it may diffuse a distance $∼0.34σ$⁠. This would mean that the force inferred from the contact size corresponds to a time-averaged version of the interparticle force. The time-averaging will act to suppress force fluctuations. This suppression is absent from simulations, where the (instantaneous) microscopic force is measured directly.

Identification of the forces associated with each contact allows us to investigate the reduced stress tensor σ [Eq. (4)]. In Fig. 4(b), we show the local anisotropy, which is the difference between the largest and smallest eigenvalues of σ. Although it may appear from visual inspection that this quantity has some spatial correlation, we have investigated such correlations and find that these are indistinguishable from the (short-ranged) density correlations expressed via the radial distribution function. We, then, plot the negative trace of the reduced stress tensor −tr(σ) that is analogous to the local pressure in Fig. 4(c). Like the number of contacts [Fig. 4(a)], this is rather heterogeneous. The trace is correlated with the number of neighbors, with higher pressure corresponding to a larger number of neighbors [Fig. 7(a)]. Here, the Pearson correlation coefficient is 0.709.

Colloidal gels have been subjected to structural analysis, particularly clusters that are local energy minima have been identified with rigidity.^8,68 Now, the local structure changes over time, leading to a larger and more complex local structure,^36,69 and at early stages like the gels of interest, here, the dominant local structure is the tetrahedron.²³ It is possible to classify the particles according to the number of local structures in which they participate, which can reveal the degree of local ordering.^70,71 Here, therefore, we count the number of tetrahedra in which each particle participates, as shown by the rendering in Fig. 4(d). Visual inspection suggests that there is some correlation between the number of tetrahedra the particles participate in, and the trace −tr(σ) [Fig. 4(c)]. This is indeed the case [Fig. 7(b)] with the correlation coefficient being 0.455.

We implement the measurement of force chains outlined in Sec. II F. In this way, we obtain the distribution of force chain lengths P(l) in our system. We emphasize that there is no reason a priori to expect that these would span the system, as is the case for granular materials in compression or under shear.⁴⁴ In fact, the majority of particles are found in force chains of a single particle. Longer force chains are rendered in Fig. 4(e). When we plot the distribution of force chain lengths l in Fig. 8, we find that in both simulation and experiment, the effect of interaction strength is weak. The force chains in experiment are rather longer. Our data are compatible with an exponential distribution, with a decay length of 3 and 3/4 particles in experiment and simulation, respectively.

Note that here we may cut some force chains at the image boundaries. Although we neglect contributions closer than a diameter d to the boundary, it is hard to remove possible boundary effects from the force chain distribution for images or the size that we acquire here. However, we may observe that in Fig. 4(e), the force chains are rather smaller than the imaging volume and, thus, we expect any boundary effects to be reasonably small, and in any case, these will tend to reduce the apparent chain length, so such boundary effects are unlikely to be the cause of the difference between the experiments and simulation that we see. Given that hydrodynamic interactions are associated with more linear structures,^23,37,38 it is tempting to suppose that these are part of the reason for the longer chains that we find in the experiments.

We have characterized the interparticle contacts in a colloidal gel of emulsion droplets. We have further investigated compressive forces between droplets related to these contacts and have semi-quantitatively benchmarked our results against computer simulations. We have fewer contacts, and particles with large numbers of contacts are not strongly correlated in space.

Turning toward the forces, we infer these from the number of pixels in the blobs in the contact image that measures the spatial distribution of the solvatochromic dye. The change in the droplet surface area due to deformation of the mesoscopic emulsion droplets incurs a high energetic cost, as the surface tension is of the order of the thermal energy for a microscopic (molecular) change in the surface area. Under the depletion forces due to the polymer, we, therefore, expect very weak deformation of the droplets. We believe that the contact volume inferred from the images of the solvatochromic dye is larger than the true contact area. Further investigations in this direction are clearly desirable, perhaps using systems with lower surface tension whose droplets would be deformed rather more,⁴¹ or so-called skinny emulsions.⁵⁷ Nevertheless, the normalized force distributions that we obtain are comparable to our simulations. The somewhat broader distribution in the simulations might be related to the softer interaction potential that we have used. The width of the interaction potential could be (somewhat) narrowed toward that assumed for the experiments to investigate if this is the cause of the difference.

We have obtained a measure for the local pressure from the reduced stress. This quantity is evaluated at the level of individual droplets and varies significantly from droplet to droplet. We attribute this variation to the heterogeneous structure of the gel and the diversity of particles’ local environments. As with the number of contacts, the stress does not show strong spatial corrections. However, it is quite well correlated with the number of neighbors and also with the local structure, as expressed by the number of tetrahedra that a droplet participates in.

Now, the system that we have considered is different from the granular materials previously considered in a number of ways. First, it is thermal; therefore, force balance only holds on average within the amorphous solid because the droplets are constantly in motion. Second, it is a gel undergoing (arrested) spinodal decomposition. This has various consequences: the network is globally under tensile, rather than compressive forces. However, our method measures compressive forces. Third, although the system is not density matched, there is no significant external field (we find no significant anisotropy in our measurements of the local stress). Indeed, in the absence of an external field, granular systems are not expected to exhibit system-spanning force chains, consistent with the fact that the force chains that we find in this thermal system are rather shorter than those encountered in granular systems under load.⁴⁴ 

Here, we have focused on a gel-forming system with effective attractive interactions driven by polymer-induced depletion. In the future, it would be interesting to apply these methods to amorphous colloidal solids formed of repulsive interactions, such as colloidal glasses. Eventually, this may enable a link to be made to earlier work on granular systems by bridging the Gardner transition.^72–74 Work in this direction is in progress.

Again, we encounter similar behavior in simulation, although the force chains are somewhat longer in our experiments, which may be related to hydrodynamic interactions in the latter that are largely neglected in the former. The effect of HI would thus be interesting to probe in the future. While granular systems with attractive interactions have been investigated, there the focus lay more toward the contacts.⁴² Given the much higher volume fraction investigated in that work, direct comparison is hard, not to mention the differences between the thermal and athermal nature of the systems. It would, nevertheless, be most attractive to explore the force chain distribution in attractive jammed materials, such as granular gels.⁷⁵ Granular systems with repulsive interactions are by their nature found at high volume fractions, and force chains typically percolate to form force networks. Nevertheless, there is some evidence for an exponential distribution in community sizes⁷⁶ in force networks, the same scaling as we find with the much shorter linear chains we measure here.

Thus, we present a colloidal model system for characterizing contacts and interdroplet forces. By considering perturbation such as shear, this system may be used to obtain a knowledge of local stress that may prove useful in understanding failure in soft solids, such as colloidal glasses and gels.

We thank Paul Bartlett, Olivier Dauchot, Jasna Brujić, and Jens Eggers for helpful discussions and Yushi Yang his valiant assistance with the TCC analysis. C.P.R. acknowledges the Royal Society for support. J.D., F.T., and C.P.R. acknowledge the European Research Council (ERC Consolidator Grant NANOPRS) for support (Project No. 617266). J.D. acknowledges Bayer AG for support. R.L.J. and C.P.R. acknowledge the EPSRC for support via Grant Nos. EP/T031247/1 and EP/T031077/1, respectively. EPSRC Grant No. code EP/H022333/1 is acknowledged for provision of the confocal microscope used in this work.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

To image the system with a confocal microscope, we employed two excitation lasers with wavelengths of 514 and 580 nm and two HyD hybrid detectors with detection wavelengths of 520–575 and 585–640 nm. These two lasers and detectors are applied to detect fluorescent signals from the bulk of the PDMS oil droplets and the contacts between the droplets, respectively. We refer to the images generated by these two channels as the droplet image and the contact image, respectively. We equalized the image histograms as a function of depth to compensate for any attenuation due to an imperfect refractive index matching between emulsion droplets and the solvent. To reduce the noise of the captured images, we applied line averaging of 32 to each frame and deconvolved the images with the Huygens software. Droplet centers were detected in the droplet image using the colloids tracking package.⁵² 

Here, the smaller lengthscale with respect to previous work with much larger droplets³⁹ necessitates a method to segregate connected contacts and determine the centers and sizes of contacts. Having obtained the droplet centers, we proceed by processing the contact images. We use two Otsu thresholds (which is a threshold based on weighted variances of intensities of pixels corresponding to features and background⁷⁷). As schematically shown in Fig. 10, we apply two Otsu thresholds to the contact images. The first distinguishes droplets (with contacts) from the solvent background. The second separates contacts (foreground) from bulk droplets (background).

To remove any contacts erroneously identified due to the residual intensity in the interior of the droplets, we use a Sobel filter to enhance the droplet edges. However, applying the Sobel filter directly to the droplet image means that the edges of each particle are not always well defined because some droplets are in contact with one another. Therefore, instead, we generate an image from the particle coordinates we have determined and apply the Sobel filter to each particle.

After thresholding images such as Fig. 1(c), we find that the contacts are frequently merged. In order to separate such connected contacts, our strategy is to add a spatial boundary to each contact. The first step is to find the weighted middle points between a reference particle and its neighbors, which are possible locations of contact centers. To determine the weighted middle points m_ij between two neighboring droplets i, j, first, m_ij needs to be located on the line connecting the centers of droplets i and j and the distances between m_ij and the two neighboring droplets b_i and b_j are proportional to particle radii a_ij, i.e., b_i/b_j = a_i/a_j. Therefore, a binary mask of the same size as the contact image is built, where the positions of weighted middle points m_ij have a value of 1 while the rest of the mask is 0.

Based on the centers of middle points, spherical blobs were created by dilating a binary kernel in three dimensions. The blobs were constructed as large as possible but without overlapping with each other. The purpose of building blobs is to contain true contacts as much as possible and build an upper boundary for the contacts to separate them from each other if they are overlapping after thresholding. Because blobs are created in between neighboring droplets (which are not necessarily in contact), so the number of blobs generated is greater than the number of true contacts.

After the initial placement of blobs [Fig. 11(a)], some are connected when we try to maximize their size as shown in Fig. 11(b). By looking at the distribution of blob volumes, it is clear that the connected blobs have noticeably larger volumes than the isolated blobs, and the binary mask with all blobs was separated into two masks: a well separated blob mask [Fig. 11(b), blob “1”] and a connected blob mask [Fig. 11(b), blob “2, 3, 4”], In the mask with connected blobs, we eroded the mask in order to separate these blobs [Fig. 11(c)]. Next, an eroded mask [Fig. 11(c)] and a non-connected blob mask [Fig. 11(b), blob “1”] were combined into a final binary mask. This mask effectively sets bounds for contacts and can be used to segregate the connected contacts. At this point, we have identified the contacts. However, we now seek to determine their size, from which we can infer the force related to each contact.

Three masks are generated in order to correctly detect the positions and sizes of contacts. The first mask [Fig. 11(a)] is the binary mask of spheres that are placed between droplets. This mask segregates some contacts that are connected after the Otsu thresholding of the contact image. It is possible that some pixels, which are located in the middle of droplets, remain after thresholding. Therefore, a second mask, which contains the edges of all droplets, is desired, in order to set constraints to contact positions. This means contacts can only be located at the edges of droplets but not inside droplets. The third mask is a thresholded contact image, which is obtained by applying the Otsu threshold to the contact image (Fig. 10). By convolving these three masks, the remaining pixels are the contacts between droplets. Each contact is then labeled with an index, and counting the number of pixels in each contact, then, gives the volume of the contact. The contact center is determined by finding the geometrical center or maximum intensity pixel in the contact.

After particle and contact tracking, both coordinates and sizes are obtained. The coordinates of droplets and contacts are p_i and c_j, respectively. The distances between each particle and contact are computed and stored in an i × j 2D matrix s_ij,

where N_p and N_c are the number of droplets and contacts, respectively. For a contact c_j, the closest two droplets p_a and p_b are detected by searching for the first two minimum values s_aj and s_bj in s_ij. These two droplets are then in contact through c_j. For each contact, we find two neighbor droplets, and, in turn, we can determine neighbor contacts for each particle, which gives the number of contacts n_c.