Influence of linear birefringence in the computation of scattering phase functions

Milos Šormaz; Tobias Stamm; Patrick Jenny

doi:10.1117/1.3503475

1 September 2010 Influence of linear birefringence in the computation of scattering phase functions

Milos Šormaz, Tobias Stamm, Patrick Jenny

Author Affiliations +

Journal of Biomedical Optics, Vol. 15, Issue 5, 055010 (September 2010). https://doi.org/10.1117/1.3503475

Abstract

Birefringent media, like biological tissues, are usually assumed to be uniaxial. For biological tissues, the influence of linear birefringence on the scattering phase function is assumed to be neglectable. In order to examine this, a numerical study of the influence of linear birefringence on the scattering phase function and the resulting backscattering Mueller matrices was performed. It is assumed that the media consist of spherical scattering particles embedded in a nonabsorbing medium, which allows us to employ the Lorenz-Mie theory. In the Monte Carlo framework, the influence of linear birefringence on the components of the electric field vector is captured through the Jones N-matrix formalism. The Lorenz-Mie theory indicates that a given linear birefringence value Δ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.⁹

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 $Δ n = n_{e} - n_{o}$ (difference between refractive indicies $n_{e}$ and $n_{o}$ along the extraordinary and ordinary axes, respectively) is small $(Δ n ⩽ 0.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 $Δ 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 $Δ 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 I (x, s, λ, t)}{d s} = - γ_{t} I (x, s, λ, t) + \frac{γ_{s}}{4 π} \int_{4 π} M (s, s^{'}) I (x, s^{'}, λ, t) d ω^{'},

where

I (x, s, λ, t)

is the Stokes vector,

x

the position in physical space,

s

the propagation direction,

λ

the wavelength, and

t

the time. The vectors

s

and

s^{'}

form a so-called scattering plane, for which the single-scattering Mueller matrix

M (s, s^{'})

is defined. Note that

M (s, s^{'})

describes the interaction between the electromagnetic wave and an isolated particle; a detailed description of the single-scattering Mueller matrix can be found in Ref. 12. Extinction and scattering coefficients are denoted by

γ_{t}

and

γ_{s}

, respectively. A single scattering event requires rotation of the Stokes vector

I = {(I, Q, U, V)}^{T}

into the scattering plane, and its multiplication with the

4 \times 4

single-scattering Mueller matrix is required. The components of the Stokes vector are defined as

2.

I = E_{l} E_{l}^{*} + E_{r} E_{r}^{*} = a_{l}^{2} + a_{r}^{2},

Q = E_{l} E_{l}^{*} - E_{r} E_{r}^{*} = a_{l}^{2} - a_{r}^{2},

U = E_{l}^{*} E_{r} + E_{l} E_{r}^{*} = 2 a_{l} a_{r} \cos δ, and

V = - i (E_{l}^{*} E_{r} - E_{l} E_{r}^{*}) = 2 a_{l} a_{r} \sin δ,

where

E_{l}

and

E_{r}

are the complex parts of the phasors

E_{l} = a_{l} \exp (- i ω t - i δ_{0})

, and

E_{r} = a_{r} \exp (- i ω t - i δ_{0} - i δ)

. Here,

E_{l}

and

E_{r}

are the parallel and perpendicular components of the electric field vector with respect to the scattering plane. The phase difference between

E_{r}

and

E_{l}

is denoted by

δ

,

a_{l} = | a_{l} |

and

a_{r} = | a_{r} |

are the amplitudes of

E_{l}

and

E_{r}

,

ω = 2 π c ∕ λ

is the frequency, and

δ_{0}

is the phase of

E_{l}

at

t = 0

(

c

is the speed of light).

The change of the electric field vector components $E_{l}$ and $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 $Δ t_{s}$ until the next scattering event has to be sampled assuming exponential distribution with expectation $τ_{s} = 1 ∕ c γ_{s}$ . The new position then becomes ${\hat{x}}^{ν + 1} = {\hat{x}}^{ν} + c {\hat{s}}^{ν} Δ t_{s}$ , where $\hat{x}$ and $\hat{s}$ are particle position and propagation direction, respectively. The superscript $ν$ denotes the state after $ν$ collisions. Consistently, the particle clock time $\hat{t}$ is updated as ${\hat{t}}^{ν + 1} = {\hat{t}}^{ν} + Δ t_{s}$ , and its weight $\hat{w}$ is decreased by the absorbed energy $Δ \hat{w} = [1 - \exp (- c γ_{a} Δ {\hat{t}}_{s})] {\hat{w}}^{ν}$ . In order to compute the scattered propagation direction ${\hat{s}}^{ν + 1}$ , which together with the incident propagation direction ${\hat{s}}^{ν}$ defines the scattering plane, the scattering angles $θ$ and $ϕ$ have to be sampled from a joint probability density function given as

Eq. 3

p (θ, ϕ) = P_{11} (θ) + P_{12} (θ) [Q \cos (2 ϕ) + U \sin (2 ϕ)] ∕ I,

where

P_{11} (θ)

and

P_{12} (θ)

are the elements of the single-scattering Mueller matrix. Then,

{\hat{E}}_{l}

and

{\hat{E}}_{r}

of the particle are changed as described in Ref. 9, whereas the new phasors of the particle, i.e.,

{\hat{E}}_{r}^{ν + 1}

and

{\hat{E}}_{l}^{ν + 1}

, point in the direction of

{\hat{s}}^{ν + 1} \times {\hat{s}}^{ν}

and

({\hat{s}}^{ν + 1} \times {\hat{s}}^{ν}) \times {\hat{s}}^{ν + 1}

, respectively. Note that not all particles experience the same number of collisions, and it is important to synchronize them. The local coordinate systems of all detected particles are first rotated into a so-called laboratory coordinate system, and then the Stokes vector is computed from the phasors

{\hat{E}}_{l}

and

{\hat{E}}_{r}

.

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 $Δ 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 $α$ between the photon’s propagation direction $s$ and the extraordinary axis $b$ , and it is given by the expression

Eq. 4

n (α) = \frac{n_{o} n_{e}}{{(n_{e}^{2} \cos^{2} α + n_{o}^{2} \sin^{2} α)}^{1 ∕ 2}} .

As the photon propagates between two successive scattering events, its polarization changes due to retardation

δ_{LB}

caused by the difference in refractive indicies

Δ n (α) = n (α) - n_{o}

, i.e.,

Eq. 5

δ_{LB} = Δ n (α) \frac{2 π d}{λ},

where

d

is the traveled distance. As described in Ref. 2, the photon’s reference frame must be rotated first such that

a_{l}

becomes parallel with the projection

b^{'}

of the extraordinary axis onto the

a_{l} - a_{r}

plane. Note that this rotation of the reference frame changes only the orientation angle of the polarization ellipse; the ellipticity stays the same. Now, the phasors

E_{l}

and

E_{r}

are expressed such that the M-matrix can be employed, i.e., the complex parts of the phasors change between two scattering events as

Eq. 6

The calculation of the M-matrix elements is provided in the next section.

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 $E_{l}$ and $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

N (g_{0}) = [\begin{matrix} n_{1} & n_{4} \\ n_{3} & n_{2} \end{matrix}],

where

g_{0} = π Δ n ∕ λ

is the retardation per unit distance. If a photon’s reference frame is chosen such that

a_{l}

is parallel to

b^{'}

, one can write

Eq. 8

N (g_{0}) = [\begin{matrix} i g_{0} & 0 \\ 0 & - i g_{0} \end{matrix}],

and the elements of the M-matrix

Eq. 9

M (g_{0}) = [\begin{matrix} m_{1} & m_{4} \\ m_{3} & m_{2} \end{matrix}],

can be computed as

10.

m_{1} = \exp (T_{N} z) [(n_{1} - n_{2}) \frac{\cosh Q_{N} z}{2 Q_{N}} + \cosh Q_{N} z)],

m_{2} = \exp (T_{N} z) [- (n_{1} - n_{2}) \frac{\cosh Q_{N} z}{2 Q_{N}} + \cosh Q_{N} z)],

m_{3} = \exp (T_{N} z) \frac{n_{3} \sinh Q_{N} z}{Q_{N}}, and

m_{4} = \exp (T_{N} z) \frac{n_{4} \sinh Q_{N} z}{Q_{N}},

where

T_{N}

and

Q_{N}

are defined as

11.

T_{N} = \frac{n_{1} + n_{2}}{2}, and

Q_{N} = [\frac{{(n_{1} - n_{2})}^{2}}{4} + n_{3} n_{4}] .

The parameter

z

in Eq. 10 denotes the distance between two successive scattering events; i.e., the distance between two collisions. In Ref. 2, it was shown how two effects like linear birefringence and optical activity can be combined in one N-matrix. Note that various polarizing effects can be modeled with the N-matrix formalism. Since we consider only linear birefringence, the elements

m_{3}

and

m_{4}

are equal to zero. Once computed, the elements of the M-matrix are applied to the complex part of the phasors

{\hat{E}}_{l}^{ν}

and

{\hat{E}}_{r}^{ν}

as

Eq. 12

(\begin{matrix} {\hat{E}}_{l}^{ν, new} \\ {\hat{E}}_{r}^{ν, new} \end{matrix}) = (\begin{matrix} m_{1} & m_{4} \\ m_{3} & m_{2} \end{matrix}) (\begin{matrix} {\hat{E}}_{l}^{ν, old} \\ {\hat{E}}_{r}^{ν, old} \end{matrix}) .

Note that the phasors

{\hat{E}}_{l}^{ν}

and

{\hat{E}}_{r}^{ν}

have to be rotated first as described in Sec. 2.2 in order to apply Eq. 12. Although the use of the N-matrix formalism is unnecessarily complex for this study, since only one polarization effect is present, i.e., linear birefringence, it is used to make the developed model more flexible and extendable for future developments.

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 $Δ 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, $Δ n = 0.001$ and $Δ 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 $Δ n$ should not be ignored.

3.1.

Validation

The Monte Carlo implementation described in the previous section was validated with experimental data.² 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 $[I = {(1, 0, 0, 1)}^{T}]$ of wavelength $λ_{vac} = 632.8 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 $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 {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 $s$ and surface normal) is $14 \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 {mm}^{2}$ ) and an acceptance angle of $20 \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 $π ∕ N$ with $N = 1000$ . All spherical particles have the same radius $r = 700 nm$ and the same refractive index $n_{sp} = 1.59$ . The background medium has an ordinary refractive index of $n_{o} = 1.393$ , and the birefringence value $Δ 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 $γ_{s} = 30 {cm}^{- 1}$ was chosen, and for the second, whose results are shown in Fig. 3 , $γ_{s}$ was set to $60 {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}$ .

Fig. 1

Simulation setup for the validation test cases.

Fig. 2

Stokes vector components for $γ_{s} = 30 {cm}^{- 1}$ detected by (a) detector 1 and (b) detector 2. Solid and dashed lines represent numerical and experimental results of Ref. 2, respectively. Dashed-dotted lines are obtained with Scatter3D. Lines labeled with circles, squares, and asterisks represent the V, U, and Q components of the Stokes vector, respectively.

Fig. 3

Stokes vector components for $γ_{s} = 60 {cm}^{- 1}$ detected by (a) detector 1 and (b) detector 2. Solid and dashed lines represent numerical and experimental results of Ref. 2, respectively. Dashed-dotted lines are obtained with Scatter3D. Lines labeled with circles, squares, and asterisks represent the V, U, and Q components of the Stokes vector, respectively.

Fig. 5

Test case 1: Scattering phase functions for same values of $r ∕ λ_{vac} \approx 1.587$ and $n_{sp} = 1.59$ , but $n_{bg}^{o} = n_{o} = 1.33$ (solid lines) and $n_{bg}^{e} = n_{e} = 1.331$ (dashed lines). Note that plots (b) and (c) are enlargements of plot (a) around the center with different magnification.

Fig. 6

Test case 1: PDFs of deflection angle $θ$ for same values of $r ∕ λ_{vac} \approx 1.587$ and $n_{sp} = 1.59$ , but $n_{bg}^{o} = n_{o} = 1.33$ (solid lines) and $n_{bg}^{e} = n_{e} = 1.331$ (dashed lines).

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 $Δ 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.

Fig. 4

Backscattering Mueller matrix for the validation test case computed with the scattering phase function based on $x \approx 4.92$ and $n_{rel} \approx 1.18$ .

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 π r n_{bg} ∕ λ_{vac}$ and the relative refractive index $n_{rel} = n_{sp} ∕ n_{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_{bg} = n_{e}$ and $n_{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 ∕ γ_{s}$ is the optical mean free path length) is considered. Here, the global coordinate system is defined by the unit vectors $e_{1}^{g}$ , $e_{3}^{g}$ , and $e_{2}^{g}$ pointing in the directions of $b$ , the upper surface normal and $e_{3}^{g} \times e_{1}^{g}$ , respectively. The local coordinate system of a particle is defined by the unit vectors $e_{1}^{p} = e_{1}^{g}$ , $e_{2}^{p} = \sqrt{1 ∕ 2} (e_{2}^{g} - e_{3}^{g})$ and $e_{3}^{p} = \sqrt{1 ∕ 2} (e_{2}^{g} + e_{3}^{g})$ . The incident laser beam intersects the medium surface at its center at an angle of $45 \deg$ and is characterized by the electric field vector $E_{i} = \cos (α) e^{- i β} e_{1}^{p} + \sin (α) e^{i β} e_{2}^{p}$ , where $α ∊ [0, π ∕ 2)$ and $β ∊ [0, π)$ are uniformly distributed random numbers. The corresponding Stokes vector of the incident laser beam (defined in the local coordinate system) is ${[1, \cos (2 α), \sin (2 α) \cos (2 β), \sin (2 α) \sin (2 β)]}^{T}$ , and the Stokes vector of the reflected particles is defined in the system with the unit vectors $- e_{1}^{g}$ , $e_{2}^{g}$ , and $- 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 ∕ λ_{vac} \approx 1.587$ and $n_{sp} = 1.59$ , but two different values for $n_{bg}$ , i.e., $n_{bg}^{e} = n_{e} = 1.331$ and $n_{bg}^{o} = n_{o} = 1.33$ , were computed. Note that this corresponds to $Δ 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 $θ$ [normalized product of scattering phase function and $\sin (θ)$ ] is plotted in Fig. 6 .

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

Fig. 7

Test case 2: Scattering phase functions for same values of $r ∕ λ_{vac} \approx 1.587$ and $n_{sp} = 1.59$ , but $n_{bg}^{o} = n_{o} = 1.33$ (solid lines) and $n_{bg}^{e} = n_{e} = 1.34$ (dashed lines). Note that plots (b) and (c) are enlargements of plot (a) around the center with different magnification.

Fig. 8

Test case 2: PDFs of deflection angle $θ$ for same values of $r ∕ λ_{vac} \approx 1.587$ and $n_{sp} = 1.59$ , but $n_{bg}^{o} = n_{o} = 1.33$ (solid lines) and $n_{bg}^{e} = n_{e} = 1.34$ (dashed lines).

The plots in Fig. 9 show the 16 components of the backscattering Mueller matrix computed with $n_{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_{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 $γ_{s}$ was the same in all directions. A quantitative comparison with the backscattering Mueller matrix based on $n_{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, $Δ n = 0.01$ was used for retardation.

Fig. 9

Backscattering Mueller matrix for test case 2 computed with $r ∕ λ_{vac} \approx 1.587$ , $n_{sp} = 1.59$ , and $n_{bg}^{o} = n_{o} = 1.33$ (birefringent medium with $Δ n = 0.01$ ).

Fig. 10

Backscattering Mueller matrix computed with $r ∕ λ_{vac} \approx 1.587$ , $n_{sp} = 1.59$ , and $n_{bg}^{o} = n_{o} = 1.33$ (isotropic medium).

Fig. 11

Normalized difference between the backscattering Mueller matrix of Fig. 9 and the one calculated with the phase function based on $n_{bg} = n_{e} = 1.34$ instead of $n_{bg} = n_{o} = 1.33$ . In both simulations, $Δ n = 0.01$ was used for the retardation. The grayscale values represent the normalized difference of the components.

Next, the difference between two backscattering Mueller matrices (one based on $n_{bg} = n_{o}$ and one based on $n_{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, $Δ 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%.

Fig. 12

Normalized difference between the backscattering Mueller matrix of Fig. 9 and the one calculated with the phase function based on $n_{bg} = n_{e} = 1.331$ instead of $n_{bg} = n_{o} = 1.33$ . In both simulations, $Δ n = 0.001$ was used to compute retardation. The grayscale values represent the normalized difference of the components.

In order to look more closely at the difference between the scattering phase functions based on $n_{bg} = n_{o}$ and $n_{bg} = n_{e}$ , the relative difference $ϵ (θ)$ computed as

Eq. 13

ϵ (θ) = \frac{| p_{n_{o}} (θ) - p_{n_{e}} (θ) |}{p_{n_{o}} (θ)},

where

p_{n_{o}} (θ)

and

p_{n_{e}} (θ)

are the scattered light intensities

[({| S_{1} |}^{2} + {| S_{2} |}^{2}) ∕ 2]

at the angle

θ

for

n_{bg} = n_{o}

and

n_{bg} = n_{e}

, respectively, is plotted in Fig. 13 . The following parameters were used for these investigations:

• Case 1: $r = 1000 nm$ , $Δ n = 0.001$ , $Δ x \approx 0.00997$ .
• Case 2: $r = 1000 nm$ , $Δ n = 0.01$ , $Δ x \approx 0.0997$ .
• Case 3: $r = 100 nm$ , $Δ n = 0.01$ , $Δ x \approx 0.00997$ .
• Case 4: $r = 10 nm$ , $Δ n = 0.01$ , $Δ x \approx 0.000997$ .

Note that in all cases,

λ_{vac} = 630 nm

and

n_{sp} = 1.59

. As expected, the average discrepancy

Eq. 14

R = \frac{\sum_{i = 1}^{N} ϵ (θ)}{N},

grows with the particle size (radius), and for example, in case 2, the average discrepancy of

R \approx 0.1080

results in a significant difference between the two corresponding backscattering Mueller matrices (see Fig. 11). On the other hand, for case 1, the average discrepancy and the difference between the backscattering Mueller matrices are much smaller (see Fig. 12). An interesting effect can be observed in Fig. 13, which shows that for the same

Δ x \approx 0.00997

but different particle radii

r = 1000 nm

(case 1) and

r = 100 nm

(case 3), the average discrepancies are

R = 0.0111

and

R = 0.0612

, respectively.

Fig. 13

Relative difference $ε (θ)$ between scattering phase functions based on $n_{bg} = n_{o}$ and $n_{bg} = n_{e}$ .

Next, scattering phase functions are calculated for two test cases with $r ∕ λ_{vac} \approx 1.587$ but $n_{sp} = 2$ (test case 5) and $n_{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 $ϵ (θ)$ on $θ$ is presented in Fig. 14 , and it can be seen that the average discrepancy $R$ does not necessarily grow with increasing relative refractive index.

Fig. 14

Relative difference $ε (θ)$ between scattering phase functions based on $n_{bg} = n_{o}$ and $n_{bg} = n_{e}$ .

4. Conclusions

For light scattering in biological tissues, the influence of linear birefringence $(Δ n ⩽ 0.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_{rel}$ . The difference in the scattered light intensity governed by phase functions based on $n_{bg} = n_{o}$ and $n_{bg} = n_{e}$ , as the Lorenz-Mie theory shows, is larger for larger scattering particles and also depends on the relative refractive index $n_{rel}$ . Thus, not just the value of linear birefringence $Δ 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_{bg} = n_{o}$ and $n_{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 $E_{l}$ and $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.

References

1.

X. Wang and L. V. Wang, “Propagation of polarized light in birefringent turbid media: a Monte Carlo study,” J. Biomed. Opt., 7 (3), 279 –290 (2002). https://doi.org/10.1117/1.1483315 1083-3668 Google Scholar

2.

M. F. G. Wood, X. Guo, and I. A. Vitkin, “Polarized light propagation in multiply scattering media exhibiting both linear birefringence and optical activity: Monte Carlo model and experimental methodology,” J. Biomed. Opt., 12 (1), 014029 (2007). https://doi.org/10.1117/1.2434980 1083-3668 Google Scholar

3.

N. Gosh, M. F. G. Wood, and I. A. Vitkin, “Mueller matrix decomposition for extraction of individual polarization parameters from complex turbid media exhibiting multiple scattering, optical activity, and linear birefringence,” J. Biomed. Opt., 13 (4), 044036 (2008). https://doi.org/10.1117/1.2960934 1083-3668 Google Scholar

4.

P. Yang, H. Wei, G. W. Kattawar, Y. X. Hu, D. M. Winker, C. A. Hostetler, and B. A. Baum, “Sensitivity of the backscattering Mueller matrix to particle shape and thermodynamic phase,” Appl. Opt., 42 (21), 4389 –4395 (2003). https://doi.org/10.1364/AO.42.004389 0003-6935 Google Scholar

5.

R. Lawless, Y. Xie, P. Yang, G. W. Kattawar, and I. Laszlo, “Polarization and effective Mueller matrix for multiple scattering of light by nonspherical ice crystals,” Opt. Express, 14 (14), 6381 –6393 (2006). https://doi.org/10.1364/OE.14.006381 1094-4087 Google Scholar

6.

M. J. Raković, G. W. Kattawar, M. Mehrübeoğlu, B. D. Cameron, L. V. Wang, S. Rastegar, and G. L. Coté, “Light backscattering polarization patterns from turbid media: theory and experiment,” Appl. Opt., 38 (15), 3399 –3408 (1999). https://doi.org/10.1364/AO.38.003399 0003-6935 Google Scholar

7.

B. Kaplan, G. Ledanois, and B. Drevillon, “Mueller matrix of dense polystyrene latex sphere suspensions: measurements and Monte Carlo simulation,” Appl. Opt., 40 (16), 2769 –2777 (2001). https://doi.org/10.1364/AO.40.002769 0003-6935 Google Scholar

8.

S. K. Gayen and R. R. Alfano, “Emerging optical biomedical imaging techniques,” Opt. Photonics News, 7 (3), 17 –22 (1996). https://doi.org/10.1364/OPN.7.3.000016 1047-6938 Google Scholar

9.

M. Xu, “Electric field Monte Carlo simulation of polarized light propagating in turbid media,” Opt. Express, 12 (26), 6530 –6539 (2004). https://doi.org/10.1364/OPEX.12.006530 1094-4087 Google Scholar

10.

N. Gosh, M. F. G. Wood, and I. A. Vitkin, “Polarimetry in turbid, birefringent, optically active media: a Monte Carlo study of Mueller matrix decomposition in the backscattering geometry,” J. Appl. Phys., 105 102023 (2009). https://doi.org/10.1063/1.3116129 0021-8979 Google Scholar

11.

R. C. Jones, “A new calculus for the treatement of optical systems. VII. Properties of the N-matrices,” J. Opt. Soc. Am. A, 38 (8), 671 –685 (1948). https://doi.org/10.1364/JOSA.38.000671 0740-3232 Google Scholar

12.

H. C. van de Hulst, Light Scattering by Small Particles, Wiley, New York (1957). Google Scholar

13.

M. Šormaz, T. Stamm, and P. Jenny, “Stochastic modeling of polarized light scattering using a Monte Carlo-based stencil method,” J. Opt. Soc. Am. A, 27 (5), 1100 –1110 (2010). https://doi.org/10.1364/JOSAA.27.001100 0740-3232 Google Scholar

14.

M. Šormaz, T. Stamm, S. Mourad, and P. Jenny, “Stochastic modeling of light scattering with fluorescence using a Monte Carlo-based multiscale approach,” J. Opt. Soc. Am. A, 26 (6), 1403 –1413 (2009). https://doi.org/10.1364/JOSAA.26.001403 0740-3232 Google Scholar

15.

M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light, Cambridge University Press, Cambridge, UK (1999). Google Scholar

16.

F. Jaillon and H. Saint-Jalmes, “Description and time reduction of a Monte Carlo code to simulate propagation of polarized light through scattering media,” Appl. Opt., 42 3290 –3296 (2003). https://doi.org/10.1364/AO.42.003290 0003-6935 Google Scholar

17.

D. Côté and I. A. Vitkin, “Robust concentration determination of optically active molecules in turbid media with validated three-dimensional polarization sensitive Monte Carlo calculations,” Opt. Express, 13 148 –163 (2005). https://doi.org/10.1364/OPEX.13.000148 1094-4087 Google Scholar

18.

P. Jenny, M. Vöge, S. Mourad, and T. Stamm, “Modeling light scattering in paper for halftone print,” 443 –447 (2006). Google Scholar

19.

P. Jenny, S. Mourad, T. Stamm, M. Vöge, and K. Simon, “Computing light statistics in heterogeneous media based on a mass weighted probability density function method,” J. Opt. Soc. Am. A, 24 (8), 2206 –2219 (2007). https://doi.org/10.1364/JOSAA.24.002206 0740-3232 Google Scholar

Citation Download Citation

Milos Šormaz, Tobias Stamm, and Patrick Jenny "Influence of linear birefringence in the computation of scattering phase functions," Journal of Biomedical Optics 15(5), 055010 (1 September 2010). https://doi.org/10.1117/1.3503475

Published: 1 September 2010

Access the abstract

JOURNAL ARTICLE
11 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 6 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Scattering

Birefringence

Particles

Backscatter

Refractive index

Monte Carlo methods

Light scattering

1.

Introduction

2.

Theory

Eq. 1

2.

2.1.

Monte Carlo Framework for Polarized Light

Eq. 3

2.2.

Linear Birefringence in a Monte Carlo Framework

Eq. 4

Eq. 5

Eq. 6

2.3.

Jones N-Matrix Formalism

Eq. 7

Eq. 8

Eq. 9

10.

11.

Eq. 12

3.

Numerical Studies

3.1.

Validation

Fig. 1

Fig. 2

Fig. 3

Fig. 5

Fig. 6

Fig. 4

3.2.

Influence of Linear Birefringence on Phase Function

Fig. 7

Fig. 8

Fig. 9

Fig. 10

Fig. 11

Fig. 12

Eq. 13

Eq. 14

Fig. 13

Fig. 14

4.

Conclusions

References

Show All Keywords

Keywords/Phrases

Search In:

Publication Years