1. Introduction

Detailed analysis of light scattered by a medium is at the core of several measurement techniques in a variety of areas ranging from atmospheric remote sensing to medical imaging. In many practical situations, the probed medium can be represented as a homogeneous horizontally infinite slab. In this case the main properties that describe the medium are its geometrical thickness $H$, its extinction and absorption coefficients $\sigma _e$ (corresponding to a mean free path $l_e$) and $\sigma _a$, and its scattering phase function $P$ (often reduced to the asymmetry parameter $g$). Generally, reflectance or transmittance is measured, under various illumination and viewing conditions, at a variety of wavelengths and spectral resolutions. In scattering media, photons follow tortuous paths before being eventually absorbed, transmitted or reflected. The diversity of possible paths can be represented by the probability distribution (hereafter called “pdf”) of total path length $l$. Investigation of the path length distribution has been a critical issue for decades [1].

The pdf of photon path length for reflected and transmitted light can be directly measured using photon time-of-flight spectroscopy [2,3] but this requires knowledge of the effective refractive index of the medium [4] and a sophisticated experimental setup that can hardly be deployed. The pdf can also be determined theoretically for certain geometries, considering successive orders of scattering [5], or by using the fact that spectral reflectance and transmittance can be written as the Laplace transform of the path length distribution [6,7], a relationship known as equivalence theorem [8]. The latter method can be applied to analytically derived expressions of reflectance and transmittance [9], or, at least in theory, to observations. However the spectral resolution and noise level required to obtain a satisfactory inversion remain very demanding for real measurements [10]. Although the full pdf contains a wealth of information on the medium [11], in most cases only some moments are accessible, for instance exploiting the Laplace transform results in the moments being related to the derivatives of the measured quantities with respect to $\sigma _a$ [12]. Hence it is worth investigating how the moments of the pdf are related to the medium characteristics. To do so, one can either assume a particular pdf and fit free parameters in order to match observations [13], or use measurements that are directly related to these moments. For instance, the first moment, namely the mean path length, is what determines overall absorption in the case of weak absorption [14].

Theoretical work has demonstrated that the mean path length of photons escaping from a non-absorbing slab illuminated by a diffuse source equals $2H$, independently of the scattering properties of the medium [15], a property that has been later on verified experimentally [16]. Similar invariance properties have been obtained for higher moments [17], but only in the case of optically thick media. More generally, the diffusion approximation has been widely used to investigate the pdf and its moments [18,19]. This approximation has proven very efficient for optically thick media but is inappropriate for optically thin media, that is, in the ballistic limit. Optically thin media are found everywhere, though, from thin cirrus clouds to human tissues, and deserve more attention.

To model the response of optically thin media to a light beam, various strategies have been explored [20]. Simple models where light can only be scattered perpendicular or parallel to the slab have been used [21,22]. It was shown that the mean path length tends towards $3H$ for normal illumination and low optical thickness, while the ballistic limit would suggest it equals $H$. This deviation was attributed to the very few photons scattered parallel to the slab, which experience such long paths that this overcompensates their scarceness. In such configuration, it was also found that directly transmitted, reflected and transmitted photons equally contribute to the mean path length. Normal illumination was also considered in [21], but with isotropic scattering in the medium. The pdf of the directly transmitted and singly-scattered (reflected or transmitted) photons was derived, to obtain a lower bound of the successive moments. It was shown that the first moment converges to the ballistic limit $H$, but that the second moment converges to something larger than $H^2$, and that the $n^{\mathrm {th}}$ moments ($n>2$) divided by $H^n$ even diverge with decreasing optical thickness $\tau$. In a complementary study [23], this divergence was further explored and it was fully explained by the long tail of the pdf for photons experiencing one or two scattering events. The latter study provides an insightful discussion on the reasons why photons scattered once and twice should equally contribute to the moments divergence.

In the present study, we extend the work of [23] to the general cases of collimated (from any direction) and diffuse illumination. The case of diffuse illumination is of particular interest because of the invariance properties that arise in such conditions. The invariance property for the mean path length has been proven for any optical thickness. In contrast, the dependence of the second moment (the variance) on optical thickness has only been studied in the case of optically thick media. For instance, its quadratic dependence has been derived theoretically in the diffusion limit [19] and observed for transmitted sunlight through opaque stratus clouds [24]. The question naturally arises how the variance behaves in the case of optically thin media. That is a central motivation for the present study, which also advances our understanding of light scattering in optically thin media.

In Section 2, the pdf for singly- and twice-scattered photons are obtained for collimated and diffuse illumination, and for directional and angularly integrated viewing configurations. Focusing on the long tails of these pdf the asymptotic behavior of the moments at low optical thickness are investigated. In Section 3, Monte Carlo ray-tracing simulations are performed to validate the analytical results, and to extend the conclusions drawn from the calculations to higher orders of scattering and anisotropic scattering (addressed theoretically in Appendix B). Some practical implications of the present study are discussed in Section 4, and we draw our conclusions in Section 5.

2. Analytical computations

Let us consider a horizontally infinite, homogeneous, purely scattering (i.e., non-absorbing) slab of thickness $H$ illuminated by a beam of photons incident at a zenith angle $\theta _0$ such that $\cos \theta _0=\mu _0$ ($s_0 = \mu _0^{-1}$). The medium is characterized by its extinction coefficient $\sigma _e$, the corresponding extinction length being denoted $l_e$. The optical thickness of the slab is defined as $\tau = H/l_e$. Let us further assume that the distribution of the photon free path $L$ between two scattering events follows a simple exponential law such that:

 

Fig. 1. Geometry of the slab, with the main quantities used in the calculations, for single (top) and double (bottom) scattering.

Download Full Size | PDF

2.1 Collimated illumination

The case of normal illumination ($\theta _0=0$) has already been treated [21,23]. Here we extend the computation to any incident angle $\theta _0$. In this case $p_0(l)$, which corresponds to directly transmitted photons, simply reads:

2.1.1 Singly scattered transmitted photons

The distribution $p_{1,t}(l)$ (equal to $0$ for $l\leq H$) can be written as:

2.1.2 Singly scattered reflected photons

The distribution $p_{1,r}(l)$ can be written as $p_{1,t}(l)$ in Eq. (5). Like previously the paths $l\geq H\left (1+s_0\right )$ (which can be obtained for any $z_0$), and those such that $l\leq H \left (1+s_0\right )$ (which can be obtained only if $z_0\leq z_\mathrm {max}$) are distinguished. Here $z_\mathrm {max}$ is such that $l=z_\mathrm {max}(1+s_0)$, i.e., $z_\mathrm {max}=\frac {l}{1+s_0}$. Defining the escape angle $\theta _r$, $\cos \theta _{r} = \frac {z_0}{l-l_0}$, and $\sin \theta _r \mathrm {d}\theta _r~=~\frac {z_0}{\left (l-z_0 s_0\right )^2}\mathrm {d}l$. $p(l|z_0)$ now reads:

 

Fig. 2. Probability distributions of total path length for singly scattered transmitted (top) and reflected (bottom) photons under collimated illumination. The black lines correspond to the analytical formulas of Eqs. (8), (9), (12), (13). The gray shades correspond to histograms obtained from Monte Carlo simulations with $10^7$ photons launched.

Download Full Size | PDF

2.2 Diffuse illumination

We now consider isotropic illumination, such that the pdf of the incident direction of the photons is given by $P(\theta _0,\phi ) = \frac {\cos \theta _0}{\pi }$. As a consequence $P(\theta _0) = 2\cos \theta _0\sin \theta _0$. In this case the pdf of $l$ for directly transmitted photons (equal to $0$ for $l\leq H$) is such that:

2.2.1 Singly scattered transmitted photons

Following Eq. (5):

2.2.2 Singly scattered reflected photons

We proceed similarly for the reflected paths to obtain the following expression:

 

Fig. 3. Same as Fig. 2 but for diffuse illumination. The black lines correspond to the analytical formulas of Eqs. (19), (23), (24).

Download Full Size | PDF

2.3 Directional observations with collimated illumination

The pdfs derived in the previous sections were obtained considering all photons escaping the slab, whatever their escape direction. Although conceptually useful, these pdfs are not applicable to standard satellite or ground observations which generally measure scattered photons at a particular viewing angle, and under collimated illumination (from the Sun or from an active source such as a lidar). Using the same approach as previously the joint probability distribution of $l$ and escape angle ($\theta _r$ or $\theta _t$), is derived for both transmitted and reflected photons. As before, this pdf is decomposed in successive orders of scatterings.

2.3.1 Reflected photons

For single scattering, the pdf simply corresponds to the probability to be scattered at the appropriate depth $z_0$, and then to escape without being scattered again:

Singly scattered photons do not correspond to very long trajectories, hence to further investigate long paths it is useful to derive $p_{2,r}(l)$. We first notice that in the case of double-scattering without $l_1=(z_0-z_1)s_1$, and $l_2=z_1 s_r$, where $s_1=\mu _1^{-1}$, and $\mu _1 = \cos \theta _1$. Also $\mathrm {d}l = \mathrm {d}z_1(s_r-s_1)$.

2.3.2 Transmitted photons

The same strategy is applied to estimate $p_{1,t}(l,\theta _r)$ and $p_{2,t}(l,\theta _r)$. This reads (for $s_t \neq s_0$):

For double scattering the strategy is the same as for reflected photons, except that the integration bounds are different. Noticing that $z_1=\frac {H s_t+z_0 s_1-l+l_0}{s_1+s_t}$, the conditions $0 \leq z_1 \leq H$ become:

 

Fig. 4. Same as Fig. 2 but for directional viewing. The black lines correspond to the analytical formulas of Eqs. (32), (36). The histograms are obtained from Monte Carlo simulations with $10^8$ photons launched.

Download Full Size | PDF

2.3.3 Specific case of nadir illumination and viewing

Although we did not derive the general expressions of $p_{2,t}(l,\theta _r)$ and $p_{2,r}(l,\theta _r)$ for low values of $l$, we provide below the expression for nadir illumination and viewing, which corresponds to backscattering measurements, for instance corresponding to monostatic lidar observations. In this case, for $l\leq 2H$:

2.4 Estimation of the moments

In this section the moments of the previously derived pdfs are estimated. The $n^{\mathrm {th}}$ moment is denoted $M^{(n)}_{i,x}$, where the lowerscript refers to the pdf of interest. The moments for collimated normal illumination ($\mu _0=1$) have already been investigated in [21,23]. The authors concluded that in the case of an optically thin slab, $M^{(n)}$ differs from the ballistic limit $H^n$ as soon as $n\geq 2$, and the normalized moments $M^{(n)}/H^n$ even diverge as $\tau ^{n-2}$ for $n\geq 3$. They attributed this to the few photons being scattered nearly parallel to the slab boundaries that would experience very long paths. These photons, making up the long tail of the pdf, would not affect the mean path but would drive the divergence for higher moments. In addition, they provided physical insight to support the fact that the moments of $p(l)$ would be dominated in such conditions by the contributions of $p_0(l)$, $p_1(l)$, and $p_2(l)$, with the contributions of photons scattered once and twice being of the same order of magnitude. Their theoretical results were supported by Monte Carlo ray tracing simulations. Here we investigate the moments of $p_0(l)$ and $p_1(l)$ in the case of collimated illumination for any incident direction, and for diffuse illumination. We also investigate the moments of $p_2(l)$ for the case of directional observations. The main focus is on the behavior of the normalized moments at low optical thickness, a normalized moment corresponding to the moment divided by the ballistic moment, i.e., $(H/\mu _0)^n$ for direct illumination and $(2H)^n$ [15] for diffuse illumination.

2.4.1 Collimated illumination

For low optical thickness the contribution of the unscattered transmitted photons to $M^{(n)}$ simply becomes $(H/\mu _0)^n$. For singly-scattered photons, using only the expressions of Eqs. (8) and (12) valid for $l \geq 2H$, the pdfs reduce to (for $\mu _0 l \gg H$):

2.4.2 Diffuse illumination

Starting from Eq. (15), $M^{(0)}_{0,t} = 2\tau ^2 \Gamma (-2,\tau )$, where we’ve used the standard notation for the incomplete $\Gamma$ function. Likewise,

2.4.3 Realistic configurations

Measuring all photons escaping from a scattering slab under collimated or diffuse illumination is not a common measurement configuration. More often, the illumination is direct, and only reflected or transmitted radiation is measured (with full angular integration or at a specific viewing angle). In such case, the relevant pdfs are those corresponding to the actually measured photons. It means that the pdfs presented in Sections 2.1, 2.2, 2.3 should be normalized by the probability of actually being reflected (respectively transmitted), in all directions or in the direction $\theta _r$ (respectively $\theta _t$). We do not derive these latter probabilities that aggregate all orders of scattering, but hypothesize that at low optical thickness, they are dominated by single scattering (excluding direct transmission measurements). As previously, the subsequent derivations are obtained in the limit $\tau \rightarrow 0$.

Directional reflection. We define the photon path pdf for reflected photons only:

 

Fig. 5. Second and third normalized moments ($\mu ^{(2)}=\mu _0^2M^{(2)}/H^2$ and $\mu ^{(3)}=\mu _0^3M^{(3)}/H^3$) of $p(l)$ for collimated illumination ($\mu _0=0.5$) and istropic scattering, as derived analytically (lines) and obtained with ray tracing simulations (symbols). Note that the moments for transmitted photons are not shown, since they would superimpose with those of reflected photons. Likewise, the second moments for twice-scattered photons are not shown because they would superimpose with those of singly-scattered photons. MC stands for Monte Carlo simulations.

Download Full Size | PDF

We now focus on $p_{2,r}(l)$. It can be shown that for $l\gg H(s_0+s_r)$,

Directional transmission. Proceeding similarly for photons transmitted in direction $\theta _t$:

Note that the independence of $M^{(1)}_{2,t}$ to $\mu _0$ and $\mu _t$ appears in the Fig. 7 of [1].

Reflection and transmission in all directions. Building on the previous derivations it is worth considering reflectance and transmittance measurements that account for photons reflected in all directions, under collimated or diffuse illumination. It implies a scaling of the pdfs derived in Sections 2.4.1 and 2.4.2 by the probability to actually be reflected or transmitted, again assumed to be dominated by single-scattering without. In this case, the scaling of $M^{(2)}$ changes. Bearing in mind that the asymptotic expressions are similar for reflected and transmitted photons it reads for collimated illumination:

For diffuse illumination:

3. Monte-Carlo ray-tracing simulations

Now that we have derived analytical expressions for the photons scattered once and twice, we compare these expressions to the results obtained with a Monte Carlo ray tracing code. This code, which is similar in essence to [25], consists in launching a number of photons at the surface of a numerical slab, and following them until they escape the slab via transmission or reflection. The distance between successive scattering events is drawn from the pdf given by Eq. (1) and the scattering angle is drawn from a Henyey–Greenstein phase function [26] described by the asymmetry parameter $g$. For diffuse illumination the incident angle is drawn from the angular distribution $2\cos \theta \sin \theta$. The code was optimized to work at low optical thickness. To this end, each photon contributes to scattering (the probability to be unscattered is taken into account for each photon), following the strategy of [21]. In addition, for diffuse illumination, more photons are sent with grazing angles to better explore the rare but long trajectories that control the behavior of $M^{(2)}$ in optically thin slabs. In practice the code is run on a supercomputer, which allows to run $30$ simulations in parallel, each using $128 \times 10^{8}$ photons ($10^{9}$ for directional observations). This amount of photons is needed to capture the divergence of the normalized moments caused by very rare events.

In Figs. 2 and 3, the shaded vertical bars show the histograms of the total paths traveled by photons in the Monte Carlo simulations. This validates the analytical formulas (Eqs. (8), (9), (12), (13)) obtained for collimated and diffuse illumination. Monte Carlo simulations are also used to obtain the pdfs for directional observations (Fig. 4), which again supports the theoretical derivations for sufficiently long paths (Eqs. (32) and (36)).

The ray tracing code is further used to investigate the behavior of the moments at low optical thickness. To this end, the moments are computed for $30$ distinct values of $\tau$, ranging from $10^{-6}$ to $1$. These simulations are first used to validate the estimations of the moments derived in Section 2. Figures 5 and 6 show the second and third moments obtained for a slab with $H=1$ m, for collimated and diffuse illumination. These simulations confirm the formulas derived analytically, in particular Eqs. (39), (40), (42), (45), (46). The same slab is used to simulate directional observations. Figure 7 shows the results for reflectance observations but the results for transmittance observations are not shown since they are redundant with those for reflectance. Again, these simulations confirm the analytical computations for singly and twice-scattered photons, namely Eqs. (50), (54), (55), (57), (58). For angularly integrated reflectance and transmittance observations, the simulations were used to confirm the assumption that singly-scattered photons dominate the $0^\mathrm {th}$ moment of $p(l)$, which validates the expressions of the moments derived for these configurations.

 

Fig. 6. Same as Fig. 5, but for diffuse illumination ($\mu ^{(2)}=M^{(2)}/4H^2$ and $\mu ^{(3)}=M^{(3)}/8H^3$). Note that the third moment for twice-scattered photons is not shown because it would superimpose with that of singly-scattered photons.

Download Full Size | PDF

Because our analytical derivations do not extend beyond singly or twice-scattered photons, the ray tracing code is used to further explore the $\tau$ dependencies for higher orders of scattering. Simulations are also performed for anisotropic scattering, i.e., $g \neq 0$ in the Henyey–Greenstein phase function. These complementary simulations are further discussed in Section 4.

4. Applications and extensions

4.1 Applications to optical medium measurements

The pdfs derived for singly- and twice-scattered photons (reflected and transmitted) suggest that such measurements (which imply time-resolved measurements) could be used to estimate $l_e$ for optically thin slabs, when low orders of scattering prevail. For instance, Eq. (13) does not depend on $H$, hence $p_{1,r}(0)$ directly gives access to $l_e$, meaning that time resolved reflectance measurements could directly provide $l_e$. In combination with an estimation of $H$ (for instance estimated from the discontinuity of the derivative of $p_{1,r}(l)$ at $l = H(1+s_0)$, or from the shape of the pdf beyond $H(1+s_0)$) this could allow to estimate the optical thickness of the slab from a single reflectance measurement. Equivalently, transmittance measurements could be used. In the case the full pdf is not measured, but its second moment is accessible (as in e.g. [10]), the obtained dependence on $\tau$ suggests that this fundamental quantity can be estimated even at very low optical thickness. More interestingly, the lower the optical thickness, the larger the variance of the photon path length distribution, suggesting an increased sensitivity to the thinnest media, contrary to what is usually experienced in remote sensing. Although the present study is mostly theoretical, it provides leads for innovative remote sensing studies of optically thin media.

4.2 Higher orders of scattering

For non directional simulations, our analytical calculations did not extend to double scattering. When accounting for double-scattering without under collimated illumination it was shown [23] that at low optical thickness $M^{(2)}_{2,t} \sim H^2/4$ (their Fig. 3), which matched very well their Monte Carlo simulations. Assuming that the pdfs of reflected and transmitted photons have similar tails for a given order of scattering, and assuming that higher-than-two orders of scattering are not significant, this suggests that for collimated illumination:

The same strategy is adopted for diffuse illumination. Contrary to collimated illumination, Monte Carlo simulations suggest that the contribution of twice-scattered photons to the second moment is half that of singly-scattered photons (i.e., $M^{(2)}_{1,r} = 2 M^{(2)}_{2,r}$). A possible explanation is that singly-scattered photons with long trajectories include those entering the slab at grazing angles (which can then be scattered in any direction) and those being scattered into the plane. For twice-scattered photons, grazing incident angles do not contribute because the probability to remain within the slab after scattering is very low. On the contrary, for photons already scattered into the plane, the probability to be scattered again before escaping is significant, without changing the overall distance. As already pointed out it seems that both contributions (grazing incident angles and scattering into the plane) are as likely, explaining the factor $1/2$ for twice-scattered photons. As for collimated illumination we can hypothesize that:

For directional observations, the calculations already encompass double scattering. Here we argue that higher orders of scattering are too rare to dominate the asymptotic behavior. Given that singly-scattered photons necessary correspond to short trajectories, it suggests that the overall behavior of the moments is driven by twice-scattered photons. This hypothesis is fully supported by the simulation results shown in Fig. 7.

 

Fig. 7. Same as Fig. 5, but for directional reflectance observations ($\mu _r=0.5$) under collimated illumination ($\mu _0=0.8$). Here $\mu ^{(2)}=M^{(2)}/H^2$ and $\mu ^{(3)}=M^{(3)}/H^3$.

Download Full Size | PDF

4.3 Anisotropic scattering

So far, all the calculations presented in this study assumed isotropic scattering. This is an ideal case, poorly representative of actual scattering media. Monte Carlo simulations using an Henyey–Greenstein phase function with $g\neq 0$ suggest that the asymptotic behavior of the moments is unchanged for collimated illumination, except for an offset corresponding to a scaling of the $\tau$ dependency. Figure 8 shows the moments computed from the simulations for $g=0$, $g=0.5$ and $g=0.8$ (for collimated illumination with $\mu _0=0.8$). In practice, this scaling corresponds to the probability to be scattered into the horizontal plane when coming from the direction $\theta _0$, namely $P(\theta _0,\pi /2)$. For diffuse illumination, the ray tracing simulations with anisotropic scattering show that the asymptotic behavior of the moments is exactly the same as for isotropic scattering. In particular the phase function has no impact on $M^{(2)}$ at low optical thickness, because the asymptotic behavior is driven by directly transmitted photons. It implies that measurements performed in these conditions would directly provide information on the optical thickness, without knowing the phase function, which is an important invariance result. In Appendix B, detailed calculations are performed to support these computational results. When considering directional reflectance or transmittance observations, the scaling includes the contribution of single scattering ($P(\theta _0,\theta _r)$), which determines the $0^{\mathrm {th}}$ moment. The angularly integrated reflectance and transmittance observations are also altered by the impact of the phase function on the $0^{\mathrm {th}}$ moment.

 

Fig. 8. Same as Fig. 5, but for $\mu _0=0.8$ and for 3 different Henyey–Greenstein phase functions ($g=0$, $g=0.5$, $g=0.8$). The lines correspond to a scaling of the paths by $P(\mu _0,\pi /2)$, where this quantity for Henyey–Greenstein phase function is computed after Eq. (76).

Download Full Size | PDF

5 Conclusions

Based on comprehensive radiative transfer calculations, we have derived the path length probability distribution functions of singly- and twice-scattered photons, reflected by or transmitted through a homogeneous, non-absorbing, isotropically scattering slab. Various illumination and viewing conditions were investigated, including collimated and diffuse illuminations, and directional reflectance and transmittance configurations. The main focus was on the asymptotic behavior of the distributions and their corresponding moments for optically thin slabs, which have been far less studied than optically thick media.

Under collimated illumination, the second moment of the photon path distribution exceeds the ballistic limit $H^2$, and the higher order normalized moments diverge as $\tau ^{2-n}$. For diffuse illumination, ballistic photons are responsible for the divergence of the second normalized moment as $\ln \tau$, which provides a means to estimate $\tau$ from the variance of the distribution. Higher order normalized moments diverge as $\tau ^{2-n}$. When only reflected or transmitted photons are considered (under directional or angularly-integrated viewing), the normalized moments diverge as $\tau ^{1-n}$, implying that the second normalized moment diverges as $\tau ^{-1}$.

Although the detailed calculations are only performed for isotropic scattering, Monte Carlo ray-tracing simulations show that the $\tau$ dependency is unchanged for anisotropic scattering, except for a scaling (related to the probability to be scattered parallel to the slab) in the case of collimated illumination. Most interestingly, under diffuse illumination the asymptotic behavior of the moments is insensitive to the phase function. Monte Carlo simulations were also used to explore the impact of higher orders of scattering, confirming that the asymptotic behavior of the moments is dominated by photons scattered once or twice, with contributions of the same order.

In summary, this mostly theoretical work sheds light on a poorly-investigated research topic. As an important application, it provides physical insight to advance the remote sensing of optically thin media. In particular, our results for collimated (and pulsed) sources associated with directional sensors could be applied to (ranged-gated) bistatic lidar systems deployed either on ground (e.g., [27,28]) or in space (e.g., [29]). By contrast, our results for diffuse illumination and angularly-integrated sensing are better suited for the analysis of laboratory measurements.