## Abstract

A two-speed photon diffusion equation is developed for light propagation in a powder bed of high volume fraction or dense particulate suspension, whereby the light speed is impacted by the refractive index difference between particles and the suspending medium. The equation is validated using Monte Carlo simulation of light propagation coupled with dynamic simulation of particle sedimentation for the non-uniform arrangement of powder particles. Frequency domain experiments at 650 nm for a 77-µm-diameter resin-powder and 50-µm-diameter lactose-powder beds as well as resin-water and lactose-ethanol suspensions confirm the scattering and absorption coefficients derived from the two-speed diffusion equation.

©2005 Optical Society of America

## 1. Introduction

Measurements of time-dependent propagation of light in particulate and powder media may enable their characterization of the multiply scattering media *in situ*. Recently, our group has focused upon employing frequency domain photon migration (FDPM) measurements and the optical diffusion equation for monitoring the blend homogeneity of active pharmaceutical ingredient (API) within pharmaceutical blending operations by relating the measured absorption coefficient of the powder bed to the API concentration [1–5].

However, the applicability of diffusion approximation in a powder bed or a dense particulate suspension may be more complicated than that in tissue or colloidal suspensions, in which the transit of photons can be approximated by an average constant velocity, *C/n̄*_{med}
, where *c* is the light speed in vacuum and *n̄*_{med}
is the average refractive index of scattering medium. For dilute systems, *n̄*_{med}*≈n*_{med}
, where *n*_{med}
is the refractive index of the suspending medium. In complex media whereby the difference of the refractive index between the scattering particles and suspending fluid can be significant, (such as the case of metal oxides, lactose, or other powders), the relative refractive index and the structure resulting from the packing of scatters may dramatically impact time-dependent measurements of multiply scattered light. As illustrated in Fig. 1, propagating photons in dense powder beds have two velocities. Outside particles, photons travel at a velocity of *c/n*_{med}
, whereas they travel within particles at velocity, *c/n*_{p}
, where *n*_{p}
is the refractive index within the particle of the powder or particulate suspension. While a one-speed diffusion equation may describe photon migration in tissue media, [6] a two-speed diffusion equation may be applicable in dense systems with disparate refractive indices.

In this contribution, we present a two-speed diffusion equation for prediction of time-dependent light propagation through powders and dense particulate suspensions and develop a simulation method, which consists of (i) dynamic simulation for generating structures of powder bed and dense particulate suspensions and (ii) Monte Carlo for tracking photon trajectories through media with generated structures. We then compare the transport coefficients predicted from both two-speed diffusion equation and Monte Carlo simulation results with experimental measurements of time-dependent light propagation in order to confirm the validity of a two-speed model.

The organization of this paper is as follows: In Section 2, we describe (i) the two-speed diffusion equation based on the hypothesis of photon equilibrium in the two-speed radiative transfer equation; (ii) dynamic simulation of powder particle sedimentation; and (iii) Monte Carlo simulation of photon migration within the simulated powder bed. We detail the experimental measurements of time-dependent light propagation made in resin and lactose powders and particulate suspensions using frequency domain photon migration (FDPM) techniques. In Section 3, we present results which (i) validate the hypothesis of photon equilibrium in two-speed diffusion equation through Monte Carlo simulations of photon density distributions; (ii) confirm the two-speed diffusion equation and the scattering coefficient predicted by Monte Carlo simulation through experimental measurements in resin and lactose powder beds; and (iii) demonstrate the impact of volume fraction upon absorption coefficients through FDPM experiments in resin-water and lactose-ethanol suspensions.

## 2. Methods and approach

#### 2.1 Two-speed photon migration diffusion equation

Under the illustrated conditions of a two-phase medium (consisting of the suspending medium and particle phases with a significant relative refractive index difference) one can represent the diffusion propagation of photons with disparate speeds in terms of a two-group model in neutron transport [7].

An expression for the density of photons, *ϕ*_{p}
, propagating at speed *c/n*_{p}
within the particle phase with particle absorption, *µ*_{a,p}
, can be written as

The density of photons, *ϕ*_{med}
, propagating at speed *c/n*_{med}
within the medium of absorption cross section, on, *µ*_{a,med}
, can be expressed from the solution of

In these expressions, the term *r*_{loss,p}*·ϕ*_{p}
represents the loss of photons from the particle phase into the suspending medium phase; and *r*_{loss,med}*·ϕ*_{med}
represents the loss of photons from the suspending medium phase into the particle phase. We assume that the photon loss comes only from the absorption processes and that photons travel between the particle and surrounding medium phases in a purely elastic scattering process. The terms $\overrightarrow{J}$
_{med}
and $\overrightarrow{J}$
_{p}
are current flux of photons within the suspending medium and particle phase, respectively.

Furthermore, we introduce a hypothesis, which will be confirmed by Monte Carlo simulation later in the contribution, that at late times when the diffusion approximation is valid (i.e., $\mid \frac{\partial \stackrel{\u0305}{J}}{\partial t}\mid \ll \frac{c}{3D}\mid \overrightarrow{J}\mid $, where *D* is a diffusion coefficient), there is a local equilibrium between photon densities within each phase such that:

where *K*_{e}
is the equilibrium constant, which depends upon the available volume fraction for photon transit inside particles and their surrounding media.

Given the equilibrium between *ϕ*_{p}
and *ϕ*_{med}
and given *Fick*’s law for multi-component diffusion,[8]

where *D*_{ij}*(i, j=med, p*) is an element of the multi-component diffusion coefficient tensor, Equations (1)–(5) can be combined to obtain

where the diffusion coefficient for the bed is given by:

and the absorption coefficient of the powder bed is given by:

Furthermore, the weak absorption coefficient in powder bed allows the approximation of the isotropic scattering coefficient as ${\mu \prime}_{s,\mathit{bed}}=\frac{1}{\left(3{D}_{\mathit{bed}}\right)}$.

Equations (6)–(8) constitute the diffusion equations for a general two-speed photon migration in powder bed or dense particulate suspension. They differ from a typical one-speed diffusion equation in that the scattering and absorption coefficients are dependent upon the refractive indices of the two phase system. In addition, the two speed equation suggests that the absorption coefficient is also dependent upon the volume fraction of particles which act as scatters.

#### 2.2 Dynamic simulation of sedimentation process

In order to simulate particle positions (or powder structure), we determine the final state through dynamic simulation as described by Yang *et al*. [9] The dynamic simulation numerically first solves the Newtonian equations of motion for up to 50,000 interacting particles, which fall and collide with each other under gravity, friction, and interaction forces, and then records the evolution of spatial particle positions and velocities with time.

The forces acting on a particle govern both translational and rotational motions within the ensemble. The Newtonian equations describing the evolution of particle positions and rotational and translational velocities within the ensemble at time *t* are:

where $\overrightarrow{v}$
_{i}
is the velocity of mass center of particle *i* and $\overrightarrow{\omega}$
_{i}
is its rotational velocity. The total force on particle *i, $\overrightarrow{F}$*_{i}
, arises from (i) gravity, *m*_{i}$\overrightarrow{g}$
; (ii) van der Waals force, $\overrightarrow{F}$
_{i,vdW}
;(iii) the normal and frictional forces due to particle contacting, $\overrightarrow{F}$
_{i,normal}
and $\overrightarrow{F}$
_{i,friction}
; and (iv) drag forces, $\overrightarrow{F}$
_{i,drag}
. $\overrightarrow{T}$
_{i}
is the total torque, which arises from $\overrightarrow{T}$
_{ij}
, friction and $\overrightarrow{T}$
_{ij,rolling}*; I*_{i}
is the inertial moment; and *N*_{i,contact}
is the number of neighboring particles that contact particle *i*. Each force in Eqs. (9) and (10) and associated parameter values used in the simulation are described in Appendix I.

At the beginning of the simulation, the vectors of particle positions are set so that particles are separated and do not contact each other. The algorithm then progresses in time in steps, *h*, of 2×10^{-8} seconds and computes particle position and translational and rotational velocities at each time step:[10]

In the dynamic simulation, periodic boundary conditions along horizontal *x* and *y* directions are used to approximate semi-infinite settling of a bed of particles. The computation required 5×10^{6} seconds for the simulation of 50,000 particles on a Pentium 4, 2.8 GHz machine with a Linux operating system. The calculation flowchart is given in Appendix II.

In order to investigate the impact of volume fraction upon bed *D*_{bed}
and *µ*_{a,bed}
, the magnitude of forces acting on the particles can be adjusted in order to generate various particle ensembles. However, it is difficult to experimentally manipulate properties which govern these forces. Consequently, in order to investigate the influence of volume fraction, we simulated fluid particulate suspension by uniformly expanding the final simulated resin and lactose powder beds at *f*_{v}
=0.64 and conducting Monte Carlo simulation of photon trajectories (as described below) on particulate arrangements at lower volume fractions. This enables comparison of Monte Carlo simulation results (described below) to experimental measurements in dense resin and lactose suspensions.

#### 2.3 Monte Carlo simulation of photon migration in particulate suspensions

The powder structure generated from dynamic simulation was employed as the scattering medium in which the statistics of time-dependent photon propagation was simulated with Monte Carlo techniques. An infinite medium of particles with diameter, *d*_{p}
, was simulated based upon the generated powder structure within a 10*d*_{p}
×10*d*_{p}
×32*d*_{p}
cubic with periodic boundary conditions along *x, y*, and *z* directions. Photon trajectories were computed using Monte Carlo technique.

In contrast to the use of *Mie* scattering for characterizing dense colloidal suspensions, where *d*_{p}*≤λ* (wavelength), [11,12] geometric optics can be employed for simulating trajectories of photon propagation when particle diameters are larger than the wavelength of near-infrared light employed [13]. Geometric optics dictates the propagation direction and scattering length in the simulated photon trajectories.

The simulation is conducted as follows: a photon package, with unity weight, is launched at a point of illumination in the suspending medium. When the photon encounters a particle, it is either reflected or refracted, depending upon (i) a randomly generated number, (ii) the relative refractive index of the particle, and (iii) the probability for reflection and refraction, which is calculated according to *Fresnel*’s formula [14].

If the photon package is reflected, it leaves the particle surface in the reflection direction and travels until it is re-directed upon encountering the next particle. The photon weight is reduced due to the absorption coefficient of the non-particle phase.

If the photon package is refracted, it transits inside the particle with an updated propagation direction and the photon weight is reduced due to the non-zero particle absorbance. Once the photon package arrives at the next particle interface, the process is repeated and the package is either reflected or refracted as described before.

The photon weight, *w*, after each *i*^{th}
scattering step is attenuated by

where the subscript *j* denotes the medium in which the step has occurred (i.e., either within the suspending, *med*, or particle, *p*, phases) and *t*_{i}
, represents the “time-of-flight” associated with the photon scattering step length.

The photon position, propagating direction, photon weight, and residing phase (i.e., whether in suspending medium or particle) are recorded with propagation time. At each recording time, a photon contributes to the photon density in suspending medium, *ϕ*_{med}
, if it is within the medium, and contributes to the photon density within particles, *ϕ*_{p}
, if it is located within the particle.

Photon trajectories are counted for propagation times less than 2.3×10^{-8} seconds. We simulate 10^{6} photon trajectories with a computational time of 4×10^{5} seconds on a Pentium 4, 2.8 GHz machine with a Linux operating system. The details of the Monte Carlo simulation of geometric optics are provided in Appendix III.

#### 2.4 Determination of scattering and absorption coefficients from Monte Carlo simulation

From the square of the mean displacement at time *t*, <*r*(*t*)>, derived from Monte Carlo statistics, the effective diffusion coefficient,*D*_{bed}
, in powder bed can be determined [15]

whereby *D*_{bed}
is a measurable quantity from time-dependent measurements, as described by Eq. (7).

The effective absorption coefficient, *µ*_{a,bed}
, in Eq. (8) for powder bed or dense particulate suspension is a function of the photon equilibrium constant, *K*_{e}
, which can be related to the ratio of mean of the total time of collected photons spent within the particle phase, <*t*_{total,p}
>, to that within the suspending medium, <*t*_{total,med}
>:

Furthermore, <*t*_{total,p}
> and <*t*_{total,med}
> can be computed from the mean scattering length within particle phase, scat, <*l*_{scat,p}
>, and the mean scattering length within the suspension medium, <*l*_{scat,med}
>,

where *N*_{scat, j}
is the number of total scattering events within suspending medium or particle phase.

In terms of *Fresnel*’s formula, [14] the reflection probability for a photon, which propagates in the particle phase and interacts with the interface, is equal to the reflection probability for a photon, which propagates in the suspending2 phase and interacts with the interface. Consequently, *N*_{scat,med}*=N*_{scat,p}
(see Appendix IV for proof). Therefore, the photon equilibrium constant and associated absorption coefficient can be expressed as:

Again, it is noteworthy that the absorption coefficient of the bed is influenced by the particle configuration as indicated by its dependence upon <*l*_{scat,med}
>, which is a function of volume fraction, *f*_{v}
, and particles diameter, *d*_{p}
. In this contribution, we us Monte Carlo simulation to investigate the dependence of <*l*_{scat,med}
> upon *f*_{v}
, which varies from 0.05 to 0.64. For a system of monodisperse particles, we provide a mathematical expression for <*l*_{scat,med}
> as a function of *f*_{v}
and *d*_{p}
.

#### 2.5 Quantitative experimental time-dependent measurements

To experimentally confirm the dynamic simulation, Monte Carlo, and two-speed diffusion equation analysis, time-dependent measurements were made using frequency domain system established in the Photon Migration Laboratories. In the following, the measured samples are described, the experimental approaches are briefly presented, and the data analysis is provided.

### 2.5.1 Samples

Measurements were conducted on resin and lactose powder beds, and in order to investigate the influence of varying volume fraction of particles, suspensions were interrogated. A uniform spatial distribution of particles was achieved by suspending the resin and lactose particles in water and liquid ethanol with use of an external agitator (H3770, Sigma-Aldrich, St. Louis, MO). The volume fraction, *f*_{v}
, of lactose-ethanol suspensions was varied from 0.17 to 0.48, and *f*_{v}
of resin-water suspension was varied from 0.20 to 0.64 by adding desired amount of resin and lactose particles into water and ethanol.

Monodisperse particles of AG 1-X8 Resin (Bio-Rad, Hercules, California), shown in the micrographs in Fig. 2(a), were employed. Resin particles have a regular spherical shape and narrow particle size distribution with a mean particle diameter of 77 µm and a standard deviation of 8 µm as determined from image analysis (Image Pro 3.0; Cybernetics, Silver Spring, MD). The polydisperse pharmaceutical lactose particles (Le-Pro^{™} Lac 455; Leprino Food, Denver, CO), shown in Fig. 2(b), have irregular shapes and a broad particle size distribution with a mean value of around 50 µm and a standard deviation of 30 µm.

The refractive index for resin particles, which are made up of polystyrene divinylbenzene, is 1.59 at 650 nm [16]. The refractive index of lactose particles in crystal form is 1.58 at 650 nm [17]. The refractive indices for purified water and liquid ethanol (A406^{P}-4; Fisher Scientific, Fair Lawn, NJ) are both 1.33 at 650 nm [18].

The volume fraction, *f*_{v}
, of the powder bed was determined by the ratio of the mass density of powder bed, *ρ*_{bed}
, to the particle mass density, *ρ*_{p}
. In a densely packed state, *ρ*_{bed}
of resin bed was determined to be 0.65 gram/cm^{3} by measuring the weight of powder bed and the associated bed volume, while *ρ*_{p}
of resin particle was determined to be 1.02 gram/cm^{3} by measuring the weight of particle bed and the associated net volume occupied by particles, which was measured by settling the resin particles in water and reading the volume increment. We found *f*_{v}
for the monodisperse resin powder bed in air to be 0.64. For polydisperse lactose powders, the same measurements gave *ρ*_{p}
=1.46gram/cm^{3}, 0.95 *ρ*_{bed}
=gram/cm^{3}, and *f*_{v}
=0.65.

### 2.5.2 Frequency domain photon migration measurements

FDPM measurements were made using the system detailed elsewhere [1]. Briefly, modulated light from a 650 nm, 35 mW laser diode was split into two beams. One beam was directed into a 1-mm-diameter optical fibers (FT-035 mm; Thorlabs, Newton, NJ) and was delivered to a reference photomultiplier tube (PMT) (H6357, Hamamatsu, Japan). The second light beam was introduced into a second optical fiber whose end was embedded within the powder bed. The light propagated through the scattering medium and was collected with a 1 mm optical fiber located (i) 5–8 mm away from the point of illumination in increments of 0.5 mm for the resin powder bed; (ii) 10–18 mm away from the point of illumination in increments of 2 mm in the lactose powder bed; and (iii) 15–22 mm away from the point of illumination in increments of 0.5 mm for lactose-water and resin-ethanol suspensions. The two PMTs were modulated at the same modulation frequencies as the laser diode (40–125 MHz), but with an additional offset frequency of 100 Hz. The resulting mixed signals from the PMTs were converted into voltage; the high frequency components were filtered by transimpedance amplifiers (70710; Oriel Co., Stratford, CT); and the heterodyned signals were acquired by LabView data acquisition software (National Instruments, Austin, TX). The phase shift (*PS*), amplitude (*AC*), and mean intensity (*DC*) of the signal from the sample PMT relative to the reference PMT were recovered for data analysis.

### 2.5.3 Analysis of phase shift, amplitude, and mean-intensity of photon density wave

FDPM measurements of phase shift (*PS*), amplitude (*AC*), and mean intensity (*DC*) for (i) the resin and lactose powder beds with a volume of 150 cm^{3} and (ii) the resin-water and lactose-ethanol suspensions with volumes of more than 400 ml were evaluated from the solution of the standard diffusion equation to provide the absorption coefficient, *µ*_{a,bed}
and scattering coefficient, *µ′*_{s,bed}
[19]. For simplicity, we denote the optical properties of the powder bed as well as the dense suspension in terms of *µ*_{a,bed}
and *µ′*_{s,bed}
.

From the experimental quantities of *PS, AC*, and *DC*, we compute a parameter, *Y*:

from which the bed scattering coefficient can be determined by

Here *ω* is the modulation frequency, *r* is the source-detector separation, and the subscript, 0, denotes the values at the smallest source-detector separation of (i) 5 mm in resin powder bed and (ii) 10 mm in lactose powder bed. Equation (21) was employed to extract the FDPM measured scattering coefficients, *µ′*_{s,bed}
, from the measurements of *Y*(*PS,AC,DC,r*) when *ω* was varied (i) from 65 to 125 MHz in the resin powder bed and (ii) from 40 to 120 MHz in lactose powder bed.

## 3. Results and discussion

#### 3.1 Powder structure from dynamic simulation of sedimentation processes

In order to illustrate the evolution of powder structure with time, the sedimentation of a powder bed, which contained 2,000 particles of 50 *p*_{d}
=50*µm* in a (5×10^{-4}
*m*).(5×10^{-4}
*m*).(2×10^{-3}
*m*) cube with periodic boundary conditions along horizontal *x* and *y* directions, is shown in Fig. 3.

The powder bed becomes more dense with the sedimentation time. At sedimentation times of 0, 0.012, 0.015, 0.024, and 0.032 seconds, the volume fractions of the powder bed, *f*_{v}
, are 0.21, 0.28, 0.32, 0.53, and 0.59, respectively. It is assumed that a mechanically balanced state is achieved when the volume fraction no longer changes with sedimentation time.

When the number of simulated lactose particles was increased to 4,000 and 50,000, the corresponding volume fraction, *f*_{v}
, grew to 0.61 and 0.64, respectively, due to the increase of vertical pressure caused by bed weight.

#### 3.2 Photon Equilibrium

Figure 4(a) is a plot of the time-dependent photon densities, *ϕ*_{med}
and *ϕ*_{p}
, statistically computed from Monte Carlo trajectories within a simulated powder bed (*f*_{v}
=0.61, 50*d*_{p}
=50*µm*, *n*_{med}
=1.0, and *n*_{p}
=1.6), which was generated based on the dynamic simulation of 4,000 particles. The squares, crosses, and triangles represent *ϕ*_{p}
and *ϕ*_{med}
at positions of 75*d*_{p}
, 165*d*_{p}
, and 255*d*_{p}
from the light source. While *ϕ*_{med}
and *ϕ*_{p}
both change with propagation time at all positions, they are related by a factor,*K*_{e}
, as hypothesized in Eq. (3), found in this case to be 2.44. To illustrate photon equilibrium, Fig. 4(b) shows *ϕ*_{p}
and (*K*_{e}
×*ϕ*_{med}
) at a series of source-detector separations for migration times varying from 6×10^{-12} to 1×10^{-8} second. Figure 4(b) indicates that at all times, there is an equilibrium between the photon density within particles and photon density in medium phase, regardless of source-detector separation in the simulated powder bed. We found that the equilibrium constant, *K*_{e}
, varies with *f*_{v}
between values of 0.05 and 0.64 (see Appendix V). These results confirm the photon equilibrium relationship between *ϕ*_{med}
and *ϕ*_{p}
, used to develop the two-speed diffusion approximation of time-dependent photon migration for characterization of powder beds and dense particulate suspensions.

The value of *K*_{e}
at each volume fraction can be computed by Eq. (18) based on the statistical quantities of <*l*_{scat,med}
> and <*l*_{scat, p}
>. Figure 5(a) illustrates the values of <*l*_{scat,med}
>/*d*_{p}
computed from Monte Carlo as a function of *f*_{v}
, which varies from 0.05 to 0.64, the relative refractive index, *n*_{p}*/n*_{med}
, which varies from 1.15 to 1.8, and particle diameter, *d*_{p}
, which varies from 20 to 50*µm*. As expected, the relative refractive index of the particle has little impact on <*l*_{scat,med}
>. The reduction of <*l*_{scat,med}
>/*d*_{p}
is consistent with the decrease of the void space between particles when the volume fraction is increased. Quantitatively, in the ranges of *f*_{v}
measured for the resin-water and lactose-ethanol suspensions, the dependence of <*l*_{scat,med}
>upon fv in Fig. 5(a) can be described empirically by a power form fit:

Figure 5(b) shows that <*l*_{scat,p}
>/*d*_{p}
increases with the increase of *n*_{p}*/n*_{med}
, which varies from 1.1 to 1.8 and that the volume fraction, which varies from 0.1 to 0.64, has little impact upon <l_{scat,p}
>*d*_{p}
. Therefore, <*l*_{scat,p}
> becomes a constant when both *d*_{p}
and *n*_{p}*/n*_{med}
are given. From Fig. 5, the mean scattering lengths were determined to be, <*l*_{scat,p}
>=0.895*d*_{p}
and <*l*_{scat,med}
>=0.575*d*_{p}
for the powder bed of *f*_{v}
=0.61 and *n*_{p}*/n*_{med}
=1.6. For this powder bed, Eq. (18) predicts 2.49 Ke=2.49, which agrees well with *K*_{e}
=2.44 obtained from the ratio of *ϕ*_{p}
to *ϕ*_{med}
computed from Monte Carlo.

#### 3.3 Measurements in powder beds for verifications of two-speed diffusion equation and simulated scattering coefficient

FDPM measurements of phase shift (*PS*), amplitude (*AC*), and mean intensity (*DC*) were made in the resin and lactose powder beds. Figure 6(a) shows the FDPM measurements for *Y*(*PS,AC,DC,r*) (i.e., Eq. (20)) versus the source-detector separation, *r*, at 65 MHz (diamonds), 95 MHz (triangles), and 125 MHz (squares) for the monodisperse resin powder bed with *d*_{p}
=77*µm* and 0.64 *f*_{v}
=0.64. Fig. 6(b) shows *Y*(*PS,AC,DC,r*) versus *r* at 40 MHz (diamonds), 80 MHz (triangles), and 120 MHz (squares) for the polydisperse lactose powder bed with a mean diameter of 50 µm and *f*_{v}
=0.65. In both cases, *Y*(*PS,AC,DC,r*) changes linearly with *r* at each modulation frequency, as predicted by Eq. (21) derived from the diffusion equation. From Eq. (21), we found the values of FDPM measured scattering coefficients to be *µ′*_{s,bed}
=47.8 cm^{-1} and *µ′*_{s,bed}
=243cm^{-1} for resin and lactose powder beds, respectively. The Monte Carlo predicted scattering coefficient, *µ′*_{s,bed}
, for a simulated monodisperse resin powder bed with particle diameter of 77*µm*, refractive index of 1.6, and volume fraction of 0.64 was found to be *µ′*_{s,bed}
=48.9 cm^{-1}. In the monodisperse resin powder bed, the scattering coefficient predicted by Monte Carlo simulation is close to the actual FDPM measured scattering coefficient and correspondingly, the predictions of *Y*(*PS,AC,DC,r*) with use of simulated *µ′*_{s,bed}
in Eq. (21) are also close to the measured values, as shown in Fig. 6(a).

For a simulated monodisperse lactose powder bed of *d*_{p}
=50*µm* and *f*_{v}
=0.64, the scattering coefficient predicted by Monte Carlo simulation was found to be *µ′*_{s,bed}
=75 cm^{-1}. Contrary to the results for the resin powder bed, the scattering coefficient predicted by Monte Carlo simulation for the monodisperse lactose powder bed of *d*_{p}
=50*µm* differs significantly from the actual FDPM measured scattering coefficient for the polydisperse powder bed with a mean diameter of *d̄*_{p.poly}
=50*µm*. The prediction of scattering coefficient in the polydisperse powder bed may require more considerations of the particle size distribution and particle shapes.

#### 3.4 Particulate suspension experiments for confirmation of the dependence of absorption coefficient upon volume fraction

Figure 7 shows the experimentally derived absorption coefficient of the suspensions, *µ*_{a,bed}
, versus volume fraction, *f*_{v}
. In the resin-water particulate suspension, *µ*_{a,bed}
decreased from 0.0193 cm^{-1} to 0.0098 cm^{-1} when *f*_{v}
was increased from 0.20 to 0.64. In the lactose-ethanol particulate suspension, *µ*_{a,bed}
decreased from 0.0143 cm^{-1} to 0.0082 cm^{-1} when *f*_{v}
was increased from 0.17 to 0.48.

The dependence of *µ*_{a,bed}
upon *f*_{v}
can be evaluated explicitly based on the combination of Eq. (19) and the empirical relationship between <*l*_{scat,med}
> and *f*_{v}
(Eq. (22)):

where *a*
_{1} and *a*
_{2} are constants defined as ${a}_{1}=\frac{{\mu}_{a,p}<{l}_{\mathit{scat},p}>}{0.395{d}_{p}}$ and ${a}_{1}=\frac{<{l}_{\mathit{scat},p}>}{0.395{d}_{p}}.\frac{{n}_{p}}{{n}_{\mathit{med}}}$.

Equation (23) was employed to fit the FDPM measured absorption coefficients in lactose-ethanol and resin-water suspensions with adjustable parameters, *µ*_{a,med}
,*a*
_{1} and *a*
_{2}. The values of *µ*_{a,med}
at 650 nm were thus found to be 0.028 and 0.022 cm^{-1} for suspending water and ethanol, respectively, and the fitting curves (solid lines) are also shown in Fig. 7. The determined *µ*_{a,med}
for suspending water is higher than the corresponding absorption coefficient for pure water reported in literature[20] possibly because of the release of chloride ions from the resin particles into the suspending water, which became slightly dark in color compared with the pure water.

Finally, Eq. (19) indicates that the particle diameter, *d*_{p}
, has little impact on the absorption coefficient of powder bed, *µ*_{a,bed}
, since both <*l*_{scat,med}
> and <*l*_{scat, p}
> are proportional to *d*_{p}
(as shown in Fig. 5.), which is then cancelled as a common factor in the denominator and the numerator of Eq. (19). The conclusion of the independence of *µ*_{a,bed}
upon *d*_{p}
is consistent with our past experimental results.^{[4]}

## 4. Conclusions

In the powder beds and dense particulate suspensions, where (i) particle size is much larger than operating near-infrared wavelength; (ii) refractive index of particle is significantly different from that of the suspending medium (i.e. air or liquid); and (iii) the volume fraction of the particles is of the same order of that of the suspending medium, we find that, (i) there is a local photon equilibrium between photon density in the suspending medium and the photon density within particles; (ii) the photon equilibrium validates the two-speed diffusion equation; and (iii) the absorption coefficient is a function of the volume fraction, but independent of particle diameter. These results suggest that for light propagation in disparate media, the diffusion approximation still applies, but that the transport coefficients can be predicted from accurate representation of a two-speed model of time-dependent light propagation.

## 5. Appendix

## 5.1 Calculation of forces on particles

(i) van der Waals force, $\overrightarrow{F}$
_{i,vdW}
[21]

## 5.2 Calculation flow-chart for dynamic simulation

* The mechanical properties are included in the table in Appendix 1.

** *See reference [29]*

## 5.3 Algorithm for Monte Carlo of geometric optics

The number of photon trajectories is set to be 10^{6} and the migration time, *t*_{end}
, for a photon in Monte Carlo is 2.3×10^{-8} second. The photon trajectories are simulated one by one. A photon package includes time, *t*, propagating direction, $\overrightarrow{v}$
, position vector, $\overrightarrow{r}$
, photon weight, *w*, and refractive index for the surrounding phase, *n*_{med}
or *n*_{p}
.

When the simulation of a photon begins, the photon is launched into the powder bed at *t*=0. The initial propagating direction, $\overrightarrow{v}$
_{0}
, at the origin, $\overrightarrow{r}$
_{0}, is randomly chosen from the isotropic distribution and the photon weight, *w*, is set to be 1. The photon package at the origin is [*t,$\overrightarrow{v}$
,$\overrightarrow{r}$
,w,n*]⃖[*0,$\overrightarrow{v}$
0,$\overrightarrow{r}$
0,1,n*_{med}
].

The particle, which the photon encounters, is identified from the photon moving direction, $\overrightarrow{v}$
_{0}, and the particle positions, which are obtained from the dynamic simulation of sedimentation process for powder bed. The time for a photon propagating from the origin, $\overrightarrow{r}$
_{0}, to the encountering point, $\overrightarrow{r}$
_{ext}
, on the external surface of the interacting particle with a speed,*c/n*_{med}
, is Δ*t*. The attenuation of photon weight, *w*, is exp(-*µ*_{a,med}
·|$\overrightarrow{r}$
_{ext}-$\overrightarrow{r}$
_{0}|). Correspondingly, the photon package at $\overrightarrow{r}$
_{ext}
, is updated [*t,$\overrightarrow{v}$
,$\overrightarrow{r}$
,w,n*]⃖[(*t+Δt*),$\overrightarrow{v}$
$\overrightarrow{r}$
_{ext}
,(*w*.exp(-*µ*_{a,med}
·|$\overrightarrow{r}$
_{ext}-$\overrightarrow{r}$
_{0}|)),*n*].

At the encountering point, $\overrightarrow{r}$
_{ext}
, the incident angle, *α*_{i}
, between the photon propagating direction, $\overrightarrow{v}$
_{0}, and the normal, $\overrightarrow{L}$
, of the external surface is employed to calculated the reflection probability in terms of *Fresnel*’s formulas[11]

where *α*_{t}
is the transmission angle determined by the Snel’s law^{[10]}

A random number, *ξ*, from the uniform distribution in the range of [0,1] is generated. The photon undergoes an external reflection if *ξ<P*. Otherwise, it undergoes a refraction from outside to inside of the particle.

If the photon undergoes an external reflection, the photon propagating direction is updated to $\overrightarrow{v}$
in terms of *Snel*’s law and the photon package becomes [*t,$\overrightarrow{v}$
,$\overrightarrow{r}$
,w,n*]⃖[*t,*$\overrightarrow{v}$
_{1},$\overrightarrow{r}$
,*w,n*].

In terms of $\overrightarrow{v}$
_{1} and the particle positions, the next interacting particle is determined and similar updating for propagating outside particles is conducted.

If the photon undergoes a refraction from outside to inside of the particle, the photon propagating direction is updated to $\overrightarrow{v}$
_{2} in terms of the *Snel*’s law. The refractive index in the photon package is changed from *n*_{med}
to *n*_{p}
since the photon transmits inside the particle. The photon package is updated [*t*,$\overrightarrow{v}$
,$\overrightarrow{r}$
,*w,n*]⃖[*t*,$\overrightarrow{v}$
_{2},$\overrightarrow{r}$
*,w,np*].

In terms of the propagating direction, $\overrightarrow{v}$
_{2}, the position of the interacting particle, and the particle radius, the ending point for the photon’s internal travel is determined as $\overrightarrow{r}$
_{int}
. The propagating time, Δ_{t}, is calculated from the distance between $\overrightarrow{r}$
_{ext}
and int $\overrightarrow{r}$_{int}
and the particle refractive index, *n*_{p}
. The attenuation of the photon weight is ext (-*µ*_{a.p}
·|$\overrightarrow{r}$_{int}-$\overrightarrow{r}$_{ext}|). At position, $\overrightarrow{r}$_{int}, the photon package is updated [*t*,$\overrightarrow{v}$,$\overrightarrow{r}$,*w,n*]⃖[(*t*+Δ*t*), $\overrightarrow{v}$$\overrightarrow{r}$
_{int}, (*w*.exp(-*µ*_{a.p}
·|$\overrightarrow{r}$_{int}-$\overrightarrow{r}$_{ext}|)), *n*].

At the interacting point, $\overrightarrow{r}$_{int}, the incident angle, *α*_{i}
, between $\overrightarrow{v}$_{2} and the normal of the particle surface, $\overrightarrow{L}$, is calculated. The reflection probability, *P*, is calculated in terms of *Fresnel*’s formulas. Then a random number, *ξ*, is generated from the uniform distribution in the range of [0,1].

The photon undergoes an internal reflection if *ξ*<*P*. Otherwise, it undergoes a refraction from inside to outside of the particle.

When the photon undergoes an internal reflection, the propagating direction becomes $\overrightarrow{v}$_{3} in terms of *Snel*’s law and the photon package is updated 3 [*t*,$\overrightarrow{v}$,$\overrightarrow{r}$,*w,n*]⃖[*t*,$\overrightarrow{v}$_{3},$\overrightarrow{r}$,*w,n*].

Similar updating for propagation within the particle is conducted.

When *ξ*≥*P*, the photon undergoes an refraction from inside to outside of the particle. The propagating direction is changed to $\overrightarrow{v}$_{4} in terms of the *Snel*’s law and the associated refractive index is changed from *n*_{p}
to *n*_{med}
. Then the photon package is updated [*t*,$\overrightarrow{v}$$\overrightarrow{r}$,*w,n*]⃖[*t*,$\overrightarrow{v}$_{4},$\overrightarrow{r}$,*w,n*_{p}
].

Similar updating for photon propagating outside particles is conducted.

In the simulation, when the photon package is updated, *t* and *t*_{end}
is compared. When *t*<*t*_{end}
, the simulation for the photon continues. Otherwise, the simulation for the photon terminates and another photon is launched into the powder bed. When the total simulated photon trajectories reaches 10^{6}, the simulation is finished.

The simulation flowchart is shown below.

where, $\overrightarrow{v}$ is the photon propagating direction; $\overrightarrow{v}$_{reflection}
is the propagating direction after reflection at an interface; and $\overrightarrow{v}$_{transmission}
is the propagating direction after refraction at an interface.

## 5.4 Monte Carlo simulation of scattering numbers

## 5.5 Photon equilibrium for powder beds with volume fraction of 0.45 and 0.38

## References and links

**1. **R. R. Shinde, G. V. Balgi, S. L. Nail, and E. M. Sevick, “Frequency-domain photon migration measurements for quantitative assessment of powder absorbance: a novel sensor of blend homogeneity,” J. Pharm. Sci. **88**, 959–966 (1999). [CrossRef] [PubMed]

**2. **T. Pan and E. M. Sevick-Muraca, “Volume of pharmaceutical powders probed by frequency-domain photon migration measurements of multiply scattered light,” Anal. Chem. **74**, 4228–4234 (2002). [CrossRef] [PubMed]

**3. **T. Pan, D. Barber, D. Coffin-Beach, Z. Sun, and E. M. Sevick-Muraca, “Measurement of low dose active pharmaceutical ingredient in a pharmaceutical blend using frequency-domain photon migration,” J. Pharm. Sci. **93**, 635–645 (2004). [CrossRef] [PubMed]

**4. **Z. Sun, S. Torrance, F. K. McNeil-Watson, and E. M. Sevick-Muraca, “Application of frequency domain photon migration to particle size analysis and monitoring of pharmaceutical powders,” Anal. Chem. **75**, 1720–1725 (2003). [CrossRef] [PubMed]

**5. **S. E. Torrance, Z. Sun, and E. M. Sevick-Muraca, “Impact of excipient particle size on measurement of active pharmaceutical ingredient absorbance in mixtures using frequency domain photon migration,” J. Pharm. Sci. **93**, 1879–1889 (2004). [CrossRef] [PubMed]

**6. **V. Venugopalan, J. S. You, and B. J. Tromberg, “Radiative transport in the diffusion approximation: an extension for highly absorbing medium and small source-detector separations,” Phys. Rev. E **58**, 2395–2407 (1998). [CrossRef]

**7. **J. J. Duderstadt and L. J. Hamilton, “*Nuclear Reactor Analysis*,” (John Wiley & Sons, 1976), *Chapter 4.*

**8. **A. F. Henry, “*Nuclear-Reactor Analysis*,” (The MIT Press, 1975), *Chapter 9.*

**9. **R. Y. Yang, R. P. Zou, and A. B. Yu, “Effect of material properties on the packing of fine particles,” J. Appl. Phys. **94**, 3025–3034 (2003). [CrossRef]

**10. **J. P. Hansen and I. R. McDonald, “Theory of Simple Liquid,” (Academic Press, 1986), *Chapter 3.*

**11. **Z. Sun and E. M. Sevick-Muraca, “Investigation of particle interactions in dense colloidal suspensions using frequency domain photon migration: bidisperse systems,” Langmuir , **18**, 1091–1097 (2002). [CrossRef]

**12. **Y. Huang, Z. Sun, and E. M. Sevick-Muraca, “Assessment of electrostatic interaction in dense colloidal suspensions with multiply scattered light,” Langmuir **18**, 2048–2053 (2002). [CrossRef]

**13. **C. F. Bohren and D. R. Huffman, “*Absorption and Scattering of Light by Small Particles*,” (John Wiley & Sons, 1983), *Chapter 7.*

**14. **L. Wang, S. L. Jacques, and L. Zheng, “MCML-Monte Carlo modeling of light transport in multi-layered tissue,” Comput. Methods Programs Biomed. **47**, 131–146 (1995). [CrossRef] [PubMed]

**15. **S. Chandrasekhar, “Stochastic problems in physics and astronomy,” Review of Modern Physics **15**, 1–89 (1943). [CrossRef]

**16. **X. Ma, J. Q. Lu, R. S. Brock, K. M. Jacobs, P. Yang, and X.-H. Hu, “Determination of complex refractive index of polystyrene microsphere from 370 to 1610 nm,” Phys. Med. Biol. **48**, 4165–4172 (2003). [CrossRef]

**17. **A. N. Winchell, “*The Optical Properties of Organic Compounds*,” (Academic Press, 1954), *Chapter 4.*

**18. **H.-J. Moon, K. An, and J.-H. Lee, “Single spatial mode selection in a layered square microcavity laser,” Appl. Phys. Lett. **82**, 2963–2965 (2003). [CrossRef]

**19. **J. B. Fishkin, P. T. C. So, A. E. Cerussi, S. Fantini, M. A. Franceschini, and E. Gratton, “Frequency-domain method for measuring spectral properties in multiple-scattering media: methemoglobin absorption spectrum in a tissuelike phantom,” Appl. Opt. **34**, 1143–1155 (1995). [CrossRef] [PubMed]

**20. **L. Kou, D. Labrie, and P. Chylek, “Refractive indices of water and ice in the 0.65- to 2.5-µm spectral range,” Appl. Opt. **32**, 3513–3540 (1993). [CrossRef]

**21. **R. Y. Yang, R. P. Zou, and A. B. Yu, “Computer simulation of the packing of fine particles,” Phys. Rev. E **62**, 3900–3908 (2000). [CrossRef]

**22. **N. V. Brilliantov, F. Spahn, and J. M. Hertzsch, “Model for collisions in granular gases,” Phys. Rev. E **53**, 5382–5392 (1996). [CrossRef]

**23. **P. A. Thompson and G. S. Grest, “Granular flow: friction and the dilatancy transition,” Phys. Rev. Lett. **67**, 1751–1754 (1991). [CrossRef] [PubMed]

**24. **L. G. Leal, “*Laminar Flow and Convective Transport Processes*,” (Butterworth-Heinemann, 1992), *Chapter 4.*

**25. **F. Podczeck, J. M. Newton, and M. B. James, “The adhesion force of micronized salmeterol xinafoate particles to pharmaceutically relevant surface materials,” J. Phys. D: Appl. Phys. **29**, 1878–1884 (1996). [CrossRef]

**26. **R. J. Roberts, R. C. Rowe, and P. York, “The relationship between Young’s modulus of elasticity of organic solids and their molecular structure,” Powder Technol. **65**, 139–146 (1991). [CrossRef]

**27. **Y. Shimada, Y. Yonezawa, H. Sunada, R. Nonaka, K. Katou, and H. Morishita, “Development of an apparatus for measuring adhesive force between fine particles,” KONA, No 20, 223–230 (2002).

**28. **K. Z. Y. Yen and T. K. Chaki, “A dynamic simulation of particle rearrangement in powder packings with realistic interactions,” J. Appl. Phys. **71**, 3164–3173 (1992). [CrossRef]

**29. **L. M. Zurk, L. Tsang, K. H. Ding, and D. P. Winebrenner, “Monte Carlo simulations of the extinction rate of densely packed spheres with clustered and nonclustered geometrics,” J. Opt. Soc. Am A **12**, 1772–1781 (1995). [CrossRef]