Effective source term in the diffusion equation for photon transport in turbid media.

The Green’s function for the diffusion equation is widely used to describe photon transport in turbid media. We have performed a series of spectroscopy experiments on a number of uniform turbid media with different optical properties (cid:126) absorption coefﬁcient in the range 0.03–0.14 cm (cid:50) 1 , reduced scattering coefﬁcient in the range 5–22 cm (cid:50) 1 (cid:33) . Our experiments have been conducted in the frequency domain, where the measured parameters are the dc intensity (cid:126) I dc (cid:33) , ac amplitude (cid:126) I ac (cid:33) , and phase (cid:126)(cid:70)(cid:33) of the light intensity wave. In an inﬁnite medium, the Green’s function predicts a linear dependence of ln (cid:126) rI dc (cid:33) and (cid:70) on the source–detector separation r . Our measurements show that the intercepts of these straight lines predicted by the Green’s function do not agree with the experimental results. To reproduce the experimental results, we have introduced an effective photon source whose spatial extent and source strength depend on the optical properties of the medium. This effective source term has no effect on the slopes of the straight lines predicted by the Green’s function at large values of r . © 1997 Optical Society of America words:


Introduction
The properties of light propagation in scattering media are important in many areas of physics and engineering such as remote sensing of the atmosphere, 1,2 studies of interstellar dust, 3 industrial production monitoring, 4 and medical diagnostics. 5 The Boltzmann transport equation describes the photon transport in the presence of scattering. In the strongly scattering regime, the Boltzmann transport equation reduces to the diffusion equation, which also describes heat conduction in solids, diffusion of neutrons in condensed media, and gas diffusion. Nearinfrared ͑NIR͒ light is strongly scattered inside most biological tissues, so the diffusion equation has been widely employed in studies focused on NIR imaging and spectroscopy of living tissue. These applications to medical diagnostics have raised much interest in recent years because of a number of promising preliminary results and because of their significant potential. 6 Optical measurements in biological tis-sues can be quantitative. In other words, they can provide absolute values of the optical coefficients that characterize the tissue. These coefficients are the absorption ͑ a ͒ and the reduced scattering ͑ s Ј͒ coefficients, both in units of inverse centimeters. It has been known for a long time that absorption properties of blood-perfused tissues can be related to hemoglobin content and saturation, 7 and recent studies have suggested that scattering properties are due to the mitochondrial compartment 8 and may be related to blood glucose concentration. 9,10 To exploit the capability of optical methods to be quantitative, it is necessary to have an accurate expression that describes light propagation in tissues. The common expressions are based on the Green's function for the diffusion equation and thereby assume that the light source is pointlike. The fact that light propagation is not diffusive near the source suggests that, to solve the diffusion equation, one should consider an effective photon source characterized by a spatial extension and a source strength that could depend on the optical properties of the medium. In this work, we have investigated the properties of the effective source ͑spatial extension and strength͒ that are consistent with the experimental results. To separate the effects of spatial extension and effective source strength on the detected signal, we have carried out a study in the frequency domain. This approach employs a light source that is intensity modulated at radio frequency ͑120 MHz in this study͒. The directly measured parameters are the average intensity ͑dc͒ and the amplitude ͑ac͒ and phase ͑⌽͒ of the light intensity wave. The phase ⌽ is independent of the source strength, so it can be employed to investigate specifically the spatial extension of the effective photon source. The source strength can then be studied with the dc or ac component of the intensity. In this study, we have considered the phase and the dc intensity, and we have employed optical fibers for both delivering and collecting the optical signal. Our measurements are conducted in a liquid turbid medium with optical properties similar to those of biological tissue in the NIR. In the theory section, Section 2, we also comment on the measurable quantity in our experiment.

A. Measurable Quantity
An optical fiber of numerical aperture N collects photons whose directions of propagation are within the solid angle 2 ͐ 0 arcsinN sin d. Within this solid angle, the photon current dI collected in d⍀ around direction ⍀ by a fiber of area ⌬A located at r is given by where v is the speed of light in the medium, u͑r, t, ⍀ ͒ is the angular photon density ͓hvu͑r, t, ⍀ ͒ is the radiance, where h is Planck's constant and is the light frequency͒, and n is the unit vector normal to the surface of the fiber tip and pointing inside the fiber. The total collected current is where the polar axis is in the direction of n . If we assume that the presence of the optical fiber does not perturb the angular photon density distribution, Eq. ͑2͒ shows that the measurable quantity is the radiance integrated over the solid angle of collection of the fiber. In the diffusion approximation, the radiance is given by the sum of an isotropic term and a smaller directional flux: where U͑r, t͒ is the photon density ͓hvU͑r, t͒ is the fluence rate͔, and D is the diffusion coefficient defined as 1͑͞3 s Ј͒, where s Ј is the reduced scattering coefficient. If we define F ϭ arcsin N ͑ F Ͻ ͞2͒ to be the acceptance angle of the detector fiber and if U ϰ ͓exp͑Ϫkr͔͒͞r, then the photon current collected by a fiber pointing toward the source ͑n ϭ r͒ is given by The ac amplitude I ac is given by the amplitude of I ͓and the dc intensity is I dc ϭ I ac ͑ ϭ 0͔͒, and the phase ⌽ is given by the argument of I. In the spherically symmetric case, ٌU is directed along r. Consequently, the flux term has its maximum value along r, whereas it is 0 in directions perpendicular to r. In other words, if the detector fiber points toward the source ͑n ϭ r͒, the radiance is and in directions and the radiance is and is proportional to the photon density. Because in the diffusion approximation the directional flux is small compared with the isotropic term, we can conclude that the detected signal is approximately proportional to the photon density U͑r, t͒ also in the case n ϭ r. This discussion is in agreement with Haskell et al., 11 although it deviates from the conclusions of other authors who consider the photon flux ͑ϪvDٌU ⅐ n ͒ ͑Ref. 12͒ or the photon current density 13 to be the measurable quantity. The photon flux describes the net flow of photons through a given surface, but an optical fiber is sensitive only to part of the solid angle, i.e., to a partial current. Liu et al. 14 have also experimentally shown that for the case n ϭ , the measured quantity is proportional to the photon density rather than to the photon flux. As a conclusive remark, we also note that in the case n ϭ the photon flux is 0, whereas the experimentally measured signal is comparable with that measured in the case n ϭ r.

B. Diffusion Equation
The diffusion equation for photon transport in a macroscopically homogeneous turbid medium is 15 where q͑r, t͒ is the source term. The frequencydomain Green's function for an infinite medium is given by 16 where k 2 ϭ ͑v a Ϫ i͒͑͞vD͒ and is the angular modulation frequency of the source intensity.
The solution to the diffusion equation ͓Eq. ͑7͔͒ for a spatially extended photon source, described by the source term q͑r͒, is given by the convolution between q͑r͒ and the Green's function: In spherical coordinates, with the polar axis along r, the convolution integral becomes The photon densities for a number of spherically symmetric photon sources q͑r͒ are reported in Appendix 1.
The anisotropic emission of an optical fiber with a large numerical aperture ͑ϳ1͒ can be modeled by an effective photon source spatially distributed within a hemisphere of radius ␥ centered at the fiber tip. In this way, we account for the directional emission in half of the solid angle ͑see Fig. 1͒. The scattering centers inside the medium act as effectively isotropic point sources, as has already been suggested, 17 so that we can apply Eq. ͑10͒ to find U hs ͑r, ͒ ͑where the subscript hs stands for hemisphere͒. The source term q͑rЈ, Ј, Ј͒ is written as where S eff ͑͒ is the total effective source strength, i.e., the number of photons per unit time emitted in the whole volume of the hemisphere. The result for U hs ͑r, ͒ at r Ͼ ␥ and ϭ 0 ͑i.e., along the fiber axis͒ obtained by the computation of the integral in Eq. ͑10͒ is We note that in the limit ␥ 3 0, U hs ͑r, ͒ reduces to S eff ͑͒U G ͑r, ͒ as it should. It is also noteworthy that in the limit r 3 ϱ, or r 2 Ͼ Ͼ ␥ 2 , U hs ͑r, ͒ becomes proportional to U G ͑r, ͒ with a constant of proportionality independent of r. The r dependence of U hs ͑r, ͒ at r 2 Ͼ Ͼ ␥ 2 is thus the same as that of the Green's function, i.e., exp͑Ϫkr͒͞r.

A. Strongly Scattering Material
The strongly scattering material employed in our study is a 9-L liquid suspension containing water, Liposyn 20%, and black India ink. Liposyn is a fat emulsion that acts as a weak absorber and a strong scatterer for NIR light. The reduced scattering coefficient of the medium is proportional to the concentration ͑v͞v͒ of Liposyn, where the proportionality factor is ϳ2.2 cm Ϫ1 ͞% ͑v͞v͒ ͓or, in terms of percentage of solids content ͑s.c.͒, ϳ11 cm Ϫ1 ͞% ͑s.c.͔͒. In our study, we have investigated reduced scattering coefficients ranging from 5 to 22 cm Ϫ1 ͑with a ϭ 0.035 cm Ϫ1 ͒. Black India ink is employed to vary the absorption coefficient in the range 0.035-0.14 cm Ϫ1 ͑with s Ј ϭ 16 cm Ϫ1 ͒. The 9-L strongly scattering suspension is placed in a rectangular container of sides 32 cm ϫ 26 cm ϫ 12 cm.

B. Frequency-Domain Spectrometer
The light source is a laser diode emitting at 780 nm ͑Sharp LT023͒. Its output, which is amplitude modulated at a frequency of 120 MHz, is coupled to a plastic optical fiber with a core diameter of 2 mm and a numerical aperture of 0.5 ͑ F ϭ 30°͒. Intensity modulation of the laser diode is accomplished by supplying a dc current of 50 mA ͑which selects a working point just above threshold͒ and a superimposed 120-MHz oscillating signal. As a result, the average emitted power is ϳ3 mW and the modulation depth is close to 100%. The radio frequency signal is provided by a frequency synthesizer ͑Marconi 2022A͒ that is phase locked to a second frequency synthesizer ͑Marconi 2022A͒. This latter synthesizer modulates the gain of the photomultiplier tube ͑PMT͒ detector ͑Hamamatsu R928͒ at a frequency of 120 MHz ϩ 400 Hz. The signals from both frequency synthesizers are amplified by radio frequency amplifiers ͑ENI Model 403LA͒. The beating between the detected signal at 120 MHz and the response function of the PMT at 120 MHz ϩ 400 Hz gives rise to a harmonic component at 400 Hz in the PMT output signal. This low-frequency component is processed to yield the frequency-domain parameters dc ͑average intensity͒, ac ͑amplitude of the intensity wave͒, and ⌽ ͑phase of the intensity wave͒ by the use of previously described methods. 18

C. Measurement Procedure
The experimental arrangement is shown in Fig. 2.
The tips of the source and the detector fibers are placed deep in the suspension such that we work in an effective infinite geometry. The two fibers face each other as shown in Fig. 1. We used two identical 9-L containers, one for media of different absorption coefficients, and one for media of different reduced scattering coefficients. We investigated seven media with a ranging from 0.035 to 0.14 cm Ϫ1 ͑and with s Ј ϭ 16 cm Ϫ1 ͒ and five media with s Ј ranging from 5 to 22 cm Ϫ1 ͑and with a ϭ 0.035 cm Ϫ1 ͒. For each medium we measured the intensity signal at sourcedetector separations ranging from 0 ͑fiber tips in con-tact͒ to 2 cm, at 0.1-cm increments. To perform the measurement in the required intensity dynamic range ͑ϳ10 5 ͒, we employed neutral density filters at the PMT detector ͑whose high voltage supply was kept constant throughout the experiment͒.

A. Determination of a and s Ј
In a previous publication, 19 we proposed a multidistance measurement protocol for the absolute determination of a and s Ј in optically dense media. This protocol consists of acquiring data at a number ͑at least two͒ of different source-detector separations, from which one can determine the slopes of the straight lines ln͑rI dc ͒, ln͑rI ac ͒, and ⌽ as a function of r. The Green's function of Eq. ͑8͒ provides expressions for these slopes that enable one to calculate a and s Ј. As we comment in Appendix 1, a spherically symmetric light source that is confined at r Յ ␥ does not change the expressions for these slopes at r Ͼ ␥. The spherically asymmetric source distribution of Eq. ͑11͒ also does not change the expression for these slopes at r 2 Ͼ Ͼ ␥ 2 . An important consequence is that a measurement of a and s Ј based on the slopes of ln͑rI dc ͒, ln͑rI ac ͒, and ⌽ at large r is not affected by the details of the source spatial distribution. We have thus determined a and s Ј of our media by applying the multidistance measurement protocol to the data collected at r between 1.5 and 2.0 cm.

B. Spatial Extension of the Effective Photon Source
Because the slopes of the straight lines discussed in Section 3 are not affected by the source spatial distribution and by the source strength, we have focused our attention on the intercepts. In particular, the phase intercept is not affected by the source strength, so that the phase intercept data give information specifically on the source spatial extent. Figure 3 shows the experimental phase intercepts as a function of a ͓Fig. 3͑a͔͒ and s Ј ͓Fig. 3͑b͔͒. We observe that the phases are relative to the phase measured at r ϭ 0, i.e., with the source and the detector fiber tips in contact. For comparison, we also report the phase intercept predicted by the integrated radiance ͓Eq. ͑4͔͒ calculated with F ϭ 30°by using U G ͑r, ͒ ͓Eq. ͑8͔͒ and U hs ͑r, ͒ ͓Eq. ͑12͔͒. In U hs ͑r, ͒ we set (13) where s0 Ј ϭ 16 cm Ϫ1 is a constant. This value of ␥, which reduces to 1͑͞3 a s0 Ј͒ 1͞2 for the dc intensity, is of the order of the diffusion length and ranges from 0.4 to 0.8 cm for the media considered in this study. The diffusion length characterizes the spatial scale of the photon density decay. It also gives an estimate of the distance from the actual source over which the diffusive regime is established. This distance is larger than the photon mean free path between effectively isotropic scattering events, which is given by 1͞ s Ј. We note here that ␥, the extension of the effective source, is a function of a and , but surprisingly not of s Ј ͓ s0 Ј in Eq. ͑13͒ is a constant͔. The constant value of s0 Ј ϭ 16 cm Ϫ1 used in Eq. ͑13͒ provides a good agreement between our model and experiment, even in the case of media that have different values of s Ј ͓see Fig. 3͑b͔͒. We also observe that U hs ͑r, ͒ with ␥ given by Eq. ͑13͒ accurately describes not only the dependence of the experimental phase intercept on a and s Ј, but also its absolute value. Fig. 2. Experimental arrangement. The laser diode ͑LD͒ emitting at 780 nm is intensity modulated at 120 MHz by a frequency synthesizer ͑Synth 1͒ by means of amplifier A1. The DC Bias selects a working point for the laser diode just above threshold. The tips of the source and the detector optical fibers are immersed in the sample to mimic an infinite geometry. The two fibers face each other, and the volume of the sample is ϳ9 L. The gain of the PMT detector is modulated at 120 MHz ϩ 400 Hz by synthesizer 2 ͑Synth 2͒ by means of amplifier A2. Neutral-density filters ͑F͒ are used to attenuate the detected signal at smaller source-detector separations. The data acquisition card in the computer and the two frequency synthesizers are phase locked ͑Synch͒.

C. Source Strength of the Effective Photon Source
Once the effective source extension has been determined by Eq. ͑13͒, the intercepts of the intensity straight lines provide information on the effective source strength. Figure 4 shows the experimental intercepts of ln͑rI dc ͒ as a function of a ͓Fig. 4͑a͔͒ and s Ј ͓Fig. 4͑b͔͒. For comparison, we also show the intercepts predicted by the integrated radiance ͓Eq. ͑4͔͒ calculated with F ϭ 30°by using U G ͑r, ͒ ͓Eq. ͑8͔͒ and U hs ͑r, ͒ ͓Eq. ͑12͔͒. In U hs ͑r, ͒ we have used the expression of Eq. ͑13͒ for ␥ and we have set where S is the actual source strength. This effective source strength is compatible with our experimental results for the intensity intercepts ͑see Fig. 4͒. In contrast to the phase intercepts, the absolute values of the experimental and the predicted intensity intercepts cannot be compared because of the unknown effect of the instrumental response. For this reason, the predicted intensity intercepts in Fig. 4 are shifted by an arbitrary offset. The fact that S eff decreases with s Ј is reasonable. In fact, the source dimension ␥ does not change with s Ј, and the number of photons emitted per unit time by the effective source is reduced by the scattering-induced absorption. The increase of S eff with a is apparently surprising. It is the net result of two competing effects determined by an increase in the absorption coefficient. The first effect is the reduction in the number of photons emitted per unit time by the effective source because of the increased absorption. The second effect is the Fig. 3. Intercept of the phase-versus-r straight line as a function of ͑a͒ a , ͑b͒ s Ј. In ͑a͒, s Ј is 16 cm Ϫ1 , whereas in ͑b͒, a is 0.035 cm Ϫ1 . The symbols are the experimental values obtained from the phase data taken at r ranging between 1.5 and 2 cm. The curves are the phase intercepts predicted by diffusion theory for a point source ͑dashed curve͒ and for the effective extended source described in the text ͑solid curve͒. reduction in the source dimensions ͑␥ ϰ 1͞ ͌ a ͒, which confines the effective source volume closer to the actual photon source. This second effect reduces the probability of photon absorption within the effective source volume. The inclusion of the following effects, ͑1͒ integrated radiance over the acceptance angle ͑Ͻ͞2͒ of the optical fiber, ͑2͒ hemispheric light source with radius ␥ given by Eq. ͑13͒, and ͑3͒ effective source strength given by Eq. ͑14͒, results in the following expression for the experimentally measured intensity at large r ͑r Ͼ Ͼ ␥͒ and ϭ 0:

Discussion and Conclusion
The experimental results shown in Figs. 3 and 4 show the inadequacy of the diffusion-equation Green's function in describing the dependence of the intercepts of dc intensity and phase on a and s Ј. Such inadequacy cannot be resolved by modeling the light source ͑a 2-mm-diameter optical fiber in our case͒ with an extended photon source of given dimensions. On the contrary, the experimental results can be reproduced by assuming that there is an effective photon source whose spatial extent and strength depend on a and s Ј. We have found that the dimensions of the effective photon source ␥ is of the order of the diffusion length, which represents the distance over which the photon migration becomes diffusive. Figure 5 shows the different r dependencies of the experimental data at r Ͻ ␥ and r Ͼ ␥ for the case a ϭ 0.035 cm Ϫ1 and s Ј ϭ 18 cm Ϫ1 . In this case, the effective source dimension ␥ is 0.71 cm. Because the diffusion-equation Green's function is accurate only at r Ͼ ␥, the different regime at r Ͻ ␥ explains the failure of the Green's function in predicting the dc and phase intercepts. This observation has important consequences in quantitative spectroscopy of turbid media. The measurement of the slopes of the dc and phase straight lines versus r provides accurate values of a and s Ј ͑Ref. 19͒. However, to measure these slopes, one needs to acquire data at multiple source-detector separations at which different volumes are sampled by the photon paths. In applications in which the sample is not macroscopically homogeneous, this fact can represent a problem. Consequently algorithms that employ a calibration measurement followed by one measurement at a single source-detector separation have been proposed. 20 A similar algorithm has been proposed in steadystate spectroscopy, 21 in which the slope and intercept of the dc straight line are used to determine a and s Ј. These algorithms are based on the diffusionequation Green's function and employ the intercepts as well as the slopes of the straight lines associated with the measured quantities. As a result of this study, these algorithms are not expected to provide accurate results, unless the calibration medium and the sample have the same optical properties.
The properties of the effective photon source ͑␥ and S eff ͒ have been empirically derived here to provide an accurate description of our experimental results. However, we have not developed a model starting from first principles, which explains the physical origin of the discrepancy between the experimental data and the predictions of diffusion theory. It is possible that the acquisition geometry ͑relative orientation of the source and the detector fibers͒ and the numerical apertures of the optical fibers play a role in our observations that has not been quantified. How- Fig. 5. Symbols represent ͑a͒ the measured phase, ͑b͒ ln͑rI dc ͒ as functions of source-detector separation r. The medium has an absorption coefficient of a ϭ 0.035 cm Ϫ1 and a reduced scattering coefficient of s Ј ϭ 18 cm Ϫ1 . Consequently the radius of the hemispheric effective photon source is ␥ ϭ 0.71 cm. The continuous curves have the slopes predicted by diffusion theory and are shifted by an offset to match the experimental data at r Ͼ ␥. At r Ͻ ␥ the experimental points deviate significantly from the linear behavior predicted by diffusion theory. ever, the effects that we observed are so large that they cannot be ignored. Other researchers have observed similar effects ͑nonzero phase intercepts in a plot of phase versus source-detector separation͒ in a geometry in which the source and detector fibers were not facing each other. 11 Their phase intercepts at 100 MHz were larger than ϩ10°, and this result could not be explained.
In conclusion, this work points out an important limitation of the Green's function solution ͑i.e., the solution for a point-like source͒ to the diffusion equation for measurements of photon transport in turbid media. More work is required for investigating the physical origin of this limitation and to assess the influence of the specific source-detector configuration used to measure the optical parameters of a turbid medium.