## Abstract

In Part I of this study [1], good agreement between experimental measurements and results from Monte Carlo simulations were obtained for the spatial intensity distribution of a laser beam propagating within a turbid environment. In this second part, the validated Monte Carlo model is used to investigate spatial and temporal effects from distinct scattering orders on image formation. The contribution of ballistic photons and the first twelve scattering orders are analyzed individually by filtering the appropriate data from simulation results. Side-scattering and forward-scattering detection geometries are investigated and compared. We demonstrate that the distribution of positions for the final scattering events is independent of particle concentration when considering a given scattering order in forward detection. From this observation, it follows that the normalized intensity distribution of each order, in both space and time, is independent of the number density of particles. As a result, the amount of transmitted information is constant for a given scattering order and is directly related to the phase function in association with the detection acceptance angle. Finally, a contrast analysis is performed in order to quantify this information at the image plane.

© 2009 Optical Society of America

## 1. Introduction

A comparison between experiment and simulation has been presented in Part I [1] for the optical radiation scattered within various turbid media. This first article demonstrated good agreement between the measurement and simulation results, effectively validating the fitness of the Monte Carlo (MC) code. Use of this model allows accurate predictions of laser radiation propagation and scattering within the intermediate single-to-multiple scattering regime, where no analytic solution can be accurately applied. The previous work also calculated probability densities for scattering orders contributing to the total forward and side-scatter intensities, which allowed extension of the Beer-Lambert relation to include a multiple-scattering term. These results motivate the current investigation, which examines the properties of individual scattering orders resolved in both space and time, and their effects on the total collected light intensity in forward and side-detection geometries.

Note that the term “*scattering order*” employed here corresponds to the number of times an individual photon packet interacts with the scattering particles prior to exiting the turbid medium. This is not to be confused with the geometrical optics term “*order of refraction*,” or “*scattering order*,” which refers to the number of times a light ray is partly reflected or refracted within a single scattering particle [2–3].

Detailed spatial distribution analysis is essential for quantifying the utility of each scattering order with regard to the process of image formation. Furthermore, the temporal signature of each scattering order is of importance for the development of modern ultra-fast laser based diagnostics [4–5]. Experimentally the unscattered light (ballistic photons) can be reliably isolated in some situations since it is not subject to the changes in polarization, phase, or direction of propagation introduced by the scattering process. A number of techniques based on spatial filtering [6], time-gating [4, 7–8], polarization gating [9–10] and coherence-based diagnostics [11] have been developed in order to discriminate ballistic photons from scattered light within scattering and/or absorbing media [12]. Nevertheless, the amount of ballistic photons can easily fall below detection limits for highly turbid media. Singly-scattered photons can be extracted temporally by detecting the first backscattered photons (using a LIDAR arrangement [5, 13–14]), and spatially by means of a crossed source/detector geometry with an “infinitely” small detection acceptance angle, as suggested in [15]. However, despite the promise of these techniques, no experimental configuration is currently able to provide a quantitative analysis of contributions from each distinct scattering order. Fortunately, this information is directly available through numerical simulation. This work presents temporal and spatial analysis of individual scattering orders, and addresses how distinct orders contribute to image formation for both the forward and side scattering.

## 2. Monte Carlo simulation for scattering orders analysis

Detailed treatments of the statistical Monte Carlo (MC) approach as well as implementations specifically related to light scattering are available in the literature [15–27]. In MC simulation, the pertinent details of light propagation, including the number of scattering events, state of polarization, coherence properties of the propagating radiation and changes of photon direction at each scattering event can be tracked individually for each photon as it traverses the scattering volume. Using this approach, the information pertaining to individual scattering orders is directly accessible since the total scattered signal is calculated by summing the detected photons. With this arrangement, scattering order information can be filtered from simulation data. These advantages make the MC method an effective tool for evaluating the contribution of individual scattering orders in image formation, where the total intensity detected equals:

In the case of low scattering, the series of scattering order intensities can be reduced to the *I*
_{(1)} term which represents the single scattering. The higher order terms of the series may be neglected or considered as perturbations to the single scattering term. In the case of multiple scattering the intensity terms have the same order of magnitude and the diffusion approximation is commonly used. Finally, in the intermediate single-to-multiple scattering regime, both single and higher order terms are considered individually and no approximation or assumption can be made [1, 15, 19–20].

The contribution from each scattering order differs depending on whether back, side, forward or the overall scattering process is considered. In parallel, a number of parameters influence the relative contribution of individual scattering orders on the total detected signal. These parameters include the optical depth (*OD* which equals the extinction cross-section *σ _{e}*, times the particle number density

*N*, times the length of the scattering volume

*l*), the detection acceptance angle (

*θ*) and the angular dependence of the scattering phase function (

_{a}*f*(

*θ*)). It has been recently demonstrated that the transition from the single to the multiple-scattering regime does not operate in the same manner for isotropic scattering as for anisotropic scattering when assuming the overall scattering process [19]. For isotropic scattering the single scattering order remains dominant until all orders become equal at high optical depths. Conversely, when considering scattering phase functions with a highly forward scattering lobe (anisotropy factor

_{s}*g*=0.82), the dominant scattering order changes from one to the next as

*OD*increases. In this case, it has been observed that the dominant scattering order and the value of

*OD*were nearly the same [19].

The work presented in this article examines the propagation of laser radiation in the intermediate scattering regime by analyzing the temporal and spatial contribution of individual scattering orders. The code employed for the simulation is identical to that used in Part I, excepting that the results are filtered at the end of the simulation for the first twelve scattering orders. The number of photons launched in the simulations is twice the amount utilized in [1], corresponding to ~6 billion photons, in order to increase statistics while detecting individual scattering orders. Similar to Part I, the scattering media employed here are cubic scattering volumes of 10 mm side length containing homogeneous solutions of monodisperse polystyrene spheres. Two types of microspheres, corresponding to *D*=1 and 5 µm, are considered and illuminated at *λ*=800 nm, resulting in particle size parameters of *x*=3.9 and 19.6 respectively. As the absorption of 800 nm light within the polystyrene sphere solution is negligible, real indices of refraction are assumed, such that *n*=1.578-0.0*i* for the scatterers and *n*=1.33-0.0i for the surrounding medium. Three optical depths, *OD*=2, 5 and 10, are used and photons are collected from both the forward and side-scattering directions with a detection acceptance angle of *θ _{a}*=8.5°. The dimensions of the scattering volume are constant and the changes in

*OD*are only due to an increase of the scatterer’s number density. The time-of-flight for each photon is directly calculated from the total photon path-length between particles considering the velocity of light in water. Due to the small size of the microspheres, the duration of photon-particle interaction is assumed to be negligible [25–27]. All photons enter the scattering medium at the same instant, resulting in a temporal profile described by a Dirac delta function. At the end of the simulation, the photon time-of-arrival curves are convolved with a realistic temporal profile for the incident pulse. These results assume a Gaussian input pulse shape with a full width at half maximum of 100 fs. Finally, the spatial profile of the laser source is modeled from the experimental image matrix recorded when the laser beam crosses a sample containing pure distilled water, similar to [1].

The results from these simulations are examined in three ways. First, the data are filtered to show the relative contribution of scattering order intensities in image formation. This approach helps to understand how experimental images, such as the ones presented in [1], are formed from complex multiple scattering processes. Simulated photons are separated for each distinct order and results of the final times-of-flight and positions in the detection plane are provided. This temporal and spatial study is performed for both the forward and side scattering detection. The second approach records the position of the final scattering events just prior to detection at the front face. The scattering events which occur within the “plane” *P* located in the center of the probed volume are considered as illustrated in Fig. 2. The last approach presents analysis of image contrast transmitted from the first four scattering orders. Here, the incident light source is spatially modulated with a bar pattern with a frequency of 2 line pairs/mm, and the case of isotropic scattering is also included to serve as a “base case” reference.

## 3. Results for the forward scattering detection

#### 3.1 Spatial analysis

Successive images of the spatial light intensity distribution representing each scattering order, from n=0 to n=12, are shown in the movies linked to Fig. 3(a) and Fig. 3(b) for the 1 and 5 *µ*m polystyrene spheres, respectively. Light is detected in the forward scattering direction and the normalized intensity profile along the vertical axis at X=5 mm is shown together with the relative value of the full width at half maximum (FWHM) of the exiting beam profile. For both particle sizes, these movies show that the FWHM increases with the value of the scattering order. This increase is more pronounced for the 1 *µ*m than for the 5 *µ*m particles. The intensity ratio, *I*
_{(n)max}/*I*
_{(incident)max}, shown in the two movies is plotted for *OD*=2, 5, and 10 respectively in Figs. 4(a), 4(b), and 4(c). These results provide a quantitative comparison of the light transmission by individual scattering orders between the two particle sizes. The shape of the curves differ slightly from the probability density per scattering order *P*
_{(n)}, presented in Part I, due to the inhomogeneous intensity distribution of the incident pulse.

At n=0, the transmission is described by the Beer-Lambert law. The effects of the large detection acceptance angle and of the scattering phase function do not influence, in this case, the relative detected intensity. However, for any n ≥ 1, it is apparent that *I*
_{(n)max}/*I*
_{(incident)max} is always larger for the 5 *µ*m than for the 1 *µ*m scattering particles. This observation demonstrates that the contributions from low and high scattering orders do not depend exclusively on *OD*, but are also closely related to the scattering phase function as well as the acceptance angle used for detection. In other words, at any given optical depth the amount of multiply scattered photons detected in the forward direction tends to be higher for larger particle sizes and detection acceptance angles (also observed in [25]). These results confirm that a modified Beer-Lambert law [1] which takes into account the intensity from the multiple light scattering is required when measuring the laser light transmission within turbid media, particularly in situations where large particles are present.

In parallel with these quantitative results, Fig. 4(d) presents the broadening of the FWHM of the normalized intensity profile as a function of scattering order. The initial laser beam intensity profile measured along the vertical axis Z, at X=5 mm equals *da*=2.55 mm. The ratio *d/da* represents the increase of the laser FWHM from its initial value *da*. At scattering order n=0, *d/da* equals 1, as the ballistics photons provide a precise representation of the incident laser beam. It is possible for this value, calculated from simulation results, to appear smaller than one, at *OD*=10, for n=0, due to the low statistics of the ballistic photons detected under these conditions. Figure 4(d) summarizes the data shown in the movies linked to Fig. 3(a) and (b), and illustrates how the laser beam is spatially broadened with increasing scattering order. A number of interesting trends can be extracted from these results. First, changes in the FWHM seem to be linear with scattering order. Second, the increase in the FWHM is mainly driven by the scattering process itself, particularly the shape of the scattering phase function. The FWHM is apparently not dependent on the *OD* of the medium or, more precisely, on the number density of scatterers. Though spatial broadening of the laser beam appears to diverge slightly at n≥6, this divergence is the result of decreasing statistics. Indeed, these results show that the spatial broadening of the laser beam remains essentially the same for all three optical depths at scattering orders n≤5. Note also that the initial inhomogeneous oval shape of the laser beam defined at n=0, tends to a smooth circular intensity profile when increasing the scattering order as a result of an increase of the FWHM in both X and Z direction. Similar effects are observed for the total light intensity, when increasing *OD*.

The most interesting observation from these initial results is that the normalized intensity distribution for a given scattering order and particle size seems to be independent of the value of the optical depth. In other words, given the dimension of the simulated volume is constant, variations in particle concentration have no apparent effect on the normalized profile of the image transmitted by a given scattering order.

#### 3.2 Temporal analysis

The transmission of an ultra-short pulse through highly scattering media has been widely investigated experimentally, especially during the late 20^{th} century [27–33]. Recently, Calba *et al*. [27] have compared numerical MC calculations to experimental results for a geometry similar to the one presented here, but at larger *OD* and for a smaller detection acceptance angle. However, the time dependence from each order has apparently not been investigated in detail in the literature [4, 19].

In this section, the successive photon times-of-arrival representing each scattering order, from n=0 to n=5, are shown in Fig. 5, together with the 10^{th} order and the total signal. These graphs provide a comparison between the 1 and 5 *µ*m polystyrene spheres for the forward scattering direction. Each curve is normalized with the maximum value of the total detected signal represented by the solid black line. The earliest photon time-of-arrival equals 44.4 ps, corresponding to the time necessary for the photons to cross 10 mm of distilled water. Results for the low optical depth, *OD*=2, are presented in Figs. 5(a) and 5(b). For the case of 1 *µ*m particles, shown in (a), the total transmitted pulse consists of a sharp peak of light intensity followed by a thin tail. This peak corresponds mostly to ballistic photons while the tail is built up from the successive contributions of the low scattering orders (1^{st}, 2^{nd} and 3^{rd} orders). For the case of 5 *µ*m particles, shown in (b), the peak of the light intensity is wider than for the 1 *µ*m particles. This is due to the temporal contribution of the low scattering orders at early arrival times which overlap with the ballistic pulse. This superposition of optical signals can be explained by the shorter time-spreading of individual scattering orders induced by the larger forward scattering lobe of the phase function for the 5 *µ*m particles at the near 0° scattering direction. At the low optical depth *OD*=2, the light intensity from scattering orders higher than 3 becomes negligible for both 1 and 5 *µ*m particles.

For intermediate optical depth, *OD*=5, presented in Figs. 5(c) and 5(d), it is observed that the total temporal profile no longer reflects the shape of the incident laser pulse. Ballistic photons continue to dominate the signal for the 1 *µ*m spheres, in contrast to the 5 *µ*m case, where the 1^{st}, 2^{nd}, 3^{rd} and 4^{th} scattering orders yield a stronger maximum peak contribution than the unscattered photons. In both cases, the contribution of the 10^{th} scattering order remains negligible compared to the total detected light intensity. The thin tail observed in (a) at *OD*=2 is now wider and longer in (c) for *OD*=5, with the appearance of the 4^{th} and 5^{th} orders.

At high optical depth, *OD*=10, the ballistic pulse remains visible for the 1 *µ*m particles, as shown in Fig. 5(e), even though the broadening of the total transmitted pulse is significant. For the 5 *µ*m particles, the peak of the total transmitted pulse shifts to later times as *OD* increases. In contrast, the peak of the 1 *µ*m curve remains at the same location for *OD*=2 and *OD*=5, corresponding to the maximum value of the ballistic light intensity.

#### 3.3 Localization of last scattering events and normalization of the spatial and time-resolved intensities within distinct scattering orders

As mentioned at the end of subsection 3.1, the normalized intensity distributions appear to be independent of the number density of particles for a given scattering order. Furthermore, it can be seen from Fig. 5 that the intensity profile peaks for individual scattering orders remain at the same location independent of *OD*. These observations suggest that the position of the last scattering events for a given order should also remain constant when varying particle concentration. This section aims to verify this hypothesis. In total, the distribution of the final scattering events prior to detection gives a spatial representation of the light source responsible for the recorded intensity distribution. One should note that different 3-dimensional distributions of final scattering events can form an identical normalized intensity pattern on the front face. However, two identical sources must, in principle, provide the same intensity pattern as well as the same normalized temporal profile.

We examine here, the location of the final scattering events within the central “plane” *P*, prior to detection on the front face (as indicated in Fig. 2). Figure 6 shows the normalized number density of final scattering events occurring within the central plane *P*, for scattering order n=4. The three rows in Fig. 6 represents *OD*’s, 2, 5, and 10, while each column displays simulation data at a given detection acceptance angle and particle size. It is clear from these results that the spatial distribution of final scattering events is identical for all three optical depths. This observation applies for both of the particle sizes (1 and 5 *µ*m,) and both of the detection acceptance angles considered (8.5° and 90°). It is important to note that this phenomenon is unique to the forward-detection geometry, and holds for all scattering orders.

The movie linked to Fig. 7(a) presents the series of final scattering events for scattering orders n=1 to 5. Results for 1 *µ*m, 5 *µ*m, and isotropic scattering cases are presented. This movie illustrates the displacement of the scattering source from one order to the next. As expected, this displacement tends to the forward direction with a side broadening effect for increasing n. At higher scattering order, the locations of the point sources tend to be closer to the exit face. However, for the isotropic scattering case, this forward displacement is not as apparent as it is for the more forward scattering phase functions.

From these last considerations, it follows that the optical paths traveled by the photons to form each “scattering order source” are also independent of particle concentration. As a result, for the forward-detection in a plane perpendicular to the incident light, each scattering order defines normalized time-resolved and spatial intensity distributions which are constant regardless the number density of particles. This observed phenomena can be analytically described such that:

where the total transmitted light intensity, *I*
_{(TOT)}, for a given number density, *N*, is made up (both in space and time) by the sum of the normalized contributions of each scattering order (independent of *N*) weighted by the appropriate probability density *P*
_{(n)} (see [1]). Figure 8 shows these normalized light intensity profiles for the odd scattering orders from n=1 to n=9. One can see that the intensity peak is shifted in time as the scattering order increases. The magnitude of this shift is directly related to the scattering phase function which follows from the particle size as shown in Fig. 7(b).

## 4. Results for the side scattering detection

#### 4.1 Spatial analysis

Successive images of the spatial light intensity distribution representing each scattering order, from n=1 to n=12, are shown in the movies linked to Fig. 9(a) and Fig. 9(b). These results correspond to simulations with 1 and 5 *µ*m polystyrene spheres, respectively. Light is detected in the side-scattering direction and the normalized intensity profile along the horizontal axis at Z=5 mm is shown.

Similar to section 3.1, these movies provide a unique view of how photons from a given scattering order spatially contribute to image formation. As expected, the movies show that the laser source penetrates the medium more effectively at lower *OD* and high scattering orders appear at greater depths within the medium than lower scattering orders. Furthermore, in contrast to the results from forward-scattering detection, the total light intensity levels remain of the same order of magnitude for the two particle size cases. As discussed earlier in section 3.1, the forward-propagating light tends to an even circular profile as scattering order increases. For the side scattering, the input light intensity is also spatially averaged or spread by scattering events in the medium. Smaller particles which exhibit more isotropic phase functions produce a stronger spreading effect, leading to a final smoothed spatial intensity distribution. Larger particles which exhibit phase functions with a narrower forward scattering produce a weaker spreading effect, but they will tend toward a similar smoothed spatial intensity if *OD* is increased to very large values. The collected intensity from side scattering is directly related to the location of the scattering events within the medium. This imparts a spatial signature to the intensity levels detected in each scattering order. In effect, each scattering order, n, is marked by the spatial signature of the scattering events from the preceding order, n-1. Consequently, effects induced by the shape of the scattering phase function become more apparent when larger scattering orders are considered. Single scattering, n=1, is unique in that photons scatter from the source directly, without the influence of spatial signatures from previous scattering. Thus, the single scattering provides a rough measure of the source penetration and exhibits an exponential decline related to *OD* due to laser light extinction. The single scattering intensity is here a direct measure of the amount of light scattered at ~90° from the scattering phase function. It is deduced from the movies linked to Figs. 9(a) and 9(b) that the maximum single scattering intensity detected on the side within 8.5° detection acceptance angle is 1.7 times larger for the 5 than for the 1 *µ*m spheres.

Figure 10 examines the data from the movies linked in Fig. 9 to show the contributions from individual scattering orders. The cross-section of the light profile along the horizontal axis at the position Z=5 mm is presented and the results are compared for 1 and 5 *µ*m particles. It is seen here that the maximum light intensity of each individual scattering order is located at various distances along the horizontal axis Y and do not superimposed upon one another. In other words, each scattering order presents a preferred area where it tends to contribute the most. The distance from “peak to peak” between successive scattering orders diminishes with increasing optical depth. It is also observed that these peaks are smoothed as n increases.

At *OD*=2, shown in 10(a) and (b), the low scattering orders contribute the most to the total detected signal. Orders higher than n=5 are of negligible importance. The single-scattering profile indicates that the source of laser light fully penetrates the medium.

At *OD*=5, shown in 10(c) and (d), the single-scattering profile indicates that source penetrates only two thirds of the scattering volume. Here, scattering orders higher than n=10 make negligible contributions to the detected signal.

At *OD*=10, shown in 10(e) and (f), all of the scattering orders included in the simulation contribute significantly to the total signal while the single-scattering profile corresponds to the lowest intensity contribution. In fact, the detection of the singly scattered light occurs in the first quarter of the cell length. After this distance, the laser beam is highly attenuated and presents a negligible contribution when compared to the multiply scattered light intensity. It is observed that the value of the highest scattering order and of the optical depth is nearly the same. More precisely, at *OD*=2, 5 and 10 the dominant scattering orders corresponds to n=2, 5 and 7 respectively. In addition, each intensity profile succeeds the next in a logical manner from one scattering order to another. This differs from the forward scattering results where the amount of ballistic light, detected at given acceptance angle, can influence the early logical succession of the series (as observed in Fig.4 (b) for the case of 1 *µ*m at *OD*=5).

It is worth mentioning that the penetration of the laser source differs from the penetration for individual scattering orders, which exhibit a dependence on the properties of the phase function. While the penetration of the source, at a given *OD*, remains constant for different particle sizes, one can observe that the maximum peak, for individual scattering orders, is always observed to be further for the 1 *µ*m than for the 5 *µ*m. This is somewhat counter-intuitive, since 5 *µ*m particles exhibit greater forward scattering, which results in better forward penetration. Nonetheless, the wide-angle forward scatter of the 1 *µ*m scattering phase function (especially the large lobe at *θ _{s}* ~64.6° [1]) results in greater penetration along the side face of the simulated volume. Finally it is important to note once more that the normalized spatial profiles depend on the number density of particles, for the side scattering, contrary to the observations made in forward scattering.

#### 4.2 Temporal analysis

The time-resolved light intensity for side scattering, contributed by each scattering order, from n=0 to n=5, is shown in Fig. 11, together with the 10^{th} order and the total signal. These graphs provide a comparison between the 1 and 5 *µ*m spheres at *OD* 2, 5 and 10. Each curve is normalized with the maximum value of the total detected signal (black solid line). From Fig. 11, the total side-scatter profile appears effectively unchanged for different *OD*’s and particle sizes. These results differ substantially from the situation observed for forward-detection, where the effects of scattering phase function and particle concentration strongly influence the shape and spreading of the total signal. For side-detection, the width of the profile is largely determined by the dimensions of the scattering medium and influenced to a lesser extent by the acceptance angle of the collection optics. The temporal and spatial profiles exhibit some similarities, suggesting a strong relation between them. It is also observed that by increasing *OD* the maximum of the curves move to a shorter time, for both particle cases. These shifts of a few picoseconds are not negligible when compared to the time scale of the forward scattering pulses. Furthermore, even if the total temporal profile can appear to be independent of *OD* and particle size, the information contained in the total signal is, in fact, strongly influenced by these parameters. This becomes apparent when one examines the relative contributions from each individual scattering order.

For *OD*=2, shown in Figs. 11(a) and 11(b), the low scattering orders contribute significantly to the total intensity. As *OD* increases to 5, as seen in Figs. 11(c) and 11(d), this contribution of low orders decreases while at the same time the contribution of the 10^{th} order begins to appear. Finally, at *OD*=10, the low scattering orders no longer contribute significantly to the total signal, which is dominated by the 10^{th} scattering order. This trend is observed for both particle sizes. In conclusion, for side-scattering detection individual orders are not well temporally separated - regardless of the particle size or concentration - and are broad in time, following the trend of the total intensity profile. This results in overlapping photon arrival times between scattering orders and suggests that a time gating approach would not extract efficiently the single light scattering intensity with such source/detector geometry.

## 5. Contrast transfer for individual scattering orders

Up to this point, the results and analysis in this article have centered on the spatially resolved intensity profiles and time propagation of individual scattering orders in media with various optical densities. This yields insight into the impact of different scattering orders on the total detected intensity and the effects of *OD*, detection parameters, and scattering phase function on each scattering order. This section now takes the same analysis a step further by introducing a spatial modulation to the input light. This makes it possible to estimate the imaging quality of the photons by calculating a contrast value from the intensity contribution. Note that, for the total transmitted signal, the effect of a turbid medium on image formation has been widely investigated in the past. Such studies have been carried out by means of computational models which determine the Modulation Transfer Function as well as for the estimation of the Point-Spread Function for a given scattering and imaging system [23, 24]. The originality of the section presented here concerns the study of contrast transfer for individual scattering orders.

In this set of simulations the input light is modulated by a square wave in the horizontal direction with a spatial frequency of 2 lp/mm. For this specific spatial frequency, the contrast transmitted to the front face of the simulated volume is evaluated for photons of different scattering orders and for the total signal detected. The contrast C of the intensity field is calculated as:

where *I*
_{max} is the maximum image intensity and *I*
_{min} is the minimum image intensity, which occurs at the centre of the bars and within the central part of the detected image. An illustration of the simulation is given in Fig. 2. The case of isotropic scattering has also been investigated for comparison purposes with the scattering process exhibited by the 1 and 5 µm particles. Figure 12 shows a movie of the resultant modulated images on the front face, for each scattering order such as 0≤n≤4.

The results of the contrast per scattering order are presented in the histogram given in Fig. 13(a) assuming a detection acceptance angle *θ _{a}*=8.5°. Scattering order zero is not plotted in this graph as that case indicates full contrast (100%) since the ballistic photons provide a faithful reproduction of the collimated incident light source being detected on the front face.

The first scattering order produces reduced contrast value, indicating that single scattering photons can also contribute to some amount of image noise for the forward-detection. It is observed that this trend increases exponentially with subsequent scattering orders. This behavior is quite intuitive since the probability that a photon is misdirected increases directly with the number of scattering events. In addition, these results clearly show that the induced scattering noise increases faster for more isotropic scattering conditions; where all scattering orders n>2 have a null contrast. For highly forward scattering phase functions, however, such as for the 5 *µ*m particles, a non-zero contrast can be still be observed even at fairly large scattering order such as n=4. This tendency to keep information after successive scattering events and to retain a positive contribution to the image is expected to increase for larger particle size parameters.

Figure 13(b) shows the value of the image contrast of the total detected intensity as a function of *OD* for the 1 *µ*m, 5 *µ*m and the isotropic scattering. Similar to Eq. (2), this final contrast, for a given number density *N*, can be calculated from the sum of the contrast per scattering order *C*
_{(n)} (given in Fig. 13 (a)) weighted by the normalized probability density *P*
_{(n)} such that:

The results of Eq. (4) are shown in dashed line in 13(b) for n≤4 only, as each scattering process *C*
_{(n)} tends here to zero at values of n>4. The best total contrast is provided for the isotropic scattering even at high *OD*. This effect is due to the greater relative concentration of ballistic light in the total detected signal due to the detection acceptance angle employed. In effect, the combination of the scattering process and a finite detection acceptance angle acts as a filter which quickly discards light that would otherwise contribute to image noise induced by multiple scattering. By means of this directional filtering, which increases as *θ _{a}* is reduced, the contrast values can remain relatively high at large

*OD*, even if the signal intensity becomes very low as a result of the same process.

However, it is also important to note that reducing the detection acceptance angle in an imaging system leads to a reduction of image resolution. Hence, in practice, it is often undesirable to decrease the detection acceptance angle. The phase functions that describe the scattering from 1 and 5 *µ*m particles both show increased forward scattering, compared to the isotropic case; however, the 5 *µ*m scatterers exhibit more strictly forward scattering [1]. This provides an interesting comparison for the total contrast:

At *OD*=2, the effect of increased forward scattering at larger angles allows the 1 *µ*m case to retain a larger relative number of lower scattering orders, yielding a greater contrast value than the 5 *µ*m case; at *OD*=5, the filtering effects of the more isotropic 1 *µ*m particles and the increased contribution of higher order photons in the 5 *µ*m case lead to similar values of contrast; at *OD*=10, the 1 *µ*m signal is overwhelmed by scattering noise as the low order photons are reduced while the contrast of the 5 *µ*m case can still be observed (~2.2%). At such OD, the scattering from the 1 *µ*m particles is too far forward to be filtered effectively by *θ _{a}* and does not forward enough to retain any positive contribution to the image at large n, resulting in a contrast that tends to zero.

## 6. Conclusion

The analysis presented here considered light propagation and image formation in turbid media for both forward and side-scattering detection by examining individual scattering order data generated from a previously validated MC model [1]. Spatial and time-resolved intensity profiles were produced for the first twelve scattering orders. Examination of the side-scattering signal reveals that although the maxima of the spatial and time-resolved profiles for each order show some degree of separation, individual scattering order profiles are broad in both space and time, resulting in significant overlap for all orders. The temporal profile of the overall signal appears similar for changes in particle concentration; however, analysis of the contribution from individual orders demonstrates that the photons collected at different *OD*’s differ significantly. Examination of the forward-scattering signal reveals that the spatial and temporal profiles of the total detected intensity are influenced by particle size and the acceptance angle of the collection optics, but the normalized intensity distribution of each order is independent of the number density of particles in both space and time. To verify this phenomenon, the distributions of final scattering events just prior to detection were examined for the first five scattering orders. These results confirm that this distribution remains constant when varying the number density of particles. Consequently, the amount of transmitted information is constant for a given scattering order regardless of particle concentration. The quality of this information was quantified by an analysis of image contrast. It has been shown that the contrast per scattering order is greater for larger scatterers which exhibit a more forward scattering phase function. Importantly, this means that image information can be extracted from multiply-scattered light, with varying degrees of fidelity for different scattering media. The analysis presented in this work, Part II in the series on light scattering, provides a unique view of the fundamental mechanisms underlying photon transport within turbid media. To our knowledge, the fact that the forward-detection profiles for each order are independent of scatterer concentration has not been discussed in the literature. This result, together with the comparative side-scatter investigation, final scattering event spatial distributions, and contrast analysis for individual scattering orders, provides important insights into the effects of multiple scattering in the intermediate scattering regime.

## Acknowledgments

This work is partly performed within the Strategic Research Centre for studies of combustion process, CECOST. The authors acknowledge the Swedish Foundation for Strategic Research (contract A3 05:183) and the European Union Large Scale Facility program (project llc001131) for their financial support.

## References and links

**1. **E. Berrocal, D. L. Sedarsky, M. E. Paciaroni, I. V. Meglinski, and M. A. Linne, “Laser light scattering in turbid media Part I: Experimental and simulated results for the spatial intensity distribution,” Opt. Express **15**, 10649–10665 (2007) [CrossRef] [PubMed]

**2. **H. C. van de Hulst, *Light scattering by small particles* (Dover, N.Y., 1981)

**3. **C. Bohren and D. Huffman, *Absorption and scattering of light by small particles* (Wiley, N.Y., 1983)

**4. **M. A. Linne, M. Paciaroni, E. Berrocal, and D. Sedarsky, “Ballistic imaging of liquid breakup processes in dense sprays,” Proc. Combust. Inst. **32**, 2147–2161 (2009) [CrossRef]

**5. **B. Kaldvee, A. Ehn, J. Bood, and M. Aldén, “Development of a picosecond lidar system for large-scale combustion diagnostics,” Appl. Opt. **48**, B65–B72 (2009) [CrossRef] [PubMed]

**6. **G. E. Anderson, F. Liu, and R. R. Alfano, “Microscope imaging through highly scattering media,” Opt. Lett. **19**, 981 (1994) [CrossRef] [PubMed]

**7. **L. Wang, X. Liang, P. Galland, P. P. Ho, and R. R. Alfano, “True scattering coefficients of turbid matter measured by early-time gating,” Opt. Lett. **20**, 913–915 (1995) [CrossRef] [PubMed]

**8. **J. C. Hebden, R. A. Kruger, and K. S. Wong, “Time resolved imaging through a highly scattering medium,” Appl. Opt. **30**, 788- (1991) [CrossRef] [PubMed]

**9. **O. Emile, F. Bretenaker, and A. Le Floch, “Rotating polarization imaging in turbid media,” Opt. Lett. **21**, 1706–1708 (1996) [CrossRef] [PubMed]

**10. **V. Sankaran, K. Schnenberger, J. T. Walsh, and D. J. Maitland, “Polarization Discrimination of Coherently Propagating Light in Turbid Media,” Appl. Opt. **38**, 4252–4261 (1999) [CrossRef]

**11. **K. A. Stetson, “Holographic fog penetration,” J. Opt. Soc. Am. **57**, 1060–1061 (1967) [CrossRef]

**12. **C. Dunsby and P. M. W. French, “Techniques for Depth-Resolved Imaging through Turbid Media including Coherence-gated Imaging,” J. Phys. D: Appl. Phys. **36**R207–R227 (2003) [CrossRef]

**13. **R. M. Measures, *Laser Remote Sensing: Fundamentals and applications* (Krieger, Florida, 1992)

**14. **M. Gai, M. Gurioli, P. Bruscaglioni, A. Ismaelli, and G. Zaccanti, “Laboratory simulations of lidar returns from clouds,” Appl. Opt. **35**, 5435–5442 (1996) [CrossRef] [PubMed]

**15. **E. Berrocal, D. Y. Churmakov, V. P. Romanov, M. C. Jermy, and I. V. Meglinski, “Crossed source/detector geometry for novel spray diagnostic: Monte Carlo and analytical results”, Appl. Opt. **44**, 2519–2529 (2005) [CrossRef] [PubMed]

**16. **I. M. Sobol, *The Monte Carlo Method* (University of Chicago, Chicago, Ill., 1974)

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

**18. **G. Zaccanti, “Monte Carlo study of light propagation in optically thick media: point source case,” Appl. Opt. **30**, 2031–2041 (1991) [CrossRef] [PubMed]

**19. **E. Berrocal, *Multiple scattering of light in optical diagnostics of dense sprays and other complex turbid media* (PhD Thesis, Cranfield University, 2006)

**20. **V. P. Romanov, D. Yu. Churmakov, E. Berrocal, and I. V. Meglinski, “Low-order light scattering in multiple scattering disperse media,” Opt. Spectros. **97**, 847–854 (2004)

**21. **I. Meglinski, M. Kirillin, V. Kuzmin, and R. Myllylä, “Simulation of polarization-sensitive optical coherence tomography images by a Monte Carlo method,” Opt. Lett. **33**, 1581–1583 (2008) [CrossRef] [PubMed]

**22. **L. R. Poole, D. D. Venable, and J. W. Campbell, “Semianalytic Monte Carlo radiative transfer model for oceanographic lidar systems,” Appl. Opt. **20**, 3653–3656 (1981) [CrossRef] [PubMed]

**23. **P. Bruscaglioni, P. Donelli, A. Ismaelli, and G. Zaccanti, “Monte Carlo calculations of the modulation transfer function of an optical system operating in a turbid medium,” Appl. Opt. **32**, 2813–2824 (1993) [CrossRef] [PubMed]

**24. **X. Gan and M. Gu, “Effective point-spread function for fast image modeling and processing in microscopic imaging through turbid media,” Opt. Lett. **24**, 741–743 (1999) [CrossRef]

**25. **C. Rozé, T. Girasole, L. Méès, G. Gréhan, L. Hespel, and A. Delfour, “Interaction between ultra short pulses and a dense scattering medium by Monte Carlo simulation: consideration of particle size effect,” Opt. Commun. **220**, 237–245, (2003) [CrossRef]

**26. **C. Calba, C. Rozé, T. Girasole, and L. Méès, “Monte Carlo simulation of the interaction between an ultra-short pulse and a strongly scattering medium: The case of large particles,” Opt. Commun. **265**, 373–382, (2006) [CrossRef]

**27. **C. Calba, L. Méès, C. Rozé, and T. Girasole, “Ultrashort pulse propagation through a strongly scattering medium: simulation and experiments,” J. Opt. Soc. Am. A **25**, 1541–1550 (2008) [CrossRef]

**28. **Y. Kuga, A. Ishimaru, and A. P. Bruckner, “Experiments on picosecond pulse propagation in a diffuse medium,” J. Opt. Soc. Am. **73**, 1812–1815 (1983) [CrossRef]

**29. **K. M. Yoo and R. R. Alfano, “Time-resolved coherent and incoherent components of forward light scattering in random media,” Opt. Lett. **15**, 320- (1990) [CrossRef] [PubMed]

**30. **G. Zaccanti, P. Bruscaglioni, A. Ismaelli, L. Carraresi, M. Gurioli, and Q. Wei, “Transmission of a pulsed thin light beam through thick turbid media: experimental results,” Appl. Opt. **31**, 2141–2147 (1992) [CrossRef] [PubMed]

**31. **Feng Liu, K. M. Yoo, and R. R. Alfano, “Transmitted photon intensity through biological tissues within various time windows,” Opt. Lett. **19**, 740–742 (1994) [CrossRef] [PubMed]

**32. **S. G. Demos and R. R. Alfano, “Temporal gating in highly scattering media by the degree of optical polarization,” Opt. Lett. **21**, 161–163 (1996) [CrossRef] [PubMed]

**33. **V. M. Podgaetsky, S. A. Tereshchenko, A. V. Smirnov, and N. S. Vorob’ev, “Bimodal temporal distribution of photons in ultrashort laser pulse passed through a turbid medium,” Opt. Commun. **180**, 217–223 (2000) [CrossRef]