We demonstrate a method for fully characterizing diffuse transport of light in a statistically anisotropic opaque material. Our technique provides a simple means of determining all parameters governing anisotropic diffusion. Anisotropy in the diffusion constant, the mean free path, and the extrapolation length are, for the first time, determined independently. These results show that the anisotropic diffusion model is effective for modeling transport in anisotropic samples, providing that the light is allowed to travel several times the transport mean free path from the source.
© 2008 Optical Society of America
There is a wide range of opaque materials that show long range anisotropic optical transport. Even though these materials randomize the information in the incoming beam due to multiple scattering, the resulting transport is anisotropic. For example, a focused beam on one side of a slab of such a material may result in an ellipsoidal pattern on the opposite side. This kind of anisotropic diffusion has been observed in, for example, aligned nematic liquid crystals [1, 2], human and animal tissue (including bone , muscle , teeth , skin , arterial walls , and blood cells under shear ), and porous semiconductors [9, 10] to name a few.
In principle, the observed anisotropy may occur for several reasons. One common situation (i.e. in many tissue and tissue-like samples) is that non-spherically symmetric scatterers in the sample are statistically correlated in their alignment resulting in an angle dependent transport mean free path [11, 12]. However even spherically symmetric or randomly oriented scatterers may result in anisotropic transport if they are dispersed in a birefringent background . In this case, the energy velocity may be anisotropic. In complex materials such as liquid crystals and strongly scattering porous semiconductors  both the energy velocity and the mean free path may, theoretically, be anisotropic .
Despite the prevalence of anisotropic opaque materials and their importance in biomedical imaging, display technology, and photonics, a number of fundamental parameters that describe anisotropic diffusion have yet to be measured. Up until now, only the anisotropy in the diffusion constant tensor has been measured. There are potentially two reasons for this fact. One reason is that the difference between the diffusion tensor and the mean free path has not been well understood. While it has been understood for some time that both the mean free path and energy velocity may be angle dependent quantities and that these two quantities combined determine the anisotropic diffusion constant tensor [13, 14], it has only recently been shown that these quantities are themselves tensors [12, 15]. The second potential reason is that it has been difficult to perform the necessary measurements to determine the mean free path with different sample orientations. The standard measurement of the mean free path is total transmission of light in a slab or wedge geometry . Proper interpretation of such measurements requires separate measurements of the angular dependence of light escaping the sample, the so-called escape function, which is used to determine the extrapolation length [17, 18]. Not all anisotropic samples easily conform to a slab or wedge geometry, especially when several orientations of the anisotropy are required, or can be easily sliced and reoriented. Furthermore, the extrapolation length is likely to itself be anisotropic. These facts make total transmission and escape function measurements difficult or impossible. Given these difficulties, it would therefore be useful to combine the proper anisotropic samples with an appropriate measurement technique for performing transmission measurements with different incoming and outgoing directions to probe each axis of the anisotropy.
Theoretically, the problem of anisotropic optical diffusion has been modeled analytically using anisotropic random walk models  and radiative transfer approaches [11, 20, 12], and numerically using Monte Carlo simulations [20, 21]. There have been several points of contention in the theoretical discussion of anisotropic diffuse light. For example, the efficacy of the diffusion model itself has been recently called into question, with the proposition that only Monte Carlo simulations of the radiative transport equations provide useful models of anisotropic multiple scattering . It has also been questioned whether the diffusion constant ratio itself can be measured by any static measurement, it being a dynamic quantity .
In this paper we fully characterize anisotropic diffusion in one sample using a relatively new imaging technique and an easily available simple model (phantom) sample. This is achieved by focusing a source beam on one side a slab of the material and imaging on a perpendicular side, a technique recently analyzed and applied to isotropic samples . By repeating this measurement in two different orientations, we are able to separately determine, for the first time, the anisotropy in the extrapolation length, the transport mean free path, and the diffusion constant. These measurements also show that the diffusion theory is an extremely useful model for characterizing details of anisotropic light diffusion. We expect our approach to be useful for characterizing anisotropic light transport in a wide variety of important materials.
In the experimental geometry, shown in Fig. 1, a continuous wave light beam is focused on one side of an opaque rectangular parallelepiped. This side will be referred to as the input side. Diffuse light scattered out of the sample from a side perpendicular to the input side creates an image. This will be referred to as the image plane. The distance of the beam to the image plane is L. The dimensions of the sample are d, w 1, and w 2 where the image plane dimensions are w 1×w 2 and the input plane dimensions are d×w 2.
where U is the energy density, D̅ is the anisotropic diffusion tensor, S is the function describing the source of diffuse light, and l̅ is the anisotropic mean free path tensor. We assume the symmetry of the sample to be such that the diffusion and mean free path tensors are diagonalized in the sample coordinate system. I.e.
To describe the source function S, we reason as follows. Due to the weak focus, the incoming light can be thought of as a narrow collimated beam. This produces a source of diffuse light that is best described by a narrow cylinder that exponentially decays along the propagation direction. In the interest of deriving a simple analytical solution, following an approach used for isotropic materials, the exponentially decaying source beam of finite width is replace by a Dirac-delta function positioned a distance of lyy inside the material . Following these assumptions and referring to the coordinate system given in Fig. 1 the source function may be defined as:
We choose boundary conditions in which the energy density goes to zero one so-called extrapolation length ex,y,z from the boundary [17, 18], where the subscripts refer to the axis perpendicular to the exiting surface. The extrapolation length extends the effective sample dimensions. Thus, for example, the distance from the delta function source to the extrapolated zero position of the energy density parallel to the x-z plane is lyy+ey.
By rescaling the parameters x′=(Dyy/Dxx)1/2 x, y′=y+ey, z′=(Dyy/Dzz)1/2(z+ez), L′=(Dyy/Dzz)1/2(L+ez), and D′=(Dxx/Dzz)1/2, the anisotropic diffusion equation can be rewritten as an isotropic diffusion equation:
with the boundary conditions U=0 at x′=±∞, y′=0,+∞, and z′=0,+∞.
We have calculated that for the anisotropic case, assuming the common simplifying conditions of l′/L′≪1, ey/L′≪1, L/d≪1, and L/w 1,2≪1, the flux through the z=0 surface may be approximated as:
I(x,y,0) may be integrated over the image plane to yield the total intensity exiting the side of the sample:
The total escaping intensity is proportional to (lyy+ey)/(L+ez). This result is similar to the result for total transmission through an optically thick infinite slab, which is a well known measurement for determining the mean free path from a series of the sample of different thicknesses . Two substantial advantages of the measurements of the mean free path described here is that they can be performed on a single sample and the sample may be anisotropic.
To fit the experimental data, Eq. 6, can be rewritten as
where β=(Dxx/Dyy)1/2 and . Thus the four independent parameters A, β,ey,and L′ define the measured image completely. After fitting an image, Itot is calculated from the four parameters analytically. The parameters Itot, β,ey,and L′ describe qualitatively different aspects of the image, namely the overall intensity, the rescaling factor in the x-y plane, the I=0 line of the function (at y=-ey, outside of the image), and the extent (size) of the image, respectively. It should be emphasized that ey can be determined directly from the fit with no assumptions about the surface reflectitivity and no knowledge of the incoming intensity, a considerable advantage over escape function measurements [26, 22].
In the experiment described, measurements were taken with two different sample orientations, one rotated 90° about the z axis with respect to the other (see Fig. 1). We refer to the orientation with the beam entering along the fast axis of diffusion as the parallel orientation and the orientation with the beam along the slow axis of diffusion as the perpendicular orientation with parameters along these beam axes subscripted with the symbols ‖ and ⊥ respectively. Thus for example, D ‖>D ⊥. The symmetry of the scatterers is such that the sample may be assumed to be uniaxial, therefore the parameters D ‖,⊥, l ‖,⊥, and e ‖,⊥ will be sufficient to fully characterize the sample.
Since the incoming beam intensity, which determines S 0, and Dzz, are the same for both orientations, the ratio of the plots of Itot vs. L+ez will yield
The limit of validity of the approximation L≪d is pushed at the largest depth measured in the experiments presented here for which L=0.45d. However the approximation of an infinite sample is nearly valid, even at these depths, affecting the value of β by several percent and the total intensity by 10%. At L=0.25d, the approximation is nearly perfect. For the measurements of total intensity in this paper it is sufficient to make a second order correction in L/d for the larger values of L, using the function f(L/d)=(4(L/d)2-1)/(4(L/d)2-3), where
The anisotropic model system was chosen taking into account several important considerations. The sample should be statistically anisotropic with a substantial degree of randomness to produce opaqueness when viewed from all directions. The sample should be uniform in structure to insure reproducible measurements on separate parts of the sample and to support the assumption of position independent parameters. Finally, the sample should be easy to obtain for use in future measurements.
One anisotropic phantom sample has been previously produced, a randomly oriented collection of parallel wax fibers in a solid cube of resin . However we wanted a more strongly scattering sample. Furthermore, we wanted no scattering from the background material, and scatterers that did not extend directly through the entire slab, to avoid issues of long range light guiding .
One sample that meets our requirements is a PET/PE (Polyester/Polyethylene) porous fiber slab from Porex Corporation. These samples are used in exacting filter applications and come in a variety of densities and fiber diameters. For our measurements, we chose a fiber diameters of 10 µm and a density of 0.37 g/cm3, as listed by the manufacturer. This corresponds to a volume fraction of polymer of roughly 30% or less, depending on the exact density of the polymer (which was unknown).
Scanning electron microscopy (SEM) photos of the samples reveal a statistically aligned fiber structure with a substantial degree of randomness (Fig. 2.) The fibers are highly interconnected, reducing the probability of light guiding through a single fiber. The SEM also reveals large air voids between the fibers suggesting a low volume fraction of fibers.
The slab dimensions used were w 1×w 2=(12×42) mm2 and (42×12) mm2 for the beam oriented parallel and perpendicular to the fibers respectively. For both ordinations, d=5 mm.
The light source used was a supercontinuum broadband laser (Fianium Femptopower 1060, 20 MHz) attenuated with neutral density filters and with a bandpass filter (Thorlabs FB650-10) producing a collimated beam at 650 nm with a wide enough band width (10 nm FWHM) to eliminate speckle when imaging. A fraction of the beam intensity was split off to a reference detector to allow for normalization for changes in input beam intensity. The beam was focused to a beam spot of ~50µm onto on the input side sample (Fig. 1). The sample was mounted on a two-axis translation stage so that several locations on the sample could be probed for each value of L, via translation along the x axis, in order to determine the measurement statistics. The value of L was adjusted by hand in 250 µm increments. The light escaping the image plane was imaged using peltier-cooled CCD (Andor iXon DV885JCS) equipped with a short range (3 cm) camera lens. By varying the the input beam intensity as well as the exposure time a wide dynamic range of measurements could be detected, allowing for a wide range of L values to be measured. Iso-intensity plots of these images are shown in Fig. 3 for orientations of the sample with the fibers parallel and perpendicular to the incoming beam.
The incoming light was linearly polarized. However the images were found to be independent of both the incoming and outgoing polarization of light for these samples. For this reason, we will not refer further to polarization effects.
Once collected, the images were analyzed using a 2d nonlinear least squares fitting routine written in Matlab. Each image was fit to Eq. 9, floating the parameters, A,β,ey,and L′, from which Itot was calculated. An examples of the data and fit for a single value of L can be seen in Fig. 3. The overlap of the fit with the data is excellent.
Repeated fitting for all values of L yielded values for the parameters for all thicknesses. This process allowed us to determine the validity of the diffusion approximation as function of L.
A second simple measurement of the mean free path was also performed. Total transmission of light in the two different orientations of the slab were performed using the same light source and an integrating sphere on the outgoing (opposite) side of the sample. We used the technique described previously on isotropic porous gallium phosphide samples  using one thickness for each orientation.
Different samples were cut from the same material to achieve an infinite slab geometry, i.e. the width of the sample (i.e. the dimensions defining the plane perpendicular to the incoming beam) was large compared to the thickness. The slab dimensions used were (12×12) mm2 with thicknesses of 5 mm and 2.5 mm for the incoming beam parallel and perpendicular to the fibers respectively.
For all images taken with values of L>0.75 mm, the data is in strong agreement with the diffusion model. Fig. 3 shows an iso-intensity plot for the data and the fit to the diffusion model for both the parallel and perpendicular orientations for one sample thickness (L=2.0 mm). Fig. 4 is a 1d slice of the same data along the y axis at the x=0 (source) position. Both of these representations convey the high quality of the fit. The residual is zero over the entire area of the fit showing only a random deviation due to surface roughness of the sample. Remarkably, the fit extends even to distances near the edge shorter than the mean free path.
There are three qualitative features in these two plots worth pointing out that indicate the effects of the three parameters D̅, l̅, and ey on the data and fit. 1) The anisotropy in the contours, lengthened for both orientations along the direction of the fibers, indicates the anisotropy in the diffusion constant tensor, i.e. D ‖>D ⊥. 2) The integrated intensity with the input beam parallel to the fibers is higher, indicating the anisotropy of the mean free path tensor along this direction, i.e. l ‖>l ⊥. 3) The data and fit extend off of the edge of the sample, into the y<0 region (i.e. the lower intensity contours in Fig. 3 are not closed but come to an abrupt end at the edge of the sample, and the function in Fig. 4(inset) clearly has a y<0 intercept) indicating the effect of the extrapolation length at the input plane, which is larger along the direction parallel to the fibers, i.e. e ‖>e ⊥.
Quantitative results for the diffusion model can be found by analyzing the fit parameters as a function of L. These are plotted in Fig. 5 for all values of L measured. Multiple points at each value of L indicate multiple measurements at different translations in the x direction. The reproducibility of these measurements for sufficiently large values of L is evident from the overlap of multiple data points.
As seen in these plots, the values of β ‖,⊥ and e ‖,⊥ plateau at thicknesses greater than L>Lmin where Lmin=0.75 mm. For values of L≤Lmin, it is likely that a non-diffuse contribution from the incoming beam is still present that is not accounted for by the diffusion model. This result suggests a useful length scale for applicability of the diffusion approximation.
The values of the parameters β ‖,⊥ and e ‖,⊥ are taken from the plots in Fig. 5a and Fig. 5b by averaging the results for L>0.75 mm. The uncertainty is the standard deviation of the data points. The values for β ‖, 1/β ‖, and β ⊥at L>0.75 are 1.84±.04, 0.54±.01 and 0.53±.02 respectively. These values are consistent, to within the experimental uncertainty, to the simple relation 1/β ‖=β ⊥, and yield a value of D ⊥/D ‖=0.28±.01. The values for the extrapolation lengths e ‖ and e ⊥ are (0.34±.04) mm and (0.10±.01)mm respectively. Note that the ratio of the extrapolation lengths is e ⊥/e ‖=0.29±.03.
The relative values of the mean free paths can be determined from Fig. 5c. The data is fit to Eq. 11 using the measured value d=5 mm and the value for ez=e ⊥=(0.10±.01) mm D ⊥/D ‖=0.28±.01 determined from Fig. 5b. This yields the ratio (l ⊥+e ⊥)/(l ‖+e ‖)=0.28±.03. Inserting the values of e ‖ and e ⊥ gives l ⊥=(0.28±0.03)l ‖±0.01 mm, which, assuming that the mean free paths are not anomolously small, gives l ⊥/l ‖=0.28 to within 10%.
The total transmission measurements yielded l ⊥+e ⊥=(0.25±0.02) mm and l ‖+e ‖=(0.93±0.04) mm, yielding (l ⊥+e ⊥)/(l ‖+e ‖)=0.27±.03. Combining these results with the results for the extrapolation lengths gives l ⊥=0.15±0.02mm l ‖=(0.59±0.06) mm and l ⊥/l ‖=(0.25±0.03). These results are fully consistent with those using the imaging approach. Furthermore they allow us to express Lmin in term of the mean free path, namely, Lmin=5l ⊥=1.3l ‖.
We have determined the anisotropy in the diffusion constant, extrapolation length, and mean free path for an anisotropic phantom system. The high quality of the fits indicates that the diffusion model describes the propagation of light in these materials well.
The ratio of the quantities determined is highly suggestive. Specifically, we have found experimentally for one type of sample that
using independent measurements for each. Therefore, for this sample,
The uncertainty on these ratios is roughly 10%. Eq. 12 suggests that the extrapolation length is proportional to the transport mean free path, as is true in the isotropic case . The isotropic result for the energy velocity in Eq. 13 is reasonable for this material, since the mean free path is large compared to the wavelength, and there is little birefringence in the material.
It has been noted previously that the values of the diffusion constant and energy velocity cannot be determined from any static (single frequency, continuous wave) measurement . Our results are consistent with this statement since we make no claim to having measured the absolute values of these quantities, but rather the dimensionless ratios.
The consistent results justify, to a degree, the assumption that the effective source can be modeled as a point source at a depth that is proportional to the anisotropic mean free path in the sample. However more detailed theory and experiment are required to fully justify this assumption.
The theory developed here assumes that the microscopic symmetry of the sample is aligned with walls of the parallelepiped section. Samples in which the internal symmetry of the sample is not aligned with the parallelepiped will show asymmetric surface images. Such cases may also, in theory, be treated with the anisotropic diffusion equation . Experimental verification of such an approach is an interesting future challenge.
Having validated our approach to characterizing transport in anisotropic materials, it would be extremely interesting to perform similar measurements on systems, such as microporous semiconductors, semiconductor nanowires, or liquid crystals, in which the energy velocity is likely not to be isotropic [2, 11].
The authors would like to acknowledge helpful discussions with Bernard Kaas and Otto Muskens. Jaime Gómez Rivas is also thanked for his careful reading of the manuscript. This work is part of the research program of the “Stichting voor Fundamenteel Onderzoek der Materie”, which is financially supported by the “Nederlandse Organisatie voor Wetenschappelijk Onderzoek”
References and links
1. M. H. Kao, K. A. Jester, A. G. Yodh, and P. J. Collings, “Observation of light diffusion and correlation transport in nematic liquid crystals,” Phys. Rev. Lett. 77, 2233–2236 (1996). [CrossRef] [PubMed]
2. D. S. Wiersma, A. Muzzi, M. Colocci, and R. Righini, “Time-resolved anisotropic multiple light scattering in nematic liquid crystals,” Phys. Rev. Lett. 83, 4321–4324 (1999). [CrossRef]
3. A. Sviridov, V. Chernomordik, M. Hassan, A. Russo, A. Eidsath, P. Smith, and A. H. Gandjbakhche, “Intensity profiles of linearly polarized light backscattered from skin and tissue-like phantoms,” J. Biomed. Opt. 10, 014012 (9 pages) (2005). [CrossRef]
4. T. Binzoni, C. Courvoisier, R. Giust, G. Tribillon, T. Gharbi, J. C. Hebden, T. S. Leung, J. Roux, and D. T. Delpy “Anisotropic photon migration in human skeletal muscle,” Phys. Med. Biol. 51, N79–N90 (2006). [CrossRef] [PubMed]
8. C. Baravian, F. Caton, J. Dillet, G. Toussaint, and P. Flaud, “Incoherent light transport in an anisotropic random medium: A probe of human erythrocyte aggregation and deformation,” Phys. Rev. E 76, 011409 (7 pages) (2007). [CrossRef]
9. P. M. Johnson, B. P. J. Bret, J. G. Rivas, J. J. Kelly, and A. Lagendijk, “Anisotropic Diffusion of Light in a Strongly Scattering Material,” Phys. Rev. Lett. 89, 243901 (4 pages) (2002). [CrossRef] [PubMed]
10. B. P. J. Bret and A. Lagendijk, “Anisotropic enhanced backscattering induced by anisotropic diffusion,” Phys. Rev. E 70, 036601 (5 pages) (2004). [CrossRef]
11. B. van Tiggelen and H. Stark, “Nematic liquid crystals as a new challenge for radiative transfer,” Rev. Mod. Phys. 72, 1017–1039 (2000). [CrossRef]
12. B. Kaas, B. van Tiggelen, and A. Lagendijk, “Anisotropy and interference in wave transport: An analytic theory,” Phys. Rev. Lett.100, 243901 (4 pages) (2008). (This reference describes angle dependent mean free path and velocity vectors which are the product of the unit vector in the direction of the wave vector and the mean free path and velocity tensors respectively.) [CrossRef]
15. L. Margerin, “Attenuation, transport and diffusion of scalar waves in textured random media,” Tectonophysics 416, 229–244 (2006). [CrossRef]
17. A. Lagendijk, R. Vreeker, and P. de Vries, “Influence of internal reflection on diffusive transport in strongly scattering media,” Phys. Lett. A 136, 81–88 (1989). [CrossRef]
20. J. Heino, S. Arridge, J. Sikora, and E. Somersalo, “Anisotropic effects in highly scattering media,” Phys. Rev. E 68, 031908 (8 pages) (2003). [CrossRef]
22. J. Taniguchi, H. Murata, and Y. Okamura, “Light diffusion model for determination of optical properties of rectangular parallelepiped highly scattering media,” Appl. Opt. 46, 2649–2655 (2007). [CrossRef] [PubMed]
23. The effect of absorption has been treated for the isotropic case in reference . For the samples studied here, the expressions given sufficiently described the data without the inclusion of absorption.
24. A. Kienle, “Light diffusion through a turbid parallelepiped,” J. Opt. Soc. Am. A 22, 1883–1888 (2005). [CrossRef]
26. M. U. Vera and D. J. Durian, “Angular distribution of diffusely transmitted light,” Phys. Rev. E 53, 3215–3224 (1996). [CrossRef]
27. J. C. Hebden, J. J. G. Guerrero, V. Chernomordik, and A. H. Gandjbakhche, “Experimental evaluation of an anisotropic scattering model of a slab geometry,” Opt. Lett. 29, 2518–2520 (2004). [CrossRef] [PubMed]