1. Introduction

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 [3], muscle [4], teeth [5], skin [6], arterial walls [7], and blood cells under shear [8]), 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 [12]. In this case, the energy velocity may be anisotropic. In complex materials such as liquid crystals and strongly scattering porous semiconductors [9] both the energy velocity and the mean free path may, theoretically, be anisotropic [11].

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 [16]. 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 [19] 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 [21]. It has also been questioned whether the diffusion constant ratio itself can be measured by any static measurement, it being a dynamic quantity [2].

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 [22]. 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.

2. Theory

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 ₁, and w ₂ where the image plane dimensions are w ₁×w ₂ and the input plane dimensions are d×w ₂.

The diffuse light is assumed to propagate through the sample in the steady state limit with no measurable absorption [23] according to the diffusion equation [11, 12]:

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.

 

Fig. 1. Cartoon of the imaging geometry. A beam is weakly focused on one (input) side of the sample. The emitted light is imaged on a side perpendicular to the input side (the image plane). The focus point is a distance L from the image plane. The focus spot size is much smaller than L, justifying a point source approximation. The sample is chosen such that w ₁,w ₂≫L, i.e. the slab is effectively semi-infinite in the image plane.

Download Full Size | PDF

and

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 l_yy inside the material [25]. 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 e_x,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 l_yy+e_y.

By rescaling the parameters x′=(D_yy/D_xx)^1/2 x, y′=y+e_y, z′=(D_yy/D_zz)^1/2(z+e_z), L′=(D_yy/D_zz)^1/2(L+e_z), and D′=(D_xx/D_zz)^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,+∞.

The flux I(x,y,0) through the z=0 surface for light governed by Eq. 5 has been solved analytically by Kienle in his treatment of isotropic media [24].

We have calculated that for the anisotropic case, assuming the common simplifying conditions of l′/L′≪1, e_y/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 (l_yy+e_y)/(L+e_z). 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 [16]. 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 β=(D_xx/D_yy)^1/2 and $A=\frac{3}{\pi }{S}_0(L+e_z)(l_{\mathrm{yy}}+e_y)D_{\mathrm{yy}}{\left(D_{\mathrm{xx}}{D}_{\mathrm{zz}}\right)}^{\frac{-1}{2}}$ . Thus the four independent parameters A, β,e_y,and L′ define the measured image completely. After fitting an image, I_tot is calculated from the four parameters analytically. The parameters I_tot, β,e_y,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=-e_y, outside of the image), and the extent (size) of the image, respectively. It should be emphasized that e_y 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 ₀, and D_zz, are the same for both orientations, the ratio of the plots of I_tot vs. L+e_z 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)²-1)/(4(L/d)²-3), where

 

Fig. 2. Porous plastic fiber sample (Porex). The scale bar to the lower left is 1 mm. The inset shows a magnified view of the same sample with a scale bar of 100 µm. From this SEM, it is clear that the fiber orientations are both highly correlated and highly connected.

Download Full Size | PDF

3. Experiment

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 [27]. 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 [27].

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/cm³, 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 ₁×w ₂=(12×42) mm² and (42×12) mm² 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,β,e_y,and L′, from which I_tot 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 [28] 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) mm² with thicknesses of 5 mm and 2.5 mm for the incoming beam parallel and perpendicular to the fibers respectively.

4. Results

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 e_y 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 _⊥.

 

Fig. 3. Iso-intensity lines for a single measured image (thin lines) and for a fit (thick lines) to Eq. 9 with L=2.0 mm for the parallel (top) and perpendicular (bottom) orientations of the beam with respect to the fibers. The fiber orientation for both plots is indicated by the double-sided arrow. The same scales are used for both axes and both plots to depict the distortion introduced by the anisotropy depending of the sample orientation. The contour lines delineate steps of 0.5, from 0.5 to 3.5 in the top plot, and from 0.2 to 2.7 in the bottom plot in units that are arbitrary in overall scale but consistent from measurement to measurement. The strong sample anisotropy in the diffusion constants is qualitatively clear. Furthermore, the extension of the data and fit over the edge (i.e. the fact that some of the low intensity contours are not closed but come to an abrupt end at the edge of the sample) indicates the effect of the extrapolation length and suggests a method for its measurement. The excellent overlap of the fit on the data indicates that the diffusion model describes the light propagation for this value of L.

Download Full Size | PDF

 

Fig. 4. Data and residuals (open symbols) and fits (solid lines), showing the intensity with the incoming beam parallel (circles) and perpendicular (triangles) to the fibers along the x=0 line. I.e. this is the intensity measured along the line parallel to the incoming beam. The inset shows the narrow range of data near the input edge of the sample (y=0), showing the negative values of the I=0 intercept of the fit. This intercept is the extrapolation length. The much larger intercept for the parallel orientation clearly indicates the larger extrapolation length at this input surface.

Download Full Size | PDF

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>L_min where L_min=0.75 mm. For values of L≤L_min, 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 e_z=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%.

 

Fig. 5. Parameters determined by fitting to the images for each value of L for the two sample orientations. The open and closed circles were measured with the incoming beam parallel and perpendicular to the sample respectively. The diffusion model becomes consistent at values of L of larger than 0.75 mm, thus the fits were performed on these data points. The average values values of β _‖,⊥ and e _‖,⊥ were computed with L>0.75 mm as shown by the horizontal line in plots (a) and (b). In the plot (c) the data was fit to Eq. 11, which is nearly linear for the values of L measured.

Download Full Size | PDF

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 L_min in term of the mean free path, namely, L_min=5l _⊥=1.3l _‖.

5. Discussion

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 [26]. 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 [2]. 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 [12]. 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].