## Abstract

Propagation of light into scattering media is a complex problem that can be modeled using statistical methods such as Monte Carlo. Few Monte Carlo programs have so far included the information regarding the status of polarization of light before and after a scattering event. Different approaches have been followed and limited numerical values have been made available to the general public. In this paper, three different ways to build a Monte Carlo program for light propagation with polarization are given. Different groups have used the first two methods; the third method is original. Comparison in between Monte Carlo runs and Adding Doubling program yielded less than 1 % error.

©2005 Optical Society of America

## 1. Introduction

Over the last fifty years polarized light travel through scattering media has been studied by many groups and by the atmospheric optics and oceanography community in particular. An exact solution of the radiative transfer for a plane-parallel atmosphere with Rayleigh scattering was derived by Chandrasekhar in 1950 [1,2]. More complicated geometries proved too complex to be solved analytically. Adding-doubling methods were developed to handle the same kind of geometry for inhomogeneous atmosphere containing randomly oriented particles [3,4].

Kattawar and Plass [5] were the first to calculate the status of polarization of the light after multiple scattering using a Monte Carlo method. They later defined the degree of polarization for homogeneous layers for different solar zenith angles and optical thickness for haze and clouds [6]. Since then many papers have been published on Monte Carlo models and their applications. Kattawar and Adams used a Monte Carlo program to calculate the Stokes vector at any location in a fully inhomogeneous atmosphere-ocean system [7]. Other groups used Monte Carlo codes to analyze changes in polarization of light pulses transmitted through turbid media [8]. Monte Carlo programs were used to measure the polarization properties of dusty spiral galaxies by Bianchi *et al.* [9,10]; a detailed description of their Monte Carlo program was included. Ambirajan and Look [11] developed a backward Monte Carlo model for slab geometry with circularly polarized incident light. Martinez and Maynard [12] wrote a Monte Carlo program for the plane slab geometry and used it to study the Faraday-effect in an optically active medium [13].

Bartel and Hielsher [14] proposed a Monte Carlo model that differs from previous programs; their program used a local coordinate system to keep track of the polarization reference frame. The coordinate system is rotated with standard rotational matrices. They calculated the diffusely back-scattered Mueller Matrix for suspension of polystyrene spheres and their numerical simulations where compared with experimental results. Cameron *et al.* showed similar images both experimentally and numerically [15]. Cote’ *et al.* [16] developed a three dimensional Monte Carlo to study optically active molecules in turbid media and Wang *et al.* [17] used Monte Carlo simulations to simulate light propagation in birefringent media.

Several authors have focused on improving the efficiency of polarized light transport simulations. Rakovic *et al.* [18] proposed a numerical method to simultaneously calculate all 16 elements of the two-dimensional Mueller matrix, and compared it favorably with their Monte Carlo program. Tynes *et al.* [19] used an effective Mueller matrix to predict the total reflectance and transmission through layered media. Kaplan *et al.* [20] measured experimentally the Mueller matrix backscattered from mono-disperse solutions of microspheres at different concentrations and scattering angles, and developed a Monte Carlo program whose speed was optimized by using the next-event point-estimator method [21]. Min Xu [22] proposed a Monte Carlo implementation that traced the electric field through a multiply scattering media in lieu of a Stokes vector; this program is unique in its ability to simulate coherent phenomenon such as speckle.

Of all the Monte Carlo programs mentioned above only Min Xu and Cote’s programs were released to the general public.

In this paper we want to illustrate three different way of implementing Monte Carlo programs that track the status of polarization of scattered light as it propagates in solutions of Mie and Rayleigh scatterers; these codes are released under the *GNU* public license and are available for download [23].

The first Monte Carlo program follows the algorithm of Chandrasekhar and Kattawar and will be called “meridian plane MC” since the photon’s polarization reference plane is always relative to the meridian plane. The second Monte Carlo program uses the algorithm described by Bartel *et al.* and will be called “Euler MC” because it relies on Euler angles used in spherical rotation problems. The third Monte Carlo program propagates a local coordinate system that is used to define the polarization state. Quaternions are used to propagate the coordinate system, hence this Monte Carlo will be termed “quaternion MC”.

## 2. Standard Monte Carlo program

Monte Carlo programs are based on a technique proposed first by Metropolis and Ulam [24]. The main idea is the use of a stochastic approach to model a physical phenomenon. In the biomedical field, Monte Carlo programs have been used to model laser tissue interactions, fluorescence and many other phenomena [25,26]. In Standard Monte Carlo the user is interested only in the distribution of absorbed, transmitter or reflected photons; in the polarized Monte Carlo other parameters must be tracked. A flow-chart that includes both Standard and Polarized Monte Carlo main steps is shown in Fig. 1. The gray boxes are necessary only in the Polarized Monte Carlo; all other steps must be included in both approaches.

A substantial difference between the polarized and standard Monte Carlo is the choice of the azimuthal and scattering angles. In the standard Monte Carlo the scattering angle *α* is chosen from a phase function, the most commonly used one is the Henyey-Greenstein phase function. Generally an azimuthal angle *ψ* is chosen uniformly between 0 and 2*π*. In the polarized Monte Carlo the scattering angle and the azimuth angle are chosen by a “rejection method” (STEP C) that will be discussed in section 9.

An important feature of all polarization codes is tracking the reference plane used to describe the light polarization. Since the status of polarization of a field is defined with respect to a reference plane, the reference plane must be defined (STEP A) and updated for each scattering event (STEP D). The different ways of tracking the reference plane account for the three different Monte Carlo programs described and are presented in detail in the next sections.

## 3. Method 1-Meridian Planes Monte Carlo

Chandrasekhar was the first to envision the scattering of a photon from one location to the next in the form of meridian and scattering planes [1]. Figure 2 shows the geometry of the problem. The photon directions, before and after scattering, are represented as points A and B respectively on the unit sphere. Each photon direction is uniquely described by two angles *θ* and *ϕ*. The first angle, *θ*, is the angle between the initial photon direction and the *Z*-*axis*. The angle *ϕ* is the angle between the meridian plane and the *X*-*Z* plane. The photon direction is also specified by a unit vector *I*_{1}
whose elements are specified by the direction cosines [u_{x}, u_{y}, u_{z}] [27]. The directions before and after scattering are called respectively *I*_{1}
and *I*_{2}
. The unit vector *I*_{1}
and the *Z*-*axis* determine a plane *COA*. The *COA* plane is the meridian plane. The incident field can be decomposed into two orthogonal components *E*
_{‖} and *E*
_{⊥}
that describe the vibration of the electrical field parallel and perpendicular to the meridian plane. The polarization relative to the new meridian plane is described by the Stokes vector *S*.

When a scattering event occurs a new meridian plane is created by the new direction of propagation * I_{2}* and the

*Z*-

*axis*, the Stokes vector must be transformed so that the polarization is properly represented with respect to the new meridian plane. Two rotations of the Stokes vector, in and out of the scattering plane, are required by this method.

## 4. Method 2-Euler Monte Carlo

The second method was developed by Bartel and Hielsher [14] and is based on the following idea: the polarization reference frame is tracked via a triplet of unit vectors that are rotated by an azimuth and scattering angle at each scattering step. The unit vectors are rotated using Euler angles by general transformation matrices [28]. In our implementation only two vectors are rotated for every scattering event * v* and

*; the third unit vector is implicitly defined by the cross product of*

**u***and*

**v***and is calculated only when a photon reaches a boundary. The advantages of this method is that the Stokes vector is rotated only once for each scattering event instead of twice as in the meridian plane method. An added advantage is that the propagation of the unit vectors is straightforward and was easier to implement.*

**u**## 5. Method 3-Quaternion Monte Carlo

The last method is based on quaternion algebra. Quaternions were used in a quantum diffusion Monte Carlo by Benoit and Clary [29] to perform calculations on molecular clusters. Richartz and Hsu [30] used quaternions to calculate the rotations and retardations of a polarized beam going through optical devices. We used quaternion algebra to track the polarization reference plane. Quaternions offer an elegant and intuitive way of handling rigid body rotations [31]. A quaternion is a 4-tuple of real numbers, and it is an element of *R*^{4}
.

A quaternion can also be defined as the sum of a scalar part *q*_{0}
and a vector part *Γ* of the form:

An angle *θ* can be associated with a quaternion, and *q* can be written as

$\mathit{q}=\mathrm{cos}\theta +\mathit{u}\mathrm{sin}\theta $

Where *u*=*q*/sin *θ* is a unit vector, and |*q*_{0}
|^{2}=cos *θ*. This construction facilitate the use of quaternions as rotational operators.

Considering a vector * T* in

*R*

^{3},

*=*

**T***it*

_{1}+

*+*

**j**t2*the quaternion operator*

**k**t_{3}*L*(

*)*

**T**can be interpreted either as a rotation of an angle *2θ* of the coordinate frame with respect to * T* about

*Γ*, or an opposite rotation of the vector

*with respect to a coordinate frame about the axis*

**T***Γ*[31].

*q** is the complex conjugate of the quaternion

*q*

*q**=*q*
_{0}-*Γ*=*q*
_{0}-* iq*1-

**j**q_{2}-

**k**q_{3}

The quaternion * q* in

*R*

^{4}can operate on the vector

*in*

**T***R*

^{3}because the latter can be represented as a quaternion

*=*

**T***t*

_{0}+

*with null scalar element*

**t***t*

_{0}. In the Monte Carlo program the rotation of the reference plane about a fixed vector is achieved with the operator

*L*(

*) where the product of quaternions*

**T***and*

**q***has the form*

**T**$\mathit{q}\mathit{T}=\left[\begin{array}{c}{q}_{0}{t}_{0}-{q}_{1}{t}_{1}-{q}_{2}{t}_{2}-{q}_{3}{t}_{3}\\ {q}_{0}{t}_{1}+{q}_{1}{t}_{0}+{q}_{2}{t}_{3}-{q}_{3}{t}_{2}\\ {q}_{0}{t}_{2}-{q}_{1}{t}_{3}-{q}_{2}{t}_{0}-{q}_{3}{t}_{1}\\ {q}_{0}{t}_{3}+{q}_{1}{t}_{2}-{q}_{2}{t}_{1}+{q}_{3}{t}_{0}\end{array}\right]$

The product of two quaternions is also a quaternion. Quaternions are associative so the * q**

*operation can be executed with two simple multiplications.*

**Tq**## 6. Polarized light Monte Carlo

We explain a Monte Carlo program in detail following the steps of the flowchart in Fig. 1. The steps for all three different methods will be described in parallel.

## 7. Launch

The first step in every Monte Carlo program is to launch one of several million photons in a media. Depending on the geometry of the beam (Gaussian or pencil beam are typical beam shapes) and the incident angle, different launching positions and photon directions cosines must be considered. When the polarization of the field is included some additional steps must be taken. The reference frame of the field must be defined and the Stokes vector describing the incident field polarization must be declared. For illumination perpendicular to the *X*-*Y* plane, the initial photon direction is defined by the following direction cosines [*u*_{x}
, *u*_{y}
, *u*_{z}
]=[0, 0, 1].

## 7.1 Launch-Meridian Planes Monte Carlo

The polarization of the field is relative to the meridian plane.

Figure 1, STEP A: The meridian plane begins with *ϕ*=0, and the reference frame is initially equal to the x-z plane. The initial Stokes vector is relative to this meridian plane.

Figure 1, STEP B: The status of polarization at launch is defined by the Stokes vector * S*=[

*I*

*Q*

*U*

*V*] specified relative to the

*X*-

*Z*plane. For example, if the user selects a Stokes vector

*=[1 1 0 0] the launched photon will be linearly polarized parallel to the*

**S***X*-

*Z*plane.

## 7.2 Launch-Euler Monte Carlo and Quaternion Monte Carlo

For both the Euler Monte Carlo and the Quaternion Monte Carlo we may select the initial reference frame.

Figure 1, STEP A: We use two unit vectors * v* and

*to define the*

**u***field reference frame. [*

**E***=[*

**v***v*

_{x},

*v*

_{y},

*v*

_{z}]=0, 1, 0],

*=[*

**u***u*

_{x},

*u*

_{y},

*u*

_{z}]=[0, 0, 1].

The unit vector * u* represents the direction of photon propagation.

Fig 1, STEP B: For any considered status of polarization of the incident field we define a Stokes vector * S*=[

*I*

*Q*

*U*

*V*] whose reference frame was established in step B; for example: if we want to model an incident field perpendicular to the reference plane we have to assume a Stokes vector

*=[1-1 0 0].*

**S**## 8. Move

The move step in polarized light Monte Carlo programs is executed as in a standard Monte Carlo program. The photon is moved a propagation distance Δs that is calculated based on pseudo random number ζ generated in the interval (0,1].

where *µ*_{t}
=*µ*_{a}
+*µ*_{s}
, *µ*_{a}
is the absorption coefficient, and *µ*_{s}
is the scattering coefficient of the media. The mean free path between every scattering and absorption event is 1/*µ*_{t}
.

The trajectory of the photon is characterized by the direction cosines [*u*_{x}
, *u*_{y}
, *u*_{z}
]. Following Witt [27] the photon position is updated to a new position [*x*
^{′}, *y*
^{′}, *z*
^{′}] with the following equations:

## 9. Drop

As in standard Monte Carlo code the absorption of light by a dye or an absorbing material present in the scattering media is tracked by giving a weight *W* to a photon and updating this weight after every absorption step according to the media albedo.

The albedo is the fractional probability of being scattered, 1-albedo is the fractional probability of being absorbed. The weight of the photon is initially equal to 1. The photon propagates through an absorbing media; after *n* steps the weight of the photon is equal to (albedo)
^{n}
. In our Monte Carlo programs we use matched boundaries so we do not consider weight reduction due to the Fresnel interfaces as in standard Monte Carlo. When the weight of the photon reaches a threshold level the photon is terminated and is considered completely absorbed. If the photon exits the media with a certain weight *W* the corresponding Stokes vector is multiplied by *W* to keep in to account the photon attenuation.

## 10. Rejection Method

A fundamental problem in every Monte Carlo program with polarization information is the choice of the angles *α*(angle of scattering) and *β* (angle of rotation into the scattering plane), Fig. 2. These angles are selected based on the phase function of the considered scatterers and a rejection method. A detailed description of how to implement a rejection method can be found in [18]. The phase function *P*(*α*,*β*) for incident light with a Stokes vector **S**
_{o}=[*I*_{o}
, *Q*_{o}
, *U*_{o}
, *V*_{o}
] is

The parameters s_{11}(*α*) and s_{12}(*α*) are two elements of the scattering matrix *M*(*α*)

The scattering matrix *M*(*α*) determines the polarization properties of the scattering element; we considered spheres. The matrix *M*(*α*) is symmetrical in the spherical case. The terms *s*_{11}
(*α*), *s*_{12}
(*α*), *s*_{34}
(*α*), and *s*_{33}
(*α*) are related to the scattered amplitudes *S*_{1}
and *S*_{2}
[32].

$${s}_{12}=\frac{1}{2}\left({\mid {S}_{2}\mid}^{2}-{\mid {S}_{1}\mid}^{2}\right)$$

$${s}_{33}=\frac{1}{2}\left({S}_{2}^{*}{S}_{1}+{S}_{2}{S}_{1}^{*}\right)$$

$${s}_{34}=-\frac{i}{2}\left({S}_{1}{S}_{2}^{*}-{S}_{2}{S}_{1}^{*}\right)$$

We omitted the angular *α* dependence. *S*_{1}
and *S*_{2}
depend from the size parameter *x*, the complex index of refraction of the particle *n*_{particle}
, the Riccati-Bessel functions, the spherical Bessel functions, and the spherical Henkel functions; a more in depth explanation of these parameters can be found in reference [32]. In our Monte Carlo programs *S*_{1}
and *S*_{2}
for every scattering angle are obtained from a Mie scattering program [32,33].

The phase function *P*(*α*,*β*) has a bivariate dependence on angles *α* and *β*. For unpolarized light where *S*_{o}
=[1, 0, 0, 0], the phase function becomes a function of just the single variable *α*,

Hopcraft *et al.* [34] graphically showed how the phase function depends on both *α* and *β* for linearly polarized light. They pointed out that the phase function is not axis-symmetric for linearly polarized light.

This is visible in Eq. (8). If *Q*_{o}
and *U*_{o}
are not both zero, then the phase function depends both on *α* and *β*. For circularly polarized light *S*_{o}
=[1 0 0 1], the phase function is axis-symmetric since there is no *β* dependence Our simulations of phase function intensities for linear and circular incidence agree with the simulations done by Hopcraft *et al.* and can be seen in Fig. 3. These simulations show the Mie phase function for linear polarized incidence and circular polarized incidence. The asymmetry of the phase function along the *Z*-*axis* for linear incidence is evident. In this Fig., the phase function shows a strong dip at 90°. A small sphere size (diameter=10 nm) was used in both graphs with index of refraction 1.59.

Figure 1, STEP C: The rejection method [35] may be used to generate random variables with a particular distribution. For functions of a single variable, two random numbers are generated. First, *P*_{rand}
is a uniform random deviate between 0 and 1. Second *α*_{rand}
is generated uniformly between 0 and π. The angle *α*_{rand}
is accepted as the new scattering angle if *P*_{rand}
≤*P*(*α*_{rand}
). If *P*_{rand}
>*P*(*α*_{rand}
) *α*_{rand}
and *P*_{rand}
are re-generated and the test is repeated.

When an angle *α*_{rand}
is accepted a similar process is implemented to calculate a new angle *β*_{rand}
. For a bivariate distribution like Eq. (7), three random numbers are generated *α*_{rand}
, *P*_{rand}
, and *β*_{rand}
. The angle *β* is uniformly distributed between 0 and 2*π*. If *P*_{rand}
≤ *P*(*α*_{rand}
, *β*_{rand}
) then both *α*_{rand}
, *ψ*_{rand}
are accepted as the new angles.

## 11. Scatter

Figure 1 STEP D. The scattering step is significantly different in each of the three programs, but it always includes at least one rotation of the Stokes vector reference plane and one multiplication to the scattering matrix *M*(*α*).

## 11.1 Scatter-Meridian Planes Monte Carlo

Three operations are necessary to scatter a photon and track the polarization; these operations are graphically described in the Fig. 4, 5, 6, and 7.

The * E* field is originally defined with respect to a meridian plane

*COA*. The field (visible in all Fig. as blue lines) can be decomposed into its parallel and perpendicular components

*E*

_{‖}and

*E*

_{⊥}. Once the scattering angle

*α*and azimuth angle

*β*have been generated with the rejection method (section 9) the Stokes vector is manipulated three times, the three manipulation will be described in the next three paragraphs.

## Rotation of the reference frame into the scattering plane

First the Stokes vector is rotated so that its reference plane is *BOA*, Fig. 5; this is the scattering plane. This rotation is necessary because the scattering matrix, that defines the elastic interaction of a photon with a sphere, is specified with respect to the frame of reference of the scattering plane. The new Stokes vector * S_{1}* defined relative to the scattering plane (

*BOA*) is found by multiplying by a rotational matrix

*R*(

*β*),

This action corresponds to a counterclockwise rotation by an angle *β* about the direction of propagation.

## Scattering of the photon at an angle α in the scattering plane

The Stokes vector at this point is referred to the scattering plane *ACB*. We can then use the scattering matrix, or Mueller matrix *M*(*α*) of Eq. (9), that is also related to this plane. The scattering matrix *M*(*α*) determines the polarization properties of the scattering sphere. The status of polarization of the field after a scattering event at an angle *α* and the relative Stokes vector can be obtained by the multiplication of the incident Stokes vector * S* to the scattering matrix

*M*(

*α*), Fig. 6.

*α*is the scattering angle between the direction of the photon before and after scattering. The angle

*α*depends on the phase function, and is obtained using the rejection method of section 9.

The direction cosines are updated to take into account the new direction. The new direction cosines will be called **û**. The new trajectory [*û*
_{x}
, *û*
_{y}
, *û*
_{z}
] is calculated based on the angle *α* and *β* and on the current trajectory [*u*_{x}
, *u*_{y}
, *u*_{z}
].

*If | u_{z}
|≈1:*

*$${\hat{u}}_{x}=\mathrm{sin}\left(\alpha \right)\mathrm{cos}\left(\beta \right)$$*

$${\hat{u}}_{y}=\mathrm{sin}\left(\alpha \right)\mathrm{cos}\left(2\beta \right)$$

$${\hat{u}}_{z}=\mathrm{cos}\left(\alpha \right)\frac{{u}_{z}}{\mid {u}_{z}\mid}$$

$${\hat{u}}_{y}=\mathrm{sin}\left(\alpha \right)\mathrm{cos}\left(2\beta \right)$$

$${\hat{u}}_{z}=\mathrm{cos}\left(\alpha \right)\frac{{u}_{z}}{\mid {u}_{z}\mid}$$

*for all other cases:*

*$${\hat{u}}_{x}=\frac{1}{\sqrt{1-{u}_{z}^{2}}}\mathrm{sin}\left(\alpha \right)\left[{u}_{x}{u}_{y}\mathrm{cos}\left(\beta \right)-{u}_{y}\mathrm{sin}\left(\beta \right)\right]+{u}_{x}\mathrm{cos}\left(\alpha \right)$$*

$${\hat{u}}_{y}=\frac{1}{\sqrt{1-{u}_{z}^{2}}}\mathrm{sin}\left(\alpha \right)\left[{u}_{x}{u}_{z}\mathrm{cos}\left(\beta \right)-{u}_{x}\mathrm{sin}\left(\beta \right)\right]+{u}_{y}\mathrm{cos}\left(\alpha \right)$$

$${\hat{u}}_{z}=\sqrt{1-{u}_{z}^{2}}\mathrm{sin}\left(\alpha \right)\mathrm{cos}\left(\beta \right)\left[{u}_{y}{u}_{z}\mathrm{cos}\left(\beta \right)-{u}_{x}\mathrm{sin}\left(\beta \right)\right]+{u}_{z}\mathrm{cos}\left(\alpha \right)$$

$${\hat{u}}_{y}=\frac{1}{\sqrt{1-{u}_{z}^{2}}}\mathrm{sin}\left(\alpha \right)\left[{u}_{x}{u}_{z}\mathrm{cos}\left(\beta \right)-{u}_{x}\mathrm{sin}\left(\beta \right)\right]+{u}_{y}\mathrm{cos}\left(\alpha \right)$$

$${\hat{u}}_{z}=\sqrt{1-{u}_{z}^{2}}\mathrm{sin}\left(\alpha \right)\mathrm{cos}\left(\beta \right)\left[{u}_{y}{u}_{z}\mathrm{cos}\left(\beta \right)-{u}_{x}\mathrm{sin}\left(\beta \right)\right]+{u}_{z}\mathrm{cos}\left(\alpha \right)$$

*Return the reference frame to a new meridian plane*

*The Stokes vector is multiplied by the rotational matrix R(-γ) so that it is referenced to the meridian plane COB (Fig. 7). The angle γ was calculated by Hovenier [36],*

*$$\mathrm{cos}\gamma =\frac{-{u}_{z}+{\hat{u}}_{z}\mathrm{cos}\alpha}{\pm \sqrt{\left(1-{\mathrm{cos}}^{2}\alpha \right)\left(1-{\hat{u}}_{z}^{2}\right)}}$$*

*where α is the scattering angle and the plus sign is taken when π<β<2π and the minus sign is taken when 0<β<π. In summary the Stokes vector S_{new} after scattering is obtained from the Stokes vectors before the scattering event S using the following matrix multiplication:*

*$${\mathit{S}}_{\mathit{new}}=\mathit{R}\left(-\gamma \right)\mathit{M}\left(\alpha \right)\mathit{R}\left(\beta \right)\mathit{S}$$*

*11.2 Scatter-Euler Monte Carlo*

*The rejection method (section 9, STEP C) establishes the scattering angle α and azimuth angle β. The Stokes vector must be rotated by an angle β so that its reference plane coincides with the scattering plane before scattering can be calculated. This is achieved with the rotational matrix R(β) of Eq. (12). The actual scattering by an angle α is done by multiplying the Stokes vector by the scattering matrix M(α), Eq. (9). After scattering the Stokes vector, the reference coordinate system must be adjusted.*

*Figure 1, STEP D: Bartel and Hielscher used a set of rotations called Euler Angles. Their implementation will be shown in Appendix 1, in this section we show what we did in our program.*

*The Stokes vector reference frame was tracked with only two vectors u and v and the rotations were implemented using the rotation matrix R_{euler}
. This formulation is often referred to as the Rodriques’ formula [28] and can be derived from the matrices used by Bartel and Hielscher. This formula is used to rotate a generic vector p about a unit vector k of an angle σ:*

*$$\mathit{p}\u2033=\mathit{p}\xb7\mathit{cos}\left(\sigma \right)+\mathit{sin}\left(\sigma \right)\xb7\left(\mathit{k}\times \mathit{p}\right)+\left[1-\mathit{cos}\left(\sigma \right)\right]\xb7\left(\mathit{k}\xb7\mathit{p}\right)\xb7\mathit{k}$$*

*in matrix format this rotation is expressed as:*

*$${\mathit{R}}_{\mathit{euler}}(\mathit{k},\sigma )=\left[\begin{array}{ccc}{k}_{x}{k}_{x}v+c& {k}_{y}{k}_{x}v-{k}_{z}s& {k}_{z}{k}_{x}v+{k}_{y}s\\ {k}_{x}{k}_{y}v+{k}_{z}s& {k}_{y}{k}_{y}v+c& {k}_{y}{k}_{z}v-{k}_{x}s\\ {k}_{x}{k}_{z}v-{k}_{y}s& {k}_{y}{k}_{z}v+{k}_{x}s& {k}_{z}{k}_{z}v+c\end{array}\right]$$*

*where k=[k_{x}
,k_{y}
,k_{z}
] is the rotational axis, c=cos(σ), s=sin(σ) and v=1-cos(σ). In the Monte Carlo program the rotations of the vector v and u occur in this order: first the unit vector v is rotated about the vector u by an angle β by multiplying v by the rotational matrix R_{euler}(u,β); u remains unchanged in this process (see Fig. 8 and 9).*

*Then u is rotated about the newly generated v by an angle α. This is done multiplying the unit vector u by the rotational matrix R_{euler}(v,α).*

*Finally the Stokes vector is adjusted for the two rotations. The Stokes vector is multiplied by the rotational matrix R(β) of Eq. (12) and then by the scattering matrix M(α) of Eq. (9). The new Stokes vector S_{new} becomes:*

*$${\mathit{S}}_{\mathit{new}}=\mathit{M}\left(\alpha \right)\mathit{R}\left(\beta \right)\mathit{S}$$*

*These steps are repeated for every scattering event.*

*11.3 Scatter-Quaternion Monte Carlo*

*In Quaternion Monte Carlo the rotations of the unit vectors u and v that identify the Stokes vector reference frame are handled with quaternion algebra.*

*STEP D: The first rotation about the vector u by an angle β is accomplished generating the quaternion qβ*

*$${\mathit{q}}_{\beta}=\beta +\mathit{u}=\beta +\mathit{i}{u}_{x}+\mathit{j}{u}_{y}+\mathit{k}{u}_{z}$$*

*and then using the quaternion operator of equation 3 on the unit vector v, *q
_{β}
vq_{β}.*

*To generate the quaternion q_{β}
we used this simple expression:*

*$${q}_{\beta 1}=\mathrm{cos}\left(\frac{\beta}{2}\right)$$*

$${q}_{\beta 2}={u}_{x}\mathrm{sin}\left(\frac{\beta}{2}\right)$$

$${q}_{\beta 3}={u}_{y}\mathrm{sin}\left(\frac{\beta}{2}\right)$$

$${q}_{\beta 4}={u}_{z}\mathrm{sin}\left(\frac{\beta}{2}\right)$$

$${q}_{\beta 2}={u}_{x}\mathrm{sin}\left(\frac{\beta}{2}\right)$$

$${q}_{\beta 3}={u}_{y}\mathrm{sin}\left(\frac{\beta}{2}\right)$$

$${q}_{\beta 4}={u}_{z}\mathrm{sin}\left(\frac{\beta}{2}\right)$$

*The second rotation is by an angle α about the axis v. This is done generating the quaternion q_{α}*

*and using the quaternion operator q*
_{α}
uq_{α}. These steps are frequently optimized for speed and computational efficiency [31]; in our program we use the simplest and more intuitive implementation and not necessarily the fastest one. The Stokes vector is also adjusted for the two rotations as in Eq. (19).*

*12. Photon life and boundaries*

*The life of a photon ends when the photon passes through a boundary or when its weight W value falls below a threshold (typically 0.001) as in standard Monte Carlo.*

*12.1 Meridian Planes Monte Carlo*

*When the photon reaches one of the boundaries and is ready to be collected by a detector, one last rotation of the Stokes vector is necessary to put the photon polarization in the reference frame of the detector.*

*STEP E: For a photon backscattered from a slab, the last rotation to the meridian plane of the detector is by an angle φ*

*The Stokes vector of the reflected photon is multiplied one final time by R(φ). For a transmitted photon the angle φ is:*

*Since Stokes vectors can be superimposed, all photons traveling in the direction of a detector can be added once they have been rotated to the detector reference frame.*

*12.2 Euler and Quaternion Monte Carlo*

*When the photon has reached one of the boundaries and is ready to be collected by a detector two final rotations of the Stokes vector are necessary to put the photon status of polarization in the detector reference frame.*

*STEP E: The first rotation is needed to return the Stokes vector to the meridian plane. To do this w is reconstructed as the cross product of v and u.*

*The angle ε needed to rotate the Stokes vector into a meridian plane is given by:*

*$$\epsilon =0\phantom{\rule{.2em}{0ex}}\mathrm{when}\phantom{\rule{.2em}{0ex}}{v}_{z}=0\phantom{\rule{.2em}{0ex}}\mathrm{and}\phantom{\rule{.2em}{0ex}}{u}_{z}=0$$*

$$\epsilon ={\mathrm{tan}}^{-1}\left(\frac{{v}_{z}}{-{w}_{z}}\right)\phantom{\rule{.2em}{0ex}}\mathrm{in}\phantom{\rule{.2em}{0ex}}\mathrm{all}\phantom{\rule{.2em}{0ex}}\mathrm{other}\phantom{\rule{.2em}{0ex}}\mathrm{cases}$$

$$\epsilon ={\mathrm{tan}}^{-1}\left(\frac{{v}_{z}}{-{w}_{z}}\right)\phantom{\rule{.2em}{0ex}}\mathrm{in}\phantom{\rule{.2em}{0ex}}\mathrm{all}\phantom{\rule{.2em}{0ex}}\mathrm{other}\phantom{\rule{.2em}{0ex}}\mathrm{cases}$$

*This rotation is about the direction of propagation of the photon, i.e. the axis u, Fig. 10.*

*After this rotation E
_{‖} is in a meridian plane, and v_{z}
=0 (Fig. 11).*

*A second rotation by an angle φ about the Z-axis, Fig. 12, will put the photon reference frame in the detector reference frame. All three reference axis w, v, u are affected by this rotation. φ is calculated as in Eq. (23) and (24)*

*$${\mathit{S}}_{\mathit{final}}=\mathit{R}\left(\phi \right)\mathit{R}\left(\epsilon \right)\mathit{S}$$*

*These steps are repeated for every scattering event. At the boundaries the last aligning rotations are the same as in Euler Monte Carlo.*

*13. Comparison with Adding Doubling program*

*In 1991 Evans designed an adding doubling code [3] that could calculate both the radiance and flux of a polarized light beam exiting the atmosphere. The atmosphere can be modeled as a disperse solution of microspheres. One important feature of Evans’ code is the way in which the flux is calculated: when a photon exits the atmosphere either in reflectance or transmission mode the values of the corresponding Stokes vector are summed without relating the polarization to a single detector reference plane. This is equivalent of having several detectors oriented to the meridian planes in different locations on the unit hemisphere. Most users use a single fixed detector to collect all the exiting photons, nevertheless we can model Evans program geometry with a simple modification to the meridian plane code by eliminating step E (Fig. 1) from the program. This prevents the last rotation to the detector meridian plane. Simulations were conducted considering a plane parallel slab of thickness L=4/µ_{s}, with absorption coefficient µ_{a}
=0, and an unpolarized incident beam (Stokes vector [1 0 0 0]) of wavelength λ=632.8 nm. The spheres index of refraction was 1.59 and the media index of refraction was 1. In the Monte Carlo program the number of photon launched was 10^{6}. The value of reflectance and transmittance through the slabs are shown in the following tables:*

*The error for the Monte Carlo runs was the standard deviation of five different runs where only the seed of the random generator was changed. The difference in between our code and Evans results is always less than 1%.*

*14. Conclusions*

*In this paper three different Monte Carlo programs that keep track of the status of polarization after every scattering event were described. The three programs are significantly different from each other yet they yield the same numerical results.*

*The programs were developed for Unix platform in C programming language. All three programs run at approximately the same speed. The three codes are available for download at [23].*

*The Euler’s angle Monte Carlo is easy to implement and intuitive, with a triplet of unit vectors following the photon reference frame after every scattering event. Euler’s rotational matrices are well known since they are often used in rigid body rotations. Moreover since for any scattering event only one rotation of the reference frame is necessary the program should in principle be faster that the meridian plane method.*

*The quaternion Monte Carlo is an elegant solution. Quaternion rotations are used more and more often in computer science to handle animations. Many improvements are available to simplify the generation of the quaternion plus the multiplicative step in order to make them faster [31,38]. Moreover Richartz and Hsu [30] demonstrated that quaternions can be used to simulate a material retardation, hence this methodology could be particularly useful in the modeling of polarized transfer in biological media, where birefringence and retardation play a large role. We intend to explore these algorithms in the future.*

*The compilation of these three models was lengthy; seemingly correct solutions were proven wrong by further testing. A key issue we identified during this process relates to the exiting routine and the backscattered reflection of photons. Many of the papers dealing with polarized Monte Carlo program use as proof of the goodness of their code the shapes of the backscattered Mueller matrix from a solution of microsphere. Since many of the elements of the Mueller matrix are axially symmetric an erroneous exiting routine could still generate apparently correct answers. The best way to test for this is to launch photons at a small angle to the normal (~30 degrees). This assures the breaking of symmetry, the backscattered pattern will be better suited for experimental comparison.*

*Another issue we encountered in the compilations of our Monte Carlo programs was related to the subdivision of the scattering function. Generally a subdivision of 180 angles is acceptable, in Henyey Greenstein based Monte Carlo programs. In our Mie theory based programs the subdivision of the scattering function had to be larger than 1000 angles to obtain accurate values of total reflectance and total transmittance.*

*Previous papers [14, 18] have discussed the correct implementation of the rejection function that chooses the scattering angle α and azimuth angle β. Some author such as Bartel and Hiescher, select first the azimuth angle uniformly between 0 and π and then enter the rejection method algorithm where the scattering angle is selected based on the original choice of β. Rakovic et al. [18] showed that the correct choice of α and β is obtained selecting both angles in the loop of the rejection method. We implemented both method and realized that the two numerical implementations are numerically very similar. Moreover the rejection method is a bottleneck for the Monte Carlo program, optimizing this section of the program would improve dramatically their efficiency.*

*Significant experimental testing was performed on all three programs. The results of these tests will be shown in the next paper [38].*

*Appendix 1*

*Starting from a frame [ e_{r}
, e_{l}
, e_{3}
] (we use the same notation as in their paper [14]) a new frame [${e}_{r}^{\u2033}$
, ${e}_{l}^{\u2033}$
, ${e}_{\mathit{3}}^{\u2033}$
] is calculated with appropriate rotational matrices. Each rotation is about an axis whose location depends upon the preceding rotations. Generally the rotational matrices that transform one reference frame to a rotated one are expressed in this form:*

*${D}_{3}\left(\delta \right)=\left[\begin{array}{ccc}\mathrm{cos}\left(\delta \right)& -\mathrm{sin}\left(\delta \right)& 0\\ \mathrm{sin}\left(\delta \right)& \mathrm{cos}\left(\delta \right)& 0\\ 0& 0& 1\end{array}\right]{D}_{l}\left(\rho \right)=\left[\begin{array}{ccc}\mathrm{cos}\left(\rho \right)& 0& \mathrm{sin}\left(\rho \right)\\ 0& 1& 0\\ -\mathrm{sin}\left(\rho \right)& 0& \mathrm{cos}\left(\rho \right)\end{array}\right]{D}_{r}\left(\omega \right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{cos}\left(\omega \right)& \mathrm{sin}\left(\omega \right)\\ 0& \mathrm{sin}\left(\omega \right)& \mathrm{cos}\left(\omega \right)\end{array}\right]$*

*where δ, ρ, and ω are three general rotation angles and the rotation sequence follows the 3-l-r Euler angle convention. The corresponding matrix multiplication occur in this order:*

*[${e}_{r}^{\u2033}$,${e}_{l}^{\u2033}$,${e}_{\mathit{3}}^{\u2033}$]= D_{3}
(δ) D_{l}
(ρ)D_{r}
(ω)[e_{r}
,e_{l}
,e_{3}
]*

*Only two rotational angles are necessary to follow the Stokes vector reference frame; these angles are α and ψ and the rotation occur about e_{r}
and e_{3}
. At launch e_{r}
=[1 0 0], e_{l}
=[0 1 0] and e_{3}
=[0 0 1], the rotated [${e}_{r}^{\u2033}$
, ${e}_{l}^{\u2033}$
, ${e}_{\mathit{3}}^{\u2033}$
] terms are:*

*${e}_{r}^{\u2033}=\left[\begin{array}{ccc}\mathrm{cos}\left(\psi \right)& -\mathrm{sin}\left(\psi \right)& 0\\ \mathrm{sin}\left(\psi \right)& \mathrm{cos}\left(\psi \right)& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{cos}\left(\alpha \right)& \mathrm{sin}\left(\alpha \right)\\ 0& \mathrm{sin}\left(\alpha \right)& \mathrm{cos}\left(\alpha \right)\end{array}\right]\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]=\left[\begin{array}{c}\mathrm{cos}\left(\psi \right)\\ \mathrm{sin}\left(\psi \right)\\ 0\end{array}\right]$*

*${e}_{l}^{\u2033}=\left[\begin{array}{ccc}\mathrm{cos}\left(\psi \right)& -\mathrm{sin}\left(\psi \right)& 0\\ \mathrm{sin}\left(\psi \right)& \mathrm{cos}\left(\psi \right)& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{cos}\left(\alpha \right)& \mathrm{sin}\left(\alpha \right)\\ 0& \mathrm{sin}\left(\alpha \right)& \mathrm{cos}\left(\alpha \right)\end{array}\right]\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]=\left[\begin{array}{c}-\mathrm{cos}\left(\psi \right)\mathrm{sin}\left(\alpha \right)\\ \mathrm{cos}\left(\psi \right)\mathrm{cos}\left(\alpha \right)\\ \mathrm{sin}\left(\alpha \right)\end{array}\right]$*

*${e}_{3}^{\u2033}=\left[\begin{array}{ccc}\mathrm{cos}\left(\psi \right)& -\mathrm{sin}\left(\psi \right)& 0\\ \mathrm{sin}\left(\psi \right)& \mathrm{cos}\left(\psi \right)& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{cos}\left(\alpha \right)& \mathrm{sin}\left(\alpha \right)\\ 0& \mathrm{sin}\left(\alpha \right)& \mathrm{cos}\left(\alpha \right)\end{array}\right]\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]=\left[\begin{array}{c}\mathrm{sin}\left(\psi \right)\mathrm{sin}\left(\alpha \right)\\ -\mathrm{cos}\left(\psi \right)\mathrm{cos}\left(\alpha \right)\\ \mathrm{cos}\left(\alpha \right)\end{array}\right]$*

*One problem with Euler angles is the occurrence of the “Gimbal Lock”. A Gimbal lock is obtained when two of the axes of rotation line up with the consequent loss of one degree of freedom. Fortunately since the D_{l}
matrix is the unity matrix the Gimbal Lock never occurs*

*References and Links*

**1. **S. Chandrasekhar, *Radiative Transfer*, (Oxford Clarendon Press; 1950).

**2. **S. Chandrasekhar and D. Elbert, “The Illumination and polarization of the sunlight sky on Rayleigh scattering,” Trans. Am. Phil. Soc.44, (1956).

**3. **K. F. Evans and G. L. Stephens, “A new polarized atmospheric radiative transfer model,” J. Quant. Spectrosc. Radiat Transfer. **46**, 413–423, (1991). [CrossRef]

**4. **S. A. Prahl, M. J. C. van Gemert, and A. J. Welch, “Determining the optical properties of turbid media by using the adding-doubling method,” Appl. Opt. **32**, 559–568, (1993). [CrossRef] [PubMed]

**5. **G. W. Kattawar and G. N. Plass “Radiance and polarization of multiple scattered light from haze and clouds,” Appl. Opt. **7**, 1519–1527, (1967). [CrossRef]

**6. **G. W. Kattawar and G. N. Plass, “Degree and direction of polarization of multiple scattered light. 1: Homogeneous cloud layers,” Appl. Opt. **11**, 2851–2865, (1972). [CrossRef] [PubMed]

**7. **G. W. Kattawar and C. N. Adams “Stokes vector calculations of the submarine light field in an atmosphere-ocean with scattering according to a Rayleigh phase matrix: effect of interface refractive index on radiance and polarization,” Limnol. Oceanogr. **34**, 1453–1472, (1989). [CrossRef]

**8. **P. Bruscaglione, G. Zaccanti, and W. Qingnong “Transmission of a pulsed polarized light beam through thick turbid media: numerical results,” Appl. Opt. **32**, 6142–6150, (1993). [CrossRef]

**9. **S. Bianchi, A. Ferrara, and C. Giovannardi, “Monte Carlo simulations of dusty spiral galaxies: extinction and polarization properties,” American Astronomical Society **465**, 137–144, (1996).

**10. **S. Bianchi*Estinzione e polarizzazione della radiazione nelle galassie a spirale*, Tesi di Laura (in Italian), 1994.

**11. **A. Ambirajan and D.C. Look, “A backward Monte Carlo study of the multiple scattering of a polarized laser beam,” J. Quant. Spectrosc. Radiat. Transfer **58**, 171–192, (1997). [CrossRef]

**12. **A. S. Martinez and R. Maynard “Polarization Statistics in Multiple Scattering of light: a Monte Carlo approach,” in *Localization and Propagation of classical waves in random and periodic structures* (Plenum Publishing Corporation New York, 1993).

**13. **A. S. Martinez*Statistique de polarization et effet Faraday en diffusion multiple de la lumiere* Ph.D. Thesis (in French and English), 1984.

**14. **S. Bartel and A. Hielsher, “Monte Carlo simulations of the diffuse backscattering Mueller matrix for highly scattering media,” Appl. Opt. **39**, 1580–1588, (2000). [CrossRef]

**15. **B. D. Cameron, M. J. Rakovic, M. Mehrubeoglu, G. Kattawar, S. Rastegar, L.-H. Wang, and G. L. Cote, “Measurement and calculation of the two-dimensional backscattering Mueller matrix of a turbid medium,” Opt. Lett. **23**, 485–487, (1998). [CrossRef]

**16. **Daniel Côté and I. Alex 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). [CrossRef] [PubMed]

**17. **X. Wang and L. V. Wang, “Propagation of polarized light in birefringent turbid media: A Monte Carlo study,” J. Biomed. Opt. **7**, 279–290, (2002). [CrossRef] [PubMed]

**18. **M. J. Rakovic, G. W. Kattawar, M. Mehrubeoglu, B. D. Cameron, L.-H. Wang, S. Rastegar, and G. L. Cote, “Light backscattering polarization patterns from turbid media: theory and experiment,” Appl. Opt. **38**, 3399–3408, (1999). [CrossRef]

**19. **H. H. Tynes, G. W. Kattawar, E. P. Zege, I. L. Katsev, A. S. Prikhach, and L. I. Chaikovskaya, “Monte Carlo and Multicomponent Approximation Methods for Vector Radiative Transfer by use of Effective Mueller Matrix Calculations,” Appl. Opt. **40**, 400–412, (2001). [CrossRef]

**20. **B. Kaplan, G. Ledanois, and B. Villon, “Mueller Matrix of Dense Polystyrene Latex Sphere Suspensions: Measurements and Monte Carlo Simulation,” Appl. Opt. **40**, 2769–2777, (2001). [CrossRef]

**21. **I. Lux and L. Koblinger, Monte Carlo Particle Transport Methods: Neutron and Photon Calculations (CRC Prss, Boca Raton, Fla., 1991).

**22. **M. Xu, “Electric field Monte Carlo simulation of polarized light propagation in turbid media”, Opt. Express **26**, 6530–6539, (2004). [CrossRef]

**23. **http://omlc.ogi.edu/software/polarization/

**24. **N. Metropolis and S. Ulam, “The Monte Carlo method,” J. Am. Stat. Assoc. **44**, 335–341, (1949). [PubMed]

**25. **S. A. Prahl, M. Keijzer, S. L. Jacques, and A. J. Welch
, “A Monte Carlo model of light propagation in tissue,” in *Dosimetry of Laser Radiation in Medicine and Biology*,
G. Mueller and D. Sliney, eds., SPIEIS 5, 102–111, (1989).

**26. **L. H Wang, S. L. Jacques, and L-Q Zheng, “MCML-Monte Carlo modeling of photon transport in multilayered tissues,” Computer Methods and Programs in Biomedicine **47**,131–146, (1995). [CrossRef] [PubMed]

**27. **A. N. Witt, “Multiple scattering in reflection nebulae I: A Monte Carlo approach,” Astophys. J. Suppl. Ser. **35**, 1–6, (1977). [CrossRef]

**28. **J. J. Craig, *Introduction to robotics. Mechanics and controls*, (Addison-Weseley Publishing Company, 1986).

**29. **D. Benoit and D. Clary “Quaternion formulation of diffusion quantum Monte Carlo for the rotation of rigid molecules in clusters,” J. Chem. Phys. **113**, 5193–5202, (2000). [CrossRef]

**30. **M. Richartz and Y. Hsu, “Analysis of elliptical polarization,” J. Opt. Soc. Am. **39**, 136–157, (1949). [CrossRef]

**31. **J. B. Kuipers, *Quaternions and rotation sequences. A primer with applications to orbits, aerospace, and virtual reality* (Princeton University Press, 1998) [PubMed]

**32. **C. Bohren and D. R. Huffman, *Absorption and scattering of light by small particles*, (Wiley Science Paperback Series,1998). [CrossRef]

**33. **http://omlc.ogi.edu/software/mie/index.html

**34. **K. I. Hopcraft, P. C. Y. Chang, J. G. Walker, and E. Jakeman
, “Properties of Polarized light-beam multiply scattered by Rayleigh medium,” in *Light Scattering from Microstructure*F. Moreno and F. Gonzales eds., Vol. 534, Lecture Notes in Physics, (Springer-Verlag, 2000), pp. 135–158. [CrossRef]

**35. **W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, *Numerical Recipes in C, the art of Scientific Computing*, (Cambridge University Press, 1992).

**36. **J. W. Hovenier, “Symmetry Relationships for Scattering of Polarized Light in a Slab of Randomly Oriented Particles,” Journal of the Atmospheric Sciences **26**, 488–499, (1968). [CrossRef]

**37. **K. Shoemake, “Animating rotation with quaternion curves,” Computer Graphics **19**, 245–254, (1985). [CrossRef]

**38. **J.C. Ramella-Roman, S.A. Prahl, and S.L. Jacques, are preparing a manuscript to be called “Three Monte Carlo programs of polarized light transport into scattering media: part II,”