Spanning the scales of granular materials through microscopic force imaging

If you walk on sand, it supports your weight. How do the disordered forces between particles in sand organize, to keep you from sinking? This simple question is surprisingly difficult to answer experimentally: measuring forces in three dimensions, between deeply buried grains, is challenging. Here we describe experiments in which we have succeeded in measuring forces inside a granular packing subject to controlled deformations. We connect the measured micro-scale forces to the macro-scale packing force response with an averaging, mean field calculation. This calculation explains how the combination of packing structure and contact deformations produce the observed nontrivial mechanical response of the packing, revealing a surprising microscopic particle deformation enhancement mechanism.

T he way that a disordered packing of macroscopic particles responds to mechanical deformations, such as compression, was investigated as early as the seventeenth century, when Stephen Hales studied the expansion of dried peas submerged in water 1 . The motivation for such studies is clear. Packings of particles surround us: coffee beans and rice; soils, embankments, many industrial processes and geophysical processes from earthquakes to landslides involve granular materials 2 , hence, the gamut of recent investigations [3][4][5] . Part of the relevant physics also applies to emulsion and foams, which consist of highly deformable frictionless particles, and to colloidal systems 6,7 , where thermal agitation is relevant. A major challenge for granular materials is relating microscopic contact and force details to global system response, such as applied stresses or strains. The link from micro-to macro-scale behaviour informs the mathematical link, known as a constitutive relation, between macroscopic strain or stress and macroscopic response. These relations appear routinely for conventional materials such as Newtonian fluids (shear stress proportional to shear strain rate) or linear elastic solids (stress proportional to strain). For granular materials, there are continuum models containing empirical constitutive relations 8 , which are widely used in the soil mechanics communities. However, they suffer from problems: the connection between grain-scale and macro-scale behaviour is still an open question, and these models are known to contain complex mathematical instabilities 9 . The lack of a wellestablished constitutive relation for granular materials is highly limiting; it is essential in applications where a continuum description is the only computationally feasible approach.
To develop models such as constitutive relations, it is essential to have direct access to microscopic experimental information that shows how granular materials respond to applied stresses or strains. The challenge is to have access to all microscopic quantities inside a granular packing, including at every contact between grains, not just for a single packing under fixed boundary conditions, but during the entire evolution of a packing under realistic deformations. Experimental probes of granular microscopic properties that provide structure and contact forces in three-dimensional (3D) packings are limited. Photoelastic methods provide this type of information in two dimensions 10 , but are hard to implement in 3D. Pioneering measurements by Brujić et al. 11 and by Zhou et al. 12 for emulsions yielded the distribution of forces at contacts, P(f), relative to the average force and information on contact networks. Tomographic techniques (such as X-ray microscopic computed tomography (microCT) 13,14 and refractive index matching 15,16 ) have been successfully implemented to measure structure and particle contacts. Extracting forces was done in microCT data 17,18 , by imposing a global force balance constraint. This constraint infers local contact forces from (measured or simulated) average particle stresses and boundary conditions, even when contact forces cannot be easily extracted. An alternative approach would be to measure forces directly, at each individual contact. Both the methods could be combined to improve accuracy. Importantly, all existing methods are limited in their ability to yield microscopic structure, particularly contact forces, over many closely spaced macroscopic state changes and with a precise control of applied strains.
In this work, we use refractive index matching tomography to provide full access to microstructure, for packings of deformable hydrogel particles that are compressed and decompressed uniaxially. The crucial aspect of our work is the unique combined ability to measure individual contact forces in vectorial detail, while straining the sample in small increments, enabling us to track the system-scale stress tensor over many small strain steps. This feature gives access to the complete micro-macro range of mechanical details of the packing, including ingredients for constitutive modelling.

Results
Description of the experiment. In a typical experiment, sketched in Fig. 1a, 514 hydrogel particles that have been saturated with fluorescent dye are contained in a Plexiglas box. The particles are roughly, but not perfectly, spherical, with a mean diameter of 2.1 cm. The box is filled with a solution of water and polyvinylpyrrolidone, such that the index of refraction of the particles is well matched to the solution. This allows optical access to a vertical laser sheet that is scanned horizontally. A camera, whose image plane is parallel to the laser sheet, records a series of images as the laser sheet is swept. To reconstruct particles and contacts, we have developed dedicated tomography algorithms that yield the force (including direction) for each individual particle contact. We discuss this method further, along with additional experimental information, in the Methods section. The particles are contained in a box with five rigid transparent walls. The sixth and topmost wall is a porous piston, which moves in the vertical direction, providing compression/decompression. A force gauge, which is in-line with the piston, measures the vertical force acting on the top layer of particles, and hence the pressure at the top boundary. The piston location and the force gauge give us the macroscopic strain and stress imposed on the system. For each scanned sequence of images, we infer the particle shapes, locations and all the contact vector forces, which give us the full microstructure (Fig. 1b). For the present system, friction is negligible and forces are nearly normal to contacts. However, the technique used here can be generalized to contact forces with friction. The density of the particles is also nearly matched to the fluid, so that the particles experience an effective gravity of about 0.01 g (g is the acceleration of gravity). The experimental protocol consists of a series of 20 uniaxial compressions, each followed by a corresponding expansion to the original boundary configuration. Each cycle consists of a compression phase imposing a strain up to 13.4% of the initial height, followed by a decompression phase returning the top plate to its original position. A full cycle is carried out in 60 quasi-static steps of 1 mm each, but the top plate does not touch the grains in the first five and last five steps. After each quasi-static step, we carry out a complete volume scan.
Microstructure and 3D force information. To obtain the contact forces from the reconstructed scanned images, we implement the following steps, which we discuss further in the Methods section. First, we identify the boundary of each particle, which is not perfectly spherical, including regions of contact with other particles. Then, from the areas of contact, we use linear elasticity, and the independently measured Young modulus of the particles (EE23 kPa), to infer the contact normal forces. The particles are nearly incompressible (volume change from uncompressed to fully compressed is o1%) so we assume that they have a Poisson ratio of E0.5. Our best estimate of the hydrogel friction coefficient is mE0.03. Hence, for these experiments, the tangential (frictional) forces are at most 3% of the corresponding normal forces and lie below experimental resolution (the minimum average contact force is hfi ¼ 10 À 2 N, see Methods). We then compute the coarse-grained continuum stress tensor 19 . We integrate the normal contribution of the stress tensor over the upper boundary to obtain the force on the top plate. We emphasize that this measurement uses only information from the micro-scale and the independently measured particle Young's modulus. The force data resulting from the stress integration and from the independent in-line force gauge are nearly identical, as we show in Fig. 2, for the last 15 of the 20 uniaxial compression/ decompression cycles. Here red and blue distinguish, respectively, the macroscopic force measurement of the gauge from the microscopically derived measurement. Although there is a small offset between the two different force measurements, the overall agreement is very good, and the results are very reproducible from cycle to cycle. The difference in the two types of measurements may be due to friction, which we neglected, as well as small differences in the properties of the hydrogels between compression and decompression 20 . These may also be responsible for part of the hysteresis. In the remainder of this paper, we use the force values derived from the images, for consistency with other microstructure measures.
With large data sets of particle-scale data, it is possible to obtain reliable statistics on all quantities of interest. We measured data from the 20 compression cycles, but we avoid the first five cycles in the analysis as they show weak transient effects. We compute statistics (microscopic averages noted by h Á i) only on grains that do not touch the walls, and on contacts between such grains, to reduce any boundary effects. We also separate the compression and the decompression motion at the same compression level to display the full mechanical trajectory. The measures presented in Fig. 3 are typically averaged over 10,400 samples for contacts and 4,100 samples for grains.
These data allow us to link microscopic quantities to similar global properties of the packing, which we discuss in the context of Fig. 3. In particular, we provide insight into the nominally power-law relation between the macro-scale force F and the total strain D: FpD b . We contrast F and D to their microscopic equivalents, hfi the average force at contacts, and hdi the average contact deformation. These last two variables are related by Hertz' law, together with the radius of curvature, r, at contacts. A key observation is that F(D) and hfi(hdi) follow substantially different functional relations 21 .
Linking microstructure and global packing response. Figure 3a shows the uncompressed and fully compressed states corresponding to one of the cycles. In Fig. 3b, we show data for the global force, F, versus the global compression D. A fit of these data to a power law yields an exponent b ¼ 2.2 ± 0.2, which is significantly larger than the microscopic (Hertzian) force law exponent of 1.5 measured for intergrain forces. This non-Hertzian packing mechanics has even been observed in much larger packings 21 , so it is not a finite size effect-the key point is that F depends on hfi as well as on additional microscopic properties. This is already partly revealed through the small hysteresis visible in the compression loop of Fig. 3b. The pressure, exerted on the top plate, can be expressed as a combination of the averaged microscopic quantities. Using a mean field argument detailed in the Methods section, we obtain: where hZi is the average number of contacts per particle, hbi is the average distance between grain centres and hji is the average packing fraction of grains within the packing (not touching the boundaries). In Fig. 3c, we show this mean field relation is well satisfied. This result gives a simple explanation of why the global force response with strain has a larger effective exponent than the particle-scale force law. Figure 3d shows data for hZi versus hji. Near jamming, for hji near 0.64, we expect hZiE6, whereas we measure Z between 4 and 5. This error is due to experimental uncertainty in distinguishing between two particles that are very close, but not in contact, versus actually in contact. By contrast, near maximal compression, the experimental error in hZi is much lower. For the largest packing fractions, hji40.74, our particles approximately organize into a crystalline lattice (see Methods), with a coexistence of hexagonal compact and cubic face-centred organization modes that implies crystal defects. Note that hji ¼ 0.74 corresponds to the packing fraction for an ordered close-packed lattice of hard spheres (FCC or HCP), although our particles are deformable. Using an overestimate for errors in all relevant quantities, including hZi, we obtain the error bars indicated by crosses for the scaling relation of Fig. 3c and equation 1. Note that the scaling relation of equation 1 is not strongly affected in absolute terms by errors in hZi, since the region where hZi is most uncertain corresponds to small values of hfi. Now that we have established the microscopic ingredients for the macroscopic mechanical response via equation 1, we can further elucidate the micro-macro link of the constitutive relation P(D) or, equivalently, F p D b . Importantly, we note that F p D b is approximately, but not strictly a power law. For example, hji, in equation 1, varies with D as hji p 1/(c-D), where the constant c depends on the packing geometry. We further see that F depends on D indirectly through hZi, so the strain-evolving topology of the microstructure matters as well-indeed, the micro-macro link is very complex. However, the fact that F is roughly a power law in D suggests that there may be other effective power-law-like relations between microscopic and macroscopic properties. We indeed find such relations. For instance, we find that hdi is well approximated by hdi p D a , with aE1.4 ± 0.1, providing we shift hdi by a small offset, hdi min , which is below the experimental resolution (Fig. 4a), and that that hfi p D g , with an empirical exponent gE1.7 ± 0.2 (Fig. 4b). It is important to note that we cannot directly link all the empirical exponents a, b and g, as hdi appears both in Hertz' law and as part of hbi in equation 1 (see the Methods section).
These observations yield the notable finding that, because a41, there is an effective mechanism through which the global strain D produces nonlinearly amplified deformations hdi at the grain level. This makes the nonlinear deformation amplification highly relevant, while also showing the deep entanglement of force and structure in the constitutive modelling of disordered packings of particles.

Discussion
We present a novel experimental technique that yields all microscopic contact and force vector information from a 3D packing of particles subject to controlled deformation. With these microscopic data, we have been able to verify a quantitative relation between the macroscopic mechanical response of the  packing and microscopic structural metrics of the packing. Furthermore, our experiment has revealed an important nonlinear enhancement of contact deformation in response to a global packing deformation, which significantly impacts the understanding and modelling of granular materials. Our new experimental approach has great potential for yielding understanding of other granular systems, including the granular response to shear, the jamming transition, particle diffusion, effects of particle shape and so on. It supplements other 3D imaging methods, such as microCT, NMR and confocal imaging, opening up a wide range of opportunities to shed new light on the poorly understood mechanics of disordered materials.

Methods
Derivation of the contact forces from geometric properties. We compute the deformation, d, at a given contact from the geometric properties of the grains involved in that contact. All grains have the same Young's modulus, so d is split equally between the two grains 22 . The radius of curvature, r, of the undeformed surface (see Fig. 5), which is different for each grain and each contact, is therefore r ¼ d þ ½d, with d the distance between the centre of the contact area and the grain centre. We observe that the force is given by F ¼ E e r e 1/2 d 3/2 and F ¼ E e d a, with E e the effective Young's modulus, 1/r e ¼ 1/r 1 þ 1/r 2 the effective contact radius and a ¼ O(A/p) the radius of the area of contact, A. Equating the two expressions for the contact force we obtain a cubic relation, which gives us d. We have independently measured the Young's modulus for our particles: EE23 kPa, see Fig. 6.
Derivation of the mean field scaling relation. The stress tensor can be computed 23 with a relation of the form s ¼ 1=V P c2V b c f c , with V an averaging volume and cAV the contacts in that volume (V is actually replaced by a smoothing kernel 19 ). For each contact c, b c is the vector between the centres of the grains and f c is the force vector along the contact normal. Neglecting friction and nonsphericity, b c and f c are nearly aligned; hence, the trace tr ðb c f c Þ % b c Á f c . The number of terms in the sum depends on the density of contacts, which is about ½ Z j, with Z the number of contacts per grain and j the grain volume fraction within the packing (for grains not touching the boundaries). The grains are nearly density matched with the surrounding fluid (density difference of E10 kg m À 3 ) so we can ignore the hydrostatic pressure. With roughly spherical grains, the isotropic pressure is very nearly proportional to tr b c f c . Hence, p p hb f Z ji. If we write any one of the quantities, q, on the right as q ¼ hqi þ e, the contact force, f, has qualitatively different statistical properties compared with b, Z and j. The contact forces are distributed broadly over 0rfrf max , while the others have relatively narrow fluctuations around a non-zero mean value. For instance, b is typically twice the roughly constant particle radius minus a typical small deformation d.
Since f is the only quantity that can change by orders of magnitude in the experiment, we use the mean values of b, Z and j in the expression for p and assume that deviations from the mean are uncorrelated, which means that pphbihZihji Á hfi. The integral of p over the top plate yields the top plate force, F, which reflects the mean field relation for p. When the top plate barely touches The surface of each grain (blue/green) is shown as solid lines, the hypothetical 'undeformed' surface of the grains is shown as dashed lines. The deformation, d, for this contact is the unknown quantity that we infer from linear elasticity. We measure the area of contact A, the position of the contact centroid C and its distance d to the center of mass G of each grain. Note that d (hence r) can vary for different contacts on the same particle, as particles are weakly non-spherical. The distance between grain centres b depends mainly on the grain polydispersity and weakly on their non-sphericity so, on average, hbiE2hdi. Our force inference algorithm uses the observed surfaces, without assuming sphericity, and provides experimental access to all the relevant quantities in order to reconstruct d. the grains, F should be null and, due to the density matching, so should hfi.
We subtract the small minimal measured values from hfi and F, which are present at the instant the plate touches the grains and which we attribute to measurement inaccuracies. This guarantees that F ¼ 0 when hfi ¼ 0. The coefficient of proportionality is the averaging volume, V, divided by the top plate area. As shown in Fig. 3c, this relation is remarkably well respected.
Crystalline structure of the packing. We analyse this crystalline structure by looking at the frequencies of the preferred directions of contacts, in each angular sector around a grain's center of mass. This gives a spherical probablility distribution, the amplitude of which in each angular sector can be represented as a surface (see Fig. 7).
Inference of 3D geometrical properties from images. For a given state of our system, we measure images of the packing cross-sections. From these we infer the geometric properties of the grains and their contacts. We use standard image denoising, local grey level renormalization, and image rescaling to compensate varying optical path lengths. We combine all the slice images into a 3D representation in terms of voxels. We detect border voxels by using a combination of statistical tests and signal processing techniques to account for physical properties of the hydrogels (for example, homogenous grey levels within grains, compensation of artifacts induced by index mismatch at grain boundaries due to their surface properties and so on). These methods are detailed in the publicly available source code. They ensure that very few false-border outliers remain. We then uniquely attribute border voxels to grains: a local tangent plane is fit at each voxel, such that the plane normal determines the direction towards a rough estimate of the grain centre at an average radius distance. These protocentres are considered to be distributed around a unique grain identifier, which we compute with 3D kernel density estimation. This process yields a cloud of border voxels uniquely attributed to each grain. We then fit an analytic shape to these borders such that for each direction u (unit vector), the surface of the grain from the center of mass in that direction is given by a function s(u) that best interpolates through all border voxels. This function s is expressed in terms of a basis of spline functions on S 2 {u}, the unit sphere. Unlike their Cartesian counterpart (for example, non-uniform rational B-splines (NURBS)), the triangular spherical b-splines 24 do not introduce any singularities on the sphere, and they are maximally isotropic. We use the convexity of the hydrogel grains as a regularizing condition, and we then perform a Levenberg-Marquardt least-square optimization to find the coefficients of s in the spline basis, so as to best fit border voxels, discarding any remaining outliers. Once we have analytic descriptions of the grain surfaces, it is easy to compute their centres of mass and any other geometric quantity, including the areas of contact.
Experimental methods. The granular material used here consists of 514 hydrogel beads 16,20 . The beads are approximately spherical and roughly monodisperse with a typical diameter of 2.1±0.1 cm. They are immersed in a waterpolyvinylpyrrolidone 360,000 MW(polyvinylpyrrolidone) solution to match the particle index of refraction to the surrounding fluid. This matching allows interior optical access, which is necessary for the refractive index-matched tomography technique used here 15 . Since the particles are almost entirely composed of water, they are also nearly density matched with the fluid; the particle density is o10 kg m À 3 greater than the fluid density. By itself, the fluid-particle system is completely transparent. To obtain contrast, we dye the particles with a hydrophobic fluorescent dye (Nile Blue 690), which can be excited with a laser sheet. This laser sheet (Lasiris SNF 635 nm, 25 mW) moves on a fast linear stage, which sequentially illuminates the sample slice-by-slice. A fast camera (Basler ava1000-120) equipped with lens and long pass filter records the fluorescent image of each slice. This produces a sequence of particle cross-sections. The camera is mounted on the same stage as the laser, as in the design of ref. 25. A complete scan typically consists of 360 slices. We confine the granular system in a rectangular box with base size 16.5 Â 16.5 cm and carry out cyclical uniaxial compression/decompression using a stage-controlled piston as described in the main text. The top piston is made from a 6 mm thick perforated sheet, which allows flow of index-matching fluid into and out of the packing. Except for gravitational gradients, the pore pressure inside the fluid is uniform and matched to the ambient atmosphere. The piston is driven by a linear stage (Newport MTM250) and controller (Newport XS4), with a step resolution of 1 mm. The force on the piston plate is measured with a force sensor (Loadstar RSB4-005M-A). The compression speed during compression/decompression is 0.1 mm s À 1 to reduce fluid-induced shear stresses on the particles and to drive the granular system completely quasi-statically. The spacing between the piston plate and the walls is small enough that particles cannot escape confinement, so the number of particles during an experiment is constant. The container height depends on the compression level, but a typical uncompressed height is about 15 cm. More experimental details can be found elsewhere 20 .