1. Introduction

Near-infrared imaging and spectroscopy (NIRI and NIRS) provide chemical components and structural properties for various kinds of random media such as biological tissue volumes [1–3] and agricultural products [4]. To interpret and understand the signal of NIRI and NIRS, theoretical and numerical examinations of photon transport for random media are indispensable. Colloidal suspensions have been widely used as liquid phantoms (e.g., Intralipid and Intralipos) for preliminary studies of biolomedical and postharvest applications of NIRI and NIRS [5–9]. Hence, high-fidelity modeling of photon transport for colloidal suspensions is crucial.

The radiative transfer theory (RTT) can describe photon transport in terms of the absorption and scattering of photons [2,10,11]. On length and time scales of the RTT, a colloidal suspension is considered as a random medium. Photon transport is roughly categorized into the ballistic and diffusive regimes from short to long scales of the RTT [12,13]. Photon absorption is quantified by the absorption coefficient in both regimes, while photon scattering in the ballistic regime is quantified by the scattering coefficient, anisotropy factor, and normalized phase function; and in the diffusive regime, it is only by the reduced scattering coefficient. As photon transport models based on the RTT, the radiative transfer equation (RTE) is valid in both regimes, while the diffusion approximation of the RTE is valid only in the diffusive regime. The absorption, scattering, and reduced scattering coefficients are determined by inverse analysis based on the photon transport model because these coefficients are parameters on the RTT-scale. The results for the coefficients by inverse analysis are strongly influenced by the model fidelity for photon transport and the numerical accuracy for calculating the photon transport model.

On the microscopic scale from the RTT-scale, the electromagnetic theory (EMT) can calculate the coefficients by modeling a colloidal suspension as a system where dielectric scatterers are suspended in a continuous background medium. Hence, the EMT-calculation allows to validate the determination procedures of the coefficients based on the RTT [14,15]. Because photon scattering is dominant over photon absorption at the near-infrared wavelength range from 600 to 1200 nm, which is used in biomedical optics, modeling of photon scattering is significant.

As the EMT in far-field, the independent and dependent scattering theories (IST and DST) have been extensively developed to calculate the scattering properties such as the scattering coefficient and the phase function in colloidal suspensions. The IST assumes no wave interactions between the scatterers, while the DST considers the wave interactions. In a dilute colloidal suspension, the scattering properties calculated from the IST coincide well with those from the inverse analysis [16]. In a dense suspension, however, the IST breaks down because the theory does not take the interference effects of the scattered fields into account. Most biological tissue volumes and agricultural products are dense media and the interference effects strongly influence the scattering properties, meaning the necessity of the use of the DST. Although various kinds of formulations for the DST have been developed, the formulation by Cartigny $et~al.$ [17] has been widely used in NIRI and NIRS [18,19], which originates from the research field of thermal engineering, because their formulation is mathematically simple by assuming a monodisperse system (single size of particles). Recently, it was reported that the scattering properties using the DST in the Cartigny’s formulation agree well with those by the inverse analysis for a dense suspension at the volume fraction up to 20% [19].

Despite the wide use of the Cartigny’s formulation, an extension of their formulation to a polydisperse system is not straightforward because of its phenomenological derivation. Most colloidal suspensions are polydisperse systems and have particle size distributions. Hence, one needs to construct the DST for a polydisperse suspension in a different formulation and to examine the influence of the polydispersity on the scattering properties. In the research field of remote sensing, the DST for a polydisperse system has been successfully developed based on the Foldy-Lax equation (FLE) and applied to snow layers [20–22]. The application of the DST based on the FLE to the research fields of NIRI and NIRS is significant because the wavelengths of the electromagnetic waves and particle sizes of the target media largely differ between the two research fields, indicating differences in the results for the scattering properties.

For modeling of photon scattering based on the RTE, the Henyey-Greenstein (HG) function [23] has been widely used as a formulation of the normalized phase function in NIRI and NIRS [24–27], because the function is a simple mathematical model to treat the highly forward scattering of photons for random media at the near-infrared wavelength range by only one parameter. However, it was reported that the shape of the HG function largely deviates from the experimental result for the phase function in colloidal suspensions, while the result calculated from the EMT coincides well with the experimental result [28]. This fact indicates that the HG function is not the high-fidelity model of the photon scattering. Hence, one needs a validation test for the HG function by comparing solutions of the RTE for the HG function with those for the phase function calculated from the EMT. Nevertheless, the validation test for the HG function has been little reported from our best of knowledge.

Our objective is to numerically examine the influences of the polydispersity and interference effects on photon scattering and photon transport for dense polydisperse colloidal suspensions at the near-infrared wavelength range by the time-dependent RTE combined with the IST or DST. In the combined model, we employ the analytical solution of the RTE, developed by Liemert and Kienle [29], by modifying parts of the normalized phase function. In the RTE-calculations, we use the results of the scattering properties calculated from the IST or DST. The other objective is to examine the validation of the HG function by comparing with the RTE-results for the phase function using the DST.

2. Models for photon transport and scattering properties

2.1 Radiative transfer theory (RTT)

2.1.1 Radiative transfer equation (RTE)

We consider photon transport in the length and time scales of the mean free path and mean time of flight as shown in the left-side of Fig. 1. Then, photon transport can be described by the RTE [10]. The time-dependent RTE in 3D random media is given by

 

Fig. 1. Schematics of photon transport in a random medium on the macroscopic scale for the RTT (left) and photon scattering in colloidal suspensions for monodisperse and polydisperse systems on the microscopic scale for the EMT (center and right).

Download Full Size | PDF

The normalized phase function $\hat {P}(\boldsymbol {\Omega } \cdot \boldsymbol {\Omega }')$ represents a probability of single scattering with the change in the direction from $\boldsymbol {\Omega }'$ into $\boldsymbol {\Omega }$ and satisfies the condition of $\int _{\mathbb {S}^{2}} d\Omega \hat {P}(\boldsymbol {\Omega } \cdot \boldsymbol {\Omega }')=1$. Also, it is often assumed that the function depends only on the scattering angle $\cos ^{-1}{(\boldsymbol {\Omega } \cdot \boldsymbol {\Omega }')}$. In NIRI and NIRS, the HG function has been widely used as a formulation of $\hat {P}(\boldsymbol {\Omega } \cdot \boldsymbol {\Omega }')$ [23]:

In the ballistic regime (short length and time scales), photon scattering is anisotropic and quantified by $\mu _s$, $g$, and $\hat {P}(\boldsymbol {\Omega } \cdot \boldsymbol {\Omega '})$. In the diffusive regime (long scales), meanwhile, it is isotropic and quantified by the reduced scattering coefficient $\mu _s'=(1-g)\mu _s$. The crossover length $r_{DA}$ from the ballistic to diffusive regimes has been evaluated as approximately $10/\mu _t'$ with $\mu _t'=\mu _a+\mu _s'$ [12,13,27].

2.2 Dependent and independent scattering theories (DST and IST)

2.2.1 Foldy-Lax equation (FLE)

As a target medium, we consider that the dielectric spherical particles are suspended in water with the total number of particles $N$. Multiple scattering of the electric fields in the colloidal suspension can be described exactly by the FLE [22,30,31]:

By the first order solution of the FLE in far-field, we derived theoretical forms of the scattering properties ($\mu _s$, $g$, and $\hat {P}$) for the DST in monodisperse and polydisperse systems of the colloidal suspensions. For the derivation, we used the spherical coordinates $(r, \theta , \phi )$ where $r$ represents the radial distance, $\theta \in [0, \pi ]$ the polar angle, and $\phi \in [0, 2\pi ]$ the azimuthal angle. We also considered the incident plane wave propagating in the $z$-direction and polarized in the $x$-direction. Then, we can adjust the scattering angle to be the same as the polar angle $\theta$, so that the normalized phase function is denoted by $\hat {P}(\theta )$. We derived theoretical forms of the scattering properties for the IST by neglecting the terms of the interference effects from those for the DST. We denote the forms of the scattering properties for the four cases of theories and systems by D-M, I-M, D-P, and I-P, e.g., D-M represents the forms using the DST for a monodisperse system.

2.2.2 Monodisperse systems

In this subsection, we provided the theoretical forms of the scattering properties using the DST and IST for a monodisperse system, denoted by D-M and I-M, before discussion for a polydisperse system. The system consists of colloidal particles with a single diameter $d$ as shown in the center of Fig. 1.

The phase function $P_{D-M}(\theta , \phi )$ for the D-M case is given as

The phase function $P_{I-M}(\theta , \phi )$ for the I-M case is given by omitting the second term of Eq. (5) because the first and second terms of Eq. (5) correspond to components of the independent scattering and interference effects, respectively. Then, the normalized phase function $\hat {P}_{I-M}(\theta )$ for the I-M case is equivalent to that using the Mie theory:

The scattering coefficient $\mu _s^{D-M}$ for the D-M case is given by the integration of $P_{D-M}(\theta , \phi )$ over the whole solid angle:

From the definition (Eq. (3)), the anisotropy factor $g_{D-M}$ for the D-M case is given as

2.2.3 Polydisperse systems

In this subsection, we discuss theoretical forms of the scattering properties using the DST and IST for a polydisperse system, denoted by D-P and I-P. The system has $N_s$-kinds of particle diameters $d_{\alpha }$ ($\alpha =1,2,\ldots , N_s$) as shown in the right-side of Fig. 1. The total number of the particles $N$ in the system is given as $\sum _{\alpha =1}^{N_s} N_{\alpha }$, and $N_{\alpha }$ is the total number of particles for the diameter $d_{\alpha }$ of the $\alpha$-class.

The phase function $P_{D-P}(\theta , \phi )$ for the D-P case is given as

The phase function $P_{I-P}(\theta , \phi )$ for the I-P case is obtained by omitting the second term of Eq. (10). From $P_{I-P}(\theta , \phi )$, the normalized phase function $\hat {P}_{I-P}(\theta )$ is obtained as

The scattering coefficient $\mu _s^{D-P}$ for the D-P case is given as

The anisotropy factors $g_{D-P}$ and $g_{I-P}$ for the D-P and I-P cases can be calculated from $\hat {P}_{D-P}(\theta )$ and $\hat {P}_{I-P}(\theta )$ (Eq. (11)) in the same way for the calculation of $g_{D-M}$ (Eq. (9)).

3. Numerical methods and conditions

3.1 Scattering properties

We consider monochromatic scattering of photons in the wavelength of $\lambda = 759$ nm, which is the same as that for the systems of NIRI and NIRS in time-domain [3,11].

In investigations of the scattering properties, we adopted the volume fraction $\rho _v$ as a controlled parameter and varied it from 0.1 to 20%. $\rho _v$ is given as $n_0 \pi d^{3}/6$ for a monodisperse system, and as $n_0 \pi d^{3}_{mean}/6$ with the (number-weighted) mean diameter $d_{mean}=\sum _{\alpha =1}^{N_s} d_{\alpha } n_{\alpha }/\sum _{\alpha =1}^{N_s} n_{\alpha }$ for a polydisperse sytem. In calculations of $H_M(\theta ;d, n_0)$ for the scattering properties, we employed the monodisperse Percus-Yevick (PY) model [35]. The PY model can describe the statistically averaged configuration of the particles in a liquid state. In calculations of $H_{\alpha \beta }(\theta )$ for the scattering properties, we employed the polydisperse PY model, which is developed by Baxter using the Wiener-Hopf factorization to the generalized Ornstein-Zernike equation [36]. For the details of the calculations of $H_{\alpha \beta }(\theta )$, please see the Ref. [22]. The monodisperse and polydisperse PY models have been extensively validated by actual measurements, molecular dynamics simulations, and Monte Carlo simulations [34,37,38]. The source codes for the numerical calculations of the scattering properties were written in MATLAB, and the calculations for the Mie theory were based on the open-sources by the Refs. [39,40].

3.2 Physical properties of colloidal suspensions in water

We numerically treated water suspensions of two kinds of spherical colloidal particles: Intralipid-10% (Kabivitrum, Stockholm, Sweden) and Psi-1.5 (Kisker Biotech, Steinfurt, Germany). We used the measurement values of the physical properties (particle size distribution and refractive index) for the colloidal suspensions from the Refs. [16,19]. Figure 2 shows the particle size distributions for Intralipid and Psi-1.5, where the fraction $f_d(d_{\alpha })$ of the distribution is normalized so that $\sum _{\alpha =1}^{N_s}f_d(d_{\alpha })=1$ with $f_d(d_{\alpha })=n_{\alpha } n_0^{-1}$. The mean diameters $d_{mean}$ for Intralipid and Psi-1.5 are 97 and 1213 nm, respectively. Intralipid has the samller mean diameter and larger deviation of the distribution than Psi-1.5. In the calculations of the scattering properties, the mean diameters were used for the monodisperse systems, while the particle size distributions were used for the polydisperse systems. Refractive indices of Intralipid, Psi-1.5, and water at $\lambda = 759$ nm are given as $n_I=1.47$, $n_P=1.44$, and $n_W=1.33$, respectively.

 

Fig. 2. Particle size distributions for the two kinds of colloidal suspensions in water: Intralipid [16] and Psi-1.5 [19]. The fraction $f_d$ at different diameters $d_{\alpha }$ is normalized so as to satisfy the condition of $\sum _{\alpha =1}^{N_s}f_d(d_{\alpha })=1$.

Download Full Size | PDF

As shown in the left-side of Fig. 2, a particle size class for Intralipid below 25 nm has not been measured in the Ref. [16] because the size classes below 25 nm have been disregarded on the electron-microscopy photographs. We preliminarily confirmed that if we additionally consider a size class below 25 nm, the results for the scattering properties are little changed.

3.3 Analytical solutions of the RTE

We employed the analytical solutions of the time-dependent RTE (Eq. (1)) for 3D infinite homogeneous media, originally developed by Liemert and Kienle [29]. Their analytical solution has been independently verified by comparisons with Monte Carlo simulations and experimental data for a colloidal suspension [41], and numerical solutions for a finite medium [42]. In the comparisons, a light source was located inside the medium to suppress the boundary effects. Their analytical solution uses the HG function as the normalized phase function and expands it by a finite series of (unassociated) Legendre polynomials. We modified their solution for applying the phase functions using the DST and IST. We decomposed the phase functions by the delta-Eddington method (dEM), which enables more accurate treatments of the highly forward-peaked shape of the phase function ($g \sim 1$) than the finite order Legendre expansion [27]. In the dEM, we calculated the expansion coefficient $C_l$ ($l=0,1,\ldots , N_C$) of $\hat {P}(\theta )$ by Legendre polynomials $P_l(\cos {\theta })$:

In the calculations of the RTE, we employed the three kinds of photon transport models as listed in Table 1: D-P model, I-P model, and HG model. We used the results of the scattering properties for the D-P and I-P cases at different volume fractions. Meanwhile, the absorption coefficient was fixed to $\mu _a=0.1\;\textrm {cm}^{-1}$, assuming that an absorber such as India ink is added to the colloidal suspension.

Table 1. Photon transport models based on the RTE with the optical properties: absorption coefficient $\mu _a$, scattering coefficient $\mu _s$, anisotropy factor $g$, and normalized phase function $\hat {P}(\theta )$.

View Table | View all tables in this article

An isotropic point source was incident at the initial time and at the origin of the $xyz$-coordinate. A temporal profile of the fluence rate $\Phi (r_{SD},t)=\int _{\mathbb {S}^{2}} d\Omega I(\boldsymbol {r}, \boldsymbol {\Omega }, t)$ was calculated at the distance $r_{SD}=|\boldsymbol {r}|$ between the source and detector. The source-detector (SD) distance $r_{SD}$ was basically set as $3.5/\mu _t'$ because we aim to examine the influences of the normalized phase function on the RTE-results in the ballistic regime, less than the crossover length $r_{DA}\sim 10/\mu _t'$. It was reported [27] that at $r_{SD} \sim 3.5/\mu _t'$, the RTE-results largely differed between the cases for the HG function (Eq. (2)) and a low order approximation of the HG function, suggesting the strong influences of a shape of the phase function on the RTE-results. Because the scattering properties largely vary with volume fractions and the ballistic regime depends on the scattering properties, the values of $r_{SD}$ also vary with volume fractions. In the calculations of $\mu _t'$ for $r_{SD}$, $\mu _s^{D-P}$ and $g_{D-P}$ were used because the D-P model is considered as the most physically accurate model among the photon transport models. However, at a quite low or high volume fraction, the value of $3.5/\mu _t'$ is out of range of the SD distance in actual measurements, because the $\mu _s'$-value is quite small or large. Here, we consider in actual measurements, the SD distance roughly ranges from 0.1 cm to 5.0 cm based on the references for time-resolved light intensity measurements [8,41]. This fact means that the RTE-results at $r_{SD}=3.5/\mu _t'$ cannot be observed in the actual measurements at a low or high volume fraction. To overcome the difficulty, we restricted the $r_{SD}$-value in the actual range from 0.1 to 5.0 cm: $r_{SD}=\max [3.5/\mu _t',\; 0.1\;\textrm {cm}]$ or $\min [3.5/\mu _t',\; 5.0\;\textrm {cm}]$. Due to the restriction, for certain cases of high volume fractions, $r_{SD}$ is comparable to $r_{DA}$, corresponding to the RTE-results in the diffusive regime.

We evaluated the differences in the $\Phi$-results between the photon transport models by $\Delta \Phi$ [%]:

Table 2 summarizes symbols for the scattering properties and the RTE, and abbreviations for the theories and systems.

Table 2. Symbols for the scattering properties and the RTE, and abbreviations for the theories and systems

View Table | View all tables in this article

4. Numerical results

4.1 Scattering properties

4.1.1 Total correlation function

Firstly, we investigated Fourier transforms of the total correlation functions $H_M(\theta )$ for a monodisperse system and $H_{\alpha \beta }(\theta )$ for a polydisperse system using the PY models in the two kinds of colloidal suspensions: Intralipid and Psi-1.5 as shown in Fig. 3. The results for $H_M(\theta )$ and $H_{\alpha \beta }(\theta )$ are ascribed to superposition of the scattered fields from the colloidal particles, reflecting the interference effects. The interference effects become enhanced when the volume fraction $\rho _v$ increases and the colloidal particles become closer to each other.

 

Fig. 3. $H_M(\theta )$ for a monodisperse system at different volume fractions $\rho _v=1$, 10 and $20\%$ (upper row) and $H_{\alpha \beta }(\theta )$ for a polydisperse system at $\rho _v=20\%$ (bottom row) as a function of the scattering angle $\theta$ in Intralipid (left column) and Psi-1.5 (right column). For monodisperse systems, the (number-weighted) mean diameters $d_{mean}$ were used: 97 nm for Intralipid and 1213 nm for Psi-1.5, respectively, calculated from the particle size distributions as shown in Fig. 2. In the figures for $H_{\alpha \beta }(\theta )$, the value of a pair of the diameters $d_{\alpha }-d_{\beta }$ is denoted as a legend.

Download Full Size | PDF

As shown in Figs. 3(a) and (b), the results for $H_M(\theta )$ in Intralipid are largely different from those in Psi-1.5 mainly because of the differences in the mean diameters. In general, $H_M(\theta )$ has a peak around the angle $\theta _p = 2 \sin ^{-1}{[\lambda /(2n_W d_{mean})]}$ due to the correlation of the local density at the nearest neighbor distance. For Intralipid, the $\theta _p$-value is evaluated as over $180^{\circ }$, so that the $H_M(\theta )$-result has no peak in the $\theta$-range. For Psi-1.5, meanwhile, the sharp peak is observed around $\theta _p = 25^{\circ }$ in $H_M(\theta )$. As the $\rho _v$-value increases, the $H_M(\theta )$-result for Intralipid decreases, while the peaks for Psi-1.5 in $H_M(\theta )$ become sharper.

As shown in Figs. 3(c) and (d), the results for $H_{\alpha \beta }(\theta )$ are different from those for $H_M(\theta )$, especially in Intralipid, meaning the influence of the polydispersity. This is because $H_{\alpha \beta }(\theta )$ depends on the diameters $d_{\alpha }$ and $d_{\beta }$ and the densities $n_{\alpha }$ and $n_{\beta }$ unlike $H_M(\theta )$.

4.1.2 Normalized phase function

We investigated the normalized phase functions for the four cases of theories and systems: $\hat {P}_{I-M}(\theta )$, $\hat {P}_{D-M}(\theta )$, $\hat {P}_{I-P}(\theta )$, and $\hat {P}_{D-P}(\theta )$ in Intralipid and Psi-1.5 as shown in Fig. 4. The formulations of the phase functions are referred to Secs. 2.2.2 and 2.2.3. $\hat {P}_{I-M}(\theta )$ and $\hat {P}_{I-P}(\theta )$ are independent of $\rho _v$, while $\hat {P}_{D-M}(\theta )$ and $\hat {P}_{D-P}(\theta )$ depend on it.

 

Fig. 4. Normalized phase functions $\hat {P}(\theta )$ for the four cases of theories and systems: the I-M and D-M cases (upper row), I-P and D-P cases (bottom row) in Intralipid (left column) and Psi-1.5 (right column) at $\rho _v=1$ and 20%. For the formulations of the phase functions, please refer to Secs. 2.2.2 and 2.2.3.

Download Full Size | PDF

All the subfigures of Fig. 4 show the nice agreement of the phase functions between the IST and DST at $\rho _v=1\%$ for each case of the systems and suspensions. At $\rho _v=20\%$, meanwhile, the phase functions using the IST ($\hat {P}_{I-M}(\theta$) or $\hat {P}_{I-P}(\theta )$) differ from those using the DST ($\hat {P}_{D-M}(\theta$) or $\hat {P}_{D-P}(\theta )$). In the angle range of $0\;<\theta \; \lesssim \;30^{\circ }$ (forward scattering components), the values of the phase functions using the DST at $\rho _v=20\%$ are reduced from those at $\rho _v=1\%$, due to the negative values of $H_M(\theta )$ and $H_{\alpha \beta }(\theta )$ in the angle range. These results suggest the less influence of the interference effects at $\rho _v=1\%$ and the strong influence at $\rho _v=20\%$.

Figures 4(a) and (c) show the differences in the phase functions for Intralipid between the monodisperse and polydisperse systems. For Psi-1.5, meanwhile, the agreement of the phase functions between the two systems is seen in Figs. 4(b) and (d). These results suggest the strong influence of the polydispersity for Intralipid and the less influence for Psi-1.5.

4.1.3 Scattering coefficient, anisotropy factor, and reduced scattering coefficient

We examined the scattering coefficients $\mu _s$, anisotropy factors $g$, and reduced scattering coefficients $\mu _s'$ for the eight cases of theories, systems, and suspensions in the range of $\rho _v$ (0.1, 1.0, 2.0, 3.0, $\cdots$, 20%) as shown in Fig. 5. The theoretical forms of $\mu _s$ and $g$ are referred to Secs. 2.2.2 and 2.2.3.

 

Fig. 5. Scattering coefficients $\mu _s$, anisotropy factors $g$, and reduced scattering coefficients $\mu _s'=(1-g)\mu _s$ at different volume fractions $\rho _v$ (0.1, 1.0, 2.0, 3.0, $\cdots$, 20%) for the four cases of theories and systems in Intralipid and Psi-1.5. The theoretical forms of $\mu _s$ and $g$ are referred to Secs. 2.2.2 and 2.2.3. The other details are the same as Fig. 4.

Download Full Size | PDF

Figures 5(a)-(d) show that the results for $\mu _s$ and $g$ using the DST (the D-M or D-P case) are smaller than those using the IST (the I-M or I-P case) at the larger $\rho _v$-values than approximately 2% for each case of the systems and suspensions. The reductions of the $\mu _s$- and $g$-values are due to the interference effects, i.e., the reduction of the forward scattering components for the phase function as shown in Fig. 4.

As shown in the results for Intralipid (Figs. 5(a) and (c)), the strong influence of the polydispersity is observed, that is, the large differences in the results of $\mu _s$ and $g$ between the monodisperse and polydisperse systems. For Psi-1.5, meanwhile, the polydispersity little influences the $\mu _s$- and $g$-results as shown in Figs. 5(b) and (d).

We calculated $\mu _s'=(1-g)\mu _s$ from the results of $\mu _s$ and $g$ as shown in Figs. 5(e) and (f). The differences in the $\mu _s'$-values between IST and DST are smaller than those in the $\mu _s$-values, especially for Psi-1.5 (Figs. 5(b) and (f)), because the $g$-value is closed to unity. These results suggest the less influence of the interference effects on photon transport in the diffusive regime, compared with that in the ballistic regime.

We discuss the influences of the polydispersity for the scattering properties by examining the scattering coefficients $\mu _s^{I-P}$ and $\mu _s^{I-M}$ for the I-P and I-M cases at $\rho _v=1\%$, because the influences of the polydispersity are observed at the low volume fraction and the discussion based on the IST is simpler than that on the DST. Figure 6 shows the distributions of $n_{\alpha } \sigma _{s, Mie}^{\alpha }$ at different diameters $d_{\alpha }$ for Intralipid and Psi-1.5 and compares with the results of $\mu _s^{I-M}=n_{0} \sigma _{s, Mie}$ at the mean diameters $d_{mean}$, indicated by the cross symbols. The values of $n_{\alpha } \sigma _{s, Mie}^{\alpha }$ are components of $\mu _s^{I-P}$, i.e., $\mu _s^{I-P}=\sum _{\alpha =1}^{N_s} n_{\alpha } \sigma _{s, Mie}^{\alpha }$. For Intralipid, the value of $d_{mean}=97$ nm is located on the left-side (smaller diameter class) in the $n_{\alpha } \sigma _{s, Mie}^{\alpha }$-distribution, and the $\mu _s^{I-M}$-value is less than the value of $n_{\alpha } \sigma _{s, Mie}^{\alpha }$ at the diameter class above 275 nm. These results lead to the difference between $\mu _s^{I-P}$ and $\mu _s^{I-M}$ for Intralipid. For Psi-1.5, meanwhile, $d_{mean}=1213$ nm is almost the same as the diameter of 1210 nm for the maximum value of $n_{\alpha } \sigma _{s, Mie}^{\alpha }$ and located at the center of the $n_{\alpha } \sigma _{s, Mie}^{\alpha }$-distribution. These results lead to the agreement between $\mu _s^{I-P}$ and $\mu _s^{I-M}$ for Psi-1.5. The difference in the results between Intralipid and Psi-1.5 is ascribed to the difference in the mean diameters or the size parameters $X=\pi n_W d/\lambda$ and variances of the particle size distributions. The values of $X$ at $d_{mean}$ are given as 0.5 for Intralipid and 6.7 for Psi-1.5, respectively. The $\sigma _{s, Mie}$-value increases as the $X$-value increases in the $X$-range of our interest [33]. For Intralipid, the $X$-value at $d_{mean}$ is less than unity and the variance of the particle size distribution is large. Hence, the larger diameter class than $d_{mean}$ strongly contributes to the $n_{\alpha } \sigma _{s, Mie}^{\alpha }$-result, although the fraction of the particle size distribution at the larger diameter class is quite small as shown in Fig. 2. For Psi-1.5, meanwhile, the $X$-value at $d_{mean}$ is larger than unity and the variance of the particle size distribution is small. Hence, the shape of the $n_{\alpha } \sigma _{s, Mie}^{\alpha }$–distribution is similar to that of the particle size distribution and $\mu _s^{I-P}$ is characterized by the mean value $\mu _s^{I-M}$.

 

Fig. 6. $n_{\alpha } \sigma _{s, Mie}^{\alpha }$ at different diameters $d_{\alpha }$ for $\rho _v=$1% in Intralipid and Psi-1.5 with the number of density $n_{\alpha }$ for the $\alpha$-class. Cross symbols represent the values of $\mu _s^{I-M}=n_0 \sigma _{s, Mie}$ at the mean diameters $d_{mean}$ with the mean number of density $n_0$.

Download Full Size | PDF

4.2 Fluence rate

4.2.1 Influences of the interference effects

We examined the influence of the interference effects on photon transport by calculating the time-dependent RTE (Eq. (1)) for Intralipid and Psi-1.5. We compared the RTE-results for the I-P and D-P models as listed in Table 1. Although not shown here, we confirmed the strong influence of the polydispersity on photon transport for Intralipid, while the little influence for Psi-1.5 as clearly suggested in Fig. 5.

As shown in Fig. 7(a), the temporal profiles of the normalized fluence rates $\hat {\Phi }(t)$ for the I-P model differ from those for the D-P model in the two suspensions, suggesting the influence of the interference effects. At $\rho _v=20$%, the values of $\mu _s$ or $\mu _s'$ are so large for the two suspensions that the crossover lengths $r_{DA}$ from the ballistic to diffusive regimes are quite short: 0.07 cm for Intralipid and 0.12 cm for Psi-1.5, respectively. Then, the source-detector distances $r_{SD}$ are comparable to the crossover lengths. Hence, the diffusion approximation is almost valid in the $\hat {\Phi }(t)$-calculations, and the differences in the $\hat {\Phi }(t)$-results between the two models come from the differences in $\mu _s'$-values.

 

Fig. 7. (a)Temporal profiles of the fluence rates $\hat {\Phi }(t)$ normalized by their peak values calculated from the RTE (Eq. (1)) based on the I-P and D-P models as listed in Table 1 for Intralipid and Psi-1.5 at $\rho _v=20$%. (b)Differences in the fluence rates $\Delta \Phi$ (Eq. (14)) between the I-P and D-P models at different $\rho _v$-values (0.1, 1.0, 2.0, $\cdots$, 20%).

Download Full Size | PDF

We calculated the differences in the fluence rates $\Delta \Phi$ between the I-P and D-P models using Eq. (14) at the different $\rho _v$-values (0.1, 1.0, 2.0, $\cdots$, 20%) as shown in Fig. 7(b). Similarly to the results of the scattering properties, the $\Delta \Phi$-value increases as the $\rho _v$-value increases for the two suspensions. The $\Delta \Phi$-values for Intralipid are larger than those for Psi-1.5, especially at a high $\rho _v$-value, because of the larger differences in the $\mu _s'$-values for Intralipid than those for Psi-1.5.

4.2.2 Validation test of the HG function

We examined the influence of the normalized phase functions on photon transport by the RTE-calculations for the HG and D-P models as listed in Table 1.

As shown in Fig. 8(a), the shapes of $\hat {P}_{HG}(\theta )$ largely deviate from those of $\hat {P}_{D-P}(\theta )$, although the $g$-values are the same for the two models. At $\rho _v=1$%, the relative differences in the phase functions between the two models are evaluated as approximately 26% for Intralipid and 306% for Psi-1.5, respectively. The differences slightly decrease to 18% and 295% when the $\rho _v$-value increases to 20%. Figure 8(b) shows the expansion coefficients $C_l$ (Eq. (13)) of $\hat {P}_{HG}(\theta )$ and $\hat {P}_{D-P}(\theta )$ by Legendre polynomials at $\rho _v=1$%. The values of $C_l$ at $l \geq 2$ differ between the two models, especially for Psi-1.5, because the shapes of the phase functions for Psi-1.5 are sharper than those for Intralipid.

 

Fig. 8. (a)HG functions $\hat {P}_{HG}(\theta )$ (Eq. (2)) and normalized phase functions $\hat {P}_{D-P}(\theta )$ for the D-P case at $\rho _v=1$% in Intralipid and Phi-1.5. (b)Expansion coefficients $C_l$ (Eq. (13)) of $\hat {P}_{HG}(\theta )$ and $\hat {P}_{D-P}(\theta )$ at different orders $l$. (c)Temporal profiles of the normalized fluence rates $\hat {\Phi }(t)$ for the HG and D-P models as listed in Table 1. (d)Differences in the fluence rates $\Delta \Phi$ (Eq. (14)) between the HG and D-P models at different $\rho _v$-values.

Download Full Size | PDF

Figure 8(c) compares the temporal profiles of the normalized fluence rates $\hat {\Phi }(t)$ for the HG model with those for the D-P model at $\rho _v=1$% in Intralipid and Psi-1.5. The peak times of $\hat {\Phi }(t)$ are almost the same between the two models, probably because the values of $\mu _s$ and $g$ are the same between the two models. Meanwhile, the $\hat {\Phi }(t)$-profiles for the D-P model decay faster than those for the HG model, especially in Psi-1.5, suggesting the influence of the differences in the $C_l$-values at $l \geq 2$ between the two models. Figure 8(d) shows the differences in the fluence rates $\Delta \Phi$ (Eq. (14)) between the HG and D-P models at different $\rho _v$-values. The $\Delta \Phi$-values are approximately 4% for Intralipid at the range of $0\;<\rho _v \lesssim 4\%$ and 12% for Psi-1.5 at the range of $0\;<\rho _v \lesssim 8\%$, respectively. These results suggest the influence of the phase functions on the RTE-results at a small $\rho _v$-value corresponding to the ballistic regime ($r_{SD}\;<\; r_{DA}$) and the limitation of the HG model. At a large $\rho _v$-value, meanwhile, the $\Delta \Phi$-values are reduced to less than 1% for the two suspensions. This is because the diffusion approximation is almost valid at a large $\rho _v$-value and the $\mu _s'$-values are the same between the two models. These results suggest that the HG model is valid as a photon transport model at a large $\rho _v$-value if the $g$-value is exactly predetermined.

5. Conclusions

We have developed photon transport models for the two kinds of polydisperse colloidal suspensions: Intralipid and Psi-1.5 by the time-dependent RTE combined with the DST or IST at the near-infrared wavelength in the range of the volume fraction from 0.1 to 20%.

We have shown the polydispersity and interference effects strongly influence the results of the scattering properties and the RTE for Intralipid, whose mean diameter is smaller and variance of the particle size distribution is larger than those for Psi-1.5. These results suggest the necessity of incorporating the polydispersity and interference effects into photon transport models for dense random media such as biological tissue volumes and agricultural products toward further improvements of NIRI and NIRS.

We have compared the RTE-results for the HG and D-P models to examine the validity of the HG function, which is conventionally used as a formulation of the normalized phase function. The RTE-results differ between the two models at a low volume fraction, especially for Psi-1.5, whose anisotropy factor is closed to unity, although the values of the scattering coefficient and anisotropy factor are the same between the two models. These results suggest the limitation of the HG function for highly forward scattering media in the ballistic regime.

A numerical calculation of the RTE combined with the IST or DST for a (more realistic) finite medium with the source and detector located at the boundary is challenging as a future work. In the accurate calculations of the RTE, the highly forward-peaked phase function should be treated appropriately based on the weighting procedure for the phase function (e.g., [25,42]).

Our developed model can directly compare the temporal profiles of the light intensity with the results of time-resolved light intensity measurements. The direct comparison enables the validation test of the photon transport model without introducing of the inverse analysis. The conventional photon transport model does not enable the direct comparison. Our works can be extended to a multi-scale modeling of photon transport from the EMT-scale to the RTT-scale in various kinds of materials such as biological tissue, agricultural products, and foods, essential in the research fields of diffuse optical imaging and spectroscopy. The extended model enables better understanding of the signals of NIRI and NIRS for the materials.