0
Research Papers: General

Influence of linear birefringence in the computation of scattering phase functions

[+] Author Affiliations
Miloš Šormaz

Swiss Federal Institute of Technology Zürich, Institute of Fluid Dynamics, Sonneggstrasse 3, 8092 Zürich, Switzerland

Tobias Stamm

Swiss Federal Laboratories for Materials Science and Technology, Laboratory for Media Technology, Überlandstrasse 129, 8600 Dübendorf, Switzerland

Patrick Jenny

Swiss Federal Institute of Technology Zürich, Institute of Fluid Dynamics, Sonneggstrasse 3, 8092 Zürich, Switzerland

J. Biomed. Opt. 15(5), 055010 (October 26, 2010). doi:10.1117/1.3503475
History: Received May 27, 2010; Revised August 31, 2010; Accepted September 01, 2010; Published October 26, 2010; Online October 26, 2010
Text Size: A A A

Open Access Open Access

* Address all correspondence to: Miloš Šormaz, Swiss Federal Institute of Technology Zürich (ETH), Institute of Fluid Dynamics, Sonneggstrasse 3, 8092 Zürich, Switzerland. Tel:41-44-632-3756; Fax:41-44-632-1147; milos.sormaz@ifd.mavt.ethz.ch

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.

Figures in this Article

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.18 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.47 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.13,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=neno (difference between refractive indicies ne and no along the extraordinary and ordinary axes, respectively) is small (Δn0.01) and, e.g., in Refs. 12, 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.

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

1dI(x,s,λ,t)ds=γtI(x,s,λ,t)+γs4π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×4 single-scattering Mueller matrix is required. The components of the Stokes vector are defined asDisplay Formula
I=ElEl*+ErEr*=al2+ar2,
Display Formula
Q=ElEl*ErEr*=al2ar2,
Display Formula
U=El*Er+ElEr*=2alarcosδ,and
Display Formula
V=i(El*ErElEr*)=2alarsinδ,
where El and Er are the complex parts of the phasors El=alexp(iωtiδ0), and Er=arexp(iωtiδ0iδ). Here, El and Er are the parallel and perpendicular components of the electric field vector with respect to the scattering plane. The phase difference between Er and El is denoted by δ, al=|al| and ar=|ar| are the amplitudes of El and Er, ω=2πcλ is the frequency, and δ0 is the phase of El at t=0 (c is the speed of light).

The change of the electric field vector components El and Er, 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 El and Er at any time for any scattering plane.

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 Δts until the next scattering event has to be sampled assuming exponential distribution with expectation τs=1cγs. The new position then becomes x̂ν+1=x̂ν+cŝνΔts, where x̂ and ŝ are particle position and propagation direction, respectively. The superscript ν denotes the state after ν collisions. Consistently, the particle clock time t̂ is updated as t̂ν+1=t̂ν+Δts, and its weight ŵ is decreased by the absorbed energy Δŵ=[1exp(cγaΔt̂s)]ŵν. In order to compute the scattered propagation direction ŝν+1, which together with the incident propagation direction ŝν defines the scattering plane, the scattering angles θ and ϕ have to be sampled from a joint probability density function given asDisplay Formula

3p(θ,ϕ)=P11(θ)+P12(θ)[Qcos(2ϕ)+Usin(2ϕ)]I,
where P11(θ) and P12(θ) are the elements of the single-scattering Mueller matrix. Then, Êl and Êr of the particle are changed as described in Ref. 9, whereas the new phasors of the particle, i.e., Êrν+1 and Êlν+1, point in the direction of ŝν+1×ŝν and (ŝν+1×ŝν)×ŝν+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 Êl and Ê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. 1617.

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=neno, where ne and no 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 expressionDisplay Formula

4n(α)=none(ne2cos2α+no2sin2α)12.
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(α)no, i.e.,Display Formula
5δLB=Δn(α)2πdλ,
where d is the traveled distance. As described in Ref. 2, the photon’s reference frame must be rotated first such that al becomes parallel with the projection b of the extraordinary axis onto the alar 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 El and Er are expressed such that the M-matrix can be employed, i.e., the complex parts of the phasors change between two scattering events as Display FormulaThe calculation of the M-matrix elements is provided in the next section.

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×2 M-matrix. Since we trace the phasors El and Er 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 asDisplay Formula

7N(g0)=[n1n4n3n2],
where g0=πΔnλ is the retardation per unit distance. If a photon’s reference frame is chosen such that al is parallel to b, one can writeDisplay Formula
8N(g0)=[ig000ig0],
and the elements of the M-matrixDisplay Formula
9M(g0)=[m1m4m3m2],
can be computed asDisplay Formula
m1=exp(TNz)[(n1n2)coshQNz2QN+cosh(QNz)],
Display Formula
m2=exp(TNz)[(n1n2)coshQNz2QN+cosh(QNz)],
Display Formula
m3=exp(TNz)n3sinhQNzQN,and
Display Formula
m4=exp(TNz)n4sinhQNzQN,
where TN and QN are defined asDisplay Formula
TN=n1+n22,and
Display Formula
QN=[(n1n2)24+n3n4].
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 m3 and m4 are equal to zero. Once computed, the elements of the M-matrix are applied to the complex part of the phasors Êlν and Êrν asDisplay Formula
12(Êlν,newÊrν,new)=(m1m4m3m2)(Êlν,oldÊrν,old).
Note that the phasors Êlν and Êrν have to be rotated first as described in Sec. 2b 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.

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 ne and no differ significantly, thus indicating that Δn should not be ignored.

Validation

The Monte Carlo implementation described in the previous section was validated with experimental data.2 The experimental and numerical results published in Figs. 56 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.8nm propagating in the positive x3-direction, entering the turbid medium at the point (2, 0.5, 0). All geometrical dimensions are in cm, and the medium size is 4×1×1 with the extraordinary axis b parallel to the x1 axis of the global coordinate system (see Fig. 1). The detectors have circular shape with an area of 1mm2, 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 14deg. 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 1mm2) and an acceptance angle of 20deg 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=700nm and the same refractive index nsp=1.59. The background medium has an ordinary refractive index of no=1.393, and the birefringence value Δn=neno varies from 0 to 1.628×105. For the computations presented here, if not mentioned otherwise, the refractive index no was used to calculate the scattering phase function; linear birefringence was considered only through retardation (as in Refs. 12). Two test cases were considered: for the first, whose results are shown in Fig. 2, a scattering coefficient of γs=30cm1 was chosen, and for the second, whose results are shown in Fig. 3, γs was set to 60cm1. The number of simulated particles was 10×106 and 30×106 for the media with lower and higher scattering coefficients, respectively. Again, note that in Ref. 2, the number of simulated particles was 108.

Figures 23 represent results “measured” with detector 1, while Figs. 23 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.1314,1819 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.

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πrnbgλvac and the relative refractive index nrel=nspnbg. While not included here, the influence of linear birefringence on the scattering phase function becomes apparent by comparing the scattering phase functions computed with nbg=ne and nbg=no. 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 20ls×20ls×4ls (ls=1γs is the optical mean free path length) is considered. Here, the global coordinate system is defined by the unit vectors e1g, e3g, and e2g pointing in the directions of b, the upper surface normal and e3g×e1g, respectively. The local coordinate system of a particle is defined by the unit vectors e1p=e1g, e2p=12(e2ge3g) and e3p=12(e2g+e3g). The incident laser beam intersects the medium surface at its center at an angle of 45deg and is characterized by the electric field vector Ei=cos(α)eiβe1p+sin(α)eiβe2p, 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 e1g, e2g, and e3g. The number of simulated particles was 108, and all reflected particles are sampled on a structured 100×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. 45.

For the first test case, phase functions for constant rλvac1.587 and nsp=1.59, but two different values for nbg, i.e., nbge=ne=1.331 and nbgo=no=1.33, were computed. Note that this corresponds to Δn=neno=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 nsp as for test case 1 were used in the second test case, Δn was increased to 0.01, resulting in nbgo=no=1.33 and nbge=ne=1.34. As can be observed in Figs. 78, the two resulting phase functions differ significantly.

Grahic Jump LocationF7 :

Test case 2: Scattering phase functions for same values of rλvac1.587 and nsp=1.59, but nbgo=no=1.33 (solid lines) and nbge=ne=1.34 (dashed lines). Note that plots (b) and (c) are enlargements of plot (a) around the center with different magnification.

Grahic Jump LocationF8 :

Test case 2: PDFs of deflection angle θ for same values of rλvac1.587 and nsp=1.59, but nbgo=no=1.33 (solid lines) and nbge=ne=1.34 (dashed lines).

The plots in Fig. 9 show the 16 components of the backscattering Mueller matrix computed with nbg=no on a 6ls×6ls 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 nbg=no, 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 nbg=ne=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.

Next, the difference between two backscattering Mueller matrices (one based on nbg=no and one based on nbg=ne 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%.

In order to look more closely at the difference between the scattering phase functions based on nbg=no and nbg=ne, the relative difference ϵ(θ) computed asDisplay Formula

13ϵ(θ)=|pno(θ)pne(θ)|pno(θ),
where pno(θ) and pne(θ) are the scattered light intensities [(|S1|2+|S2|2)2] at the angle θ for nbg=no and nbg=ne, respectively, is plotted in Fig. 13. The following parameters were used for these investigations:

  • Case 1: r=1000nm, Δn=0.001, Δx0.00997.
  • Case 2: r=1000nm, Δn=0.01, Δx0.0997.
  • Case 3: r=100nm, Δn=0.01, Δx0.00997.
  • Case 4: r=10nm, Δn=0.01, Δx0.000997.

Note that in all cases, λvac=630nm and nsp=1.59. As expected, the average discrepancyDisplay Formula

14R=i=1Nϵ(θ)N,
grows with the particle size (radius), and for example, in case 2, the average discrepancy of R0.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 Δx0.00997 but different particle radii r=1000nm (case 1) and r=100nm (case 3), the average discrepancies are R=0.0111 and R=0.0612, respectively.

Grahic Jump LocationF1 :

Simulation setup for the validation test cases.

Grahic Jump LocationF2 :

Stokes vector components for γs=30cm1 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.

Grahic Jump LocationF3 :

Stokes vector components for γs=60cm1 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.

Grahic Jump LocationF4 :

Backscattering Mueller matrix for the validation test case computed with the scattering phase function based on x4.92 and nrel1.18.

Grahic Jump LocationF5 :

Test case 1: Scattering phase functions for same values of rλvac1.587 and nsp=1.59, but nbgo=no=1.33 (solid lines) and nbge=ne=1.331 (dashed lines). Note that plots (b) and (c) are enlargements of plot (a) around the center with different magnification.

Grahic Jump LocationF6 :

Test case 1: PDFs of deflection angle θ for same values of rλvac1.587 and nsp=1.59, but nbgo=no=1.33 (solid lines) and nbge=ne=1.331 (dashed lines).

Grahic Jump LocationF9 :

Backscattering Mueller matrix for test case 2 computed with rλvac1.587, nsp=1.59, and nbgo=no=1.33 (birefringent medium with Δn=0.01).

Grahic Jump LocationF10 :

Backscattering Mueller matrix computed with rλvac1.587, nsp=1.59, and nbgo=no=1.33 (isotropic medium).

Grahic Jump LocationF11 :

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

Grahic Jump LocationF12 :

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

Grahic Jump LocationF13 :

Relative difference ε(θ) between scattering phase functions based on nbg=no and nbg=ne.

Next, scattering phase functions are calculated for two test cases with rλvac1.587 but nsp=2 (test case 5) and nsp=2.5 (test case 6). The refractive indices of the host medium in ordinary and extraordinary directions are no=1.33 and ne=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.

Grahic Jump LocationF14 :

Relative difference ε(θ) between scattering phase functions based on nbg=no and nbg=ne.

For light scattering in biological tissues, the influence of linear birefringence (Δn0.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 nrel. The difference in the scattered light intensity governed by phase functions based on nbg=no and nbg=ne, as the Lorenz-Mie theory shows, is larger for larger scattering particles and also depends on the relative refractive index nrel. 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 nbg=no and nbg=ne. 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 El and Er, a 2×2 M-matrix can be directly applied. The implementation was validated with results from an established reference before numerical studies were conducted.

Wang  X., and Wang  L. V., “ Propagation of polarized light in birefringent turbid media: a Monte Carlo study. ,” J. Biomed. Opt..  1083-3668 7, (3 ), 279–290  ((2002)).
Wood  M. F. G., , Guo  X., , and Vitkin  I. A., “ Polarized light propagation in multiply scattering media exhibiting both linear birefringence and optical activity: Monte Carlo model and experimental methodology. ,” J. Biomed. Opt..  1083-3668 12, (1 ), 014029  ((2007)).
Gosh  N., , Wood  M. F. G., , and Vitkin  I. A., “ Mueller matrix decomposition for extraction of individual polarization parameters from complex turbid media exhibiting multiple scattering, optical activity, and linear birefringence. ,” J. Biomed. Opt..  1083-3668 13, (4 ), 044036  ((2008)).
Yang  P., , Wei  H., , Kattawar  G. W., , Hu  Y. X., , Winker  D. M., , Hostetler  C. A., , and Baum  B. A., “ Sensitivity of the backscattering Mueller matrix to particle shape and thermodynamic phase. ,” Appl. Opt..  0003-6935 42, (21 ), 4389–4395  ((2003)).
Lawless  R., , Xie  Y., , Yang  P., , Kattawar  G. W., , and Laszlo  I., “ Polarization and effective Mueller matrix for multiple scattering of light by nonspherical ice crystals. ,” Opt. Express.  1094-4087 14, (14 ), 6381–6393  ((2006)).
Raković  M. J., , Kattawar  G. W., , Mehrübeoğlu  M., , Cameron  B. D., , Wang  L. V., , Rastegar  S., , and Coté  G. L., “ Light backscattering polarization patterns from turbid media: theory and experiment. ,” Appl. Opt..  0003-6935 38, (15 ), 3399–3408  ((1999)).
Kaplan  B., , Ledanois  G., , and Drevillon  B., “ Mueller matrix of dense polystyrene latex sphere suspensions: measurements and Monte Carlo simulation. ,” Appl. Opt..  0003-6935 40, (16 ), 2769–2777  ((2001)).
Gayen  S. K., and Alfano  R. R., “ Emerging optical biomedical imaging techniques. ,” Opt. Photonics News.  1047-6938 7, (3 ), 17–22  ((1996)).
Xu  M., “ Electric field Monte Carlo simulation of polarized light propagating in turbid media. ,” Opt. Express.  1094-4087 12, (26 ), 6530–6539  ((2004)).
Gosh  N., , Wood  M. F. G., , and Vitkin  I. A., “ Polarimetry in turbid, birefringent, optically active media: a Monte Carlo study of Mueller matrix decomposition in the backscattering geometry. ,” J. Appl. Phys..  0021-8979 105, , 102023  ((2009)).
Jones  R. C., “ A new calculus for the treatement of optical systems. VII. Properties of the N-matrices. ,” J. Opt. Soc. Am. A.  0740-3232 38, (8 ), 671–685  ((1948)).
van de Hulst  H. C.,  Light Scattering by Small Particles. ,  Wiley ,  New York  ((1957)).
Šormaz  M., , Stamm  T., , and Jenny  P., “ Stochastic modeling of polarized light scattering using a Monte Carlo-based stencil method. ,” J. Opt. Soc. Am. A.  0740-3232 27, (5 ), 1100–1110  ((2010)).
Šormaz  M., , Stamm  T., , Mourad  S., , and Jenny  P., “ Stochastic modeling of light scattering with fluorescence using a Monte Carlo-based multiscale approach. ,” J. Opt. Soc. Am. A.  0740-3232 26, (6 ), 1403–1413  ((2009)).
Born  M., and Wolf  E.,  Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light. ,  Cambridge University Press ,  Cambridge, UK  ((1999)).
Jaillon  F., and Saint-Jalmes  H., “ Description and time reduction of a Monte Carlo code to simulate propagation of polarized light through scattering media. ,” Appl. Opt..  0003-6935 42, , 3290–3296  ((2003)).
Côté  D., and Vitkin  I. A., “ Robust concentration determination of optically active molecules in turbid media with validated three-dimensional polarization sensitive Monte Carlo calculations. ,” Opt. Express.  1094-4087 13, , 148–163  ((2005)).
Jenny  P., , Vöge  M., , Mourad  S., , and Stamm  T., “ Modeling light scattering in paper for halftone print. ,” in  Proc. CGIV 2006-Third European Conference on Color in Graphics, Imaging and Vision. ,  IS&T ,  Leeds, UK , pp. 443–447  ((2006)).
Jenny  P., , Mourad  S., , Stamm  T., , Vöge  M., , and Simon  K., “ Computing light statistics in heterogeneous media based on a mass weighted probability density function method. ,” J. Opt. Soc. Am. A.  0740-3232 24, (8 ), 2206–2219  ((2007)).
© 2010 Society of Photo-Optical Instrumentation Engineers

Citation

Miloš Šormaz ; Tobias Stamm and Patrick Jenny
"Influence of linear birefringence in the computation of scattering phase functions", J. Biomed. Opt. 15(5), 055010 (October 26, 2010). ; http://dx.doi.org/10.1117/1.3503475


Access This Article

Figures

Grahic Jump LocationF7 :

Test case 2: Scattering phase functions for same values of rλvac1.587 and nsp=1.59, but nbgo=no=1.33 (solid lines) and nbge=ne=1.34 (dashed lines). Note that plots (b) and (c) are enlargements of plot (a) around the center with different magnification.

Grahic Jump LocationF8 :

Test case 2: PDFs of deflection angle θ for same values of rλvac1.587 and nsp=1.59, but nbgo=no=1.33 (solid lines) and nbge=ne=1.34 (dashed lines).

Grahic Jump LocationF1 :

Simulation setup for the validation test cases.

Grahic Jump LocationF2 :

Stokes vector components for γs=30cm1 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.

Grahic Jump LocationF3 :

Stokes vector components for γs=60cm1 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.

Grahic Jump LocationF4 :

Backscattering Mueller matrix for the validation test case computed with the scattering phase function based on x4.92 and nrel1.18.

Grahic Jump LocationF5 :

Test case 1: Scattering phase functions for same values of rλvac1.587 and nsp=1.59, but nbgo=no=1.33 (solid lines) and nbge=ne=1.331 (dashed lines). Note that plots (b) and (c) are enlargements of plot (a) around the center with different magnification.

Grahic Jump LocationF6 :

Test case 1: PDFs of deflection angle θ for same values of rλvac1.587 and nsp=1.59, but nbgo=no=1.33 (solid lines) and nbge=ne=1.331 (dashed lines).

Grahic Jump LocationF9 :

Backscattering Mueller matrix for test case 2 computed with rλvac1.587, nsp=1.59, and nbgo=no=1.33 (birefringent medium with Δn=0.01).

Grahic Jump LocationF10 :

Backscattering Mueller matrix computed with rλvac1.587, nsp=1.59, and nbgo=no=1.33 (isotropic medium).

Grahic Jump LocationF11 :

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

Grahic Jump LocationF12 :

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

Grahic Jump LocationF13 :

Relative difference ε(θ) between scattering phase functions based on nbg=no and nbg=ne.

Grahic Jump LocationF14 :

Relative difference ε(θ) between scattering phase functions based on nbg=no and nbg=ne.

Tables

References

Wang  X., and Wang  L. V., “ Propagation of polarized light in birefringent turbid media: a Monte Carlo study. ,” J. Biomed. Opt..  1083-3668 7, (3 ), 279–290  ((2002)).
Wood  M. F. G., , Guo  X., , and Vitkin  I. A., “ Polarized light propagation in multiply scattering media exhibiting both linear birefringence and optical activity: Monte Carlo model and experimental methodology. ,” J. Biomed. Opt..  1083-3668 12, (1 ), 014029  ((2007)).
Gosh  N., , Wood  M. F. G., , and Vitkin  I. A., “ Mueller matrix decomposition for extraction of individual polarization parameters from complex turbid media exhibiting multiple scattering, optical activity, and linear birefringence. ,” J. Biomed. Opt..  1083-3668 13, (4 ), 044036  ((2008)).
Yang  P., , Wei  H., , Kattawar  G. W., , Hu  Y. X., , Winker  D. M., , Hostetler  C. A., , and Baum  B. A., “ Sensitivity of the backscattering Mueller matrix to particle shape and thermodynamic phase. ,” Appl. Opt..  0003-6935 42, (21 ), 4389–4395  ((2003)).
Lawless  R., , Xie  Y., , Yang  P., , Kattawar  G. W., , and Laszlo  I., “ Polarization and effective Mueller matrix for multiple scattering of light by nonspherical ice crystals. ,” Opt. Express.  1094-4087 14, (14 ), 6381–6393  ((2006)).
Raković  M. J., , Kattawar  G. W., , Mehrübeoğlu  M., , Cameron  B. D., , Wang  L. V., , Rastegar  S., , and Coté  G. L., “ Light backscattering polarization patterns from turbid media: theory and experiment. ,” Appl. Opt..  0003-6935 38, (15 ), 3399–3408  ((1999)).
Kaplan  B., , Ledanois  G., , and Drevillon  B., “ Mueller matrix of dense polystyrene latex sphere suspensions: measurements and Monte Carlo simulation. ,” Appl. Opt..  0003-6935 40, (16 ), 2769–2777  ((2001)).
Gayen  S. K., and Alfano  R. R., “ Emerging optical biomedical imaging techniques. ,” Opt. Photonics News.  1047-6938 7, (3 ), 17–22  ((1996)).
Xu  M., “ Electric field Monte Carlo simulation of polarized light propagating in turbid media. ,” Opt. Express.  1094-4087 12, (26 ), 6530–6539  ((2004)).
Gosh  N., , Wood  M. F. G., , and Vitkin  I. A., “ Polarimetry in turbid, birefringent, optically active media: a Monte Carlo study of Mueller matrix decomposition in the backscattering geometry. ,” J. Appl. Phys..  0021-8979 105, , 102023  ((2009)).
Jones  R. C., “ A new calculus for the treatement of optical systems. VII. Properties of the N-matrices. ,” J. Opt. Soc. Am. A.  0740-3232 38, (8 ), 671–685  ((1948)).
van de Hulst  H. C.,  Light Scattering by Small Particles. ,  Wiley ,  New York  ((1957)).
Šormaz  M., , Stamm  T., , and Jenny  P., “ Stochastic modeling of polarized light scattering using a Monte Carlo-based stencil method. ,” J. Opt. Soc. Am. A.  0740-3232 27, (5 ), 1100–1110  ((2010)).
Šormaz  M., , Stamm  T., , Mourad  S., , and Jenny  P., “ Stochastic modeling of light scattering with fluorescence using a Monte Carlo-based multiscale approach. ,” J. Opt. Soc. Am. A.  0740-3232 26, (6 ), 1403–1413  ((2009)).
Born  M., and Wolf  E.,  Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light. ,  Cambridge University Press ,  Cambridge, UK  ((1999)).
Jaillon  F., and Saint-Jalmes  H., “ Description and time reduction of a Monte Carlo code to simulate propagation of polarized light through scattering media. ,” Appl. Opt..  0003-6935 42, , 3290–3296  ((2003)).
Côté  D., and Vitkin  I. A., “ Robust concentration determination of optically active molecules in turbid media with validated three-dimensional polarization sensitive Monte Carlo calculations. ,” Opt. Express.  1094-4087 13, , 148–163  ((2005)).
Jenny  P., , Vöge  M., , Mourad  S., , and Stamm  T., “ Modeling light scattering in paper for halftone print. ,” in  Proc. CGIV 2006-Third European Conference on Color in Graphics, Imaging and Vision. ,  IS&T ,  Leeds, UK , pp. 443–447  ((2006)).
Jenny  P., , Mourad  S., , Stamm  T., , Vöge  M., , and Simon  K., “ Computing light statistics in heterogeneous media based on a mass weighted probability density function method. ,” J. Opt. Soc. Am. A.  0740-3232 24, (8 ), 2206–2219  ((2007)).

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging & repositioning the boxes below.

Related Book Chapters

Topic Collections

Advertisement

Buy this article ($18 for members, $25 for non-members).
Sign In