Gigahertz photon density waves in a turbid medium: Theory and experiments

The predictions of the frequency-domain standard diffusion equation (SDE) model for light propagation in an infinite turbid medium diverge from the more complete P& approxI'matron to the linear Boltzmann transport equation at intensity modulation frequencies greater than several hundred MHz. The P& approximation is based on keeping only the terms 1=0 and 1= 1 in the expansion of the angular photon density in spherical harmonics, and the nomenclature P& approximation is used since the spherical harmonics of order 1=1 can be written in terms of the first order Legendre polynomial, which is traditionally represented by the symbol P&. Frequency-domain data acquired in a quasi-infinite turbid medium at modulation frequencies ranging from 0. 38 to 3. 2 GHz using a superheterodyning microwave detection system were analyzed using expressions derived from both the P& aproximation equation and the SDE. This analysis shows that the P& approximation provides a more accurate description of the data over this range of modulation frequencies. Some researchers have claimed that the P& approximation predicts that a light pulse should propagate with an average speed of c/&3 in a thick turbid medium. However, an examination of the Green's function that we obtained from the frequency-domain P, approximation model indicates that a photon density wave phase velocity of c/&3 is only asymptotically approached in a regime where the light intensity modulation frequency aproaches infinity. The Fourier transform of this frequency-domain result shows that in the time domain, the P & approximation predicts that only the leading edge of the pulse {i. e. , the photons arriving at the detector at the earliest time) approaches a speed of c /+3.

The predictions of the frequency-domain standard diffusion equation (SDE) model for light propagation in an infinite turbid medium diverge from the more complete P& approxI'matron to the linear Boltzmann transport equation at intensity modulation frequencies greater than several hundred MHz.
The P& approximation is based on keeping only the terms 1=0 and 1=1 in the expansion of the angular photon density in spherical harmonics, and the nomenclature P& approximation is used since the spherical harmonics of order 1=1 can be written in terms of the first order Legendre polynomial, which is traditionally represented by the symbol P&. Frequency-domain data acquired in a quasi-infinite turbid medium at modulation frequencies ranging from 0.38 to 3.2 GHz using a superheterodyning microwave detection system were analyzed using expressions derived from both the P& aproximation equation and the SDE. This analysis shows that the P& approximation provides a more accurate description of the data over this range of modulation frequencies. Some researchers have claimed that the P& approximation predicts that a light pulse should propagate with an average speed of c/&3 in a thick turbid medium. However, an examination of the Green's function that we obtained from the frequency-domain P, approximation model indicates that a photon density wave phase velocity of c/&3 is only asymptotically

I. INTRODUCTION
Frequency-domain methods, in which the light source intensity is modulated at radio frequency, have been successfully applied to spectroscopy studies of turbid media [1 -4]. In macroscopically homogeneous media, simple expressions for the absorption and reduced scattering coefficients were obtained in the strongly scattering regime, where the photon mean-free path is much shorter than the source-detector separation [5,6]. These expressions are based on the standard di6'usion equation (SDEj, which is valid when the photon current density J is related to the photon density Uby Fick's law [7,8]. This condition is satisfied in turbid media with relatively large scattering coefficients and at modulation frequencies up to hundreds of MHz. The refraction and diffraction of the diffuse photon density waves predicted by the SDE has been experimentally studied [9,10]. At GHz modulation frequencies, the SDE is expected to break down [11,12]. In this paper we study this GHz frequency region theoretically and experimentally. Theoretically of the first order Legendre polynomial, which is traditionally represented by the symbol P&. Expansion of the angular photon density in terms of the first I. +1 spherical harmonics is generally referred to as a PL approximation [7]. Experimentally, we have employed a frequency down-conversion technique. This technique extends the useful bandwidth of our microchannel plate detector to about 4 GHz. In the application of frequency-domain methods to the optical study of biological tissue, there are a number of reasons for considering the use of modulation frequencies in the GHz region. In spectroscopy, intrinsic tissue inhomogeneities can affect the intensity and the phase shift measurement to different extents. In fact, the inhuence of the optical properties of a specific spatial region can be substantially different in the intensity and phase measurements [13]. The use of only one parameter, be it the intensity, the demodulation, or the phase shift information, might thus be preferable due to instrument configuration considerations. To be effective in recovering the tissue optical parameters, though, measurements in a modulation frequency range extending from a few MHz up to a 1063-651X/96/53(3)/2307(13)/$10. 00 2307 1996 The American Physical Society FISHKIN, FANTINI, vandeVEN, AND GRATTON few GHz may be required. In this paper we discuss the range of optical parameters that can be recovered from the measurement of a single quantity, i.e. , the phase of the photon density wave in a range of modulation frequencies. Another reason for doing measurements at GHz modulation frequencies is to perform localized tissue spectroscopy. The use of higher Inodulation frequencies reduces the probed volume of tissue [12,14]. In imaging, it has been shown that measurements at higher modulation frequencies provide better resolution than measurements at lower modulation frequencies [14 -16].
In this work, we also discuss the limiting phase velocity of the photon density wave predicted by the frequency-domain solution of the P, approximation to the Boltzmann transport equation. We describe how this limiting phase velocity relates to the leading edge of pulse propagation in a turbid medium.

II. THEORY
In the frequency domain, the source term S' '(r, t) in Eq. (1) is given by S' '(r, t) = 6(r)S I 1+ 2 exp[i (cot + E) ] ], where 6(r) is a three-dimensional Dirac delta function located at the origin, S is the source strength (in photons per second). 3 is the modulation of the source, i =v' -1, ru is the angular modulation frequency of the source, and c. is an arbitrary phase.
Using Eq. (2) and the above mentioned assumption that 1 )) 3p, D (which means that p, ' )) p, ), we write S' '(r, t)=qo(r, t)+ -3DV' q, (r, t) . c Bt U(r, t) is the density of photons (with units of photons/volume), r is position, t is time, c is the speed of a photon in the transporting medium (i.e. , water in our experiments), 1 D=-3(p, +p, ') (2) is the diffusion coefficient (with units of distance), is the reduced scattering coefficient (with units of inverse distance), g is the average of the cosine of the scattering angle after a photon collision, and p, and p"respectively, are the inverse of the mean-free path for photon absorption and the inverse of the mean-free path for elastic scattering collisions (p, and p, have units of inverse distance). S' '(r, t) is the photon source, with qo(r, t) being the isotropic term in a first-order spherical harmonics expansion of the source term in the Boltzmann transport equation, and q, (r, t) is the first moment in this expansion which describes the dipole-like anisotropy of the source.
Equation (1) reasonably approximates the Boltzmann transport equation when r is far from sources and boundaries, and p, /(p, +p, ) is close to unity, i.e. , p, »p, [7,8]. c (p, +p, ') (Neglect of p, in the diffusion coefficient D in the SDE approximation has been previously noted by Furutsu and Yamada) [17]. We have now reached a critical point in our discussion.
In the P] approximation to the Boltzmann transport equation as represented by Eq. (1), the assumption is that p, »p, . Since a condition for the validity of the SIDE is 1 &&3p, a, the further assumption that p, ' &)p is necessary in order to reduce the I'] approximation equation to the SDE, as represented by Eq. (5). The P, approximation equation thus involves an assumption about the scattering coeKcient p"whereas the SDE involves an additional assumption about the reduced scattering coe%cient p, '. p, ' is typically an order of magnitude less than p, in biological tissue, since g is of the order of 0.9 [18]. In these conditions of highly forward scattering (g =0. 9), Fig. 1(a) shows the regions of the p, 'p, plane where one can apply the I'] approximation and the SDE. In Fig. 1(a), the conditions of applicability of I'i and SDE (p, )) p, and p. , ' »p, "respectively) are written p, & 10p, and p, & 10p"respectively. The choice of the factor 10 is arbitrary, so that Fig. 1(a) should not be interpreted rigorously. Nevertheless, it gives an idea of the values of the optical coefficients p, and p, ' for which the I'& and SI3E approximations to the Boltzmann transport equation (BTE) are applicable. We observe that in the case of isotropic scattering (g =0), the regions of the p, 'p, plane where I'] and SDE are applicable coincide. The other condition for the validity of the SDE, namely 1))3~a/c, requires that the modulation frequency~/2~be much smaller than the eA'ective scattering rate c(p, +p, ') divided by 2' This conditio. n is shown in Fig. 1(b), where we plot the product rp, ' as a function of r~/2~. The condition for the validity of the SDE (i.e. , 1))3coD/c) is written -, )3~D/c in Fig. 1(b).
The reason for having rp', on they axis of Fig. 1(b) is that we can thus plot (using a single line) the condition r ))1/p, for the validity of the I'& approximation. This condition (written as rp, & 10 or rp, ') 1 by considering again g =0.9) limits the applicability of the P, approximation to the region above the line rp, '= I in Fig. 1(b).
We will discuss this condition for the validity of the P, approximation below. Here we stress again that the factor 10 used to plot the conditions ip )3~D /c aIld rp, )10 is an arbitrary factor used to interpret the actual conditions 1))3coD/c and rp, ))l, respectively. For this reason, Fig. 1 Fig. 1(a). Figure l(b) tells us that the P, approximation is applicable for source-detector separations down to about a few mm (or about 1 cm if we consider that g can be as small as 0.6 in tissues [18]), and the SDE is applicable for modulation frequencies up to about 1 G-Hz.

B. Limiting value of the phase velocity
The first two terms in Eq. (1) have the form of a wave equation, and the first term in Eq. (1) has led to claims that the P, approximation equation predicts that a light pulse in a turbid medium propagates with an average speed of c/&3, c being the speed of a photon in the medium surrounding the scattering particles [19 -21]. We  [14]. Equation (7) yields expressions for the frequency-domain quantities measured at r relative to the corresponding quantities measured at ro, namely, the steady-state photon density (dc), the amplitude of the photon density oscillations (ac), and the phase shift of the photon density wave (4&). The relative quantities are given by ' The relative demodulation of the photon density wave is given by 53 The phase velocity V of a photon density wave in the P& approximation may be derived from Eq. (7): The ratio V /(c/&3) is platted versus modulation frequency in Fig. 2, where V is given by Eq. (15), c = 2.26 X 10' cm/s (the speed of light in water), p, '=4. 30 cm ', and p, =0.008 cm '. One can see from Zaccanti et aI. performed time-domain measurements of light pulse propagation in turbid media with optical parameters comparable to those given in Fig. 3 of this work [23]. They have compared their experimental results to the predictions of the SDE, and have stated in the conclusion of their paper that "the shapes of the pulses obtained with the difFusion approximation are similar to the measured pulses but significantly wider" [23]. Note that in Fig. 3 of our paper the pulse shape predicted by the SDE is considerably broader than the pulse shape predicted by the Pi equation. This result indicates that the P & equation may provide a better description than the SDE of the pulse propagation data presented by Zaccanti et al. We assume that the time-domain Green's function solution of the P& equation adequately describes the transport of a photon in a turbid medium when the photon undergoes a sufficient number of collisions in its transit from the source to the detector. This means that the source/detector separation r should be su%ciently large such that r )) 1/(p, +p, ) (or r ))1/p"since p, )) p, ), i.e. , r is much greater than the mean-free path of a photon collision event. In these conditions, the transmitted light pulse is adequately described by the Pi approximation in the time domain, and should thus also be adequately described by the P, approximation in the frequency domain at all modulation frequencies (the frequency-domain and time-domain expressions are related by the Fourier transformation).
Once the sourcedetector separation is sufficiently large with respect to 1/p" there is not an upper limit on the modulation frequencies at which the P, approximation equation is ade-quate. Yoo, Liu, and Alfano have determined the minimum value of the product rp, ' is of the order of 10 for adequate description of light pulse transport within the SDE [24]. Further studies are required to determine the minimum value of the product rp, for which light transport is adequately described within the P, approxi- 10" In Figs. 4 -7, r =1.8 cm, ro=1. 1 cm, and c =2.26X10' cm/s. Figure 4 shows plots of relative phase shifts versus modulation frequency, calculated for three different reduced scattering coefficients and two different absorption coefficients with both the SDE expression given by Fishkin and Gratton [14] and the I', approximation expression [Eq. (13)]. Figure 5 is derived from Fig. 4 (as well as from one set of curves at p, '=25. 80 cm ' not shown in Fig. 4). Figure 5 shows the difference in phase between the P, approximation model and the SDE model at four different reduced scattering coefficients and two different absorption coefficients for the medium. Since typical uncertainties in a measurement of N"& are of the order of 0.2, at modulation frequencies above 500 MHz we should be able to distinguish between the two models from a measurement of 4"& obtained from a medium with p, '  To experimentally study the deviations between the I', approximation equation and the SDE predictions, we used a superheterodyning microwave detection system to detect photon density waves generated in a turbid medium by a mode-locked laser. The superheterodyning technique that we employ extends the useful bandwidth of our microchannel plate detector to about 4 GHz (this frequency range is determined by the mixerrf amplifier combination used in our instrument). We describe our frequency-domain instrument in detail in Sec. III.

AND METHOD
A. Description of the experimental apparatus l. Laser system A block diagram of the frequency-domain spectrophotometer is shown in Fig. 8. The light source used in our measurements is a mode-locked neodymium-yttriumaluminum-garnet laser (Nd YAG, Coherent Antares 76-S, Santa Clara, CA) which produces a train of equally spaced light pulses with a repetition rate of 76.20 MHz and with its output frequency doubled to a wavelength of 532 nm. Each pulse is -150 ps at full-width, halfmaximum. Fourier analysis of the 76.20-MHz light pulse train yields a series of harmonic intensity-modulation frequencies with a spacing of 76.20 MHz between each harmonic [25]. A master oscillator is used as a reference for the frequency synthesizer that drives the neodymiurn-YAG laser mode locker. The same master oscillator is also used as an external reference for the microwave synthesizer (Hewlett-Packard, model 8341B), i.e. , LOCAL OSCIL 1 in Fig. 8. A personal computer controls all instrument operations, which includes positioning of the ends of optical fibers F, and Fd in the turbid medium.
The 76.20-MHz light pulse train (attenuated to -50-mW average power) is split into two parts via a beam splitter, as shown in Fig. 8. One of the beams, which contains the majority of the power of the original beam, is focused via a lens onto optical fiber F"which conveys the laser light to the turbid medium under study. The other beam is directed onto a fast photodiode (35-ps rise time photodiode from Antel Optronics Inc. , model AR-S2), which serves as the reference channel for our frequency-domain measurements [26,27].  Fig. 9.
All modern frequency-domain instruments use the cross-correlation techniques introduced by Spencer and Weber [29]. The basic element of a frequency-domain instrument is the light source, which should have a power spectrum in the modulation frequency range of the instrument's operation. Using a narrow-pulsed, modelocked laser, the power spectrum extends to about 50 -100 GHz [30]. The light detector generally limits the bandwidth of a given instrument. Obtaining the maximum bandwidth from a given detector is a crucial point.
For microchannel plate photomultiplier detectors (MCP PMT), the anode bandwidth is essentially the same as the cathode bandwidth, and the signal can be handled directly at the MCP PMT output [31].
The scheme of our electronic superheterodyning detection system is shown in Fig. 9. The output of the MCP PMT is separated into its ac and dc components. Special care is used in this part of the circuit. Since this is the only microwave frequency part of the signal processing, all components (capacitor, cables, connectors) must be rated for up to 4-GHz operation. We have verified, using a Hewlett-Packard spectrum analyzer (model 8569A), that the frequency response of this part of the circuit extends to about 4 6Hz. After the Inicrowave mixer M1 (Anzac, model MDC123), a frequency component of the signal is at low frequency (100 kHz). It can then be handled using standard electronic components. It should be noted that the average dc current from our MCP PMT must not exceed 200 nA, i.e. , approximately 2 V after the dc current-to-voltage converter (again 10 V/A). About the same signal intensity with respect to the dc signal is available at each one of the harmonic frequencies over a spectrum which extends up to the limit imposed by the detector. However, the ac part of the signal has a 50-0 impedance due to the mixer input characteristics, which typically results in a 10-pV signal. These signal levels are very low, and a direct conversion to the 40-Hz region is difficult, since the ubiquitous 60-Hz line signal dominates.
This difficulty is one of the reasons for using the intermediate conversion frequency, where amplification can be obtained in a low noise frequency region. The mixer has a conversion loss of about 7 dB (a factor of -5), which reduces the signal output from the mixer to about 2 pV. The voltage level at the rf input of the mixer was fixed at about +10 dBm. This value was found to be optimal for the Anzac mixer. Some harmonics are produced at this high rf signal level, but we found that filtering of the second harmonic of the microwave frequency carrier is excellent and the linearity of the mixer output is also good.

B. Measurement protocol in the turbid medium
The turbid medium under study was a macroscopically uniform 68-1 mixture of a 1.3 skim milk to water ratio.
The skim milkwater mixture was held in a 76 1 glass tank of dimensions 60X30X42 cm . The geometrical configuration of the detector optical fiber Fd with respect to the source optical fiber F, was such that the aperture of fiber F, directly opposed the aperture of fiber Fd in the turbid medium (Fig. 8). The ends of optical fibers F, and  Table, New Hyde Park, NY), with the position reproducibility in r to within 10 pm). At a given modulation frequency cu/2m, the phase shift 4 and demodulation M quantities measured at a distance r from the source were respectively compared to the N and M quantities measured at the r0=1. 1-cm source-detector separation. This enables us to obtain the relative quantities @"& and M"&, defined in Eqs. (11) - (14). As noted by Fishkin et al. [6], the measurement of the relative demodulation and phase shift quantities (i.e. , M"i and N"i) at a given value of co/2n has the following advantage: terms which are dependent on the source emission properties are eliminated, as are the spectral response factors of the superheterodyning phase-sensitive detection system. When fitting our frequency-domain data (@", and M"i) to Eqs. (13) and (14), respectively, to obtain the medium absorption and reduced scattering coefficients (p, and p, ') at the light wavelength 532 nm, we assume that n =1.33 for the mul- mixture of a 1:3 skim milk:water ratio. For clarity, we only show data acquired at four of the 14 modulation frequencies used in acquiring our data. The solid lines are linear fits to these data. The modulation frequencies at which we acquired our data ranged from 381.0 to 3200.4 MHz. The uncertainty in the data points is smaller than the size of the symbols representing the data. Figure 10 shows plots of the measured relative phase shift versus source-detector separation r in the 68-1 1:3 skim milk to water medium. All phase values were measured relative to the phase at the ro=1. 1 cm sourcedetector separation. For clarity, only four of the above mentioned modulation frequencies used in our experiments are shown in this plot. The solid lines are linear fits to these data. The linearity of the measured phase shifts with respect to r in Fig. 10 is consistent with the predictions of Eq. (13). Figure 11 shows plots of the natural logarithm of the measured relative demodulation versus r. As in Fig. 10, we only show four of the 14 modulation frequencies. The solid lines are linear fits to these data. The linearity of the data with respect to r at the three lowest modulation frequencies in Fig. 11 is consistent with the predictions of Eqs. (11), (12), and (14). detector are immersed in a 68-1 mixture of a 1:3skim milk:water ratio. For clarity, we only show data acquired at four of the 14 modulation frequencies used in acquiring our data. The solid lines are linear fits to these data. The modulation frequencies at which we acquired our data ranged from 381.0 to 3200.4 MHz. The uncertainty in the data points is smaller than the size of the symbols representing the data.
The data acquired at 3200.4 MHz in Fig. 11 are not perfectly linear with respect to r. Rather, they show a slight downward curvature. This result is inconsistent with the predictions of both the P& approximation equation and the SDE.  Table I, where the correlated uncertainties given in this table were calculated using 0.98 probability limits for the confidence interval). These fits were obtained using the commercial package GLoBALs UNLtMITED (Laboratory for Fluorescence Dynamics, Dept. of Physics, University of Illinois at Urbana-Champaign) [34]. The residues of the SDE fit of Fig. 12(a) and the P, approximation equation fit of Fig.   13(a) are shown in Figs. 12(b) and 13(b), respectively.
The residues of the SDE fit show strong systematic trends, and most of these residues are beyond the range of uncertainty of our measurements (i.e., b, C&"t=+0. 2' and bM"&=+0.004). This indicates that the SDE does not provide an accurate description of the data over the range of modulation frequencies studied. The residues of the P& approximation equation fit show much less of a systematic trend, and are for the most part within the range of uncertainty of our measurements.
To further verify the relative accuracy of the P, approximation equation as compared with the SDE, we have also extracted the optical parameters p, ' and p, TABLE I. Optical parameters of the 68-l 1:3 skim milk:water solution obtained by fitting the frequency-domain data obtained from this medium to the I'& equations and the SDE. All of the data fits were performed using 0.98 probability limits for the confidence interval. Table I Table I Table I). The fitted parameter p, ' is relatively uncorrelated with the fitted parameter p, in all of the P, approximation fits to our data. The p, ' and p, values obtained from the SDE fits to our data are more correlated than those obtained from the Pi approximation equation fits.  Table I). ( Fig. 4(a), p, is of the order of 0.01 cm '. The maximum curvature in the phase shift with respect to modulation frequency in Fig. 4(a) occurs roughly between 20 and 500 MHz, whereas the phase shift in this figure is relatively linear at modulation frequencies greater than 500 MHz and at modulation frequencies less than 20 MHz. In Fig. 4(b), p, is of the order of 0.1 cm '. The maximum curvature in the phase shift with respect to modulation frequency in Fig. 4 . Smaller values of p, ' yield larger deviations between the two models at a given modulation frequency. Figure 1(b) shows an rp, ' versus.
re/2m diagram of the regions of validity of the SDE and Pi approximations.
At higher modulation frequencies and relatively small reduced scattering coefficients, the P, approximation should be used for the description of the propagation of light in turbid media. However, for relatively large reduced scattering coefficients and/or low frequencies, the SDE approximation can be adequate for the description of light propagation in turbid media. We note that at very small values of the scattering coefficient, neither approximation is adequate and higher order approximations to the Boltzmann transport equation may be required to describe photon migration in the lower scattering regime adequately. It is also important to realize, as shown in Fig. 1(b), that the distance between the source and detector is another important factor for the applicability of the P& approximation. The smaller the scattering coefficient, the larger the source-detector separation should be for the applicability of the P, approximation.
The frequency-domain P & approximation equation predicts that as the modulation frequency approaches, infinity, the phase velocity of the photon density wave V asymptotically approaches c/&3 (see Fig. 2). We have plotted the Fourier transform of our frequency-domain solution to the P& approximation equation in Fig. 3 using the same optical parameters that generated Fig. 2 [22]. Examination of these time-domain calculations shows that at a distance r from the light source, the earliest arriving photons generated by the pulsed source 5' '(r, t)=5(r)6(t) arrive at a time t;"=r/(c/v'3) according to the P& approximation equation. This result contradicts previous claims [19 -21] Table I) indicates that the Pi approximation equation gives a more accurate and self consistent description of our data than the SDE. However, there is still a systematic trend in the residues obtained from the P& approximation equation fits, which may indicate that the Pi approximation equation does not provide a complete description of our data. Boas et al. [35] have derived higher order frequency-domain P3 approximation expressions which may provide a better representation of our observations than the P3 approximation expressions.
The conditions for the applicability of higher order approximations to the Boltzmann transport equation for describing photon migration in turbid media is an important problem which may be of use to the community of researchers who use dift'using light to image turbid media.
In this paper, we have studied the physical meaning of the P, approximation. The mathematical limit of the higher order Pz approximations to the Boltzmann transport equation has been established, but the physical picture of the process described by the addition of higher order terms has not been investigated, as yet. In any case, our studies of the g surface calculated from the phase only indicate that choosing an appropriate modulation frequency range is critical for extracting an accurate value of the absorption of the medium. The range of absorption coefficients that can be well determined by a measurement of the phase scales with the modulation frequency range utilized: the larger the absorption coefficient, the higher the optimal modulation frequency range for determining the absorption coefticient of the medium.

A.CKNOWLEDGMENTS
These experiments and analyses of the data produced were performed at the Laboratory for Fluorescence Dynamics (LFD) in the Department of Physics at the University of Illinois at Urbana-Champaign (UIUC). The LFD and this work is supported by the National Institutes of Health (RR03155 and CA57032) and by UIUC.