Numerical Study of Near-Infrared Light Propagation in Aqueous Alumina Suspensions Using the Steady-State Radiative Transfer Equation and Dependent Scattering Theory

Abstract

Understanding light propagation in liquid phantoms, such as colloidal suspensions, involves fundamental research of near-infrared optical imaging and spectroscopy for biological tissues. Our objective is to numerically investigate light propagation in the alumina colloidal suspensions with the mean alumina particle diameter of 55 nm at the volume fraction range 1–20%. We calculated the light scattering properties using the dependent scattering theory (DST) on a length scale comparable to the optical wavelength. We calculated the steady-state radiative transfer and photon diffusion equations (RTE and PDE) using the DST results based on the finite difference method in a length scale of the mean free path. The DST calculations showed that the scattering and reduced scattering coefficients become more prominent at a higher volume fraction. The anisotropy factor is almost zero at all the volume fractions, meaning the scattering is isotropic. The comparative study of the RTE with the PDE showed that the diffusion approximation holds at the internal region with all the volume fractions and the boundary region with the volume fraction higher than 1%. Our findings suggest the usefulness of the PDE as a light propagation model for the alumina suspensions rather than the RTE, which provides accurate but complicated computation.

1. Introduction

Near-infrared imaging and spectroscopy offer chemical components (e.g., hemoglobin concentration) and structural properties (e.g., cellular density) for biological tissues and agricultural products [1,2]. Among the techniques, near-infrared optical tomography (NIROT), conventionally called diffuse optical tomography, has the potential to image a deep region of tissue volumes [3,4], such as the brain in the human head and the trachea in the human neck. For the development of NIROT, elucidation of light propagation in a target medium is indispensable because NIROT employs a theoretical model to describe light propagation, and image qualities by NIROT depend on the accuracy of the light propagation model [5].

In the near-infrared wavelength range, the scattering process of light propagation is dominant over the absorption process. The radiative transfer theory (RTT) can describe light propagation in a target medium on a length scale of the mean free path between light and the medium. As the RTT, the radiative transfer equation (RTE) and photon diffusion equation (PDE) have been widely used. The RTE can accurately describe light propagation, but it requires complicated computation and high computational loads. The PDE, the diffusion approximation to the RTE, provides fast calculation with a simple numerical treatment but is invalid at several conditions, such as low-scattering medium [5]. A target medium is regarded as a scattering medium (continuous medium) with optical properties on the length scale. The optical properties (e.g., scattering and absorption coefficients) specify the strengths of light scattering and absorption. The optical properties have been evaluated by inverse analysis as parameters for biological tissues and agricultural products [6,7,8]. The inverse analysis optimizes a difference in light intensity between measured and calculated results using the RTT and incorporates multiple computations of light propagation. This fact suggests that an efficient model of light propagation is desirable in addition to accuracy. Because the validation of the results by inverse analysis is quite difficult for biological tissues and agricultural products, the validation test is necessary for tissue phantoms [9].

As liquid tissue phantoms, colloidal suspensions have been widely used. Colloidal particles correspond to scatterers, and the scattering properties vary with the volume fractions of the particles. The electromagnetic theory (EMT) can calculate the scattering properties for colloidal suspensions on a length scale comparable to the optical wavelength. As the EMT, the independent scattering theory (IST) and dependent scattering theory (DST) have been developed. The IST treats no interaction of the electric fields scattered by particles, while the DST treats the interference of the fields in the far-field [10,11]. The scattering properties using the IST agree well with those by the inverse analysis at a lower volume fraction than approximately 5%, while the results using the DST do at a high volume fraction up to approximately 20% [12,13,14]. For a dilute suspension (low volume fraction), it is easy to adjust values of the scattering properties because of the linear dependence of the scattering properties on the volume fraction. Meanwhile, the adjustment is difficult for a dense suspension (high volume fraction) due to a complicated dependence of the scattering properties on the volume fraction. Hence, clarifying the complicated dependence is significant as fundamental research of optical imaging techniques. Many researchers have examined the scattering properties for dense suspensions by the DST [10,13,14,15,16,17,18,19]. Meanwhile, research for the dependence of light propagation on the volume fraction has been less reported to the best of our knowledge. In the research field of microwave remote sensing, Tsang and co-workers have examined radiative transfer for snow layers using the RTE combined with the DST (RTE-DST) [11]. Recently, some authors have examined the near-infrared light propagation by the RTE-DST for dense colloidal suspensions based on the analytical solutions of the time-dependent RTE [19]. However, the dependence of light propagation on the volume fraction has not been fully clarified. In particular, a range of the volume fraction where the diffusion approximation holds, is unclear, although the length and time scales for the diffusion approximation have been extensively discussed [1,20,21,22]. Because the volume fraction is a control parameter for phantom experiments, the information of the volume fraction range is helpful.

In this study, we numerically calculated the steady-state RTE-DST for aqueous alumina suspensions, widely used as a liquid phantom [9], at different volume fractions up to 20%. The numerical schemes for the RTE-DST have been little discussed in previous researches. We investigated the volume fraction range where the diffusion approximation holds by comparing numerical results using the PDE combined with the DST (PDE-DST). The RTE-DST employs the scattering coefficient and phase function calculated from the DST, while the PDE-DST does the reduced scattering coefficient only.

2. Theoretical Models of Light Propagation and Light Scattering

2.1. Radiative Transfer Equation (RTE)

The RTE, a linear Boltzmann equation, can describe light propagation [23] on a length scale of the mean free path between light and a medium. On the length scale, a colloidal suspension is considered as a scattering medium (continuous medium) with the optical properties as shown in Figure 1a. In this study, we considered homogeneous scattering media. For a 3D medium, the steady-state RTE is given by

$$\left[\mathbf{\Omega }·\nabla +{\mu }_a\left(\mathit{r}\right)+{\mu }_s\left(\mathit{r}\right)\right]I(\mathit{r},\mathbf{\Omega })={\mu }_s\left(\mathit{r}\right){\int }_{4\pi }d{\mathbf{\Omega }}^{\prime }p(\mathbf{\Omega }·{\mathbf{\Omega }}^{\prime })I(\mathit{r},{\mathbf{\Omega }}^{\prime })+q(\mathit{r},\mathbf{\Omega }),$$

(1)

where $I(\mathit{r},\mathbf{\Omega })$ in $\mathrm{J}{\mathrm{cm}}^{-2}{\mathrm{str}}^{-1}$ represents the light intensity as a function of spatial position vector $\mathit{r}=(x,y,z)\in {\mathbb{R}}^3$ in $\mathrm{cm}$; and angular direction (unit direction vector) $\mathbf{\Omega }=({\mathbf{\Omega }}_x,{\mathbf{\Omega }}_y,{\mathbf{\Omega }}_y)$ in str. ${\mu }_a\left(\mathit{r}\right)$ and ${\mu }_s\left(\mathit{r}\right)$ in ${\mathrm{cm}}^{-1}$ are the absorption and scattering coefficients, respectively. $p(\mathbf{\Omega }·{\mathbf{\Omega }}^{\prime })$ in ${\mathrm{str}}^{-1}$ is the normalized phase with ${\mathbf{\Omega }}^{\prime }$ and $\mathbf{\Omega }$ denoting the incident and scattering directions, respectively; ${\int }_{4\pi }d{\mathbf{\Omega }}^{\prime }$ represents the integration over the whole solid angle; and $q(\mathit{r},\mathbf{\Omega })$ in $\mathrm{J}{\mathrm{cm}}^{-3}{\mathrm{str}}^{-1}$ is a source function. The anisotropy factor g characterizes the anisotropy of light scattering and is defined as ${\int }_{4\pi }d{\mathbf{\Omega }}^{\prime }p(\mathbf{\Omega }·{\mathbf{\Omega }}^{\prime })cos\mathbf{\Omega }·{\mathbf{\Omega }}^{\prime }$. The steady-state light propagation model based on the RTE has been widely discussed in various research fields, including photoacoustic imaging [24,25] and photocatalyst [26,27,28]. The light propagation model based on the RTE has been validated by light reflectance measurements for tissue phantoms [29].

We consider the refractive-index mismatched boundary condition (RMBC) at the boundary [21,29]. The boundary condition treats the light reflection and refraction at a medium-air interface induced by a difference in the refractive indices between the medium and air as shown in Figure 1b. Fresnel’s law and Snell’s law give the reflectivity and transmissivity coefficients.

2.2. Numerical Schemes and Conditions for the RTE

We numerically calculated the RTE based on the finite difference method. We employed the third order upwind scheme for spatial discretization and discrete ordinate method (DOM) for angular discretization [30,31,32]. As discrete angular directions of light intensity, we used the level symmetric even (LSE) quadrature set [33] with the total number of the directions, $N_{\mathbf{\Omega }}=48$. Several researchers have reported the superiority of the LSE quadrature set [34,35]. The discrete phase function is renormalized by Liu’s method [36] for accurate calculations of the scattering integral (first term in the right-hand side of Equation (1)). Liu’s method works well in the current condition where the anisotropy factor g is small. For details of the upwind scheme and DOM, see references [30,31,32], respectively. Although the references [31,32] focus on a case of highly forward peaked scattering rather than a case of weakly forward peaked scattering, they summarize the DOM with Liu’s method. We numerically solved the discrete forms of the RTE iteratively by Bi-CGSTAB method [37] based on Eigen 3 [38]. We set the tolerance error to ${10}^{-16}$. We numerically treated the RMBC based on Klose’s procedure [29]; we operated the matrix for the RMBC to the results for the light intensity vector after iterations. We preliminarily confirmed that results when operating the RMBC matrix at each iteration were almost the same as the results in the current procedure. However, this procedure requires more computational time. The spatial step size was uniformly given as $\mathrm{\Delta }r$ = 0.025 cm. Our preliminary study showed that the RTE results with a smaller $\mathrm{\Delta }r$-value were almost the same as those with $\mathrm{\Delta }r$ = 0.025 cm. Our in-house source code was written in the C++ programming language, and all the matrices were compressed to vectors in the compressed row storage format. We implemented parallel CPU programming with 128 threads work station (AMD EPYC 7452) by OpenMP. The computational time at each volume fraction was approximately 1 h.

We considered a homogeneous cubic medium with a length of 2.5 cm and an isotropic point source at the origin. In the RTE calculations, we employed the DST results for the scattering coefficient ${\mu }_s$ and normalized phase function $p\left(\theta \right)$ at each volume fraction. For the refractive index of the alumina suspension, we used the following equation: $n_{cs}\left(\eta \right)=\eta n_p+(1-\eta )n_w$ [39], where $\eta$ is the volume fraction; and $n_p$ and $n_w$ represent the refractive indices for alumina particles and water, respectively. We referred values of $n_p$ and $n_w$ as $1.768$ and $1.332$, respectively [40], with the optical wavelength of 600 nm. As shown in Figure 2a, the $n_{cs}$-value increases linearly as the volume fraction increases. According to the change of $n_{cs}$, the reflectivity coefficient for the RMBC changes with the scattering angle as shown in Figure 2b. For the absorption coefficient, we set it to 0.1 ${\mathrm{cm}}^{-1}$ because we focus on the scattering process of light propagation.

2.3. Photon Diffusion Equation (PDE)

The steady-state PDE is obtained by the diffusion approximation to the RTE [41]:

$$\left[-\nabla ·D\left(\mathit{r}\right)\nabla +{\mu }_a\left(\mathit{r}\right)\right]\mathrm{\Phi }\left(\mathit{r}\right)=q_{DE}\left(\mathit{r}\right),$$

(2)

where the fluence rate $\mathrm{\Phi }\left(\mathit{r}\right)$ is defined by ${\int }_{4\pi }d\mathbf{\Omega }I(\mathit{r},\mathbf{\Omega })$; the diffusion coefficient $D\left(\mathit{r}\right)$ is by ${\left[3{\mu }_s^{\prime }\left(\mathit{r}\right)\right]}^{-1}$ with the reduced scattering coefficient ${\mu }_s^{\prime }\left(\mathit{r}\right)=(1-g){\mu }_s\left(\mathit{r}\right)$; and the isotropic source $q_{DE}\left(\mathit{r}\right)$ is by ${\int }_{4\pi }d\mathbf{\Omega }q(\mathit{r},\mathbf{\Omega })$.

We implemented the Robin boundary condition, which is the diffusion approximation of the RMBC [42]:

$$\mathrm{\Phi }\left(\mathit{r}\right)+D\left(\mathit{r}\right)\gamma \left(n\right){\mathbf{e}}_b·\nabla \mathrm{\Phi }\left(\mathit{r}\right)=0,$$

(3)

where $\gamma \left(n\right)$ represents the coefficient for the diffusive reflection [43] and ${\mathit{e}}_b$ represents the outward normal unit vector at the boundary.

2.4. Numerical Schemes and Conditions for the PDE

We numerically calculated the PDE based on the finite difference method. We employed the central difference scheme for spatial discretization. In the PDE calculations, we employed the results of the reduced scattering coefficient ${\mu }_s^{\prime }$ using the DST at each volume fraction. The other numerical schemes and conditions were the same as those for the RTE (in-house C++ code, CPU parallel computing, etc.). The computational time at each volume fraction was approximately 20 min.

2.5. Analytical Solutions for the RTE and PDE

We employed analytical solutions of the steady-state RTE and PDE for a 3D infinite slab to verify our numerical results and investigate the diffusion approximation. We considered the slab with two z-planes ($z=$ 0.0 cm and 2.5 cm with a width $d_{slab}=2.5$ cm) under the Robin boundary condition (Equation (3)). We focused on calculating points inside the medium on the line perpendicular to the z-planes.

The analytical form of the fluence rate for the PDE is given by [44]

$${\mathrm{\Phi }}_{PDE,slab}(z,\mathit{l})=\frac{1}{4\pi D}\sum _{s=-N_s}^{N_s}\left[\frac{1}{r_s^{+}}{e}^{{\mu }_e{r}_s^{+}}-\frac{1}{r_s^{-}}{e}^{{\mu }_e{r}_s^{-}}\right],$$

(4)

where $\mathit{l}=(x,y)$, $N_s$ is the number of the summation with s ($s=0,\pm 1,\pm 2,\cdots ,\pm N_s$); $r_s^{+}$ is a distance between the s-th positive source and calculating position ($z,\mathit{l}$); $r_s^{-}$ is a distance between the s-th negative source and calculating position; and ${\mu }_e=\sqrt{{\mu }_a/D}$. $r_s^{+}$ and $r_s^{-}$ are given by

$$\begin{array}{ccc}\hfill r_s^{+}& =& \sqrt{{(z-z_s^{+})}^2+l^2},r_s^{-}=\sqrt{{(z-z_s^{-})}^2+l^2},\hfill \\ \hfill z_s^{+}& =& 2s(2z_b+d_{slab})+z_0,z_s^{-}=2s(2z_b+d_{slab})-2z_b-z_0,\hfill \end{array}$$

(5)

where $z_b=\gamma \left(n\right)D$ is an extrapolated length and $z_0=\min [\mathsf{\Delta }r,1/{\mu }_s^{\prime }]$ is a small positive valued.

The analytical solution (Equation (4)) is constructed by summing analytical solutions for an infinite medium over the different distances to treat the boundary condition based on the method of images. We considered the analytical form for the RTE in the same way as that for the PDE

$${\mathrm{\Phi }}_{RTE,slab}(z,\mathit{l})=\sum _{s=-N_s}^{N_s}\left[{\mathrm{\Phi }}_{RTE,inf}\left(r_s^{+}\right)-{\mathrm{\Phi }}_{RTE,inf}\left(r_s^{-}\right)\right],$$

(6)

where ${\mathrm{\Phi }}_{RTE,inf}$ is the analytical solution of the RTE for an infinite medium, derived by Liemert and Kienle [45]. Equation (6) treats the boundary condition under the diffusion approximation, while ${\mathrm{\Phi }}_{RTE,inf}$ is the exact solution of the RTE. We modified the open-source code of ${\mathrm{\Phi }}_{RTE,inf}$ for using the phase function calculated from the DST instead of the Henyey–Greenstein (HG) phase function, originally used.

2.6. Dependent Scattering Theory (DST)

The DST can calculate the scattering properties (${\mu }_s$, g, ${\mu }_s^{\prime }$, $p\left(\theta \right)$) for dense colloidal suspensions [10,19] on a length scale comparable to the optical wavelength. On the length scale, we considered an aqueous alumina suspension as a system consisting of multi-sized spherical particles in the background medium (water). The system corresponds to the infinitesimal volume of the scattering medium in Figure 1a. The DST is based on the first-order solution of the Foldy–Lax equation [46,47], which is equivalent to the Maxwell equation [11]. The DST describes well the experimental results of the scattering properties [13,14,17].

We denote a particle diameter class by $\alpha$ or $\beta$ and the total number of the diameter classes by $N_d$. The scattering coefficient ${\mu }_{s,D}$ for the DST is given by

$${\mu }_{s,D}={\int }_0^{\pi }d\theta {\int }_0^{2\pi }d\varphi sin\theta \sum _{\alpha =1}^{N_d}\sum _{\beta =1}^{N_d}\sqrt{n_{\alpha }{n}_{\beta }}{\mathit{F}}_{\alpha }^{Mie}(\theta ,\varphi )·{\mathit{F}}_{\beta }^{Mie*}(\theta ,\varphi )S_{\alpha \beta }\left(\theta \right),$$

(7)

where $\theta$ is the polar angle and set to the scattering angle ${cos}^{-1}(\mathbf{\Omega }·{\mathbf{\Omega }}^{\prime })$; $\varphi$ is the azimuthal angle; $n_{\alpha }$ or $n_{\beta }$ is the number density for the $\alpha$-class or the $\beta$-class ($\{\alpha ,\beta \}=\{1,2,\cdots ,N_d\}$); ${\mathit{F}}_{\alpha }^{Mie}(\theta ,\varphi )$ is the scattering amplitude vector using the Mie theory [48]; ${\mathit{F}}_{\beta }^{Mie*}(\theta ,\varphi )$ is the complex conjugate of ${\mathit{F}}_{\beta }^{Mie}(\theta ,\varphi )$; and $S_{\alpha \beta }\left(\theta \right)$ is the partial static structure factor. See the reference [19] for the formulations of the other scattering properties and the numerical schemes.

2.7. Numerical Conditions of the DST

We considered aqueous alumina suspensions as a size-polydisperse colloidal suspension. We used the particle size distribution of an alumina suspension after sonification measured by Mahmoud et al. [49]. Figure 3 shows the fraction of particle size with the mean diameter $d_{mean}$ of 55 nm and $N_d=10$. We varied volume fractions $\eta [-]$ of the alumina particles from 0.01 to 0.20 with the optical wavelength of 600 nm.

3. Numerical Results

3.1. Scattering Properties

Figure 4 shows the scattering properties (${\mu }_s$, ${\mu }_s^{\prime }$, g) calculated from the DST for the aqueous alumina suspensions at different volume fractions. As shown in Figure 4a, at a lower volume fraction less than approximately 5%, the scattering and reduced scattering coefficients depend on the volume fraction linearly. Meanwhile, at a higher volume fraction, the curvilinear dependence of the coefficients is observed. As shown in Figure 4b, the anisotropy factor is almost zero, although the value decreases as the volume fraction increases. The g-result indicates the light scattering is almost isotropic for the suspensions. This is because the particle diameter (mean diameter $d_{mean}$ of 55 nm) is much smaller than the optical wavelength $\lambda$ of 600 nm; i.e., the size parameter $x_s=2\pi n_{cs}{d}_{mean}/\lambda$ is much smaller than unity. This fact suggests that the Rayleigh theory holds. Because of the g-result, the reduced scattering coefficient values are almost the same as those of the scattering coefficient.

In Figure 5, we plot the normalized phase functions $p\left(\theta \right)$ calculated from the DST with the scattering angle $\theta$. The results using the DST almost agree with those using the Rayleigh theory ($3[1+{cos}^2\theta ]/\left(16\pi \right)$) because of the small $x_s$-value. We also plot the HG phase functions [50] with the same g-values as those using the DST because the function has been widely used. The HG phase functions deviate from the results using the DST, although the deviations are small in the current conditions.

3.2. Fluence Rate

We investigated the fluence rates at different positions on the line parallel to the z-axis with $(x,y)=(0,0)$ inside the medium (blue dash line in Figure 1a). Here, the distances between source and calculating positions (SC distances) range from 0.05 to 2.50 cm at equal intervals of 0.05 cm.

Before discussion of the RTE results, we investigated the PDE results. In Figure 6a,b, we compared the numerical solutions (NSs) of the PDE with the analytical solutions (ASs) of the PDE for the infinite slab (Equation (4)) with the boundaries at $z=0.0$ cm and $2.5$ cm for different volume fractions of 1% and 20%. The results are normalized by the values at the SC distance of 1.25 cm. The numerical results nicely agree with the analytical solutions at all the SC distances, verifying our numerical calculations of the PDE. The agreement suggests that the boundaries except those at $z=0.0$ cm and $2.5$ cm less influence the numerical results on the line inside the medium, because the NSs are for the cubic medium while the ASs for the infinite slab. In Figure 6a,b, we plot the AS for the semi-infinite medium at the $z=0$ boundary, obtained from Equation (4) with $s=0$. We observe differences in the ASs between the infinite slab and semi-infinite medium in the SC distances from 2.3 cm to 2.5 cm due to the boundary effects at $z=2.5$ cm. We plot the relative differences (RDs) between the NS and AS for the infinite slab at each figure’s bottom. The small RDs are observed near the boundaries of $z=0.0$ cm and $2.5$ cm.

We evaluated the differences in the normalized fluence rates between the NS and AS for the infinite slab by the mean absolute percentage difference (MAPD):

$$D_{\mathrm{\Phi }}(N,P;A,P)=\mathrm{Mean}\left[\left|\frac{{\widehat{\mathrm{\Phi }}}_{AS,PDE}^c-{\widehat{\mathrm{\Phi }}}_{NS,PDE}^c}{{\widehat{\mathrm{\Phi }}}_{NS,PDE}^c}\right|\right],$$

(8)

where ${\widehat{\mathrm{\Phi }}}_{AS,PDE}^c$ and ${\widehat{\mathrm{\Phi }}}_{NS,PDE}^c$ represent the AS and NS at the c-th calculating point, respectively. As shown in Figure 6c, the MAPDs are less than 3% at all the volume fractions.

Next, we investigated the diffusion approximation at the boundary by comparing the numerical results for the RTE with the AS (Equation (6)) for the infinite slab at different volume fractions. The AS of the RTE (Equation (6)) is based on the diffusion approximation to the boundary region. Meanwhile, the numerical calculations of the RTE treat the RMBC. As shown in Figure 7a,b, the numerical result nicely agrees with the AS except the boundary region (the SC distances from 2.3 cm to 2.5 cm) at the volume fraction of 1%. The agreements suggest the verification of our numerical calculations of the RTE. The difference at the boundary region ($z=2.5$ cm) in Figure 7a is probably due to the fact that the diffusion approximation is invalid at the volume fraction of 1%. We plot the AS of the RTE for the semi-infinite medium at the $z=0$ boundary obtained from Equation (6) with $s=0$. Similar to the AS of the PDE, we observe the difference between the infinite slab and semi-infinite medium in the SC distance region from 2.3 cm to 2.5 cm. As shown in Figure 7c, the MAPDs, $D_{\mathrm{\Phi }}(N,R;A,R)$, between the NS of the RTE and AS of the RTE for the infinite slab are less than 5% except for the case for the volume fraction of 1%. The small MAPD-values indicate that the diffusion approximation is valid at the boundary at the volume fraction range higher than 1%.

Finally, we investigated the diffusion approximation inside and at the boundary by comparing the numerical results for the RTE with those for the PDE. Figure 8a shows that at the volume fraction of 1%, the numerical results agree with each other except the boundary regions at $z=0.0$ cm and 2.5 cm. The behavers of the RD in the numerical calculations between the RTE and PDE are similar to those between the NS and AS of the RTE (Figure 7a). These results suggest that the diffusion approximation is valid inside the medium even at the lower volume fraction of 1%. As shown in Figure 8b, at the volume fraction of 20%, the numerical results agree well with each other, suggesting the diffusion approximation is valid at all the SC distances. At the $z=2.5$ boundary, the RTE results slightly differ from the PDE results at the two-volume fractions, probably because of the difference in the boundary conditions between the RTE and PDE (RMBC and Robin BC). Figure 8c shows the MAPDs, $D_{\mathrm{\Phi }}(N,R;N,P)$, between the RTE and PDE are less than approximately 5%, meaning the diffusion approximation holds at all the SC distances with the volume fractions higher than 1%. Constant values of the MAPD (∼2%) in the volume fractions from 7% to 20% are probably because the spatial discretization is the same for the RTE and PDE. From Figure 8a,b, the MAPD-values mostly come from the differences at the boundary regions ($z=0.0$ cm and 2.5 cm).

From the results in Figure 7c and Figure 8c, we conclude the diffusion approximation holds at the internal regions with the whole volume fraction range 1–20% and all the SC distances with the volume fraction higher than 1%.

4. Discussion

The DST results showed the ${\mu }_s^{\prime }$-values range from 11 to 70 ${\mathrm{cm}}^{-1}$ for the alumina suspension at the volume fraction range 1–20% in Figure 4. These values are within the range for biological tissue volumes [51]. Meanwhile, the g-values are almost zero, different from a typical value of biological tissue volumes ($g\lesssim 0.8$) [51]. These results suggest that the alumina suspension is an appropriate tissue phantom modeling the reduced scattering coefficient rather than the scattering coefficient and anisotropy factor. Our comparative study of the RTE calculations with the PDE calculations showed the diffusion approximation holds for light propagation, meaning that ${\mu }_s^{\prime }$ primarily characterizes light propagation. This result supports our previous suggestion for the applicability of the alumina suspension.

According to the previous research works of time-dependent light propagation [20,21,22], the characteristic length where the diffusion approximation holds has been evaluated at approximately $10/{\mu }_t^{\prime }$ with ${\mu }_t^{\prime }={\mu }_s^{\prime }+{\mu }_a$. In our current conditions, the characteristic lengths are given as 0.89 cm for the volume fraction of 1% and 0.15 cm for 20%, respectively. As shown in Figure 8a,b, the RTE results agree with the PDE results inside the distance of the characteristic length from the boundaries. Hence, the current results for steady-state light propagation are consistent with those for time-dependent light propagation. Our numerical study suggests that the PDE-DST is useful as a light propagation model for the alumina suspensions rather than the RTE-DST because the diffusion approximation is almost valid and the PDE calculations are faster than the RTE calculations.

The numerical investigation of light propagation in the alumina suspension (liquid phantom) is essential for developing near-infrared spectroscopy and imaging toward their practical application, such as non-invasive cancer detection in the human body. Moreover, the RTE-DST and PDE-DST, developed here in the alumina suspensions, can be applied to other colloidal suspensions and research fields. For example, a titanium dioxide suspension is one of the photocatalytic systems. Because the evaluation of light propagation in a photocatalytic system is significant, the investigation in this paper probably provides valuable information in the research field of photocatalysis.

Based on the following two reasons, we concluded that PDE-DST is beneficial over the RTE-DST in the alumina suspension at the volume fraction range of 1–20%, meaning the diffusion approximation is valid. The other is because the computational loads of the PDE-DST are smaller, and the numerical treatment is simpler. The RTE-DST provides accurate computation of light propagation but requires complicated treatments of the phase function using the DST and RMBC. As shown in Figure 8a,b, the difference between the numerical results mostly comes from the differences at the boundary regions ($z=0.0$ cm and 2.5 cm). Hence, developing an efficient and accurate treatment of RMBC to a boundary condition of the PDE is crucial for future work.

In future work, examining light propagation for the colloidal suspensions with larger particle sizes will be challenging. The scattering properties become larger with an increasing particle size. The large ${\mu }_s$-values indicate the diffusion approximation probably holds at a broader condition, such as the volume fraction compared to the suspension with a smaller particle size. Meanwhile, this result requires a finer spatial step size for the numerical calculation of the RTE-DST. In the suspension with the larger particle size, the g-values become close to unity, meaning highly forward scattering. The g-results prevent the diffusion approximation, where light scattering is isotropic. Moreover, the g-results require accurate and efficient treatments of the highly forward peaked phase function for the numerical calculation of the RTE-DST.