Spatial blurring in laser speckle imaging in inhomogeneous turbid media.

Laser speckle imaging (LSI) has developed into a versatile tool to image dynamical processes in turbid media, such as subcutaneous blood perfusion and heterogeneous dynamics in soft materials. Spatially resolved information about local dynamics is obtained by measuring time-dependent correlation functions of multiply scattered light. Due to the diffusive nature of photons in highly scattering media, the measured signal is a convolution of the local dynamics in the material and the spatial distribution of photons. This spatial averaging inevitably leads to a loss of resolution, which must be taken into account for a correct interpretation of LSI measurements. In this paper we derive analytical expressions to quantify the effects of spatial blurring in backscatter LSI for materials with heterogeneous dynamics. Using the diffusion approximation, we calculate the photon density distribution for a semi-infinite material, and we predict the effect of dynamic heterogeneity on the measured correlation function. We verify our theoretical expressions using random walk simulations. Our results show that LSI measurements in dynamically heterogeneous materials should be interpreted with caution, especially when only a single wavelength and correlation time are used to obtain the dynamical map.

Introduction

Many disordered materials, such as rubbers, plastics, foams, suspensions, and biological tissue strongly scatter light, which makes a direct observation of structure and dynamics, for example with light microscopy, very difficult. However, detailed information about the dynamical processes in such non-transparent materials can be obtained by analyzing the temporal fluctuations of light that is multiply scattered in the material. This forms the basis of diffusing-wave spectroscopy (DWS)[1-3], a powerful technique to study Brownian motion or flow in concentrated colloidal systems[4-6] or to extract viscoelastic moduli of soft materials from the dynamics of embedded tracer particles[7]. While DWS has been applied mostly to uniform samples, it can also be used as an imaging tool for characterizing spatially resolved dynamics in heterogeneous turbid media[8-10]. A particularly useful realization of this is laser speckle imaging (LSI). In LSI, the sample is illuminated with a plane wave of coherent light. Photons that enter the sample undergo many scattering events before leaving the sample again and reaching the camera. Each camera pixel receives many photons that have scattered from different positions in the sample; the resulting path length differences create a random interference pattern that is known as a speckle pattern. Movement of the scattering particles causes temporal fluctuations in the speckle pattern, which can be quantified by calculating for each speckle the intensity autocorrelation function g 2(τ) = 〈I(0)I(τ)〉/〈I〉2, with I(τ) the intensity at time τ. The decay rate of this autocorrelation function is directly related to the dynamics of the scatterers in the material. Because the photons that reach the detector have scattered many times inside the sample, decorrelation already occurs when the scatterers have moved only a fraction of the wavelength. This makes it possible to measure particle displacements as small as a few nanometers with LSI[1].

Because of its high sensitivity, its relative cost-effectiveness, and its non-invasive nature, LSI has become an attractive tool in bio-medical imaging[11-16] and materials science[17-22]. However, obtaining quantitative information from LSI measurements remains difficult, because the decay of the autocorrelation function depends not only on the internal dynamics of the material, but also on the distribution of photon paths through the sample. This means that an accurate model is needed for the propagation of photons in the turbid medium. Usually, it is assumed that the propagation of light in the material can be approximated as a diffusion process, so that the path followed by an individual photon can be described as a random walk[23-25]. The distribution of path lengths can then be calculated by solving the diffusion equation. For samples in which the particle dynamics does not depend on position, this yields a direct relation between the intensity autocorrelation function and the mean-square displacement of the particles 〈Δr 2(τ)〉[2,4]. However, for dynamically inhomogeneous samples, the information carried by the speckle pattern depends on which regions the photons have probed. Since each photon takes a different path through the sample, LSI measurements performed on inhomogeneous samples actually represent averages, resulting from a convolution of dynamic processes in the sample and the spatial probability distribution of photons in the sample[26,27]. This leads to spatial blurring and loss of resolution and thus limits the ability of LSI to detect dynamic heterogeneities[28]. For a correct interpretation of LSI measurements on dynamically heterogeneous samples, it is therefore necessary to quantify the effects of spatial averaging.

The typical distance travelled by photons before leaving the sample in a backscatter experiment is a few times the transport mean free path l * [20], which is the average step length of the photon random walk. One would therefore expect the penetration depth of light into the sample and the extent of spatial blurring to be also a few times l * [29,30]. However, it is clear that the effect of blurring must also depend strongly on the correlation time τ, since the initial decay of the correlation function (at small τ) is much more sensitive to long photon paths than the long-time decay[1]. Hence, to fully appreciate how spatial averaging affects the measured LSI signal at different time scales, a more detailed calculation, taking all possible photon paths into account, is needed. Such an analysis has been performed for a multilayer medium, in which the sample is heterogeneous only in the z-direction (perpendicular to the surface)[10,31,32]. Here, we extend these findings to dynamically heterogeneous materials with arbitrary spatial distributions of particle dynamics. Using the diffusion approximation, we calculate the photon density distribution in the material and we derive an analytical expression that relates the autocorrelation function measured in a particular location to the distribution of mean-square displacements in the material. We then use this result to evaluate the lateral resolution of LSI, and we show how this depends on the correlation time. We verify our results using random walk simulations.

LSI theory

We consider LSI in the backscattering geometry for samples that are much thicker than the penetration depth of light in the material, so that the material can be considered as a semi-infinite half-space. We assume that no photon absorption takes place. Photons enter the material at a particular location and then undergo a sequence of scattering events before leaving the sample again at a point where they are detected by a camera (Fig. 1). The transport of photons in the material is characterized by the mean free path l, which is the average distance between two scattering events; it depends on the number density ρ and the scattering cross-section σ of the particles, l = 1/ρσ. The transport mean free path l * is the distance over which the direction of light becomes randomized; it is related to the mean free path by l * = l/〈1 − cosθ〉, where θ is the scattering angle and where the average is taken over the scattering form factor of the particle. For very small particles the scattering is isotropic and l * ≈ l, while for larger particles the scattering is peaked in the forward direction so that l * > l. The transport mean free path can be determined experimentally[33,34] or it can be calculated from Mie theory[23]. In the Rayleigh limit, for particles that are small compared to the wavelength of light, l * can be several tens to hundreds of micrometers large and increases very strongly with increasing wavelength, as .

Figure 1

A photon trajectory: after entering the sample, a photon performs a random walk with mean free path l through the material before exiting the sample at . The signal detected at is the result of all the random walks ending at . The random walk is assumed to start at the first scattering event, . Movement of the particles in a time τ leads to a change in the phase of the light wave, Δϕ (τ) = ∑ ⋅ Δ (τ) with  =   −  the scattering vector transfer of scattering event i.

When the total path length of a photon is much larger than l *, the photon trajectory can be described as a random walk with average step size l *. The total electric field collected at a particular location is the superposition of the fields of all photon paths p ending at :where E is the amplitude of the electric field from path p and ϕ (t) its phase at time t. The electric field correlation function g 1(τ), which is related to the experimentally measured intensity correlation function g 2(τ) by the Siegert relation, g 2(τ) = 1 + βg 1(τ)2 with β an experimental constant of order unity, can then be written aswhere is the fraction of the scattered intensity at coming from path p. Terms with p ≠ p′ do not contribute, because light waves belonging to different paths are uncorrelated. The phase shift Δϕ (τ) = ϕ (τ) − ϕ (0) is due to movement of the N scatterers along path p and can be written as[1,4] where  =   −  is the wave vector transfer of scattering event i, with a magnitude that is determined by the scattering angle as q = 2k 0sin(θ/2), and with Δ (τ) =  (τ) −  (0) the displacement of scatterer i in the time interval τ. For Brownian motion, the Δ (τ) are random Gaussian variables and the average in Eq. (2) is readily obtained[1] to give:with . Here, we will assume that l * does not depend on the position in the material; in other words, that the material is dynamically heterogeneous, but optically homogeneous. This is a reasonable approximation for many applications in materials science; for example, spatial variations in crosslink density in a polymer material will lead to variations in the local rigidity (and therefore in the local dynamics), without causing large variations in the refractive index. To simplify the analysis further, we pass to a continuum limit and replace the summation over all paths by an integral over all path lengths s = N l:where P( , s) is the relative contribution of paths of length s to the intensity at , s min the minimum path length (which is of order l *), and Γ( , τ, s) the weighted average mean square displacement of scatterers along all paths of length s that end at . This can be written aswhere 〈Δr(, τ)2〉 is the mean square displacement of particles at position and ρ(; , s) is the normalized photon density at position for diffusion paths of length s ending at [27]. For a homogeneous medium 〈Δr(, τ)2〉 does not depend on , so that Γ( , τ, s) = 〈Δr(τ)2〉 = 6Dτ with D the diffusion coefficient of the particles. However, for a medium with an inhomogeneous diffusivity, Γ( , τ, s) depends on what parts of the sample have been probed by the light. In this case, both the path length distribution P( , s) and the spatial density distribution of photon paths ρ(; , s) are needed to interpret the measured g 1( , τ) in terms of the spatiotemporal dynamics in the sample.

Photon density distribution in scattering media

To derive expressions for P( , s) and ρ (; ), we start by defining G(, n; , , s) as the probability for a random walk of length s, starting from a point on the surface of the material and ending at , to pass through point after a distance n. Since this random walk consists of two random walks, first from to in n steps and then from to in s − n steps, we can write:where the propagator G( , , n) gives the relative probability of paths from point to point with path-length n. The probability of all photon paths of length s to a point iswhere J( ) is the intensity of the incident light at . We consider here the case of homogeneous illumination, so that J( ) is constant along the surface z = 0. With this, we can writefor the probability of all walks of length s ending at to pass through r after a distance n. The spatial density distribution of photon paths ρ(; , s) is then obtained by averaging G(, n; , s) over all steps[27].

The distribution of path lengths can be calculated as:

Within the random walk approximation that we adopt here, the propagator G( , , s) is the solution of the diffusion equation[1,2].where the diffusion coefficient of photons in the sample can be expressed as D p = vl */3, with v the speed of light in the material. Since the path length s = vt, we can also write this aswith initial condition G( , , 0) = δ( −  ). Since the incident light becomes diffuse at a depth z ≈ l *, we take the beginning of the random walk in the sample to be at a distance z 0 = l * inside the sample. The boundary condition can be specified by requiring that for s 0 the net flux into the sample is zero[1]. It has been shown that this is almost equivalent to forcing G( , , s) to become zero a small distance outside the sample, at the extrapolated boundary condition z = −z , with z  ≈ 0.7l * [23,35]. The solution of the diffusion equation with this boundary condition (and constant l *) can be obtained using the method of images[36]:

Using this result with  =   = (x , y , 0), we can obtain the path length distribution from Eqs (11) and (8). Carrying out the integrations over and s gives (with s min = 0):

Since we consider homogeneous illumination and homogeneous optical properties of the material, the path length distribution does not depend on the location of the detector, so that we have dropped the dependence on . The path length distribution P(s) has a maximum around s ≈ l *, while for long paths, , , as found previously[1].

Next, we calculate the spatial photon density ρ(; , s) from Eq. (10). The convolution integral is most conveniently calculated by using the properties of the Laplace transform[27]:where denotes the Laplace transform with respect to s. Using the Laplace transform of Eq. (14), we findwith

Using Eqs (5), (6), (15), and (17), we can directly relate the field correlation function measured at location , g 1( , τ), to the spatial distribution of particle mean-square displacements in the material, 〈Δr(, τ)2〉.

Results

Intensity distribution

The goal of an LSI experiment is to generate a spatially resolved image of the dynamics in the sample. The intensity of multiply scattered light that leaves the sample at z = 0 is recorded with a camera; from the intensity fluctuations the correlation function can be obtained for each position in the z = 0 plane. Since the photons that reach have taken different paths through the sample, the correlation function is a weighted average of the dynamics in the region probed by these photons. The spatial distribution of photons arriving at can be calculated asand represents the relative contribution of particles in the region around to the signal measured at . As shown in Fig. 2a, most photons probe a region within a distance l * from the detector; nevertheless, a significant fraction of the photons also explores regions farther out in the sample. We have also compared our prediction for the intensity distribution with the distribution obtained from random walk simulations (Fig. 2b). For distances from the detector that are significantly larger than l * the agreement is very good; however, for small distances, there are clear differences, which are due to the fact that short paths are not well described by the diffusion approximation.

Figure 2

(a) Photon intensity distribution as a function of and z, obtained from Eq. (18), with s min = 1.3l *. (b) Intensity distribution obtained with random walk simulations. (c) Spatial density distribution ρ(z; s) for photon paths of various path lengths s/l * = 5, 10, 50, and 100, and overall photon density I = ∫P(s)ρ(z; s)ds compared to random walk simulations (symbols).

To investigate how this spatial distribution of photons affects the LSI signal in a dynamically heterogeneous material, we consider two different examples of materials that contain a layer of thickness d in which the diffusion coefficient of the particles D 1 differs from that in the rest of the material D 0 (Fig. 3). In case A the layer is parallel to the surface of the material, and positioned between z = d and z = 2d. This situation allows us to investigate how sensitive the LSI technique is to dynamic regions that are hidden at some depth below the surface of the material. A similar geometry has been considered previously using a slightly different approach, based on solving the correlation transport equation[31]. Here, we present an alternative expression for the autocorrelation function for this case, which we use to study the effect of τ and the ratio l */d on the measured autocorrelation function. We then consider case B, in which the layer is perpendicular to the surface of the material. This situation, which has not been studied quantitatively before, allows us to assess the effect of resolutional blurring in the lateral direction.

Figure 3

Two dynamically heterogeneous materials containing a layer of thickness d in which the diffusion coefficient of particles D 1 differs from that in the rest of the material D 0.

Dynamic heterogeneity in the z-direction (case A)

When the material is only heterogeneous in the z-direction, but not in the xy directions, the correlation function does not depend on the position of the detector and Eq. (6) can be simplified to , where ρ(z; s) is the spatial distribution of photons in the z-direction, obtained by integrating ρ(; , s) over x and y:with z 3 = z + 2z . For long paths, this can be approximated as

As shown in Fig. 2c, the photon density for diffusion paths of a given length s has a maximum around a depth z max ≈ (sl */6)1/2. Particles at this depth contribute most to the decorrelation of these paths. The total photon density, obtained by integrating over all path lengths, I(z) = ∫P(s)ρ(z; s)ds, is highest close to the surface of the material and decays as for larger z.

In Fig. 4a, we compare the field correlation function g 1(τ) for a homogeneous sample (D 1 = D 0) with that in a sample in which the diffusion coefficient in the layer D 1 is either a factor of 10 higher or lower than that in the rest of the material D 0, for l * = 2d/3. Clearly, the decorrelation is slowed down by a slowly diffusing layer, while it is accelerated by a rapidly diffusing layer. For a homogeneous sample, the correlation function is given by[1,4] with γ = (z 0 + z )/l * ≈ 1.7 and the time for particles in the medium to diffuse a distance . It follows that a plot of ln[g 1(τ)] versus τ 1/2 should give a straight line with slope . This is indeed what we find (orange data in Fig. 4b). Here, the short correlation times, corresponding to fast decorrelation, originate from the long paths with many scattering events, while the long correlation times originate from the short paths with only few scattering events. For a material with a layer of different diffusivity (D 1 ≠ D 0), we find that the initial decay follows that of the homogeneous sample. This initial decorrelation is due to the very long paths, that mostly sample the region z > 2d. For intermediate correlation times, we find a faster decay for D 1 D 0 and a slower decay for D 1 < D 0, with a slope in Fig. 4b that differs roughly by a factor . This different decay is due to the paths with a length on the order of 6d 2/l * that sample the diffusive layer. At longer correlation times, where we mainly see the short photon paths that sample the region z < d, the slope returns to that of the homogeneous sample.

Figure 4

LSI results for samples with a layer positioned between z = d and z = 2d in which the diffusion coefficient D 1 differs from that in the rest of the material D 0 (case A in Fig. 3). Lines represent theoretical results and symbols random walk simulations. (a) Field correlation functions g 1(τ) as a function of the renormalized time τ/τ 0 with , for l * = 2d/3 and different ratios D 1/D 0. (b) ln[g 1(τ)] versus τ 1/2 for the same parameters. (c) Correlation function g 1(τ) for τ/τ 0 = 6 × 10−3, 3 × 10−2, and 6 × 10−2 and D 1 = 10D 0 as a function of l */d. (d) Apparent mean square displacements obtained from the g 1(τ) values in (c) using Eq. (21).

The extent to which photons probe different depths in the sample depends on the transport mean free path l * (and thus on the wavelength). We therefore expect the measured dynamics in a material with a z-dependent diffusion coefficient to depend strongly on the value of l *. To investigate this, we calculate the correlation function for the case with D 1 = 10D 0 for different ratios l */d. The largest effect of the diffusive layer is observed when l * is on the order of the layer depth d (Fig. 4c). When l * is much smaller, the photons do not penetrate deep enough in the sample to reach the layer, while for much larger l * the photons travel much deeper in the sample and only spend a small fraction of the time in the layer. We also find that for longer correlation times τ, the optimum l * for which the effect of the diffusive layer is highest gradually shifts to higher values. Again, the reason for this is that longer decay times correspond to paths with fewer scattering events, so that a larger l * is needed for these paths to reach the diffusive layer. These findings are in agreement with previously obtained results that showed that the correlation time at which the largest effect of the diffusive layer is observed decreases with increasing depth d of the layer[31]. Usually, in the analysis of LSI data the correlation functions g 1( , τ) are converted to mean square displacements using Eq. (21), implicitly assuming that the material is homogeneous in the z-direction. Doing this, one obtains an apparent mean square displacement, which is a weighted average over the z-range probed by the photons. As shown in Fig. 4d, this apparent diffusion coefficient lies between D 0 and D 1, depending on the ratio l */d and the correlation time τ, and reaches at most 30 to 40 percent of D 1.

Dynamic heterogeneity in the lateral direction (case B)

When the material is homogeneous in the z-direction, but inhomogeneous in the xy-plane, the measured correlation function depends on the position of the detector and is given by the photon density ρ(x − x , y − y ; s), which is obtained from ρ(; , s) by integrating over z:with the exponential integral function. The lateral spreading of photons can be quantified by the second moment of this distribution,with r − r  = ((x − x )2 + (y − y )2)1/2. If the material is heterogeneous only in one direction (as in Fig. 3B), we can also integrate over y, to obtainwhich has a second moment 〈(x − x )2〉 = sl */3.

To see how this lateral spreading of photons influences the resolution of LSI in dynamically heterogeneous materials, we consider a material that contains a layer of width d perpendicular to the imaging plane in which the diffusion coefficient is 10 times higher than that in the rest of the material (case B in Fig. 3). We calculate the correlation function as a function of the detector position x using Eqs (15) and (24), and convert this into an apparent mean square displacement using Eq. (21). Figure 5 shows the apparent mean square displacements for two different correlation times and various l *. It is clear that the actual mean square displacement is more accurately followed for small l * values. As expected, when l * > d, the photons explore a region that is much larger than the width of the layer, leading to smoothing of the profile and a strongly reduced imaging contrast. However, even when l * is ten times smaller than d, the smoothing is still very significant for short correlation times (Fig. 5a), and the apparent mean square displacement measured in the layer is almost two times smaller than the actual mean square displacement. For longer correlation times, the blurring is significantly smaller (Fig. 5b); this is explained by the fact that these longer correlation times correspond to shorter paths with fewer scattering events, which sample a smaller region of the sample (Eq. (23)), therefore causing less blurring.

Figure 5

Apparent mean square displacement for samples with a layer of higher diffusivity with D 1 = 10D 0, placed between x = −d/2 and x = d/2 (case B in Fig. 3), for l */d = 0.1, 0.3, 1, and 3, and for τ/τ 0 = 3 × 10−4 (a) and 6 × 10−2 (b). Lines represent theoretical results and symbols random walk simulations. The dashed black line represents the real mean square displacements.

Discussion

We have obtained theoretical expressions for the path length distribution and spatial density distribution of photons in a strongly scattering, semi-infinite medium under plane wave illumination. These expressions allowed us to calculate the dynamic correlation function that would be measured in a backscatter LSI experiment for dynamically heterogeneous materials. We have applied these results to two simple geometries of a layered medium with different diffusivities, but our expressions can be used for any spatial distribution of diffusion coefficients. Moreover, the effects of directional flows with spatially varying shear rates can be included straightforwardly, by adding a term to Eq. (6), which accounts for the decorrelation due to spatially inhomogeneous flows[5,27].

Our model makes use of the diffusion approximation for describing the transport of photons in the medium, which is known to be accurate for sufficiently long paths, . It is not obvious that this condition holds for backscattering LSI, in which short paths contribute significantly to the signal; for example, according to Eq. (15) the path length distribution has a maximum for s ≈ l *. This is aggravated by our neglect of polarization effects, which are important for short photon paths for which depolarization does not yet occur[33]. To suppress the contribution of these short paths, we have assumed that the random walk starts at a depth z ≈ l * and we have imposed a minimum path length s min ≈ l *. Likewise, in the simulations, we only considered photons that were scattered at least twice. This corresponds to the experimental situation, where singly scattered photons are usually excluded from the analysis by using a polarization filter. With these conditions, we find remarkably good agreement between the theoretical results (with s min = 1.3l *) and the random walk simulations (Figs 2, 4, and 5). The agreement is especially good for short correlation times τ, which correspond to the long diffusion paths, with many scattering events. For longer τ, corresponding to shorter paths, the diffusion approximation loses its validity, leading to larger differences between the theoretical results and the simulations (see e.g. Figs 4c,d and 5b). Nevertheless, even in this case there is still very good qualitative agreement.

Our results show that the interpretation of LSI measurements in dynamically heterogeneous samples should be done with care. The apparent diffusion coefficient measured at a particular position can depend strongly on the transport free path l * and the correlation time τ. For larger l * and for shorter τ, the decorrelation is due to longer photon paths, which probe regions that are deeper inside the material, while for small l * and large τ the decorrelation is due to short paths that probe the regions near the surface (Fig. 4). For the same reason, the resolution of LSI in the lateral direction is smaller for large l * and small τ due to the spatial averaging that is inherent for the longer photon paths. For short correlation times, significant spatial blurring occurs even on length scales that are more than ten times larger than l * (Fig. 5). The outcome of an LSI measurement in a dynamically heterogeneous material thus depends strongly on the value of l * and τ. For a reliable interpretation, it is therefore recommended to perform measurements for a range of correlation times, and preferably also for different l *. We believe that our theory will be useful for analyzing the results of such measurements and for estimating the effects of spatial averaging, and will therefore contribute to an improved accuracy of LSI measurements.

Methods: random walk simulations

To validate the expressions obtained for the photon density distribution and the autocorrelation function, we use random walk simulations[37]. We collect statistics for 106 photons, which are launched one at a time in the +z-direction at z = 0, and allowed to perform a random walk until they leave the sample again at z = 0. The step length is sampled from a Poisson distribution with mean l *, and we assume isotropic scattering, so that the direction of each step is random. For each step we calculate the transfer wave vector , and we record the accumulated phase shift , with 〈Δr( , τ)2〉 = 6D( )τ the mean square displacement of particles at the location of scattering event i. The field correlation function g 1(τ) is obtained from this by averaging over all random walks[37]. Since we are dealing with multiple scattering, we only consider trajectories with at least 2 scattering events.