*n*has a bigger impact on the scattering phase function for large particles. This conclusion is further supported by Monte Carlo simulations, where the phase function was calculated based on the refractive index once in the ordinary direction and once in the extraordinary one. For large particles, comparisons of the resulting backscattering Mueller matrices show significant differences even for small Δ

*n*values.

## 1.

## Introduction

Light scattering and the modeling thereof plays an important role for a wide range of scientific topics and engineering applications ranging from atmospheric research to medicine. It has been demonstrated that certain properties of a turbid medium can be determined from measurements of reflected polarized light.^{1, 2, 3, 4, 5, 6, 7, 8} This kind of diagnostics, however, is based on reverse engineering, which relies on accurate and efficient simulation tools. Two major theories exist to describe such phenomena. The more rigorous approach (analytical theory) is based on Maxwell’s electromagnetic equations, but due to its mathematical complexity and associated high computational cost, the simpler and computationally more efficient transport theory is preferred for many applications. This is based on a vector Boltzmann equation for the Stokes vector, which determines light intensity and polarization state as a function of position, propagation direction, and time. An established approach to solve this equation is the Monte Carlo method, where a large number of computational particles with individual properties are employed to represent the Stokes vector distribution. One possibility is to assign a Stokes vector, position, and propagation direction to each particle and use the so-called Stokes-Mueller formalism to evolve these properties.^{4, 5, 6, 7} In another equivalent approach, which is adopted in this paper, the evolution of complex electric field vectors is computed for each particle.^{9}

Recently, there is a large interest in applying polarized light to examine properties of biological tissues.^{1, 2, 3, 10} Due to anisotropy of the tissue structure and due to the presence of chiral molecules (e.g., glucose and proteins), effects such as linear birefringence and optical activity have to be incorporated into the Monte Carlo methods. The usual assumption regarding linear birefringence is that a tissue is uniaxial, which means that there exists only one axis of anisotropy. The main axis along which the refractive index differs from those along the other main directions is called the extraordinary axis. Note that the ordinary axes along which the refractive index is uniform are perpendicular to the extraordinary axis. In most biological tissues, the birefringence value
$\Delta n={n}_{e}-{n}_{o}$
(difference between refractive indicies
${n}_{e}$
and
${n}_{o}$
along the extraordinary and ordinary axes, respectively) is small
$(\Delta n\u2a7d0.01)$
and, e.g., in Refs. 1, 2, where birefringence was modeled, the assumption that it has no significant influence on the scattering phase function was adopted. However, the effect of linear birefringence on the components of the electric field vector between scattering events is captured by the Jones N-matrix formalism as described in Refs. 2, 11 (retardation). In this paper, we want to examine the impact of
$\Delta n$
on the backscattering Mueller matrices. Assuming spherical scattering particles embedded in a nonabsorbing, homogeneous, isotropic medium allows us to employ the Lorenz-Mie theory for the computation of the phase functions. While it is not straightforward to also take anisotropic background media into account, it is shown how the scattering phase function alters if the refractive index along the extraordinary instead of the ordinary axis is used for its computation. From theory, it also follows that the same value of
$\Delta n$
has a stronger impact on the scattering phase function as the radius of the scattering particles grows.

In Sec. 2, the governing equations and some basic theory are presented together with the method for modeling linear birefringence. In Sec. 3, validation and computational studies are presented, and last, conclusions are given in Sec. 4.

## 2.

## Theory

In the transport theory for polarized light, one has to consider a vector Boltzmann equation, which reads

## Eq. 1

$$\frac{d\mathbf{I}(\mathbf{x},\mathbf{s},\lambda ,t)}{ds}=-{\gamma}_{t}\mathbf{I}(\mathbf{x},\mathbf{s},\lambda ,t)+\frac{{\gamma}_{s}}{4\pi}{\int}_{4\pi}\mathbf{M}(\mathbf{s},{\mathbf{s}}^{\prime})\mathbf{I}(\mathbf{x},{\mathbf{s}}^{\prime},\lambda ,t)\mathrm{d}{\omega}^{\prime},$$## 2.

The change of the electric field vector components ${\mathbf{E}}_{l}$ and ${\mathbf{E}}_{r}$ , which are in fact the Jones vector parameters, has to be computed at every scattering event; a complete description of the method can been found in Refs. 9, 13. Note, however, that the Stokes vector can be computed from ${E}_{l}$ and ${E}_{r}$ at any time for any scattering plane.

## 2.1.

### Monte Carlo Framework for Polarized Light

If the vector Boltzmann Eq. 1 is solved with a Monte Carlo method, a large number of particles has to be employed. To evolve a particle, first the random time $\Delta {t}_{s}$ until the next scattering event has to be sampled assuming exponential distribution with expectation ${\tau}_{s}=1\u2215c{\gamma}_{s}$ . The new position then becomes ${\widehat{\mathit{x}}}^{\nu +1}={\widehat{\mathit{x}}}^{\nu}+c{\widehat{\mathit{s}}}^{\nu}\Delta {t}_{s}$ , where $\widehat{\mathit{x}}$ and $\widehat{\mathit{s}}$ are particle position and propagation direction, respectively. The superscript $\nu $ denotes the state after $\nu $ collisions. Consistently, the particle clock time $\widehat{t}$ is updated as ${\widehat{t}}^{\nu +1}={\widehat{t}}^{\nu}+\Delta {t}_{s}$ , and its weight $\widehat{w}$ is decreased by the absorbed energy $\Delta \widehat{w}=[1-\mathrm{exp}(-c{\gamma}_{a}\Delta {\widehat{t}}_{s})]{\widehat{w}}^{\nu}$ . In order to compute the scattered propagation direction ${\widehat{\mathit{s}}}^{\nu +1}$ , which together with the incident propagation direction ${\widehat{\mathbf{s}}}^{\nu}$ defines the scattering plane, the scattering angles $\theta $ and $\varphi $ have to be sampled from a joint probability density function given as

## Eq. 3

$$p(\theta ,\varphi )={P}_{11}\left(\theta \right)+{P}_{12}\left(\theta \right)[Q\phantom{\rule{0.2em}{0ex}}\mathrm{cos}\left(2\varphi \right)+U\phantom{\rule{0.2em}{0ex}}\mathrm{sin}\left(2\varphi \right)]\u2215I,$$Last, in order to correctly capture effects at the interface between two media with different refractive indicies, Fresnel’s and Snell’s laws were implemented. An implementation of surface effects in the case of unpolarized light, where the scalar Boltzmann equation is solved, is described in Ref. 14. A nice description of the physics of light crossing an interface can be found in Refs. 15. Implementation in the Monte Carlo framework was already described in Refs. 16, 17.

## 2.2.

### Linear Birefringence in a Monte Carlo Framework

Media with anisotropic structure cause light to experience linear birefringence. Here, uniaxial media with a single axis of anisotropy are considered. The birefringence magnitude is defined by $\Delta n={n}_{e}-{n}_{o}$ , where ${n}_{e}$ and ${n}_{o}$ are the refractive indicies along the extraordinary and ordinary axes, respectively. The Monte Carlo method was extended to include the effect of linear birefringence through the Jones N-matrix formalism. Here, we assume that linear birefringence has no effect on reflection and transmission of light at the medium interfaces and focus on its impact on the scattering phase function.

The refractive index is a function of the angle $\alpha $ between the photon’s propagation direction $\mathbf{s}$ and the extraordinary axis $\mathbf{b}$ , and it is given by the expression

## Eq. 4

$$n\left(\alpha \right)=\frac{{n}_{o}{n}_{e}}{{({n}_{e}^{2}\phantom{\rule{0.2em}{0ex}}{\mathrm{cos}}^{2}\phantom{\rule{0.2em}{0ex}}\alpha +{n}_{o}^{2}\phantom{\rule{0.2em}{0ex}}{\mathrm{sin}}^{2}\phantom{\rule{0.2em}{0ex}}\alpha )}^{1\u22152}}.$$## 2.3.

### Jones N-Matrix Formalism

Details about this method can be found in Refs. 11, but for completeness here, the most important part is described. The effect of linear birefringence over an infinitely small optical path segment can be described by the so-called differential N-matrix and integration along the trajectory between two scattering events leads to the $2\times 2$ M-matrix. Since we trace the phasors ${\mathbf{E}}_{l}$ and ${\mathbf{E}}_{r}$ instead of the Stokes vector, the M-matrix can directly be applied. In our study, the depolarization effect occurs only due to multiple scattering, i.e., no depolarization occurs between scattering events.

The N-matrix for linear birefringence is given as

## Eq. 7

$$\mathbf{N}\left({g}_{0}\right)=\left[\begin{array}{cc}{n}_{1}& {n}_{4}\\ {n}_{3}& {n}_{2}\end{array}\right],$$## Eq. 8

$$\mathbf{N}\left({g}_{0}\right)=\left[\begin{array}{cc}i{g}_{0}& 0\\ 0& -i{g}_{0}\end{array}\right],$$## Eq. 9

$$\mathbf{M}\left({g}_{0}\right)=\left[\begin{array}{cc}{m}_{1}& {m}_{4}\\ {m}_{3}& {m}_{2}\end{array}\right],$$## 10.

## 11.

## Eq. 12

$$\left(\begin{array}{c}{\widehat{E}}_{l}^{\nu ,\mathit{\text{new}}}\\ {\widehat{E}}_{r}^{\nu ,\mathit{\text{new}}}\end{array}\right)=\left(\begin{array}{cc}{m}_{1}& {m}_{4}\\ {m}_{3}& {m}_{2}\end{array}\right)\left(\begin{array}{c}{\widehat{E}}_{l}^{\nu ,\mathit{\text{old}}}\\ {\widehat{E}}_{r}^{\nu ,\mathit{\text{old}}}\end{array}\right).$$## 3.

## Numerical Studies

Before performing more relevant numerical investigations, models and implementations of the linear birefringence and surface effects (Fresnel and Snell laws) are validated. Then, scattering phase functions are shown for particles with different radii and refractive indices and for different refractive indices of the background medium, reflecting the difference $\Delta n$ between the refractive indices along extraordinary and ordinary axes. For many applications, it is of interest to estimate properties of scattering media by comparing measured and calculated Mueller matrices, which involves reverse engineering. Therefore, it is essential that the Monte Carlo method properly takes into account all relevant physical effects. So far, however, linear birefringence has typically been neglected. In the following, we want to demonstrate that this simplification may be questionable. For this paper, $\Delta n=0.001$ and $\Delta n=0.01$ are considered, and it is shown, based on the Lorenz-Mie theory, that the scattering phase functions based on ${n}_{e}$ and ${n}_{o}$ differ significantly, thus indicating that $\Delta n$ should not be ignored.

## 3.1.

### Validation

The Monte Carlo implementation described in the previous section was validated with experimental data.^{2} The experimental and numerical results published in Figs. 5 and 6 of Ref. 2 were kindly provided by the first author thereof. The simulation setup is shown here in Fig. 1
. A turbid medium is illuminated with circularly polarized laser light
$[\mathbf{I}={(1,0,0,1)}^{T}]$
of wavelength
${\lambda}_{\mathit{vac}}=632.8\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
propagating in the positive
${x}_{3}$
-direction, entering the turbid medium at the point (2, 0.5, 0). All geometrical dimensions are in cm, and the medium size is
$4\times 1\times 1$
with the extraordinary axis
$\mathbf{b}$
parallel to the
${x}_{1}$
axis of the global coordinate system (see Fig. 1). The detectors have circular shape with an area of
$1\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{2}$
, and their centers are at (2, 0.5, 1) (detector 1) and at (2, 0, 0.5) (detector 2). The acceptance angle (angle between photon’s propagation direction
$\mathbf{s}$
and surface normal) is
$14\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$
. Note that in Ref. 2, the same detectors were used in the experiments, but for the numerical studies, detectors with a slightly different geometry (rectangular area with a surface area of
$1\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{2}$
) and an acceptance angle of
$20\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$
were considered. The turbid medium, surrounded by air with refractive index of 1, is composed of spherical scattering particles embedded in a nonabsorbing host medium, which allows us to employ the Lorenz-Mie theory. The resolution to precompute and tabulate the scattering phase function is
$\pi \u2215N$
with
$N=1000$
. All spherical particles have the same radius
$r=700\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
and the same refractive index
${n}_{\mathit{sp}}=1.59$
. The background medium has an ordinary refractive index of
${n}_{o}=1.393$
, and the birefringence value
$\Delta n={n}_{e}-{n}_{o}$
varies from 0 to
$1.628\times {10}^{-5}$
. For the computations presented here, if not mentioned otherwise, the refractive index
${n}_{o}$
was used to calculate the scattering phase function; linear birefringence was considered only through retardation (as in Refs. 1, 2). Two test cases were considered: for the first, whose results are shown in Fig. 2
, a scattering coefficient of
${\gamma}_{s}=30\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
was chosen, and for the second, whose results are shown in Fig. 3
,
${\gamma}_{s}$
was set to
$60\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
. The number of simulated particles was
$10\times {10}^{6}$
and
$30\times {10}^{6}$
for the media with lower and higher scattering coefficients, respectively. Again, note that in Ref. 2, the number of simulated particles was
${10}^{8}$
.

Figures 2 and 3 represent results “measured” with detector 1, while Figs. 2 and 3 show results for detector 2. In both figures, curves labeled with circles, squares, and asterisks represent changes of the V, U, and Q components, respectively, of the detected Stokes vectors as functions of
$\Delta n$
. Solid and dashed lines show numerical and experimental results published in Ref. 2, while dashed-dotted lines represent results obtained with our method, i.e., with the code Scatter3D.^{13, 14, 18, 19} It can be observed that the results obtained with Scatter3D are in the same range as the numerical results published in Ref. 2. The small difference between the numerical simulations may be attributed mainly to the slightly different detector geometries.

For further validation, the backscattering Mueller matrix computed with Scatter3D (see Fig. 4 ) was compared with the one presented in Fig. 4 of Ref. 1, and very good agreement was achieved.

## 3.2.

### Influence of Linear Birefringence on Phase Function

For the numerical studies presented in this and the following sections, the scattering phase functions were computed according to the Lorenz-Mie theory, whose input consists of the particle size parameter $x=2\pi r{n}_{\mathit{bg}}\u2215{\lambda}_{\mathit{vac}}$ and the relative refractive index ${n}_{\mathit{rel}}={n}_{\mathit{sp}}\u2215{n}_{\mathit{bg}}$ . While not included here, the influence of linear birefringence on the scattering phase function becomes apparent by comparing the scattering phase functions computed with ${n}_{\mathit{bg}}={n}_{e}$ and ${n}_{\mathit{bg}}={n}_{o}$ . Note, however, that this comparison delivers only an indication of the error due to neglecting this influence of linear birefringence, which can be further highlighted by comparing the corresponding backscattering Mueller matrices. Therefore, a domain of size $20{l}_{s}\times 20{l}_{s}\times 4{l}_{s}$ ( ${l}_{s}=1\u2215{\gamma}_{s}$ is the optical mean free path length) is considered. Here, the global coordinate system is defined by the unit vectors ${\mathbf{e}}_{1}^{g}$ , ${\mathbf{e}}_{3}^{g}$ , and ${\mathbf{e}}_{2}^{g}$ pointing in the directions of $\mathbf{b}$ , the upper surface normal and ${\mathbf{e}}_{3}^{g}\times {\mathbf{e}}_{1}^{g}$ , respectively. The local coordinate system of a particle is defined by the unit vectors ${\mathbf{e}}_{1}^{p}={\mathbf{e}}_{1}^{g}$ , ${\mathbf{e}}_{2}^{p}=\sqrt{1\u22152}({\mathbf{e}}_{2}^{g}-{\mathbf{e}}_{3}^{g})$ and ${\mathbf{e}}_{3}^{p}=\sqrt{1\u22152}({\mathbf{e}}_{2}^{g}+{\mathbf{e}}_{3}^{g})$ . The incident laser beam intersects the medium surface at its center at an angle of $45\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$ and is characterized by the electric field vector ${\mathbf{E}}_{i}=\mathrm{cos}\left(\alpha \right){e}^{-i\beta}{\mathbf{e}}_{1}^{p}+\mathrm{sin}\left(\alpha \right){e}^{i\beta}{\mathbf{e}}_{2}^{p}$ , where $\alpha \u220a[0,\pi \u22152)$ and $\beta \u220a[0,\pi )$ are uniformly distributed random numbers. The corresponding Stokes vector of the incident laser beam (defined in the local coordinate system) is ${[1,\mathrm{cos}\left(2\alpha \right),\mathrm{sin}\left(2\alpha \right)\mathrm{cos}\left(2\beta \right),\mathrm{sin}\left(2\alpha \right)\mathrm{sin}\left(2\beta \right)]}^{T}$ , and the Stokes vector of the reflected particles is defined in the system with the unit vectors $-{\mathbf{e}}_{1}^{g}$ , ${\mathbf{e}}_{2}^{g}$ , and $-{\mathbf{e}}_{3}^{g}$ . The number of simulated particles was ${10}^{8}$ , and all reflected particles are sampled on a structured $100\times 100$ grid at the surface of the turbid medium. Studies also reveal, as theoretically expected, that the influence of linear birefringence is more prominent for larger scattering particles and shorter wavelengths. Note that in all simulations, the influence of linear birefringence on reflection/transmission at the substrate surface was neglected. For retardation, however, it was considered according to Eqs. 4, 5.

For the first test case, phase functions for constant $r\u2215{\lambda}_{\mathit{vac}}\approx 1.587$ and ${n}_{\mathit{sp}}=1.59$ , but two different values for ${n}_{\mathit{bg}}$ , i.e., ${n}_{\mathit{bg}}^{e}={n}_{e}=1.331$ and ${n}_{\mathit{bg}}^{o}={n}_{o}=1.33$ , were computed. Note that this corresponds to $\Delta n={n}_{e}-{n}_{o}=0.001$ . It can be observed in Fig. 5 that the difference between the two phase functions is negligible. The probability density function of the deflection angle $\theta $ [normalized product of scattering phase function and $\mathrm{sin}\left(\theta \right)$ ] is plotted in Fig. 6 .

While the same $r$ , ${\lambda}_{\mathit{vac}}$ , and ${n}_{\mathit{sp}}$ as for test case 1 were used in the second test case, $\Delta n$ was increased to 0.01, resulting in ${n}_{\mathit{bg}}^{o}={n}_{o}=1.33$ and ${n}_{\mathit{bg}}^{e}={n}_{e}=1.34$ . As can be observed in Figs. 7 and 8 , the two resulting phase functions differ significantly.

The plots in Fig. 9 show the 16 components of the backscattering Mueller matrix computed with ${n}_{\mathit{bg}}={n}_{o}$ on a $6{l}_{s}\times 6{l}_{s}$ patch on the surface centered around the incident point. Note that the values of each matrix element were normalized by the maximum value of the (1,1) component. The scattering medium is the one used in test case 2, and more details are provided in Ref. 13. As can be seen by comparing Fig. 9 with Fig. 10 , where no retardation is considered (isotropic medium) but with same ${n}_{\mathit{bg}}={n}_{o}$ , linear birefringence results in more diffuse patterns for all components except for (1,1), (1,2), (2,1), and (2,2). Note that the scattering coefficient ${\gamma}_{s}$ was the same in all directions. A quantitative comparison with the backscattering Mueller matrix based on ${n}_{\mathit{bg}}={n}_{e}=1.34$ for the phase function computation is presented in Fig. 11 , where the normalized difference is shown for each of the 16 components. In the calculations of both backscattering Mueller matrices, $\Delta n=0.01$ was used for retardation.

Next, the difference between two backscattering Mueller matrices (one based on ${n}_{\mathit{bg}}={n}_{o}$ and one based on ${n}_{\mathit{bg}}={n}_{e}$ for the computation of the scattering phase function) for scattering media having the properties of test case 1 is shown in Fig. 12 . In both simulations, $\Delta n=0.001$ was used for retardation. As can be observed, the difference for the elements (1,2) and (2,1) is less than 10%.

In order to look more closely at the difference between the scattering phase functions based on ${n}_{\mathit{bg}}={n}_{o}$ and ${n}_{\mathit{bg}}={n}_{e}$ , the relative difference $\u03f5\left(\theta \right)$ computed as

## Eq. 13

$$\u03f5\left(\theta \right)=\frac{|{p}_{{n}_{o}}\left(\theta \right)-{p}_{{n}_{e}}\left(\theta \right)|}{{p}_{{n}_{o}}\left(\theta \right)},$$• Case 1: $r=1000\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ , $\Delta n=0.001$ , $\Delta x\approx 0.00997$ .

• Case 2: $r=1000\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ , $\Delta n=0.01$ , $\Delta x\approx 0.0997$ .

• Case 3: $r=100\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ , $\Delta n=0.01$ , $\Delta x\approx 0.00997$ .

• Case 4: $r=10\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ , $\Delta n=0.01$ , $\Delta x\approx 0.000997$ .

Next, scattering phase functions are calculated for two test cases with $r\u2215{\lambda}_{\mathit{vac}}\approx 1.587$ but ${n}_{\mathit{sp}}=2$ (test case 5) and ${n}_{\mathit{sp}}=2.5$ (test case 6). The refractive indices of the host medium in ordinary and extraordinary directions are ${n}_{o}=1.33$ and ${n}_{e}=1.34$ , respectively. The dependence of $\u03f5\left(\theta \right)$ on $\theta $ is presented in Fig. 14 , and it can be seen that the average discrepancy $R$ does not necessarily grow with increasing relative refractive index.

## 4.

## Conclusions

For light scattering in biological tissues, the influence of linear birefringence $(\Delta n\u2a7d0.01)$ on the scattering phase function is usually neglected, which was the motivation for the performed study. It can be concluded that the influence of linear birefringence in scattering phase function computations depends on the particle size parameter $x$ and the relative refractive index ${n}_{\mathit{rel}}$ . The difference in the scattered light intensity governed by phase functions based on ${n}_{\mathit{bg}}={n}_{o}$ and ${n}_{\mathit{bg}}={n}_{e}$ , as the Lorenz-Mie theory shows, is larger for larger scattering particles and also depends on the relative refractive index ${n}_{\mathit{rel}}$ . Thus, not just the value of linear birefringence $\Delta n$ has to be considered in order to neglect its influence on scattering, but instead the comparison between two scattering phase functions has to be made. All this results in a variable so-called average discrepancy $R$ , which reflects both dependencies and quantifies the difference between two scattering phase functions computed with ${n}_{\mathit{bg}}={n}_{o}$ and ${n}_{\mathit{bg}}={n}_{e}$ . Note that the modeling of the effect of linear birefringence on a scattering phase function would require an extension of the Lorenz-Mie theory for an anisotropic refractive index of the background medium. The present study is just an indication of the problem and shows errors produced if one neglects the influence of linear birefringence on scattering. In this study, the proposed boundary value for the average discrepancy $R$ is 0.01. In all calculations of backscattering Mueller matrices, an influence of linear birefringence on reflection/transmission at the sample boundaries was neglected.

Birefringence is incorporated in our Monte Carlo framework through the Jones N-matrix formalism, and since we deal with the components of the electric field vector ${\mathbf{E}}_{l}$ and ${\mathbf{E}}_{r}$ , a $2\times 2$ M-matrix can be directly applied. The implementation was validated with results from an established reference before numerical studies were conducted.