Simultaneous measurement of the microscopic dynamics and the mesoscopic displacement field in soft systems by speckle imaging

The constituents of soft matter systems such as colloidal suspensions, emulsions, polymers, and biological tissues undergo microscopic random motion, due to thermal energy. They may also experience drift motion correlated over mesoscopic or macroscopic length scales, \textit{e.g.} in response to an internal or applied stress or during flow. We present a new method for measuring simultaneously both the microscopic motion and the mesoscopic or macroscopic drift. The method is based on the analysis of spatio-temporal cross-correlation functions of speckle patterns taken in an imaging configuration. The method is tested on a translating Brownian suspension and a sheared colloidal glass.


Introduction
Rough surfaces and scattering media generate a characteristic speckle pattern [1] when illuminated by coherent light, e.g. from a laser. By analyzing a time sequence of speckle patterns, valuable information can be retrieved on the sample evolution. Broadly speaking, one may distinguish between "static" speckle patterns generated by solid objects and "dynamic" speckles formed by soft matter systems (e.g. colloidal suspensions, emulsions, polymer solutions, biological tissues), whose components undergo Brownian motion, thereby continuously reconfiguring the scattered speckle pattern. In the former case, relevant to metrology and interferometry [2], a rigid displacement or a long wave-length deformation is often measured, e.g. in response to vibrations, applied load, or a change of temperature; speckle patterns are recorded onto a 2D detector such as a CCD or CMOS camera, using an imaging optics. By contrast, the microscopic (e.g. Brownian) dynamics of soft systems is quantified by g 2 (τ) − 1, the autocorrelation function of the temporal fluctuations of the scattered intensity. In these dynamic light scattering measurements (DLS, a.k.a. photon correlation spectroscopy [3]) a point-like detector (e.g. a phototube) is placed in the far field, where it collects light within a few speckles.
Recent developments have made the distinction between these two research fields increasingly fuzzy. On the one hand, imaging geometries similar to those for static speckles have been used to detect motion, e.g. in vascular flow ( [4,5] and references therein). Motion is typically quantified by inspecting the contrast of the speckle pattern, < I 2 > / < I > 2 , where I is the local intensity and · · · is an average over a small region centered around the point of interest. The method is based on the fact that speckles are blurred and the contrast reduced in those regions where significant motion occurs [6]. Rigid motion has been measured also by laser speckle velocimetry [7,8], by applying spatial cross-correlation methods to speckle images, an approach similar to particle imaging velocimetry [9] and image correlation velocimetry [10] often used in fluid mechanics. On the other hand, CCD and CMOS cameras are now routinely used as detectors for DLS, especially for samples exhibiting slow dynamics, i.e. speckle fluctuations on time scales from a fraction of second up to several hours. In the first implementations of these so-called multispeckle approaches [11,12], the far field detection scheme of traditional DLS was used and the intensity correlation function was calculated from where I p (t) is the intensity of the p−th pixel at time t and < · · · > p and < · · · > t indicate averages over pixels and time, respectively. In the far field detection scheme, each pixel receives light issued from the whole scattering volume and all pixels are associated to nearly the same scattering vector q = 4πnλ −1 sin θ /2, where n is the solvent refractive index, λ the in-vacuo laser wavelength and θ the scattering angle. As in traditional DLS, the decay time of g 2 (τ) − 1 is related to the time it takes a scatterer to move (relative to the other scatterers) over a distance ∼ q −1 [3]. A further step in bridging the gap between DLS and speckle imaging methods is represented by photon correlation imaging (PCI) [13]. In PCI, DLS data are obtained by analyzing time series of speckle patterns acquired using a 2D detector and a low-magnification imaging optics. Similarly to conventional imaging, a given area of the detector corresponds to a well-defined region in the sample. Unlike conventional imaging, however, the image is formed using only light scattered in a narrow range of scattering vectors q . This allows one to calculate a spatially-resolved version of Eq. (1), where pixel averages are performed on small regions of the detector, thereby providing information on the local dynamics. This method has been applied to systems whose microscopic dynamics are significantly heterogeneous in space, such as glassy or jammed soft matter [13,14], and it has been extended to highly turbid media such as drying coatings [15] and granular systems [16]. As highlighted by this short overview, previous works have focussed either on the measurement of the rigid displacement of a set of scatterers, regardless of any relative motion between them, or, conversely, on their relative motion due to the microscopic dynamics, regardless of any average drift component. Here, we present a new method that allows one to quantify in a single measurement the contribution of each of these phenomena to the evolution of speckle patterns formed in the imaging geometry [17]. Figure 2, which will be discussed later, shows the essence of the method in a glimpse: the spatial cross-correlation of two speckle images taken at a time lag τ exhibits a peak, whose position and height yield the sample rigid shift and its internal dynamics over the time τ, respectively. Applications of this method to gels submitted to gravitational [18] or internal [19] stress have been presented in previous publications, but the method itself was not discussed there. In this paper, we provide a detailed description of the algorithm used to implement the method, addressing in particular the challenges inherent to PCI experiments, i.e. the reduced size of the speckles and the need to process a very large num-ber of images in a reasonably short time. Finally, we demonstrate our method by testing it on a model system, a suspension of Brownian particles contained in a cell displaced by a motor, and by measuring the velocity profile and the microscopic dynamics of a sheared colloidal glass.

Sub-pixel Digital Imaging Correlation algorithm
The first step in our method is to find the local rigid displacement with sub-pixel resolution. To this end, we use a cross-correlation technique inspired by particle imaging velocimetry [9] and image correlation velocimetry [10]. A time series of images of the sample is taken, using a PCI setup. Each image is divided into a grid of regions of interest (ROIs), corresponding to square regions of side L in the sample. Under a rigid drift, the speckle pattern in a ROI at time t will appear in a shifted position in a successive image taken at time t + τ. If the scatterers undergo relative motion, in addition to a rigid shift, the shifted ROI will not be an identical copy of the original one. Still, the displacement field can be estimated by calculating by what amount a ROI of the second image has to be back-shifted in order to maximize its resemblance with the corresponding ROI of the first image. In practice, for a given ROI the shift along the horizontal and vertical directions, ∆x and ∆y, is determined by searching the maximum of corr[I, J], the spatial crosscorrelation of the intensity of the two images, defined by with In the above equations, I r,c is the intensity at time t of the pixel at row r and column c, J r,c is the intensity at the same location but at time t + τ, and k and l are the shifts expressed in number of rows and columns, respectively. Here and in the following, double sums over r and c extend over all rows and columns for which the terms of the sum exist, in this case the N pixels of the overlap region between the full image and the shifted ROI. For computational efficiency, covar[J, I] is usually calculated in Fourier space [20]. Note that covar[J, I] → 0 far from the peak, where I and J are uncorrelated. The position (k, l) of the global maximum of corr[J, I] yields the desired displacement along the direction of columns (x axis) and rows (y axis), ∆x = l and ∆y = k respectively, with pixel resolution. Several schemes have been proposed to improve this resolution, e.g. by calculating the centroid of corr[J, I], or by fitting its peak to a 2-dimensional analytical function such as a Gaussian. While both methods work well for broad, circularly symmetric peaks, they tend to be less robust when the peak is sharp or it has an asymmetric shape. The shape of the peak is determined by the spatial autocorrelation of the intensity pattern; for our speckle images, it depends on the shape and size of the speckles, which may not be symmetrical, depending on the shape of the illuminated sample volume and the imaging optics [1]. Moreover, the peak usually extends over just a few pixels, because one minimizes the speckle size in order to maximize the information content in the image. To overcome the limitations inherent to peakbased schemes, we use an alternative approach based on a least-square method that allows us to obtain the displacement field with a typical resolution of a few hundredths of a pixel, with no requirements on the shape or broadness of the peak and without using any fitting function. The first step of our shift-finding algorithm is the same as in standard PIV methods: corr[J, I] is calculated and the pixel-resolved shift, (k, l), is determined from the position of its global maximum. The next step consists in the refinement of such displacement with sub-pixel resolution. For the sake of simplicity, let us first assume that J is simply a shifted version of I, with displacement (∆x = l + δ x, ∆y = k + δ y), with |δ x| < 1, |δ y| < 1. The intensity J r,c may then be expressed as a weighted average over a suitable set of pixels of the intensity of the image I. In principle, an infinite number of terms are needed to obtain J r,c , if δ x and δ y are non-integer [21]. In practice, linear interpolation is usually sufficient to reconstruct J to a good approximation, thereby greatly simplifying the calculation. Using linear interpolation, the intensity J r,c at a given pixel is expressed as the weighted sum over the (at most) four pixels of I that partially overlap with that pixel, as exemplified in Fig. 1a. Thus, The coefficients d are the overlap areas shaded in Fig. 1a: and similarly for the other terms. The term ε r,c has been added to account for the fact that in general J will not be an exact (albeit shifted) replica of I, because of experimental noise and due to any evolution of the speckle pattern, which in our case is due to the microscopic dynamics of the scatterers. For distinct speckles, these fluctuations are uncorrelated [1,3]; we thus treat ε as a noise term and determine the displacement (∆x, ∆y) as the rigid shift that minimizes, in a least-squares sense, the difference between J and the linear combination of I in the r.h.s. of Eq.(5). More specifically, we search for a set of four coefficients a = {a 1 , a 2 , a 3 , a 4 } that minimizes the cost function χ 2 defined as where explicit expressions relating a to the sub-pixel shift and k i , l i to the pixel-resolved shift will be provided in the following. Note that in principle nine coefficients are required, instead of the four introduced here, since the direction of the shift is not known a priori. However, we expect that only up to four of them differ significantly from zero. In order to speed up the determination of a, we calculate the centroid of the peak of corr[J, I] in order to predetermine the direction of the shift (e.g. top-left, top-right etc.), so that only four coefficients have to be computed. This is done by calculating in which of the four quadrants labeled by A, ..., D in Fig. 1b lays the center of mass of the crosscorrelation peak, as explained in detail in Sec. 7, where we also provide explicit expressions for k i , l i . Note that the centroid algorithm is only used to determine the direction of the shift, not its subpixel magnitude, thus avoiding the limitations recalled above for sharp or non-symmetric peaks. The optimum shift is obtained by imposing ∂ χ 2 /∂ a i = 0, i = 1, .., 4. By exchanging the order of the sums in Eq. (6) one recognizes that a is the solution of b = M · a, with: Any standard method is suitable to solve the above set of equations; in our implementation, we use singular value decomposition [20]. It should be emphasized that b is known, since covar[J, I] has been already calculated to estimate the pixel-resolved displacement (see Eqs. (2,3)). Therefore, the only extracost required for calculating the displacement with subpixel resolution is the computation of covar[I, I] and the solution of the set of linear equations b = M · a. This moderate computational extracost is due to the fact that a linear interpolation scheme has been adopted in Eq. (6): higher-order interpolations, although more precise, would lead to a much more complex, non-linear minimization problem. Finally, we note that in a typical multispeckle DLS experiment, one calculates the correlation functions for a given starting time (i.e. a given image I) and several time delays τ (i.e. several images J). Therefore, the computational cost for calculating covar[I, I] is shared between several lags, further increasing the efficiency of the algorithm. Once a is computed, the shifts along the x and y directions are calculated with subpixel precision according to (see Eqs. (27) and (29) in Sec. 7 for the definition of k 1 and l 1 ).

Dynamic Light Scattering: corrections to g 2 (τ) − 1 for drifting samples
In order to quantify the internal dynamics, one needs to compute the intensity correlation function g 2 − 1 between image J and a shifted version of I, so as to avoid any artifact due to the rigid shift of the speckles. Denoting by I ′ the image I shifted by (∆x, ∆y), the (un-normalized) intensity correlation function corrected for the shift contribution is The shifted image I ′ may be constructed using an interpolation method. Tests on real speckle images show that linear interpolation, although suitable for determining the shift with good accuracy, is not precise enough to reconstruct a shifted version of I suitable for the calculation of g 2 − 1. Higher-order interpolation schemes are thus required. As shown in Ref.
[21], image shifting by interpolation is equivalent to convolving the original image with a suitable kernel: where we have assumed for simplicity that the kernel is symmetrical and separable, i.e. that 2D h(x, y) = h(x)h(y). Unfortunately, in our case this approach would be too time-consuming, because it requires a convolution operation, Eq. (12), in addition to the calculation of the correlation function, Eq. (11). We introduce below an alternative method that leads to a much faster algorithm, where the only computational cost is that of evaluating the kernel for a few points, with no need for the calculation of convolution and correlation functions. It is convenient to consider only nonnegative fractional shifts δ x, δ y, given by with j x = floor(∆x), i y = floor(∆y), where floor(x) is the largest integer ≤ x. As we shall discuss it below, the choice of the kernel is not crucial; a good choice is a truncated, windowed sinc function with an even number, M, of supporting points: a property that will be of use in the following. Using Eqs. (13,14), the convolution product (12) may be rewritten as By replacing the r.h.s. of Eq. (18) in Eq. (11) and by exchanging the order of the sums, we obtain: Finally, by recalling the definition of the covariance between J and I, Eq. (3), and using the fact that the kernel is DC-constant, Eq. (17), one obtains or, equivalently, where I = N −1 ∑ r,c I r,c and similarly for J. Equation (21) is the central result of our method. It shows that the intensity correlation function corrected for the shift contribution can be simply obtained as a linear combination of a few terms of covar[J, I], weighted by the kernel. Since covar[J, I] has already been calculated to determine the shift, the extra cost is essentially limited to the evaluation of M 2 values of the kernel, which is typically negligible. Finally, we note that covar[J, I] vanishes on the length scale of the speckle size as its argument departs from (i y , j x ), which is close to the location of the peak of the covariance. Hence, it is sufficient to take M on the order of a few speckle sizes (in units of pixels), because in Eq. (21) the contribution of the kernel for larger lags would be multiplied by a vanishingly small quantity. For example, we find that for images with a speckle size of about 5 pixels, the correction is virtually independent of M for M ≥ 8.

Experimental tests
We test our method on two samples: a suspension of Brownian particles loaded in a cell displaced by a motor, and a colloidal glass to which a shear deformation is applied. The Brownian sample is a dispersion of polystyrene microspheres (radius R = 0.265 µm) in an aqueous solution of fructose at 75% weight fraction. The particle volume fraction is 10 −5 and the sample is kept at a temperature T = 9 • C. The setup is described in [23]; we use the imaging geometry shown in Fig. SM1 a) of [13], where an image of the sample is formed onto a CCD detector using light scattered at θ = 90 deg, corresponding to a scattering vector q = 2.46 µm −1 . The field of view is 1820 × 364µm 2 and images are acquired at a rate of 10 Hz. The sample cell is attached to a motor that can impose a drift in the y direction at a controlled speed, v y = 10µm s −1 . Figures 2 a)-c) show the spatial crosscorrelation calculated applying Eq. (2) to pairs of speckle images taken while displacing the sample, for three different time lags. For τ = 0 s, Eq. (2) yields the spatial autocorrelation of the speckle pattern: accordingly, a sharp peak of height one and centered at ∆x = ∆y = 0 is visible, whose FWHM ≈ 2.9 pixels provides the speckle size. As the lag increases, the peak position drifts in the y direction, due to the translation motion imposed by the motor. Additionally, its height decreases, due to the relative motion of the Brownian particles that reconfigure the speckle pattern. Figure 2 d) shows a cut of the crosscorrelation peak along the ∆x = 0 direction, for four values of τ. From this plot, it is clear that if g 2 (τ) − 1 was to be computed from a purely temporal crosscorrelation, as in Eq. 1, one would observe a spurious, fast decay, essentially due to the rigid shift only. This would correspond to follow corr[I(t), I(t + τ)] at ∆x = ∆y = 0, as a function of τ. By contrast, if the relative motion of the Brownian particles is to be obtained, one has to measure the height of the peak as it drifts, as in the method proposed here. The inset of Fig. 2 d) shows the same data, plotted as a function of distance along y with respect to the (subpixel) peak position. It is worth noting that the peak width remains constant, in contrast to what suggested (but not demonstrated, to our knowledge) in patent literature [24], where it was proposed that the peak would broaden with τ as a result of the internal motion of the scatterers. Thus, the relevant parameter for extracting the relative motion is indeed the peak height, not its width.
In figure 3 a), we show the displacement of the speckle pattern as a function of τ, obtained from the sub-pixel peak position averaged over 200 pairs of images (i.e. 20 s), taken while translating the sample. The data are very well fitted by a linear law (red line) as expected for motion at constant speed, thus indicating that our algorithm captures very well the drift   Fig. 2. a)-c) spatial cross-correlation between speckle images generated by a diluted Brownian suspension that is translated along the y direction by a motor. The speckle patterns are recorded on a CCD using an imaging collection optics (see text for more details). As the delay τ between pair of images is increased, the peak position shifts to larger y and its height decreases, due to the relative motion of the Brownian particles. d) cut of the crosscorrelation along the ∆x = 0 line. Curves are labeled by τ. Inset: same data, replotted as a function of spatial shift with respect to the peak position. component of the speckle pattern, from a fraction of a pixel up to tens of pixels. The error bars are the standard deviation of the displacement over the measurement time, σ y . For τ ≤ 4.2 s, σ y /∆y < 4% and the error bars are smaller than the symbol size, indicating that the detection of the peak position is very reliable, even when the peak height is as low as 0.1 or the displacement is just a fraction of a pixel. For τ = 5.6 s, the error bar is significantly larger, because the peak height becomes comparable to the noise level and the peak position can be hardly resolved. Beyond τ = 5.6 s, no peak can be reliably found, thus preventing the displacement to be measured. The slope of the linear fit to the data is 3.30 ± 0.03 pixel s −1 .
Recalling that the nominal speed of the motor is v y = 10µm s −1 , this implies that 1 pixel corresponds to 3.03±0.03 µm in the sample, in excellent agreement with 3.15±0.15 µm/pixel as obtained from the magnification of the imaging system, evaluated using geometrical optics. Figure 3 b) shows the intensity correlation function g 2 (τ) − 1, averaged over 20 s. If the sample is kept at rest during the measurement (black squares), the intensity correlation function exhibits an exponential decay, as expected for diluted Brownian suspensions [3], as better seen in the inset that shows the same data in a semilog plot. When the sample is translated at constant speed, the uncorrected g 2 − 1 decays on a much shorter time scale (blue circles) and its shape departs from a simple exponential. Clearly, no information on the microscopic dynamics can be obtained from the uncorrected data. The red crosses are the data corrected according to Eq. (21): for τ ≤ 4.2 s the corrected g 2 − 1 is very close to that measured for the stationary sample, thereby demonstrating the effectiveness of our correction scheme. For larger lags, the corrected data tend to overestimate g 2 − 1: this is consistent with the fact that the displacement cannot be reliably measured, as discussed in relation to fig. 3 a). Indeed, in this case the peak-finding algorithm spuriously interprets the highest level in the noisy base line of corr[I(t), I(t + τ)] as the (higher-than-expected) degree of correlation. Having validated our method on a model sample whose displacement is well controlled, we test it on a more realistic experimental situation, a sheared colloidal glass for which the velocity field is not uniform over the field of view. The sample is a dense suspension of hard-sphere-like colloidal particles, a widely-studied model system for the glass transition [26]. The particles have radius ≈ 100 nm (as in [25]) and volume fraction ϕ ≈ 0.6. The setup is similar to that for the Brownian sample, but here the sample is kept in a square cell of section 10 × 10 mm 2 , in which a glass bead of diameter D = 5 mm is inserted. The bead is attached to a motor that displaces it in the y (vertical) direction, parallel to the cell wall, at a speed v y = 0.1 µm s −1 .
The minimum gap e between the wall and the bead surface is 1280 ± 90 µm. For the small displacements studied here (≤ 850 µm) and given that e << D, the deformation is close to a simple shear. The sample is illuminated by a laser sheet in the vertical (x, y) plane, of thickness ≈ 100 µm. We image a region of size 710 × 530 µm 2 using light scattered at θ = 90 deg, corresponding to q = 20.6 µm −1 . To obtain space-resolved information on the mesoscopic displacement and the microscopic dynamics, we run our algorithm on ten ROIs of size 31 × 264 µm 2 regularly spaced at a growing distance x from the moving wall. Figure 4a shows the velocity profiles close to the bead wall, for two times t after starting shearing the sample. In this representation, the slope of the data is the local shear rateγ. The solid line shows the velocity profile expected for homogeneous shear, corresponding to an average shear rate across the whole gap ofγ = 7.9 × 10 −5 s −1 . It is clear that already at t = 120 sγ is non-uniform across the gap, with a highly-sheared band close to the moving surface (x ≥ −164 µm,γ ≈ 1.4 × 10 −4 s −1 ), followed by a low shear region (x ≤ −220 µm, γ ≈ 6.0 × 10 −5 s −1 ). Similar shear banding has been reported for other colloidal glasses [27].
Interestingly, shear banding is seen to evolve with time. At t = 820 s, the shear rate for the high-and low-shear bands is comparable to that at t = 120 s (γ ≈ 1.3 × 10 −4 s −1 anḋ γ ≈ 6.1 × 10 −5 s −1 , respectively), but the boundary between the two zones has moved from x = −190 µm to x = −240 µm. Additionally, the occurrence of a marked drop of v y close to the moving surface suggests slipping. This behavior is reminiscent of that reported for a variety of jammed or glassy soft materials, see e.g. [28], which exhibit complex spatio-temporal shear patterns. Figure 4b) shows the intensity correlation function measured for the ROI at the position indicated by the arrow in a). For the unsheared sample (crosses) no dynamics is observed on time scales up to 20s, about 2000 times the Brownian time for the same particles in the diluted regime. This is consistent with the notion that the microscopic dynamics of a sample at rest is slowed down by orders of magnitude on approaching the glass transition. The open symbols show the uncorrected g 2 − 1: a fast decay is observed, essentially due to the translation of the speckle pattern due to the imposed shear. Once corrected, the data still show a decay of g 2 − 1 (albeit a slower one), thus indicating that particles move with respect to each other, in addition to be advected by the shear deformation. We emphasize that the corrected g 2 − 1 is sensitive to the component of the particle displacement along the direction of q, which lays in the horizontal plane, perpendicular to the shear direction. Therefore, the decay of g 2 − 1 is not due to the affine component of the particle displacement along y, but rather to irreversible rearrangements associated with flow in glassy systems [29]. Interestingly, we find that the decay of g 2 − 1 is faster at t = 820 s, when both the localγ and its gradient are larger. This suggest a direct relation between (local) shear rate and plastic rearrangements, as proposed for granular materials [16] and emulsions [30].

Conclusions
We have introduced a method to obtain the mesoscopic displacement field and the microscopic dynamics in soft materials where the constituents undergo both a drift motion and a relative displacement. The algorithm proposed here is optimized for the typical features of speckle images in PCI experiments, where a small speckle size is highly desirable to maximize the spatial resolution and the statistics of the measurement. The algorithm is highly efficient in that the correction of g 2 − 1 does not requires any significant computational extracost, besides that necessary to determine the displacement field. The method has been successfully tested on a Brownian suspension and a colloidal glass. Although similar information may be in principle obtained using confocal or optical microscopy, our method allows one to investigate samples that are difficult or impossible to visualize in real space, such as the very small particles of our colloidal glass. A generalization to speckle patterns obtained under partially coherent illumination, such as in recent microscopy developments [31, 32] is also possible [33]. The method presented here should be particularly valuable for soft materials where slow dynamics is coupled to the effects of an external stress, as in rheological experiments or in samples submitted to an external field such as gravity [18], or in disordered jammed materials, where internal stress is known to play a major role [19,34] in the sample dynamics.