1.

Introduction

When light propagates through a diffusive medium, it experiences absorption and scattering. If the scattering coefficient is considerably larger than the absorption coefficient, the medium can be regarded as diffusive, and light propagation obeys a diffusion process.¹ However, diffusion is meaningful only beyond transport mean free path (TMFP). For media that are considerably shorter than the TMFP, the media are almost ballistic.² Hence, this length scale (TMFP) is extremely important in an optical imaging system and in practice determines the upper limit beyond which no ballistic imaging of the media can be reconstructed. Penetration depths of imaging technologies, which are based on ballistic imaging, such as optical coherence tomography,^3–5 are limited by this length scale. For biological tissues (of almost any kind—skin, brain, liver, and so on) that are thicker than 2 mm, the ballistic component of the light [in the visible or near-infrared (IR) regimes] is negligible.^6,7 In order to investigate a thick medium, most optical imaging methods are using diffusion-based techniques, such as photon density waves and inverse-scattering solutions of the diffusion equation.^8–10

Since the TMFP,^11,12 which is the reciprocal of the reduced scattering coefficient, is the diffusion length scale, it was assumed that the length of the transition from the ballistic regime to the diffusive one should also be equal approximately to this length scale. However, experiments and simulations show disagreement regarding this point.^13–18 The transition length scale varies from one experiment (or simulation) to the next.^15–18 Recent experiments reveal that the transition length depends not only on the scattering coefficient but on the collecting angle as well. Moreover, it was shown both experimentally and theoretically that the transition occurs within a much shorter distance, which is approximately the reciprocal of the scattering coefficient ${\mu }_{\mathrm{s}}$ (instead of the reduced scattering coefficient ${\mu }_{\mathrm{s}}^{\prime }$).^19,20

It should be stressed that the determination and classification of the transition point are not merely an academic issue. In their pioneering works, Yoo and Alfano²¹ have been able to differentiate between the ballistic and diffusive photons with streak camera, and with an ultrafast optical shutter, the same group created a ballistic image of a 3.5-mm thick human tissue.²² This technology was developed in different directions,^15,23–33 and recently, it was demonstrated that by applying the conclusions of Refs. 19 and 20 to this technology, a ballistic image can be reconstructed even when the thickness of the medium is increased substantially.³⁴

The mathematical model that was utilized in Ref. 19 was based on a superposition of the Beer–Lambert term, which describes the ballistic domain, and the diffusion equation solution, which describes the diffusive regime. Since, unlike the ballistic light, the diffusive one is scattered in every direction, then the latter term is multiplied by the collecting angle of the detector. As was later demonstrated experimentally,^19,20 despite its simplicity, this model anticipated the transition from the ballistic regime to the diffusive one with high accuracy.

This result is quite unexpected due to the following reasons: (1) the optimized model works with great accuracy even when the absorption coefficient is relatively large and (2) the model is based only on two terms (ballistic and diffusive) and ignores other types of transport, such as quasiballistic ones.^35,36

The simplicity of this model is very appealing since the generic solution requires solving the radiative transfer equation,^17,37,38 which is too complicated for an analytical solution. The equation can be solved numerically with a Monte Carlo (MC) simulation; however, it is time-consuming. Even faster techniques, such as the discrete ordinates and the adding–doubling methods are still computationally intensive and are, therefore, used primarily for layered, or quasi-one-dimensional (1-D), media (see, e.g., Refs. 39 and 40).

Moreover, the simple analytical solution, namely Ref. 19, teaches about the main parameters that affect the transition point, which, as was shown in Ref. 34, has a practical implication in ballistic imaging, namely, imaging an object hidden in 1 cm of chicken breast.

It is the object of this paper to investigate the validity of this model for different absorption values by investigating the transition from the ballistic regime to the diffusive one with much higher accuracy. To achieve the high accuracy, an MC simulation^41–51 was utilized in an MATLAB computational platform. This simulation allows, except for comparison to experiment, a very flexible method to investigate the transition with the required accuracy. It should be emphasized that since the MATLAB program was written in a parallel form, the simulation time was reduced by more than two orders of magnitude.

The results reveal that the model does agree with the simulations and is valid for higher absorption coefficient with a simple power-law prefactor.

2.

System and its Model

The model of the experiment is shown in Fig. 1. Since this paper is not an experimental research and the original experiment is described in details in Refs. 19 and 20, only the main points will be listed below. The light source was a laser with a wavelength 840 nm, the medium was an intralipid solution with a scattering coefficient of approximately ${\mu }_{\mathrm{s}}\cong 140{\mathrm{cm}}^{-1}$, and varying values of absorption coefficients ${\mu }_{\mathrm{a}}\cong 0.04$, 2.4, 4.8, and $9.6{\mathrm{cm}}^{-1}$, which were constructed by varying the India ink concentration in the solution (as can be seen from the experimental results in the ballistic regime, the additional ink has a negligible effect on the scattering coefficient). The light source was placed near the sample, so the beam remained parallel at the entrance to the medium. Two pupils were placed on both sides of the medium to keep the beam at the same width of the detector. The pupils’ diameter was smaller than the width of the sample to keep the quasi-1-D approximation valid. In the experiment, the setup was oriented vertically (and not horizontally as in Fig. 1), so the width of the sample ($L$) can be varied by changing the solution level in the glass (see Ref. 20 for details).

Fig. 1

System schematic. $L$ is the slab thickness and $L_0$ is the distance between the slab and the detector.

Clearly, if the medium’s width $L$ is short enough, then most scattered photons would be blocked by the barrier and mainly ballistic photons would be detected. However, if the medium is wider than the transition length, the number of ballistic photons is reduced dramatically, and as a consequence, most of the detected particles experience scattering.

3.

Theoretical Model

The original model, which was developed in Ref. 19, is based on the premises that the detected radiation consists of only two terms: ballistic and diffusive. The ballistic term is governed by the Beer–Lambert’s law,^14–18,37 i.e., the intensity decreases exponentially

The diffusion term is a degenerated 1-D solution of the three-dimensional diffusion equation⁸ (see Ref. 52 for a solution of the diffusion equation for a slab geometry).

Since the light source is continuous, the stationary diffusion equation for the photon density $\rho (\mathbf{r})$ can be utilized⁵³

When the beam’s cross section is larger than the medium’s thickness, there is degeneracy in the transversal coordinate and the beam decays approximately exponentially in the propagation direction ($z$). Since in the diffusion approximation the medium is isotropic, the local density is proportional to the local intensity, i.e.,

However, while the ballistic photons, which have survived the medium, suffer no additional losses, the diffusive ones are scattered in all directions and only a fraction eventually reach the detector. This process is shown in Fig. 2.

Fig. 2

Diffused photons detection.

Equation (3) is a good approximation of the diffusion equation at the end of the medium [inside the medium an additional exponentially increasing term should be introduced, i.e., $I_{\text{diffusive}}(z)/I_{\text{diffusive}}(0)\cong \exp (-{\mu }_{\mathrm{eff}}z)+C\exp (+{\mu }_{\mathrm{eff}}z)$; however, since this term is multiplied by an exponentially small coefficient, i.e., $C\sim \exp (-2{\mu }_{\mathrm{eff}}L)$, eventually, at the end of the medium, i.e., at $z\cong L$, they both have the same exponentially decaying dependence on $z$ (for elaboration, see, e.g., Ref. 54)]. The effects of the boundary can be neglected in the model, since the attenuation due to the boundary’s reflectivity is below the experimental accuracy.

If the collecting angle of the detector is $\delta \mathrm{\Omega }$, then the diffusive term should be multiplied by the factor $\delta \mathrm{\Omega }/4\pi$, which in our system can be approximated by $(\delta \mathrm{\Omega }/4\pi )\cong (d^2/4\pi L_0^2)$, where $d^2$ is the rectangular cross section of the detector (the fact that the prefactor is proportional to $L_0^{-2}$ is consistent with the result of Ref. 55).

Hence, since the model consists of both the ballistic and diffusive terms, it can be approximated by

4.

Monte Carlo Simulation

The simulation algorithm was similar to Ref. 18, and its main features are presented here for completeness purposes.

It should be stressed, however, that within the terminology of Refs. 5657.–58, the simulation included contributions from all types of particles, i.e., all orders of scattering. There is no distinction between multiply scattered and low-order scattered photons.

To keep the Beer–Lambert’s law, the probability density of the distances between scatterings ($s$) is

The probability density of the cosine of the elevation angle $\theta \in [0,\pi ]$ obeys the Henyey–Greenstein^38,59 phase function

Fig. 3

Photon propagation in a diffusive slab.

Eventually, each trajectory that hits the detector surface is multiplied by $\exp (-{\mu }_{\mathrm{a}}\sum _{i=1}^N{s}_i)$, where $N$ is the number of scattering events, to account for the absorption.

By multiplying the arriving photons by their absorption attenuation, the number of simulating photons can be reduced substantially, since each particle in the simulation simulates a group of particles. It was found that 20 million such groups of particles are sufficient to simulate with good accuracy (higher than the experimental one) a range of 10 orders of magnitude (eight of which were measured).

5.

Comparison Between Experiments, Simulations, and the Mathematical Model

In Fig. 4, the MC simulation was compared to the experimental results of Ref. 19 and to the theoretical model [Eq. (4)].

Fig. 4

Comparison of the MATLAB (blue circle), experiment (red crosses), and theoretical results (magenta dashed line). The scattering coefficient of the experiment, simulations, and the theoretical model is $140{\mathrm{cm}}^{-1}$ and the absorption coefficient is $9.6{\mathrm{cm}}^{-1}$. The anisotropic coefficient is 0.9. The beam, the barrier, and the detector size are $0.15\mathrm{cm}\times 0.15\mathrm{cm}$ and the distance from the slab to the detector is 30 cm. The $x$-axis is the slab thickness and the $y$-axis is the relative intensity decay in logarithmic scale.

In Fig. 4, there is a good agreement between the simulation and the experimental results. However, there is a disagreement with the theoretical model. The exponential decay agrees with the theoretical model both in the ballistic regime and the diffusive one; however, there is a disagreement on the transition point, and as a consequence, there is a downward shift in the diffusive part of the plot. This discrepancy was well known in Ref. 20 but its source was unknown. In Sec. 6, we will quantify this discrepancy.

6.

Introducing a Correction Prefactor to the Diffusive Term

Since the discrepancy occurs at the transition, it can be corrected by an additional prefactor ($a$) on the diffusive term (since there is no change in the slope of the graph). To quantify this term, an MC simulation was carried on the same geometrical system ($d=0.15\mathrm{cm}$, $L_0=30\mathrm{cm}$, and $g=0.9$) but with different scattering and absorption coefficients. In every scenario (specific ${\mu }_{\mathrm{s}}$ and ${\mu }_{\mathrm{a}}$), the simulation intensity data were used to calculate the prefactor $a$ that fits the following equation:

Table 1

The a prefactor for different scenarios.

In Fig. 5, the data are plotted on a single graph. The data can be fitted to a power-law

Fig. 5

The prefactor $a$ versus the absorption coefficient in a log–log plot, the scattering coefficient of the blue squares is $100{\mathrm{cm}}^{-1}$, of the green diamonds is $120{\mathrm{cm}}^{-1}$, and of the red circles is $140{\mathrm{cm}}^{-1}$. The solid line stands for Eq. (12).

These results suggest an interesting conclusion that the analytical expression can be utilized even when the validity of the diffusion approximation is questionable. The prefactor $a$ allows using the model for absorption coefficient (at least) as large as ${\mu }_{\mathrm{a}}\cong {\mu }_{\mathrm{s}}/10$.

Therefore, Eq. (3) can be refined to

Figures 6 and 7 present the effect of the additional prefactor $a$. In Fig. 6, the prefactor is absent (or equal to 1), and the discrepancy is evident. On the other hand, in Fig. 7, the prefactor is present, i.e., the data are compared with Eq. (13), and the agreement is excellent.

Fig. 6

Simulation and theoretical results for a scattering coefficient of $140{\mathrm{cm}}^{-1}$. The absorption coefficient of the blue circles (simulation) and blue line (theoretical) is $0.09{\mathrm{cm}}^{-1}$, of the red X’s (simulation) and the red line (theoretical) is $0.8{\mathrm{cm}}^{-1}$, of the green triangles (simulation) and the green line (theoretical) is $2{\mathrm{cm}}^{-1}$, and of the pink squares (simulation) and the pink line (theoretical) is $5{\mathrm{cm}}^{-1}$. The anisotropic coefficient is 0.9. The beam, the barrier, and the detector size are $1\mathrm{cm}\times 1\mathrm{cm}$ and the distance from the slab to the detector is 30 cm.

Fig. 7

The same simulation (as Fig. 6) results with the prefactor $a$.

In Figs. 8 and 9, the relative error is plotted versus the sample’s width with and without the prefactor. The relative error is defined $E\equiv 2|I_{\mathrm{S}}-I_{\mathrm{M}}|/|I_{\mathrm{S}}+I_{\mathrm{M}}|$, where $I_{\mathrm{S}}$ and $I_{\mathrm{M}}$ are the intensities of the simulation and the model [Eq. (13)], respectively. Despite the great improvement of the prefactor, there is still some error at the transition area, which may indicate the need for a third transport term (see, e.g., Ref. 61). This discrepancy is relatively large when there is a sharp transition between the ballistic and the diffusive regimes, i.e., when the absorption is low. However, for many applications, the model’s accuracy is surprisingly good.

Fig. 8

The relative error $E$ between the simulation data and the model without the prefactor for the coefficients of Fig. 6.

Fig. 9

The relative error $E$ between the simulation data and the model with the prefactor for the coefficients of Fig. 7.

7.

Summary

An MC simulation was conducted to investigate the model presented in Refs. 19 and 20. The main purpose was to validate the premises of Ref. 19 that the transition from the ballistic to the diffusive regimes can be described with high accuracy by only two mathematical terms: ballistic and diffusive, namely $I(z)=I(0)[\exp (-{\mu }_{\mathrm{t}}z)+(\delta \mathrm{\Omega }/4\pi )\exp (-{\mu }_{\mathrm{eff}}z)]$. The simulation, as well as the experimental results, indicates that this is indeed a good model. Moreover, it was found that this model is valid even for relatively large absorption coefficients in which the validity of the diffusive approximation is dubious. To fix the model for larger absorption coefficients, a power-law prefactor should be added to the second term, namely $a={({\mu }_{\mathrm{a}}/{\mu }_{\mathrm{a}}^0)}^{1/3}$, where ${\mu }_{\mathrm{a}}^0=0.007{\mathrm{cm}}^{-1}$. This prefactor teaches that the original model is valid provided the scattering coefficient is 20,000 times larger than the absorption coefficient. It should be stressed that the conclusions of the research are limited to the regime of the simulation, namely, to the propagation of visible and near-IR light in biological tissues. A more generic claim requires an extension of the research.