Photon-density wave correlation spectroscopy detects large-scale fluctuations in turbid media

We study the ﬂuctuations in the photon-density wave parameters (cid:64) average intensity (cid:126) dc (cid:33) , modulation amplitude, and phase (cid:35) caused by macroscopic ﬂuctuations in the optical properties of turbid media. We present both a theoretical analysis based on diffusion theory and its experimental veriﬁcation on a strongly scattering suspension containing absorbing particles (cid:126) 1–1.6 mm effective diameter (cid:33) in turbulent motion. The photon-density waves are induced by the laser diode output (cid:126) 750 nm (cid:33) , which is intensity-modulated at 110 MHz. The dc, amplitude, and phase are acquired with an acquisition time per data point of 8 ms, which corresponds to a frequency bandwidth of 62.5 Hz. We have found that in the presence of the absorbing particles, the dc and phase average values and power spectra are in good agreement with our theoretical predictions. We have veriﬁed that our instrument can extend the measured frequency band up to the kHz region, which is appropriate for the study of ﬂuctuations of optical parameters in biological tissues. (cid:64) S1063-651X (cid:126) 98 (cid:33) 01708-5 (cid:35)


I. INTRODUCTION
In recent years the correlation spectroscopy technique was used to study optically thick media that exhibit a high degree of multiple scattering ͓1-4͔. This technique, called diffusing wave spectroscopy ͑DWS͒, relates scattered light fluctuations to the motion of scattering particles. DWS was proposed as a tool for the study of microscopical particle motion ͑over a fraction of the optical wavelength͒ in multiply scattering media such as colloids and microemulsions.
Biological tissues typically display strong light scattering, and therefore can be important subjects for DWS applications. However, in addition to the small-scale fluctuations due to the motion of the microscopical scatterers, biological activity in tissues may produce large-scale spatial and temporal fluctuations of tissue optical properties in the nearinfrared band. It was demonstrated that such fluctuations in the brain may result, for example, from hemodynamics ͓5,6͔ and neuronal activity ͓7͔. In particular, it was found that visual stimulation induces local changes in the optical properties of the human brain visual area ͓7͔. It was also noted that near-infrared tissue spectroscopy has an advantage over other techniques for studying neuronal processes in that it potentially combines good temporal resolution ͑10-50 ms͒ with a spatial resolution of the order of 5 mm ͓5,7͔. Therefore, the analysis of large-scale fluctuations in near-infrared photon migration may provide a useful tool for the study of functional processes.
In this paper we present a theoretical and experimental study of a system exhibiting large-scale (ϳmm) optical parameter fluctuations. The idea is to induce photon-density waves in a turbid medium using an intensity-modulated light source. The phase velocity and the attenuation of these photon density waves depend on optical properties of the medium and on the modulation frequency ͓8-11͔. The localized areas having optical properties differing from those of the background medium affect the photon-density waves by absorption and refraction processes. Similarly to the DWS method, we relate the fluctuations in the photon-density wave parameters to the optical fluctuations in turbid media. The differences from DWS are that ͑i͒ in our case the measurable parameters are the average intensity ͑dc͒, the modulation amplitude ͑ac͒, and the phase (⌽) of the photon density wave ͑having a frequency of ϳ100 MHz͒, while in DWS the only measurable value is the light intensity, and ͑ii͒ in our case the process under study is the large-scale (ϳ mm͒ local optical property change, while in DWS it is the microscopic particle displacement (ϳ nm͒.
The purpose of our investigation is to show the potential of this approach for in vivo studies of biological systems, and to determine the instrumental capabilities. The experimental part of our work is based on frequency-domain measurements of the diffuse photon-density wave propagating in a Liposyn suspension with fluctuating optical properties. The spatial and temporal fluctuations of the scattering and absorption coefficients are caused by the motion of absorbing particles having a size of 1.0-1.6 mm. The frequencydomain parameters analyzed in this study are the dc and phase of the photon-density wave. Similarly to the standard DWS approach, we use the diffusion approximation of light transport in a turbid medium ͓8-11͔ to express the statistical characteristics of frequency domain parameters ͑average values and autocorrelation functions͒ in terms of those of the particle motion.
The paper is organized as follows. In Sec. II we describe our experimental apparatus and methods. In Sec. III we derive the expressions for the mean values and autocorrelation functions of the frequency-domain parameters based on a simplified statistical model for the motion of the absorbing particles. In Sec. IV we compare the experimental measurements with the theoretical predictions. The discussion of the results is presented in Sec. V.

II. EXPERIMENTAL SETUP AND METHOD
The central part of our experimental setup ͑Fig. 1͒ consists of a container ͑a 1-L beaker͒ filled with an aqueous *FAX: ͑217͒ 244-7187.
Liposyn suspension ͑Abbott Laboratories͒, in which small black rubber particles undergo turbulent motion caused by a magnetic-stick stirrer. The average velocity of particle motion is controlled by adjustment of the stirrer rotation frequency. In each experiment we use particles of approximately equal size, either 1 mm or 1.6 mm effective diameter.
The 750 nm light emitted by the laser diode (ϳ2 mW average power, Sharp LT030MD͒ is guided to the medium through a multimode silica optical fiber ͑source fiber͒ having a core diameter of 600 m and a numerical aperture of 0.39 ͑Thorlabs, FT-600-EMT͒. A glass fiber bundle ͑detector fiber, internal diameter: 3.2 mm, numerical aperture: 0.56, Oriel, Model No. 77527͒ collects the scattered light and conducts it to the photomultiplier tube ͑PMT, Hamamatsu R928͒. To prevent significant light blockage by particles coming very close to the fiber ends, both fibers are placed in glass tubes such that the ends are approximately in the center of the hemisphere at the bottom of the tubes, which have a diameter of 1.3 cm. The tubes are filled with the same medium as the entire container. The distance between the source and detector fibers ranges from 1.4 to 2.4 cm, which corresponds to a minimum distance between the tube walls ranging from 0.1 to 1.1 cm.
The laser current and the PMT voltage are modulated by the synchronized synthesizers at 110 and 110.005 MHz respectively. The signal from the PMT at the difference frequency ͑also called the cross-correlation frequency͒ is applied to the input of an interface card ͑ISS A2D, ISS Inc., Champaign, IL͒ for an IBM-PC computer, where the data processing is performed to obtain dc, ac, and ⌽ ͓12͔. Figure  2 shows a sample of the dc and phase temporal evolution obtained with an acquisition time of 8 ms. The left part of the FIG. 2. Time traces of the dc and phase signals before and after the stirrer is turned on. The average particle radius is 0.79 mm; the particle concentration is 0.3 cm Ϫ3 ; the magnet rotation frequency is about 2.5 Hz; the acquisition time is 8 ms per data point.
FIG. 3. Phase and dc normalized autocorrelation functions obtained from the experimental power spectra using the discrete inverse Fourier transform. Experimental conditions are the same as in Fig. 2, and the total time per sweep is 8 sec.
FIG. 1. Schematic of the experimental setup, PMT is the photomultiplier tube. The Synthesizer 3 is used as a clock for the sampling synchronization.
PRE 58 plot corresponds to the situation when the mixer is off, so that the particles rest on the bottom of the beaker, and the right part ͑after tϷ2.5 s͒ corresponds to the fluctuation regime, when the stirrer is on. One can see that when the mixer is on, both the phase and the dc exhibit fluctuations over some average value, which is different from that when the mixer is off. The time series of the dc and phase values were further processed to obtain the time-averaged values of the frequency domain parameters and their fluctuation power spectra. For the power spectrum analysis two characteristics of the data acquisition are important, namely the sampling rate ͑which is inverse of the acquisition time͒ and the total length of the time sweep. The upper limit of the sampling rate is determined by the cross-correlation frequency ͑5 kHz for our instrument͒ and is the higher the larger the frequency is. This is because the frequency-domain parameter evaluation employs the fast Fourier transform of a cross-correlation waveform, in general obtained by averaging over a number of acquired waveforms ͓12͔. For the measurements reported in this work, the acquisition time per point was 8 ms, which corresponds to 40 periods of the cross-correlation wave. All results were verified to be in good quantitative agreement with those obtained at different acquisition times from 4 to 16 ms per point. An acquisition time of 4 ͑16͒ ms corresponds to a frequency bandwidth of 125 ͑31.25͒ Hz in the power spectrum. The total time of a sweep is restricted only by the capacity of the computer memory. In our measurements it was 8.2 s, which corresponds to a frequency reso- FIG. 4. View of the beaker setup from above. Coordinate axes x and z show only the spatial directions and are not attached to any fixed Cartesian coordinate system. The vertical y axis is perpendicular to the plane of the figure.  Figure 3 shows typical dc and phase normalized autocorrelation functions ͑autocorrelation functions normalized to unity at zero time͒ obtained from the experimentally acquired power spectra by means of the inverse discrete Fourier transform. Our technique also allows us to obtain the absorption and reduced scattering coefficients of the medium from the plots of the dc and phase vs source-detector distance in the absence of particles. The details on the corresponding technique may be found in ͓13,14͔. In all the measurements reported here the reduced scattering coefficient s Ј of the Lyposin suspension was 9.8Ϯ0.1 cm Ϫ1 . The absorption coefficient a was always 0.042Ϯ0.001 cm Ϫ1 , apart from the experiment where we varied a ͑data set 3 in Table I͒. Figure 4 shows a close up view of the beaker setup and explains the geometry-related notations used below.

A. The model
We describe the signal fluctuations due to the particle motion using the diffusion approximation to the theory of light propagation in turbid media ͓8-11͔. In this approximation the light energy density as a function of time and space coordinates ͑the so-called ''photon density''͒ obeys the general diffusion equation and the particular boundary conditions specified by the light absorption and scattering at the boundaries. We assume that the rubber particles contribute to the signal change by absorbing photons from all possible paths passing through a given particle. Since the particles used in the experiment are small, one can neglect the contributions from those photon paths that include more than one particle. In this approximation the solution to the diffusion equation may be written as where UϭU dc ϩU ac is the sum of dc and ac photon densities created by the light source in the absence of particles, ␦u j ϭ␦u dc (rj គ )ϩ␦u ac (r j ,) is the sum of dc and ac contributions from the particle having coordinate r j , and N p is the total number of particles in the beaker. The simplest analytical form of ␦u j is that for a spherical particle at position r, for which it may be expressed as the sum of the dc and ac terms both having the form where and are the distances from the particle to the source and the detector, respectively ͑see Fig. 4͒. Note that as a result of the cylindrical symmetry with respect to the z axis passing through the source and detector ͑Fig. 4͒, and are the only independent coordinates entering Eq. ͑2͒, since the angle may be expressed in their terms. In Eq. ͑2͒, h l (1) (x) are the spherical Hankel functions of the first kind, Y l,m (,) are spherical harmonics, A l are the numerical constants defined by the boundary conditions ͓15͔, v is the speed of light in the medium, and S() is the source power at angular frequency ͑for the dc, ϭ0). In Eq. ͑2͒, k is the diffusion wave number,

͑3͒
where a and s Ј are, respectively, the absorption and the reduced scattering coefficients of the medium.
The constants A ł in Eq. ͑2͒ depend on the size of the object and on its scattering and absorbing properties ͓15͔. Using the expressions for A ł derived in ͓15͔ one can check that for perfectly absorbing objects having a diameter smaller than about 5 mm, the magnitudes of A ł monotonically decrease by several orders of magnitude for each successive l. Keeping only the zero-order term we have

͑5͒
In Eq. ͑5͒, a is the radius of the spherical object. Note that if one takes typical values for near-infrared a and s Ј in biological tissues ͑e.g., a ϭ0.1 cm Ϫ1 , s Јϭ8 cm Ϫ1 ), where k Ϸ1cm Ϫ1 , and since aӶ1 cm, we can set cos(k a)ϭ1 and sin(k a)ϭk a, so we obtain the approximation This is unlike the case of a nonperfect absorber, for which the zero-order term in Eq. ͑2͒ is proportional to a 3 . Note that the modulation frequency does not enter Eq. ͑6͒, so that in this approximation the coefficients A for the dc and ac parts of the signal are equal.

B. Average values of dc and phase fluctuations
If one considers a homogeneous spatial distribution of absorbing particles, and by assuming that all particles are spheres of the same radius a, the ensemble average of the dc fluctuation ͗␦U dc ͘ϭ͚͗ j ␦u(r j ),0͘ is given by the volume integral of the right-hand side of Eq. ͑4͒ at ϭ0, multiplied by the particle concentration n. Taking into account that in terms of the variables and the volume element is (2/R)dd ͑where R is the source-detector separation͒, and that the source and the detector are protected from the particles by tubes of radius ⌬, one obtains for the average dc change the expression where ␣ϭͱ3 a s Ј, or using the approximations in Eq. ͑6͒,

͑8͒
For comparison with the experimental measurements it is more convenient to use the relative average dc fluctuation ͗␦U dc ͘U dc . Assuming an infinite medium, for which U dc ϭS(0)(3 s Јe Ϫ␣R /4vR), one gets Note that the right-hand side of Eq. ͑9͒ is a linear function of R.
To estimate the average phase fluctuation we write the ac part of the signal in the form where the factored term is the ac signal with no particles. It follows from Eq. ͑10͒ that the tangent of the phase change ␦⌽ due to the particles is tan͑␦⌽ ͒ϭ Note that for particles less than a few millimeters in diameter the absolute value of A is much smaller then R. Let us also take into account the fact that the values of j and j in Eq.
͑10͒ are of the order of R/2 and greater than 0.65 cm ͑the test tube radius͒. For this reason the second term in the denominator of Eq. ͑11͒ may be neglected. This is valid for particles of small enough size if their concentration is not too large. Since in our measurements the instantaneous phase deviation values from the value in the absence of particles lie within a few degrees ͑see Fig. 2͒, one can also replace tan(␦⌽) by ␦⌽, so that finally the phase fluctuation may be expressed as the sum over the single-particle contributions ␦(r j ) as This approximation leads to the expression for the average phase fluctuation: from which one can see that ͗␦⌽͘ is a linear function of the source-detector distance R. Using the approximation ͑6͒ for A , Eq. ͑14͒ becomes

C. Autocorrelation functions
If one assumes a homogeneous particle distribution in the medium and the independence of the motion of each particle from the others, the dc autocorrelation function K dc () may be given by the expression ͓16͔ K dc ͑ ͒ϭn ͵ ͵ ␦u͑r 1 ,0͒␦u͑r 2 ,0͒W ͑r 1 →r 2 ͒d 3 r 1 d 3 r 2 ,

͑16͒
where ␦u(r,0) is the single-particle contribution to the dc change given by Eq. ͑4͒, r 1 and r 2 are the coordinates of a particle at time moments t and tϩ, respectively, W (r 1 →r 2 ) is the probability density for the particle displacement from point r 1 at time t to point r 2 at time tϩ, and n is the concentration of particles. The particular form of W (r 1 →r 2 ) is determined by the nature of the particle motion and is not known for the given environment. We assume that the vector of a particle displacement r 2 Ϫr 1 may be represented as a sum of the regular periodical component rЈ() describing the rotation in the horizontal plane around the beaker vertical axis, and a random component rЉ. For this model the probability density may be written as ͓17͔ and the integration over r 2 in Eq. ͑16͒ reduces to that over rЉ. Since the detector sensitivity to the presence of a particle is restricted to distances of about ␣ Ϫ1 Ϸ1 cm, we assume for simplicity that rЈ() evolves along the x axis from Ϫϱ to ϩϱ at each rotation period ͑see Fig. 3͒. Using these assumptions and the approximation of the function ␦u(r,0) by the series of Laguerre functions ͑see the Appendix͒ one can get the following expression for the dc normalized autocorrelation function ͓16͔: where d() is the x component of the vector rЈ(), which describes the regular periodic motion and Ϫ1/2 is the characteristic decay length of the function ␦u(r,0) along the x and y directions in the first order of the Laguerre function approximation ͓see Eq. ͑A2͔͒, and erf(s)ϭ(2/ͱ)͐ 0 s e Ϫt 2 dt is the error function.
Since Eq. ͑12͒ gives the instantaneous phase fluctuation in terms of the sum over single-particle contributions ͓Eq. ͑13͔͒ as is done for the dc, the correlation function for the phase may be estimated using the general expression of the same form as Eq. ͑16͒ for the dc: K phase ͑ ͒ϭn ͵ ͵ ␦͑r 1 ͒␦͑r 2 ͒W ͑r 1 →r 2 ͒d 3 r 1 d 3 r 2 ,

͑19͒
where W was defined above. In the Appendix, using the first-order approximation of function (r) by the series of associated Laguerre functions, we show that the normalized autocorrelation function for the phase has the form where ␤ Ϫ1/2 is the decay length for the function ␦(r) along the x and y directions ͓see Eq. ͑A5͔͒.
To discuss Eqs. ͑18͒ and ͑20͒ let us first note that both formulas have the exponential factors identical apart from the replacement of ␤ by in Eq. ͑20͒. These factors give the damped-pulse form ͑see Fig. 3͒ to the curves of autocorrelation functions, while the terms in the denominators provide the correlation decay for large . If one sets dϭ0 in Eqs. ͑18͒ and ͑20͒, the corresponding curves decay monotonically. The origin of the peaks in the autocorrelation function is the periodic component of the particle motion, while the decay is due to the random diffusion of particles. Note that enters the denominator of the exponent argument as a multiple of D x , so that this diffusion constant is responsible for the broadening of successive correlation peaks. At small ͓Ӷ(␤D x,y ) Ϫ1 ,(D x,y ) Ϫ1 ͔ the evolution of both the dc and phase curves is mostly affected by the exponential terms. Particularly, if one assumes a linear dependence of d on , these terms provide a Gaussian decay with the characteristic decay constants proportional to and ␤ for the dc and phase autocorrelation functions, respectively. Typical values of that provide the best approximation of ͐ 0 ϱ ␦u 2 (,0)d ͑see the Appendix͒ with only the zero order term in Eq. ͑A2͒ are about 3.6 cm Ϫ2 , while the optimum values of ␤ are somewhat less ͑about 3.3 cm Ϫ2 ). In other words, when a particle recedes from the probe, the probe feels the particle in the phase longer than in the dc. This is why, according to Eqs. ͑18͒ and ͑20͒, the autocorrelation function for the phase decays slower than that for the dc. For large ͓ӷ(␤D x,y ) Ϫ1 ,(D x,y ) Ϫ1 ͔ one can neglect the unity in the denominators of the right-hand sides of both Eqs. ͑18͒ and ͑20͒. Also, all terms of order less than six in the numerator polynomial in Eq. ͑20͒ can be neglected. Consequently, at large both curves decay as Ϫ1 . The factor erf͓RϪ2⌬/ͱD z ͔ in Eqs. ͑18͒ and ͑20͒ describes the correlation decay due to the diffusion along the z axis. This decay is slower than those introduced by the exponential and denominator terms because of the probe weak sensitivity to particle displacement in the z direction. One should note that Eqs. ͑18͒ and ͑20͒ essentially represent the normalized autocorrelation functions of the dc and phase responses for a single-particle random perturbation. Their proportionality to the dc and phase correlation functions, where the coefficients linearly depend on the particle concentration, is the consequence of the assumption of the statistical independence of the motion of each particle from the others. For the phase correlation function the linear dependence on n is also a result of the approximation made to obtain Eq. ͑12͒ from Eq. ͑11͒.

A. Fit of average values
Equations ͑9͒ and ͑15͒ derived in Sec. III may be compared with the frequency domain parameter time-averaged data obtained using the instrumentation described in Sec. II. For a discrete set of M source-detector distances R i we obtained three data sets for the time-averaged dc and phase values y i, j . These data sets correspond to measurements made at various ͑i͒ concentrations n j of large particles with an effective radius aϭ0.79 mm; ͑ii͒ concentrations n j of small particles with an effective radius aϭ0.50 mm; ͑iii͒ absorption coefficients a j for the constant concentration n ϭ0.6 cm Ϫ3 of small particles. The effective radius of the particles was estimated from the measurement of the volume occupied by a known large number ͑1800͒ of particles. The smaller particles were produced from the larger ones by cutting each into four pieces, so that the average radius ratio equals 4 1/3 . To vary the absorption coefficient of the medium we employed various concentrations of black India ink. The data were fitted by curves Y a,q j (R) given by Eq. ͑9͒ for the dc or Eq. ͑15͒ for the phase. The parameter a was imposed to have the same value for a given data set and corresponds to the particle effective radius. N values q 1 , . . . ,q N of the parameter q ͑where N is the number of curves in each data set͒ vary for different curves within each data set and correspond to n j for the first and second data sets and to a j for the third one. The fitting algorithm is based on the minimization of the 2 function ͓18͔, ͑where i, j are the standard deviation errors in the measurements y i, j ) with respect to the values of the parameters a,q 1 , . . . ,q N . The parameters providing the best fit are reported in Table I, while Figs. 5͑a͒-5͑c͒ illustrate the good-ness of the fit. From Figs. 5͑a͒-5͑c͒ and Table I one can see that Eqs. ͑9͒ and ͑15͒ fit the data with parameters that are in good agreement with the expected ones. A small negative curvature exhibited by the sets of data points may be related to the nonhomogeneity of the particle spatial distribution, which is larger near the beaker wall and smaller near its center. Since in our measurements the fiber nearest to the wall was fixed and the source-detector distance was changed by moving the other fiber towards the beaker center, the effective concentration of the particles was decreasing with distance. Note that while the fitted radii in the first two data sets are consistently smaller than the expected ones, the ratio of the large to small fitted radii is 1.7, close to the expected value of 4 1/3 Ϸ1.6. Figure 6 shows the experimental power spectra of the phase fluctuations obtained for three different mixing magnet rotation rates. One can see that each curve has a maximum peak at a frequency that increases with the magnet rotation rate. There are also peaks at the harmonics of these frequencies. This confirms that the particle motion may be considered as the superposition of the periodic and random components assumed for the derivation of Eqs. ͑18͒ and ͑20͒. The fact that the maximum peak frequency for each curve in Fig.  6 is more than two times smaller than the corresponding frequency of the magnet rotation is explained by the spirallike character of particle motion, such that particles return to the probe after traveling along a number of vertically stretched loops.

B. Fit of power spectra
Another assumption used to obtain the dc and phase autocorrelation functions is the statistical independence of the particle motion. To check its validity, power spectra were obtained for different concentrations of particles. From Fig.  7͑a͒ one can see that most of the spectrum increases homogeneously with concentration and approximately doubles when the concentration doubles, in agreement with the assumption made. To check this point quantitatively we plot the integral of the dc power spectrum as a function of particle concentration ͓Fig. 7͑b͔͒. Figure 7͑b͒ shows that the integrated dc power spectrum increases somewhat faster than linearly, and the divergence from a straight line becomes more significant at large concentrations ͑note that the inset, FIG. 5. ͑a͒ Average values of the phase change for different concentrations of large particles ͑0.79 mm effective radius͒; ͑b͒ average relative dc change for different concentrations of small particles ͑0.50 mm effective radius͒; ͑c͒ average phase change for different absorption coefficients of the turbid medium. The points represent the experimental data as a function of source-detector distance, while the lines are the best theoretical fits.
FIG. 6. Phase power spectra acquired at three successively increasing magnet rotation frequencies. The arrows indicate the frequencies corresponding to the peaks. The experimental conditions are the same as in Fig. 2, apart from the effective particle radius ͑0.50 mm͒. corresponding to very small particle concentrations, shows an almost perfect linearity͒. The curve deviation from the linearity may be related to the power spectrum peak at the average particle rotation frequency, which increases faster than the rest of the spectrum and becomes more pronounced at large concentrations ͓Fig. 7͑a͔͒. One can attribute this fast growth of the resonant peak to the correlations among particles, which are not included in the theoretical model. The real picture of particle motion may be more complicated than the one assumed in our theoretical model, and its detailed investigation is beyond the goals of the present paper. However, some of the correlation function characteristics predicted by the model are observed experimentally. Figure 3 shows the normalized autocorrelation functions for the dc and phase obtained from the power spectra by means of the inverse discrete Fourier transform. The behavior of these curves is in qualitative agreement with the theoretical results: ͑i͒ both curves exhibit damped pulsations with a period comparable to that of the mixer rotation rate, and ͑ii͒ the dc autocorrelation decays faster than that of the phase. To verify experimentally Eqs. ͑18͒ and ͑20͒, we fitted the experimental power spectra to the theoretical curves by adjusting the values of the diffusion constants D x ,D y ,D z . Figure 8 shows fits for the dc and phase spectra, where the theoretical curves are obtained from the discrete Fourier transform of the autocorrelation functions given by Eqs. ͑18͒ and ͑20͒. The function d() in Eqs. ͑18͒ and ͑20͒, representing the regular periodic component of particle motion, is simulated by the function L tan(⍀), where Lϭ4 cm is the radius of a horizontal circle centered at the beaker vertical axis and passing through the middle of the line joining the source and detector fibers ͑Fig. 4͒, and ⍀ is the rotation rate. Such a choice of d() provides a kind of motion in which the particle periodically recedes to infinity and returns back to the probe. To obtain the theoretical curves shown by the continuous lines in Fig. 8, all the parameters entering Eqs. ͑18͒ and ͑20͒ ͑apart from the diffusion constants͒ are set equal to their experimental values. The diffusion constant values providing the best fit are D x ϭD y ϭD z ϭ6 cm 2 /s. One can see that the theoretical curves are in good quantitative agreement with the experimental data.

V. DISCUSSION AND CONCLUSION
Data acquired with a sampling rate of 125 Hz show that our instrument is appropriate to investigate optical fluctuations in the frequency band of 0-62.5 Hz ͑processes on a time scale of 16 msec or slower͒, and over a spatial size as small as 1 mm. Reference ͓7͔ has presented measurements of temporal changes in the human occipital cortex optical parameters during visual stimulation. Since the characteristic time of the observed process is about 100 ms, the authors assign their observation to the neuronal activity. As reported by Gratton et al. ͓7͔, optical signals corresponding to specific FIG. 7. ͑a͒ dc power spectra for different particle concentrations. From bottom to top, the curves correspond to particle concentrations of 0.008, 0.016, 0.032, 0.064, 0.128, and 0.256 cm Ϫ3 . ͑b͒ Integrated dc power spectrum vs particle concentration. The inset shows an enlarged view of the boxed part of the curve in the main plot.
FIG. 8. The dc ͑a͒ and phase ͑b͒ power spectra. The points represent the experimental data, while the lines are the best fits provided by Eq. ͑18͒ ͓panel ͑a͔͒ and Eq. ͑20͒ ͓panel ͑b͔͒. The data are obtained in the following experimental conditions: particle effective radius and concentration are, respectively, 0.5 mm and 0.3 cm Ϫ3 , magnet rotation frequency is 4 Hz. In the theoretical curves the effective particle rotation frequency is 1.5 Hz. Although the frequency band investigated extends up to 62.5 Hz, the frequency axis is limited to 20 Hz because the particle motion does not contribute to the power spectrum at larger frequencies.
brain functions are very weak, having amplitude comparable to that of the underlying noise due to biological activity and the instrument noise. The instrument noise is statistically independent of the specific process under study and subtracting its spectrum from the power spectrum of the process can filter it out. Following this approach, we performed a measurement to obtain the dc and phase fluctuation power spectra using only one cross-correlation period per data point. The minimum period corresponded to a time resolution of 0.2 ms and a maximum sampling frequency of 5 kHz. In our test we collected data for the same total amount of time, but using different sampling rates. These power spectra and those obtained at a sampling interval of 8 ms essentially coincide in the frequency band 0-62.5 Hz. This test demonstrates that the maximum effective sampling rate of our instrument extends up to 5 kHz.
Theoretical calculations confirmed by measurements reveal that both phase and dc power spectra are proportional to the concentration of particles and to the fourth power of their effective radius. For particles having an effective diameter of 1 mm at a concentration of 0.3 cm Ϫ3 we measured a signalto-noise ratio ͑SNR͒ of 22 dB for the phase power spectrum at the maximum peak frequency. Using the above theoretical result to scale this value, we found that the phase SNR becomes 0 dB for particles with an effective diameter of 0.6 mm at a concentration of 0.3 cm Ϫ3 . Alternatively, it decreases to 0 dB for particles with an effective diameter of 1 mm and at a concentration of 0.03 cm Ϫ3 . The dc SNR's are typically much higher. For the same conditions at which the phase power spectrum SNR is 22 dB, the dc maximum peak SNR is 41 dB. These results give an indication of what is the smallest measurable size of particles at a given concentration ͑and in our experimental conditions͒.
Our measurements show that the phase and the dc of photon-density waves propagating through a turbid medium with fluctuating optical properties have quantitatively different statistical properties. Namely, the phase correlation function decays slower than that of the dc. This effect reflects the well-established property of the dc and phase signals to explore different spatial regions. The dc signal is very sensitive to the region close to the source and detector fibers. Instead the phase signal responds more uniformly to the region between the source and the detector.
The statistical model we have used to describe the particle motion corresponds to the so-called model of homogeneous stationary turbulence ͓17͔. We realize that this model is very rough and does not include many features of the real particle motion. However, the goal of the present work is not the investigation of the properties of the particle motion in the flow, but the demonstration of the method to study fluctuations of the optical parameters in turbid media. The values of the particle diffusion constants ͑about 6 cm 2 /s) obtained from the experimental power spectra are consistent with the observed average speed of particles ͑about 5-10 cm/s͒, which confirms the validity of the method.
We note that although the system we have studied is far from being a perfect model of any real tissue, some of its features are similar to those observed in tissues ͓7͔. For instance, the occurrence of areas having optical properties differing from those of the background medium, the temporal fluctuations in the local optical properties, and the presence of a periodic component in the process ͑provided in tissues by the heart beat͒.
In conclusion, using the diffusion approximation for light transport in turbid media we have obtained analytical expressions for the mean changes in the dc and phase due to the motion of absorbing particles. Our measurements show that these expressions are in good agreement with the experimental data. In our case, the measurement of macroscopic fluctuation of optical parameters in turbid media provides information about the underlying physical processes and their time characteristics.
The results and methods presented in this work may be used to develop a high temporal resolution technique for the study of functional processes related to blood flow and oxygenation, and to neuronal activity.

ACKNOWLEDGMENT
This work was supported by NIH Grant Nos. RR03155 and CA57032.

APPENDIX
To perform the integration in Eq. ͑16͒ we replace the function ␦u(r,0) by a suitable approximation. Due to the axial symmetry, this function depends only on two spatial coordinates, which may be chosen as the distance from the particle to the source-detector line ͑see Fig. 4͒ and the longitudinal coordinate z along this line. Using Eq. ͑4͒ one can check that the change in ␦u(r,0)ϭ␦u(,z,0) with z is only significant in the vicinity of the source or the detector. Since the source and detector are confined by the tubes, we let z ϭR/2, so that Eq. ͑4͒ becomes ␦u͑,R/2,0͒Ϸ A 4 exp͓Ϫ2␣ͱ͑R/2͒ 2 ϩ 2 ͔ ͑ R/2͒ 2 ϩ 2 . ͑A1͒ The next step is to approximate Eq. ͑A1͒ by the weighted series of N Laguerre polynomials L n () ͓19͔: A 4 ͩ ͚ nϭ1 N h n L n ͑ ͒ ͪ exp͑Ϫ/2͒, ͑A2͒ where ϭ 2 and h n ϭ ͵ 0 ϱ exp͓Ϫ2␣ͱ͑R/2͒ 2 ϩϪ/2͔ ͑ R/2͒ 2 ϩ L n ͑ ͒d.

͑A3͒
The substitution of Eq. ͑A1͒ with the series ͑A2͒ is advantageous to estimate the correlation function for the following reasons. First, one can vary the parameter , which specifies the ''width'' of the Laguerre functions, such that the first term in Eq. ͑A2͒ agrees with Eq. ͑A1͒ to within 10%. Second, in the limit of large the double integral in Eq. ͑16͒ turns into the square of an integral of the form ͐ 0 ϱ ␦u(,0)d. The approximation of a function by the series of weighted orthogonal polynomials provides the minimization of the error in the estimation of the integral of the function with respect to all possible polynomial approximations of the same order. Therefore, the series ͑A2͒ gives a good approximation for the correlation function at large . In addition, the value of the parameter may be chosen to minimize not the error in the function approximation itself, but the error in the estimation of ͐ 0 ϱ ␦u 2 (,0)d, which is the principal term entering the expression for the dc autocorrelation function at ϭ0. In this way the parameter naturally becomes dependent on the physical characteristics of the system, such as the source-detector distance R and the absorption and reduced scattering coefficients a and s Ј .
Keeping a large enough number of terms in Eq. ͑A2͒ one may calculate the correlation function to any desired accuracy. However, since we have no detailed information about the statistics of the particle motion, we can analyze the situation only qualitatively. For this reason we keep only the first term in the series ͑A2͒ to derive the expression for the dc autocorrelation function given by Eq. ͑18͒.

͑A4͒
The next step for calculating the integral in Eq. ͑19͒ is the substitution of the right-hand side of Eq. ͑A4͒ by a weighted series of orthogonal polynomials. The only difference from the dc case is that a quicker convergence to the approximated function may be obtained using the set of associated Laguerre polynomials of the first order L 1 n (␤), where ␤ is the characteristic width of the function in the right-hand side of Eq. ͑A4͒. This function is approximated by the series Rͩ ͚ n g n L 1 n ͑ ␤͒ͪ ͑␤͒exp͑ Ϫ␤/2͒, ͑A5͒ where the coefficients are g n ϭ␤⌫͑nϩ1 ͒/͓⌫͑ 1ϩnϩ1 ͔͒ 3 ͵ 0 ϱ Im"A exp͕2ik͓ͱ͑R/2͒ 2 ϩϪR/2͔͖… ͑ R/2͒ 2 ϩ exp͑Ϫ␤/2͒L 1 n ͑ ␤͒d, ͑A6͒ and ␤ is chosen to minimize the error in the calculation of the integral of ␦ 2 (x). Keeping only the first term in the series ͑A5͒, one obtains the autocorrelation function for the phase given in Eq. ͑20͒.