*g*), but has minimal dependence on higher moments of the phase function. An empirical expression, having a similar form as the double scattering approximation for LEBS, is found to accurately predict the average penetration depth in the multiple scattering regime. The expression is shown to be accurate for a broad range of experimentally relevant optical properties and spatial coherence lengths.

## 1.

## Introduction

The depth-selective characterization of tissue has many applications for biomedical research. One important example is the early detection of epithelial cancers, which account for more than 95% of all cancer cases. Epithelial layers are typically less than several hundred micrometers thick, thus requiring the penetration depth of the signal to be significantly smaller than the transport mean free path of light. Specialized optical methods must be employed in order to accomplish this high degree of depth-selectivity. Furthermore, the epithelial layer thickness can vary in different organs and locations within an organ. Therefore, a technique should ideally possess the ability to optimize the probed penetration depth to best suit the disease being studied. Proposed mechanisms for rejecting light outside of the tissue layer of interest have included spatial gating,^{1} polarization gating,^{2, 3} directional gating,^{4} and coherence gating.^{5, 6, 7}

Low-coherence enhanced backscattering (LEBS) is a coherence gating technique that extends conventional enhanced backscattering (EBS) by using broadband partial spatial coherence illumination combined with wavelength-resolved detection. LEBS allows for depth-selective spectroscopy and overcomes all of the limitations of conventional EBS measurements from tissue.^{8} Additionally, the instrument is simpler to implement than other interferometer-based techniques and a single measurement can obtain scattering information about a broad range of penetration depths.^{9} The depth-selective characteristics of LEBS have resulted in potential noninvasive cancer screening approaches for colorectal and pancreatic cancers.^{10, 11} However, the penetration depth of the LEBS signal in these commonly utilized experimental scenarios is not simple to estimate. Experimental measurements of the penetration depth can be burdened with the difficulties of measuring samples that are tens of micrometers thick. Monte Carlo simulations can be used to solve the equation of radiative transfer and obtain a simulation of LEBS.^{12, 13} Although accurate, numerical simulations of LEBS are typically time consuming and a closed form expression that can accurately predict the LEBS penetration depth for a large variety of optical properties and experimental scenarios is of great value for the purposes of designing and validating LEBS instruments.

Previous work has resulted in a double-scattering model of the penetration depth of LEBS which is accurate for spatial coherence lengths (*L*
_{sc}) that are much smaller than the scattering mean free path (*l*
_{s}).^{12} This model is informative in demonstrating how the LEBS penetration depth is determined by the coherence volume, however, the range of accuracy for the presented expressions is limited to *L*
_{sc} values that are atypical for LEBS experiments in tissue.^{10, 11} LEBS experiments have utilized *L*
_{sc} values ranging from 17 to 140 μm,^{9} while the double scattering approximation was demonstrated to be accurate for radial separations of only 5 μm for tissue-relevant optical properties. Although the most probable penetration depth was demonstrated to have approximate agreement between the double scattering model and numerical results for *L*
_{sc} < *l*
_{s}, the translation between the most probable penetration depth and more broadly relevant average penetration depth is not straightforward. In this work, we extend this concept to multiple scattering by adjusting the relative proportion of the lateral and axial contributions to the coherence volume. The accuracy of the presented expression is demonstrated with comparisons to Monte Carlo simulations utilizing the Mie, Henyey–Greenstein, and more general phase function obtained from the Whittle–Matérn correlation function. The results are also presented for cases of co-polarized and unpolarized light illumination and collection.

## 2.

## Methods

A publicly available light Monte Carlo code was modified to store the maximum penetration depth of each backscattered ray.^{14} As described in earlier work, the code was implemented to simulate the LEBS peak by obtaining the Fourier transform of the radial reflectance distribution.^{15, 16} One implementation of the code was used to simulate Mie theory scattering and included polarization tracking. In this case, simulations were performed with sphere sizes having the anisotropy factor (*g*) ranging from 0.06 and 0.96. A second implementation of the code was used to evaluate the effect of varying types of phase functions on the LEBS penetration depth. In this case, the scalar-wave version of a general phase function that is based on the Whittle–Matérn correlation function was implemented.^{17, 18}

## 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} F(\theta) = \frac{{2(kl_c)^2 (m - 1)}}{{1 - [ {1 + (2kl_c)^2 }]^{1 - m} }} \cdot \frac{1}{{\{ {1 + [ {2kl_c \sin ({\theta / {2)}}}]^2 }\}^m }},\nonumber\\ \end{eqnarray}\end{document} $$\begin{array}{c}\hfill F\left(\theta \right)=\frac{2{\left(k{l}_{c}\right)}^{2}(m-1)}{1-{\left[1+{\left(2k{l}_{c}\right)}^{2}\right]}^{1-m}}\xb7\frac{1}{{\left\{1+{\left[2k{l}_{c}\mathrm{sin}(\theta /2)\right]}^{2}\right\}}^{m}},\end{array}$$*k*is the wavenumber, and the two parameters

*l*

_{c}and

*m*relate to the correlation length and functional form of the refractive index correlation function. An

*m*< 1.5 corresponds to a power law correlation function (this is the mass fractal regime of light scattering with the fractal dimension

*D*

_{mf}= 2

*m*), an

*m*of 2 indicates an exponential correlation function, and as

*m*approaches large values, a Gaussian correlation function is approached.

^{17}Any intermediate cases are also included as intermediate values of

*m*. When

*m*is 1.5, Eq. 1 is the well known Henyey–Greenstein phase function, in which case, [TeX:] $g = \left\langle {\cos \theta } \right\rangle = 1 - [\sqrt {1 + 4(kl_c)^2 } - 1/2(kl_c)^2]$ $g=\u2329\mathrm{cos}\theta \u232a=1-[\sqrt{1+4{\left(k{l}_{c}\right)}^{2}}-1/2{\left(k{l}_{c}\right)}^{2}]$ . For other values of

*m*,

## 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} 1 - g = \frac{1}{{2\left({m - 2} \right)\left({kl_c } \right)^2 }} - \frac{{2\left({m - 1} \right)}}{{\left({m - 2} \right)\{ {[ {1 + 4\left({kl_c } \right)^2 }]^{m - 1} - 1}\}}}.\nonumber\\ \end{eqnarray}\end{document} $$\begin{array}{c}\hfill 1-g=\frac{1}{2\left(m-2\right){\left(k{l}_{c}\right)}^{2}}-\frac{2\left(m-1\right)}{\left(m-2\right)\left\{{\left[1+4{\left(k{l}_{c}\right)}^{2}\right]}^{m-1}-1\right\}}.\end{array}$$^{17, 18}For our purposes, Eq. 1 will simply be used as a convenient and flexible two parameter phase function that can have varying first and higher order moments.

Due to the scaling property of Monte Carlo simulations, the output of a single simulation with a given phase function can be rescaled to obtain any value of *l*
_{s} by simply rescaling all of the output lengths (e.g., depths, radii, and bin sizes).^{19} A series of LEBS peak enhancement factors were calculated from the zero value of the Fourier transform:
[TeX:]
$E = I_{{\rm LEBS}} (0) = {{\int}\!\!{\int}} p(\vec r)c(\vec r)d \vec r$
$E={I}_{\mathrm{LEBS}}\left(0\right)=\int \int p\left(\stackrel{\u20d7}{r}\right)c\left(\stackrel{\u20d7}{r}\right)d\stackrel{\u20d7}{r}$
, where
[TeX:]
$p(\vec r)$
$p\left(\stackrel{\u20d7}{r}\right)$
is the two-dimensional backscattering probability distribution as a function of the exit radius
[TeX:]
$\vec r$
$\stackrel{\u20d7}{r}$
, and
[TeX:]
$c(\vec r)$
$c\left(\stackrel{\u20d7}{r}\right)$
is the spatial coherence of illumination.^{15, 16} For a given optical property, a saturation curve was constructed for *E* by limiting the maximum depth from which rays were reflected. In other words, the thickness of the medium was modified via post-processing. The normalized derivative of the saturation curve yielded the probability distribution as a function of depth, *p*(*z*). The average penetration depth was then calculated according to the first moment: *D*
_{P} = ∫*z* · *p*(*z*)*dz*.

## 3.

## Results

Under the double scattering approximation, the penetration depth has been shown to follow the form
[TeX:]
$D_{mp} = a\left({l_s^* } \right)^b \left({L_{\rm sc} } \right)^{1 - b}$
${D}_{mp}=a{\left({l}_{s}^{*}\right)}^{b}{\left({L}_{\mathrm{sc}}\right)}^{1-b}$
, with *D*
_{mp} being the most probable penetration depth,
[TeX:]
$l_s^* = {{l_s } \mathord{\left/ {\left({1 - g} \right)}}$
${l}_{s}^{*}={l}_{s}/\left(1-g\right)$
, and *b* being approximately 1/3.^{12} This penetration depth can be physically interpreted as being related to the cube-root of the effective coherence volume as characterized by a cylinder with a radius that is proportional to *L*
_{sc} and a height that is proportional to *l*
_{s} or
[TeX:]
$l_s^*$
${l}_{s}^{*}$
. Here, we build upon this concept by accounting for the changes in the dimensions of the coherence volume due to multiple scattering. Thus, we utilize the same form for the average penetration depth under multiple scattering but modify the definitions of *a* and *b* based on the empirical fitting of Monte Carlo data:

## 3

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {{D_P } \mathord{\left/ {L_{\rm sc} }} = a\left({{{l_s ^* } \mathord{\left/ {L_{\rm sc} }}} \right)^b, \end{equation}\end{document} $${D}_{P}/{L}_{\mathrm{sc}}=a{\left({l}_{s}^{*}/{L}_{\mathrm{sc}}\right)}^{b},$$## 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} a &=& a_1\left({1 - g} \right)^{a_2 }, \nonumber\\ b &=& b_0 + b_1 \left({1 - g} \right)^{b_2 }. \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill a& =& {a}_{1}{\left(1-g\right)}^{{a}_{2}},\hfill \\ \hfill b& =& {b}_{0}+{b}_{1}{\left(1-g\right)}^{{b}_{2}}.\hfill \end{array}$$*a*and the power

*b*depend on the anisotropy factor

*g*, as described by Eq. 4 and the constants in Table 1. The least square fit was optimized to be accurate within the range of 2 < [TeX:] $l_s^*$ ${l}_{s}^{*}$ /

*L*

_{sc}< 33. For most soft tissues where

*g*is approximately 0.9,

*a*and

*b*can be approximated as constants (see Table 1 for values). These values depend on whether or not the measurement configuration utilizes a co-polarized configuration or omits the polarizers, as approximated with the scalar wave case.

## Table 1

Values of constants from Eq. 2.

Constant | Unpolarized | Co-polarized |
---|---|---|

∼a (g = 0.9) | 1.415 | 1.329 |

∼b (g = 0.9) | 0.4773 | 0.4620 |

a_{1} | 1.071 | 0.9479 |

a_{2} | –0.1208 | –0.1468 |

b_{0} | 0.3000 | 0.3182 |

b_{1} | 0.5101 | 0.5523 |

b_{2} | 0.4589 | 0.5842 |

It is important to recognize that Eqs. 3, 4 include the first moment of the phase function (*g*) but do not account for higher moments. However, earlier numerical studies have shown that the angular distribution of the LEBS peak (e.g., the angular peak width) has a significant dependence on higher moments of the phase function.^{15} For example, the Mie phase function with a *g* of 0.9 results in a different LEBS peak width than the Henyey–Greenstein phase function with an identical value of *g*. The effect of higher moments on the LEBS penetration depth can be evaluated with the scalar Whittle–Matérn phase function described by Eq. 1. The first moment of the phase function, *g*, is obtained from Eq. 2. We can modify higher moments of the phase function by varying *m* for a fixed value of *g*.^{18, 20} This results in phase functions that have a different shape but an identical anisotropy factor. Figure 1a is a comparison of penetration depths from three phase functions with varying shapes but an identical anisotropy factor (*g* = 0.9). The penetration depths from the Henyey–Greenstein phase function and phase functions with *m* values of 1.3 and 1.7 are nearly identical. Figure 1b is a plot of the penetration depth as a function of *m* for
[TeX:]
$l_s^*$
${l}_{s}^{*}$
= 1000 μm and *g* = 0.9. The average deviation from the model for the entire range of *m* values shown is 4.5%. Figure 1c shows the dependence of the penetration depth on
[TeX:]
$l_s^*$
${l}_{s}^{*}$
for several values of *g* and an *L*
_{sc} value 60 μm based on Monte Carlo simulations that utilize the Henyey–Greenstein phase function. For small values of
[TeX:]
$l_s^*$
${l}_{s}^{*}$
, the penetration depth does not have a strong dependence on *g*. This is because the coherence area is large enough to capture multiply scattered diffuse light that only contains information about the transport mean free path of the medium. On the other hand, when
[TeX:]
$l_s^*$
${l}_{s}^{*}$
> *L*
_{sc}, there is a clear dependence of the penetration depth on *g* due to the selection of light paths that enter and exit the medium at small separation distances. The model (solid lines) closely predicts the LEBS penetration depth for the entire range of properties shown. The results in Fig. 1 indicate that the first moment of the phase function (*g*) has the most significant effect on the LEBS penetration depth, with higher order moments having a minor contribution. It is therefore justifiable to neglect higher order moments of the phase function in model estimates of the LEBS penetration depth, as this amount of error is acceptable for most applications.

Tissue phantom experiments of LEBS most often use polystyrene microsphere suspensions, and are more accurately modeled with the Mie theory phase function than the Henyey–Greenstein or Whittle–Matérn phase functions. Furthermore, Mie theory-based simulations allow for the possibility of accounting for the polarization state of backscattered light. Figure 2 compares simulations that utilize the Mie theory phase function with predictions of the model. Unpolarized light simulations using the Mie phase function are compared to simulations that utilize the Henyey–Greenstein phase function in Figs. 2a and 2b. Figure 2a shows a plot of the LEBS penetration depth as a function of *g* for unpolarized light simulations that utilize the Mie theory phase function, Henyey–Greenstein phase function simulations, and the model. The anisotropy factor for the Mie theory phase function is determined by selecting the diameter of the particle, while higher order moments of the phase function are unconstrained. Higher values of *g* correspond to larger diameter spheres in the Mie theory-based simulations shown. Figure 2b shows the same phase function comparison for the exponent on
[TeX:]
$l_s^*$
${l}_{s}^{*}$
*/L*
_{sc} from Eq. 4. This exponent can be interpreted as the relative radius to depth ratio of the coherence volume within the scattering medium. It is important to point out that the exponent, *b*, has a significant dependence on *g* and this dependence is well-modeled with Eq. 4. These results again suggest that the effect of the type of phase function plays a minor role in determining the penetration depth. Figure 2b shows that *b* decreases with increasing *g* suggesting that the spatial coherence has a larger influence in determining the penetration depth for large values of *g.* Figure 2c shows the agreement between the model and the often used co-polarized illumination/collection geometry, as simulated with the Mie theory phase function.

A demonstration of the trends that influence the penetration depth of the LEBS signal is illustrated in Fig. 3. The penetration depth increases with increasing
[TeX:]
$l_s^*$
${l}_{s}^{*}$
, as shown in Fig. 3a. The LEBS penetration depth is approximately proportional to *l*
_{s}
^{b}, where *b* < 1, in contrast to conventional EBS which has a penetration depth that is directly proportional to *l*
_{s}. Figure 3b shows the dependence of the penetration depth on *g*. Interestingly, the penetration depth decreases for higher values of the anisotropy factor, an effect that is likely due to the increased number of scattering events that accompany a decrease in *l*
_{s} (higher *g* and constant
[TeX:]
$l_s^*$
${l}_{s}^{*}$
results in smaller *l*
_{s}). The influence of a smaller value of *l*
_{s} outweighs the forward-directed property of the phase function. Figure 3c demonstrates the possibility of adjusting the penetration depth by selecting an *L*
_{sc}. The penetration depth is proportional to *L*
_{sc}
^{1-b}. Again, *b* approaches 1 in the limiting case where *L*
_{sc} approaches ∞ (the conventional EBS regime). The errors between the simulation and model are plotted at the bottom of each panel. In all of the cases the error is less than 10%.

## 4.

## Conclusions

In summary, a simple model for the penetration depth of LEBS that is valid for tissue relevant properties and previously used experimental instruments was presented. It was found that the penetration depth of light contributing to the LEBS peak is primarily dependent on
[TeX:]
$l_s^*$
${l}_{s}^{*}$
, *L*
_{sc}, and *g* with a small effect of higher order moments of the phase function. A generalized model with a similar form to the double scattering approximation was shown to agree with Monte Carlo simulations for both the scalar wave and polarized light case. The maximum error of the model is estimated to be approximately 10%, being sufficient for either system design or approximate data interpretation of acquired tissue measurements because most tissue layers have a relatively large thickness relative to their depth location. Additionally, the penetration depth distribution from LEBS is a broad function^{12} and a 10% error is small relative to the width of this distribution. Thus, we anticipate that the presented model will aid in LEBS system design and allow for a better understanding of the depth location of cancerous and precancerous alterations that have previously been observed in human tissue by providing simple and accurate estimations of the average penetration depth of the LEBS signal.

## Acknowledgments

This work was supported by the National Institute of Health through Grant Nos. R01CA128641 and R01 EB003682.