1. Introduction

The present work constitutes a continuation of the comparative study initiated in [1] and aimed at inspecting the (dis)similarities between the Stokes–Mueller and the Jones formalisms in the framework of Monte Carlo (MC) simulations. Whereas the first part of our investigation dealt with theoretical considerations, this second segment provides an empirical support, whereby the results yielded by a series of numerical experiments are examined. We conducted light propagation simulations with two different and independent “polarized MC” programs: (i) The first, “Program U,” uses the Stokes–Mueller notation, implemented in Ulm—likewise subject to a comparison with Maxwell’s theory [2])—and (ii) the second, “Program B,” employs the Jones formalism, developed in Bern—based on the concepts detailed by Rička and Frenz [3,4] and comprehensively described elsewhere [5–7].

The purpose of this study is twofold. On the one hand we demonstrate that both models yield identical simulation results if the coherent properties of multiply scattered light are disregarded, no matter whether the probing light is polarized or unpolarized. On the other hand, this work serves as a “validation” for both programs. This validation bears a particular significance since existing polarized MC programs principally rest on qualitative comparisons with real experiments. Quantitative comparisons, be it with experiments or other simulation programs, are difficult to perform [8], and to our knowledge, a comparison as the one presented here has not been undertaken thus far. The major difficulty stems from trying to “reconcile differences caused by conventions” [9] involved in the modeling of light polarization. The ellipsometry community has attempted to circumvent this issue by introducing the so-called Nebraska conventions [10]. The authors of the latter have acknowledged that “arbitrary choices have to be made at many points. Some of these choices, although of profound consequence, are not recognizable from the results, and the comparison of data from different sources may be impossible.” Hence, we pursue here the effort launched in [1] by putting a special emphasis on the conventions adopted by both programs. Moreover, our simulated dataset is made available online [11], for we deem it to constitute a valuable reference for other groups who develop MC programs.

This paper is organized as follows. In Section 2, we present the polarimetric imaging experiments simulated by both MC programs. For simplicity, our study is confined to semi-infinite samples composed of uncorrelated, optically inactive spherical polystyrene particles. Nonetheless, we distinguish two kinds of illumination beams: fully polarized ones and totally unpolarized ones. In Section 3, we proceed with the description of the software, i.e., how the polarized photon paths are generated within the two different programs. In Sections 4 and 5, we present the simulation results in the form of output images and tabulated basic statistics relative to the photon paths, before discussing the similarities/differences that can be observed, especially in the presence of an unpolarized light source.

2. Simulated Polarimetric Imaging Experiments

A. Principle

Polarimetric experiments such as the one sketched in Fig. 1(a) have a long history. Analogous setups were developed almost two centuries ago and were commonly employed to assess the quality of sugar production [12]. Nowadays, polarimetric experiments based on imaging find a large number of applications, ranging from astrophysics and x ray physics to ellipsometry and biomedical optics [13]. In the latter case, the purpose is to develop new diagnostics procedures for tissues located beneath the skin surface [14–17]. It is believed that changes in the polarization altering properties of tissues can be valuable to reveal their pathology. Standard imaging polarimetric experiments consist in illuminating a sample with a given input polarization state, ${|e〉}_i$, and then recording the spatial distribution of the backscattered light with different output polarization analyzers, $〈f|$ [see Fig. 1(c)]. The purpose is to map out the correlations between the various input and output states, so that a thorough understanding of the sample’s polarization altering properties can be achieved. The mapping of the correlations can be visualized in various forms, though the backscattering imaging Perrin–Mueller (PM) matrix has quickly been established as the primary means of visualization [18–20]. This matrix [see Fig. 1(d)], which is composed of $4\times 4$ images and which can be obtained from at least 16 independent measurements [21], constitutes not only a convenient extension for MC programs of the Stokes–Mueller type, but also a suitable way to represent synthetically the mapping [22]. This imaging PM matrix, as we define it in Eqs. (12)–(14) in [1], is experimentally obtained in three steps.

 

Fig. 1. (a) Schematic drawing of the simulated polarimetric imaging experiment designed to record the spatial distribution of the light backscattered from a sample. The sample—a suspension of polystyrene spheres in water—is illuminated by a polarized or unpolarized laser beam. (Note that in our simulations, the illumination is done normal to the sample surface, and not with an oblique incidence as suggested here.) The backscattered light goes through the detection setup which is composed of a polarization analyzer and imaging optics, e.g., a CCD array coupled with a lens. This is implemented as shown in (b): a photon path that reaches the sample’s backscattering surface within the detection area is repositioned on the detector plane at $z=0$ and is recorded by the imaging array, whose normal is given by ${\widehat{\mathbf{n}}}_D=(0,0,-1)$. Such a path is drawn in red (lighter color); the photon exits the sample with an angle, ${\theta }_{\text{out}}$, with respect to ${\widehat{\mathbf{n}}}_D$. The Cartesian lab frame is denoted $({\widehat{\mathbf{x}}}_L,{\widehat{\mathbf{y}}}_L,{\widehat{\mathbf{z}}}_L)$ and the coordinate system attached to the imager array is $({\widehat{\mathbf{x}}}_D,{\widehat{\mathbf{y}}}_D,{\widehat{\mathbf{z}}}_D)$. (c) Displays an example of a polarization sensitive intensity distribution image, $J_{〈o|i〉}$, recorded by Program B (see also Fig. 3). (d) Shows an imaging PM matrix obtained from $6\times 6J_{〈o|i〉}$ measurements.

Download Full Size | PDF

Within the scope of our comparative study, we simulated the above-mentioned polarimetric imaging experiments with both Program U, using steps 2 and 3, and Program B, using all steps from 1 to 3. Note that, although these three steps can easily be reproduced in any polarized MC program, simulations employing the Stokes–Mueller formalism tend to overlook the initial analysis procedure and record directly the output Stokes vectors of the backscattered photon paths, whereby the interpretation of the polarimetric images becomes even more difficult (see Sections 4 and 5). In the following paragraphs, we detail the simulation parameters that were used to model the examined samples, as well as the illumination/detection schemes.

B. Samples

One of the intentions of the present work is to demonstrate that MC simulations, both of Stokes–Mueller type and Jones type, are capable of treating depolarization effects. As already stated in [1], depolarization effects can be the result of two distinct kinds of randomness: the randomness of the “quantum world” and that due to thermal fluctuations. As the latter has not been included heretofore into MC simulations, we choose to focus on the former, where the depolarization manifests itself upon detection by the incoherent addition of the randomly sampled photon paths.

For this purpose, we simulated light propagation in uniformly distributed suspensions of nonabsorbing polystyrene spheres (refractive index, $n_{\text{sphere}}=1.59$) immersed in water (refractive index, $n_{\text{water}}=1.33$). We thereby ensure the easy reproducibility of our simulations by other groups, and discard the kind of randomness we are not interested in. Such spheres obey the Mie scattering model [24] and are inherently nondepolarizing, since scattering by a polystyrene sphere is a deterministic process [25] that respects, for instance, the nondepolarizing criterion proposed by Gil and Barnabeu [26].

The modeled samples are of dimensions $x_{\text{sample}}=y_{\text{sample}}=3000\mathrm{mm}$ and $z_{\text{sample}}=10\mathrm{mm}$, so that their interface in the lab’s $x_L,y_L$-plane extends indefinitely, while they maintain a finite thickness along $z_L$. These laterally infinite slabs are assumed to be surrounded by air ($n_{\text{air}}=1$) at $z=0\mathrm{mm}$ and $z=z_{\text{sample}}=10\mathrm{mm}$. We covered both the “Rayleigh” and “Mie” scattering regimes, in that we distinguished nine different samples, divided into three distinct categories according to the particles’ size parameter (see Table 1 and Fig. 2). For each category, we varied the scattering coefficient, ${\mu }_s$, between $0.1{\mathrm{mm}}^{-1}$ and $5{\mathrm{mm}}^{-1}$, in order to account for both single and multiple scattering.

 

Fig. 2. Polar plots of probability densities for the scattering angle, $\theta$ (see Section 3). (a) Samples of type §1 (radius $a=50\mathrm{nm}$, size parameter, 0.50); the dipole-like scattering profile is characteristic of Rayleigh-type scatterers, (b) samples of type §2 (radius $a=200\mathrm{nm}$ size parameter 1.98), and (c) samples of type §3 (radius $a=800\mathrm{nm}$, size parameter, 7.94); the forward scattering behavior is typical of the Mie regime.

Download Full Size | PDF

Table 1. List of the 9 Samples Modeled in the Simulations, Together with Their Respective Specifications

View Table | View all tables in this article

Precision errors are given with a 95% confidence interval (and rounded for readability), as calculated with Program B.

C. Illumination

The samples are illuminated at the middle of their surface at $z=0\mathrm{mm}$ by means of an infinitesimally narrow, collimated “pencil beam” ($\lambda =633\mathrm{nm}$) normal to that surface. Naturally, this illumination model, together with the previously described semi-infinite samples, are nothing but theoretical constructs that we choose to employ, once again, for the sake of reproducibility.

We considered $n_s=6+1$ different input polarization states (see the Table in Appendix A). The first six states $\{(|LX〉,|LY〉),(|L+〉,|L-〉),(|C-〉,|C+〉)\}$ correspond to three fundamental orthogonal pairs of fully polarized states and were used to construct the samples’ imaging PM matrix. The seventh and last input, $|O〉$, coincides with unpolarized illumination, which constitutes a case of particular interest in the frame of our comparative investigation, as the modeling of an unpolarized light source with the Jones formalism has seldom been studied.

In Program U, unpolarized illumination is produced by using the pre-averaged Stokes vector, ${\mathrm{s⃗}}_O=\left(\begin{array}{cccc}1& 0& 0& 0\end{array}\right)$. In Program B, a random ensemble is created by “mixed sampling,” which is a befitting description of natural unpolarized light. More specifically, a random initial polarization state, $|O〉$, is assigned to each photon path as follows:

D. Imaging

Here, we are interested in the intensity/polarization distribution in the backscattering geometry, as shown in Fig. 1(b): backscattered photon paths are recorded by a square-shaped CCD ($150\times 150$ pixels centered on the illumination point, with a pixel size of $0.025\times 0.025\mathrm{mm}$ and a numerical aperture, $\text{N.A.}=1$), after undergoing subsequent polarization analysis. The implementation is detailed in Section 3.D.

3. Software Description

In this section, we describe in further detail the generation/detection of polarized photon paths with both programs. The most general specifications relative to both programs are listed in Table 5 in Appendix B. As the generation of photon paths by means of polarized MC simulations has already been subject to numerous publications (e.g., [6,27,28]), we confine our description to parameters relevant to our investigation and stress the conventions employed by both programs. (One point that is not addressed here is the comparison of both programs’ computational performances.)

The prime idea of the simulations consists in constructing independent photon paths by generating, for each path, successive path segments, where each segment, $i$, is defined by: a start/end position, ${\mathrm{r⃗}}_i/{\mathrm{r⃗}}_{i+1}$, a propagation direction, ${\widehat{k}}_i$, and a polarization state, ${|e〉}_i$. Each photon path is constructed with the following steps.

The pivotal step here is the realization of the scattering event, which defines a new path segment with a new propagation direction, ${\widehat{k}}_{i+1}$, and a new corresponding polarization state, ${|e〉}_{i+1}$. This is achieved by sampling two Euler angles. $0\le \phi <2\pi$ and $0\le \theta \le \pi$, the former being the “azimuthal angle” and the latter being referred to as the “scattering angle.” However, as explained below, the conventions relative to the rotations and the sampling methods differ for both programs.

It is worthwhile mentioning that, as concise as the Fresnel formulas for reflection/refraction may be, their implementation into a polarized MC program is not a straightforward task, despite the abundant literature available (e.g., [29–32]). Nevertheless, an elaborate explanation of such a task is beyond the scope of this paper: we will only hint at the importance of remaining consistent with the used conventions.

Likewise, the two programs differ in the scope of purpose. Program U is designed especially for the investigation of light propagation in semi-infinite planar layered samples. On the other hand, Program B is a general purpose program. It combines a ray tracing algorithm with a simple construction tool, which allows to simulate light propagation in almost arbitrarily structured samples. This feature is not directly exploited in the present work, but it is extremely valuable for debugging purposes; individual elements of the algorithm, such as for example Fresnel’s interface laws, can be individually tested in simple yet nontrivial geometries.

A. Rotation Conventions

The Euler angles and the conventions with which they are used have been thoroughly described by Goldstein et al. [33]. The first angle, $\phi$, unambiguously refers to a rotation of the photon’s local system around its $z$ axis. (Note that all the rotations are performed clockwise when looking into the direction pointed by the corresponding axis-vector.) The second angle, $\theta$, can correspond to a rotation around the local $x$ axis (this is the so-called “$x$ convention,” which gives the most common definition of the scattering angle [34]) or to a rotation around the local $y$ axis (“$y$ convention”).

Both programs employ the right-hand rule (or right-handed coordinate system), yet with different conventions with respect to the definition of the scattering angle, $\theta$: Program B employs the “$x$ convention,” whereas Program U adopts the “$y$ convention.”

B. Scattering Geometry and Amplitude Scattering Matrix

Program U, like most programs of the Stokes–Mueller type, is based on the conventions defined by Bohren and Huffman [24]. This concerns not only the definition of circular polarization (not given here), but also the definition of the amplitude scattering matrix, $S_U(\theta )$:

C. Sampling of the Angles

As scattering from a polystyrene sphere implies azimuthal symmetry, the $\phi$-angle is uniformly distributed between 0 and $2\pi$, while the Mie scattering law gives the distribution of the $\theta$-angles. Both programs sample the scattering angle, $\theta$, from 7200 tabulated angle values, prepared prior to the simulation runs. More precisely, Program B samples the $\theta$-angle with the inversion technique and the $\phi$-angle with the rejection method, as described by Jaillon and Saint-Jalmes [35]. As to Program U, it samples both angles with the rejection method, analogously to [19,27,36].

D. Detection of the Photon Paths

Every time a photon path exits the sample at $z=0$ within the detection area, the pixel that has been hit needs to be incremented accordingly in the imager arrays [see Fig. 1(b)]. This entails calculating the probability that a photon gets detected after polarization filtering. As discussed in [1], the modeling of polarization filtering constitutes a delicate subject matter. Here, we adopt the ray optical approximation; that is, Program B and Program U apply the formulas of Eqs. (17) and (19) in [1], respectively.

In the Jones formalism, the polarization filtering is performed explicitly. At each simulation run (where a given initial state, ${|e〉}_i$, is assigned to all photon paths), Program B sequentially records six CCD arrays, each of which is coupled with one polarization filter (these six analyzers, $〈f|$, coincide with the three orthogonal pairs we have previously defined $\{(〈LX|,〈LY|),(〈L+|,〈L-|),(〈C-|,〈C+|)\}$).

On the other hand, Program U carries out the polarization filtering implicitly, as this is inherent to the “measurable” Stokes vector. In other terms, every time a photon path is to be detected in Program U, the four components of the associated Stokes vector are directly recorded into the four imaging arrays designed for this purpose.

4. Simulation Results Obtained with Both Programs

The remainder of this paper will concentrate on presenting and discussing the simulation results yielded by both programs. For each sample listed in Table 1 we performed two series of $n_s=6+1$ simulations, first by neglecting the refractive index mismatch between the sample and the surrounding air, second by taking these “Fresnel interfaces” into account. Each simulation run corresponds to one polarization state, ${|e〉}_i$, of the illuminating beam (see Appendix A) and generates:

Our results are presented in three consecutive steps that lead to the imaging PM matrix, as explained in Section 2.A.

A. Basic Statistics and Intensity Distributions—“Primary Data”

The first set of results is reported in Table 2. For all nine samples, we give the fraction of photon paths that get backscattered at the illumination surface (reflectance) and the fraction of photon paths that get detected by the CCD arrays. Given the geometry of our problem, where lateral infinity is assumed, these values do not depend on the input state, ${|e〉}_i$. These elementary statistics exhibit perfect agreement between both programs.

Table 2. Simulated Reflectance Values for the Different Samples ^a

View Table | View all tables in this article

In addition to the data listed in Table 2, an example of the resulting $(6+1)\times 6$ intensity distribution images, $J_{〈o|i〉}$, obtained with Program B is presented in Fig. 3. These images constitute the “primary data” that would be recorded in a real experiment. The center of each image coincides with the illumination point. Without delving too much into the reading of such backscattering polarimetric images, a few basic observations can be made. For instance, the symmetric nature of the examined sample is self-evident, or the co- (e.g., $J_{〈LX|LX〉}$) and cross-polarization (e.g., $J_{〈LX|LY〉}$) patterns are finely discernible. Another trivial observation is that the copolarization patters obtained with the $|L+〉/|L-〉$ inputs are identical to those yielded by the $|LX〉/|LY〉$ inputs, when rotated by $\pm 45°$. This elemental symmetry feature becomes even more discernible in the output imaging Stokes vectors (see the images in Section 4.B).

 

Fig. 3. “Primary data” obtained with Program B for the sample §3c (where the Fresnel interfaces were neglected); these $(6+1)\times 6$ images constitute the intensity distributions, $J_{〈o|i〉}$, as they would be recorded in a real experiment. Each column, $i$, corresponds to one of the $n_s=6+1$ input states, ${|e〉}_i$; each row, $o$, corresponds to a polarization analyzer, $〈f|$. $J_{〈o|i〉}^a=J_{〈o|i〉}/\text{total}(J_{〈o|i〉})$, “max” stands for the overall maximum value, and $\gamma =3.3$. Images such as $J_{〈LX|LX〉}$ and $J_{〈LX|LY〉}$ reveal co- and cross-polarization properties, respectively.

Download Full Size | PDF

B. Output Imaging Stokes Vectors

Following the intensity distributions, $J_{〈o|i〉}$, we present the output Stokes vectors, ${\mathbf{S}}_i^{\text{out}}$, obtained with both programs. However, rather than showing the entirety of our results (which can be accessed online [11]), we display a selection of relevant examples: Two distinct series of output Stokes vectors calculated from samples §1c and §3c are given in Figs. 4 and 5, respectively. It should be pointed out that we chose not to present results obtained for samples where the interfaces of Fresnel type were included, since, in our examples, the inclusion of index mismatches does not affect the qualitative patterns observed in the various elements.

 

Fig. 4. Top: $n_s=6+1$ output Stokes vectors calculated with Program B for the highly scattering sample §1c, composed of smaller polystyrene spheres. Like in Fig. 3, each column, $i$, corresponds to one of the input states, ${|e〉}_i$. The first elements, which represent intensities, are displayed in a similar fashion to the distributions of Fig. 3 (thus, $J_{〈o|i〉}^b$ here is analogous to $J_{〈o|i〉}^a$). Bottom: Same series of ${\mathbf{S}}_i^{\text{out}}$ obtained with Program U. The helicity flip that occurs near the illumination point is visible in the $V^{\text{out}}$ elements of the vectors. ${\mathbf{S}}_{C-}^{\text{out}}$ and ${\mathbf{S}}_{C+}^{\text{out}}$. Note also that the $U^{\text{out}}$ elements of the ${\mathbf{S}}_{L+}^{\text{out}}/{\mathbf{S}}_{L-}^{\text{out}}$ vectors are identical to the $Q^{\text{out}}$ elements of the ${\mathbf{S}}_{LX}^{\text{out}}/{\mathbf{S}}_{LY}^{\text{out}}$ vectors, when rotated by $\pm 45°$. This is a basic symmetry feature, useful for an assessment of the consistency of the data evaluation.

Download Full Size | PDF

 

Fig. 5. Top: $n_s=6+1$ output Stokes vectors calculated with Program B from the intensity distributions of Fig. 3, i.e., for the highly scattering sample §3c, composed of larger polystyrene spheres. Bottom: same series of ${\mathbf{S}}_i^{\text{out}}$ obtained with Program U. Stokes’ equivalence theorem [40] can be verified, e.g., ${\mathbf{S}}_O^{\text{out}}={\mathbf{S}}_{LX}^{\text{out}}+{\mathbf{S}}_{LY}^{\text{out}}$.

Download Full Size | PDF

For each sample, we generated a set of $n_s=6+1$ output Stokes vectors, each of which corresponds to an illumination polarization state, ${|e〉}_i$. These ${\mathbf{S}}_i^{\text{out}}$ vectors lie at the forefront of our study as the analysis of these images allows not only for the qualitative and quantitative comparison of both programs’ results, but also for the verification of the proper use of the conventions on both sides. This will be further developed in the discussion; nonetheless, an initial examination of the images establishes the perfect qualitative agreement between the results yielded by both programs, for all of the input states. An immediate conclusion that arises from the viewing the output Stokes vectors generated for the three different orthogonal pairs is again the symmetric nature of our sample models.

Note that the visual representation of these images is nontrivial; the data needs to be presented in an intelligible yet relevant manner, so as to accentuate the physical significance of the Stokes vectors’ four components. To this end, we shall once more refer to the definitions given in [1]. The first element, $I^{\text{out}}$, designates the “total intensity” of the detected light after passing through the analyzers. For this reason, the elements, $I^{\text{out}}$, are displayed analogously to the intensity distributions, $J_{〈o|i〉}$. The three remaining elements, $Q^{\text{out}}$, $U^{\text{out}}$, and $V^{\text{out}}$, represent polarization distributions and are normalized by $I^{\text{out}}$ (some authors might even assimilate these distributions to “degrees of polarization” [37]).

Here, the color coding is such that, in each image, red and blue designate two orthogonal states. Purple stands for depolarized regions, where an even and incoherent mixture of both states is present. For example, in the elements, $Q^{\text{out}}$, regions colored in red indicate the prominence of $|LX〉$ states, whereas blue is assigned to regions where $|LY〉$ states are prevalent. Similarly, in the $V^{\text{out}}$ elements, red and blue correspond to regions where $|C+〉$ and $|C-〉$ states are dominant, respectively.

Thus, several assertions regarding the $V^{\text{out}}$ elements can be made. Looking at the vectors, ${\mathbf{S}}_{C-}^{\text{out}}$ and ${\mathbf{S}}_{C+}^{\text{out}}$, we can observe the helicity flip that occurs in the direct vicinity of the illumination point by single backscattering (this is clearly visible in Fig. 4, but might be difficult to observe in Fig. 5; see Fig. 6 for a magnified version of the elements). Yet, further from the illumination point, Figs. 4 and 5 display different behaviors. In the sample §1c, the multiple scattering by smaller spheres leads to depolarization, while in the sample §3c, the larger spheres tend to preserve the helicity through multiple forward scattering [38,39]. Another peculiar remark concerns the $V^{\text{out}}$ elements obtained with linear illumination states: the polystyrene spheres in the sample §3c have a size parameter large enough to distinctly generate circular states from the scattering of linear ones, which is not the case for the sample §1c.

 

Fig. 6. We show here magnified versions of the $V^{\text{out}}$ elements from the vectors ${\mathbf{S}}_{C-}^{\text{out}}$ (left) and ${\mathbf{S}}_{C+}^{\text{out}}$ (center) in Fig. 5, as well as the $m_{44}$ element from the imaging PM matrix at the right column of Fig. 7. These images were obtained with the strongly scattering sample §3c, composed of larger polystyrene spheres. The center of the images coincide with the light detected after single backscattering, i.e., with a helicity flip. The regions far from the center correspond to the light detected after multiply scattering; the forward scattering behavior of the spheres preserves the helicity of the input state.

Download Full Size | PDF

Furthermore, it can be verified that the output Stokes vector for unpolarized input, ${\mathbf{S}}_O^{\text{out}}$, respects Stokes’ equivalence theorem [40], insofar as it is equal to the sum of one of the orthogonal pairs, e.g., ${\mathbf{S}}_O^{\text{out}}={\mathbf{S}}_{LX}^{\text{out}}+{\mathbf{S}}_{LY}^{\text{out}}$.

C. Sample’s Backscattering Imaging Perrin–Mueller Matrix

In this subsequent and final step, we expose the imaging PM matrices calculated from the series of output Stokes vectors shown previously (see Fig. 7). Even though these imaging PM matrices do not provide additional information, per se, with respect to the comparison of both programs, it might be of some interest to dwell upon them.

 

Fig. 7. Two sets of imaging PM matrices obtained with Program B (top) and Program U (bottom). The matrices on the left were calculated from the Stokes vectors of Fig. 4, while the matrices on the right were calculated from the data in Fig. 5. The first element, $m_{11}$, is distinguished from the others; it represents the same intensity distribution as the element $I_{\text{out}}$ of $S_O^{\text{out}}$ (and $J_{〈o|i〉}^c$ here is analogous to $J_{〈o|i〉}^a$ and $J_{〈o|i〉}^b$). The remaining elements, $m_{ij}$, are normalized with respect to $m_{11}$, so as to represent polarization distributions.

Download Full Size | PDF

One especially important point concerns the symmetrical relations within the PM matrix, as it is difficult to find concordance on the (anti)symmetry considerations in the literature. In accordance with our sample properties, these relations can be retrieved from Perrin’s or van de Hulst’s works [41,42]: $m_{21}=m_{12}$, $m_{31}=-m_{13}$, $m_{32}=-m_{23}$, $m_{42}=m_{24}$, $m_{41}=m_{14}$, and $m_{43}=-m_{34}$ (where $m_{ij}$ corresponds to the element of the PM matrix at the row $i$ and column $j$).

Analogously to the output Stokes vectors, the representation of the imaging PM matrix plays a significant role. The very first element, $m_{11}$, represents an intensity distribution, akin to the $I^{\text{out}}$ element of the output Stokes vector obtained with unpolarized input. Actually, it is evident from Eq. (5) that the first column, ${\mathrm{M⃗}}_1$, is nothing but an output Stokes vector that coincides with ${\mathbf{S}}_O^{\text{out}}$. The first row of the PM matrix equally seems to carry a particular meaning, albeit one that is difficult to interpret. Besides, the last element, $m_{44}$, is rather revealing, as it can yield information on the size of the scatterers constitutive of the sample, as explained in Section 4.B and shown in Fig. 6 [43].

More generally, the decoding of the imaging PM matrix is an arduous task, and has been subject to numerous investigations (e.g., [44–47]). There is perhaps a legitimate need to introduce alternative and more judicious encodings of the data collected from experiments; likewise, the primary data, $J_{〈o|i〉}$, and output Stokes vectors, ${\mathbf{S}}_i^{\text{out}}$, can be further exploited when characterizing samples. (Actually, the output Stokes vectors and the imaging PM matrix are nothing but two different visualizations of the information contained in such primary data, albeit with a different and compressed encoding.)

5. Discussion

A. Unpolarized Illumination: Mixed Sampling versus the Pre-averaged Stokes Vector

Our early observations established the perfect qualitative agreement between both programs’ results, this for the entire set of numerical experiments we have conducted. Closer inspection reveals that in the case where the samples are illuminated with pure polarization states, Programs B and U even exhibit the same statistics: This demonstrates, as expected, both formalisms’ equal ability to treat not only polarized light propagation, but also depolarization engendered by the incoherent summation of randomly sampled photon paths.

Howbeit, second order statistical differences appear when unpolarized illumination is simulated, particularly in the $V^{\text{out}}$ element of the corresponding output Stokes vector $S_O^{\text{out}}$. To illustrate this in the case of the samples §3a and §3c, we plotted in Fig. 8 cross sections of the $V^{\text{out}}$ elements obtained with Programs B and U for linear, $|LX〉$ (a), circular right, $|C-〉$ (b), and unpolarized $|O〉$ (c) inputs, respectively. For instance, for the sample §3a, it is apparent that despite both datasets displaying comparable mean values, the standard deviation ${\sigma }_B$ of the data from Program B is 60 times higher than that of Program U, ${\sigma }_U$. The evolution of the ratio, ${\sigma }_B/{\sigma }_U$, with different types of samples is shown in Table 3 for the four elements $(I_{\text{out}},Q_{\text{out}},U_{\text{out}},V_{\text{out}})$ of ${\mathbf{S}}_O^{\text{out}}$.

 

Fig. 8. Cross sections of the $V^{\text{out}}$ elements obtained with Programs B (in black) and U (in blue), for the samples §3a (left column) and §3c (right column). The three different rows correspond from top to bottom to (a) linear, $|LX〉$, (b) circular right, $|C-〉$, and (c) unpolarized, $|O〉$, inputs. While the fluctuations are alike for inputs that are pure polarization states, second order statistical differences are visible for unpolarized input (third row), especially in the single scattering regime (left column).

Download Full Size | PDF

Table 3. Ratio ${\mathsf{\sigma }}_{\mathsf{B}}/{\mathsf{\sigma }}_{\mathsf{U}}$ Given for All Elements of ${\mathsf{S}}_{\mathsf{O}}^{\mathsf{\text{out}}}$, Where the Input ${\mathsf{s⃗}}_{\mathsf{O}}=\mathsf{(}\mathsf{1}\mathsf{0}\mathsf{0}\mathsf{0})$ in Program U

View Table | View all tables in this article

The following essential conclusions can be made:

It is noteworthy that the differences predominantly affect the element $V^{\text{out}}$, which is associated with circular polarization. In Program B, the fluctuations in $V^{\text{out}}$ suggest an even mixture of right-handed and left-handed states that cancel each other out and lead to a mean value of zero, whereas in Program U, the statistics of the same $V^{\text{out}}$ element expose the quasi nonexistence of these circular states. This is explained by the shape of the PM matrix ${\mathbf{M}}_{\mathrm{iso}}$ of a spherical and isotropic scatterer [24,41]:

Multiplying ${\mathrm{s⃗}}_O=(\begin{array}{cccc}1& 0& 0& 0\end{array})$ by this matrix cannot produce circular states. As such, the Stokes–Mueller formalism favours linear states upon the single scattering of the pre-averaged unpolarized input state. The smaller the scatterers (i.e., dipole like), the more this phenomenon is accentuated, since $m_{34}=0$ [48]. The larger the scatterers (i.e., in the Mie regime), the more likely they are to create circular states, as ascertained in Fig. 5 (see Section 4.B for the explanations), and the smaller the differences become.

Further, it can be seen in Table 3 that with larger spheres, where the forward scattering preserves the circular states (see Section 4.B), the second order statistical differences between the two programs also involves the output Stokes vector’s $Q_{\text{out}}$ and $U_{\text{out}}$ elements. These differences are emphasized when Fresnel interfaces are taken into account within the simulations.

We verified that all the differences between both programs disappear upon implementation of mixed sampling in lieu of ${\mathrm{s⃗}}_O$ in Program U. Assuredly, mixed-sampling is a better suited model to describe natural unpolarized light; and fluctuations similar to those observed in Program B are bound to happen in nature. However, such fluctuations would be quite difficult to measure in a real experiment and in essence, the differences in the second order statistics can be acknowledged as being inconsequential in most cases (with perhaps some exceptions, e.g., in optically active samples).

Yet, we considered the origin of these differences to be an intriguing question that lies beyond a mere academic consideration. Hence, in order to better fathom the differences between both formalisms, we pursued our investigation by comparing the results obtained with mixed sampling and orthogonal sampling when modeling the unpolarized light source [1].

B. Unpolarized Illumination: Mixed Sampling versus Orthogonal Sampling

As already explained in [1], an alternative to the mixed sampling for the modeling of an unpolarized light source would be to apply Stokes’ equivalence theorem, i.e., to sample the photon paths’ initial polarization state from one orthogonal pair. Thus, we repeated our simulations with Program B by replacing the mixed sampling with orthogonal sampling from the following two pairs: linear $(|LX〉,|LY〉)$ and circular $(|C-〉,|C+〉)$. We compared the newly obtained output Stokes vectors, $S_O^{\text{out}}$, with the ones yielded by mixed sampling. The different sampling methods yield quasi-identical results for the $I^{\text{out}},Q^{\text{out}}$, and $U^{\text{out}}$ elements. Yet, peculiarly, the $V^{\text{out}}$ elements are anew affected by noticeable second-order statistical differences. Such differences are illustrated in Fig. 9(a), where we plotted cross sections of the said elements for the sample §3a. Possibly, the elements, $Q^{\text{out}}$, and $U^{\text{out}}$, generally remain unaffected because they correspond to linear states that can be expressed with the circular pair, in which case differentiating the linear/circular orthogonal samplings from the mixed one becomes irrelevant. The opposite does not hold true, i.e., circular states cannot be expressed with the linear pairs; therefore, the differences are more visible in the $V^{\text{out}}$ element.

 

Fig. 9. Cross sections of the $S_O^{\text{out}}$ vector’s $V^{\text{out}}$ element obtained with the sample §3a. Each curve corresponds to a different sampling method used for the modeling of the unpolarized light source: mixed sampling in black, orthogonal sampling from the linear pair $(|LX〉,|LY〉)$ in orange, and orthogonal sampling from the circular pair $(|C-〉,|C+〉)$ in green. In (a) the “classical” detection method is employed (as detailed in Section 3.D), whereas in (b) the modified detection scheme (where photons paths do not simultaneously pass through the different analyzers) was adopted.

Download Full Size | PDF

The statistics of the $V^{\text{out}}$ element generated with sampling from the linear pair are similar to those obtained from ${\mathrm{s⃗}}_O$ in Program U. Au contraire, the statistics which correspond to sampling from the circular pair follow those obtained from mixed sampling. Our observations matched those made when comparing mixed sampling to the pre-averaged Stokes vector, ${\mathrm{s⃗}}_O$, in Section 5.A. That is, the differences are most noticeable for samples where the backscattered light is detected only after few scattering events. Yet, before deliberating further over the meaning of these findings, we shall report a more striking observation: in the case of the sample §3a, even the recorded patterns are dissimilar on the $V^{\text{out}}$ element for the various sampling techniques (see Fig. 10).

 

Fig. 10. $S_O^{\text{out}}$ vector’s $V^{\text{out}}$ element obtained with the sample §3a. Each image corresponds to a different sampling method used for the modeling of the unpolarized light source. Left: orthogonal sampling from the linear pair $(|LX〉,|LY〉)$. Center: orthogonal sampling from the circular pair $(|C-〉,|C+〉)$. Right: mixed sampling. The “classical” detection method (as detailed in Section 3.D) leads to correlations in the resulting output Stokes vectors’ $V^{\text{out}}$ elements, as can be seen with the difference in patterns.

Download Full Size | PDF

After some consideration, we came to the conclusion that these strongly contrasting behaviors in statistics and patterns originate from a fundamental assumption, which is part of the detection procedure and that has not been mentioned so far. As explained in Section 3.D, each photon path that reaches the detector is simultaneously sent, as it were, through the different analyzers (again, this is carried out explicitly in Program B and implicitly in Program U). In other words, the same photon path delivers information through the different polarization analyzers: certainly, this is rather unrealistic (“a photon cannot be divided into parts!”) and leads to correlations in the output Stokes vector elements.

These correlations are especially pronounced in the case of photon paths detected after single scattering, and mainly concern the $V^{\text{out}}$ element. While such a simplification of the detection process would be inconceivable in the quantum world, it does not bear profound consequences in polarized MC programs, which might explain the lack of discussions on the topic.

Our hypothesis was corroborated by once again repeating the simulations with Program B, this time by sending each detected photon path through a single analyzer. Analogous to the treatment of reflection/refraction from a Fresnel interface, where a random number is generated to determine whether the photon path is reflected or refracted, we resorted for each detected path to random numbers in order to determine which analyzer is to be employed amidst each pair.

Our subsequent results have demonstrated that the differences in patterns vanish and that the second order statistical differences disappear. As an example, we show in Fig. 9(b) the same $V^{\text{out}}$ element cross-sections as in Fig. 9(a), but obtained with the modified detection scheme, again for the sample §3a. The standard deviation values ($\sigma \simeq 0.71$) are quasi-identical for all three curves.

We can now relate the conclusions above to the Stokes–Mueller formalism, and concede the following.

A final question entitled to be answered is whether the sampling schemes differ in terms of convergence. For this purpose, we have repeated the simulations by illuminating the sample §3c with the different sampling methods. By keeping the modified detection scheme, we performed successive simulation runs, where for each run we gradually increased the number of simulated photon paths and recorded the mean value over a small area ($10\times 10$ pixels) near the center of the obtained $V^{\text{out}}$ element. The results are shown in Fig. 11. We recorded for each model the resulting root-mean-square (RMS) value in regions where enough statistics were achieved (i.e., with a number of photon paths exceeding 40 M). With the mixed sampling, the orthogonal sampling from the linear pair $(|LX〉,|LY〉)$, and the orthogonal sampling from the circular pair $(|C-〉,|C+〉)$ the RMS values are 0.024, 0.022, and 0.018, respectively. The results for the mixed sampling differ from that of the orthogonal sampling from the circular pair by almost a factor of $\sqrt{2}$.

 

Fig. 11. Convergence of the mean value over a small area ($10\times 10$ pixels) near the center of the $S_O^{\text{out}}$ vector’s $V^{\text{out}}$ element, obtained with the sample §3c. Each curve corresponds to a different sampling method used for the modeling of the unpolarized light source: mixed sampling in black, orthogonal sampling from the linear pair $(|LX〉,|LY〉)$ in orange, and orthogonal sampling from the circular pair $(|C-〉,|C+〉)$ in green. The inset displays a magnified view of the curves in the region where the number of simulated photon paths exceeds 40 M.

Download Full Size | PDF

6. Conclusion

By comparing the simulation results yielded by Programs B and U, we addressed the two points raised in the introduction, namely, the validation of both programs and the elucidation of the ubiquitous Stokes–Jones duality. It is remarkable that the two programs produce identical results and exhibit the same statistics, this despite having dissimilar algorithms. The set of data generated by this comparative study will hopefully constitute a valuable asset to ascertain the trustworthiness of other polarized MC programs, in a quantitative yet simple fashion. Forasmuch as we omitted comparing the computational efficiency of both programs, measuring the respective performances of the two formalisms would be an interesting subject of further investigation.

During the course of our work, we have experienced first hand the critical role of the conventions involved in the modeling of polarimetric experiments. These conventions, whose weight cannot be stressed enough, should be unequivocally defined, not only for the easy comparison of the data, but also for the proper interpretation of the results. This is best asserted by Holm: “Almost everyone, even old-hands, has experienced that frustrating feeling of trying to reconcile differences caused by conventions. However, most authors (for whatever reasons) sidestep the convention issue completely” [9].

Likewise, we have demonstrated that Jones’ original writings, which state that his formalism cannot treat depolarization effects, become obsolete in the context of polarized MC simulations. The stochastic nature of these simulations allows both the Jones and the Stokes–Mueller formalisms for the accurate modeling of depolarization effects. We have purposefully avoided engaging in depolarization effects that arise from thermal fluctuations. Nevertheless, the latter can be modeled with both formalisms, similarly to unpolarized light sources. Depending on the simulations’ intent and the nature and geometry of the examined samples, the appropriate method (pre-averaging or mixed-/orthogonal-sampling) can judiciously be chosen, bearing in mind that the convergence and second order statistics of the results may differ from one method to another.

Meanwhile, we shall reiterate the findings of the first paper. On the one hand, compared with the Stokes vector, the Jones vector forms a relatively tangible physical entity that can be manipulated in a rather intuitive manner, thereby reducing the risk of implementation errors in a MC program. On the other hand, the polarization analysis that intervenes in the imaging process should not be overlooked. This constitutes a central step, where the approximations employed should be distinctly specified, as they might affect the outcome of the simulations, especially in the case of single scattering.

Last, but not least, we anticipate the accomplishment of more sophisticated numerical experiments (with intricate geometrical structures, optically active and/or anisotropic materials, etc.) based on the premise of reliable simulation programs. The examination of light polarization in general and circular states in particular are bound to unveil some of nature’s subtle phenomena [50].

Appendix A: List of Input States Used in the Simulations: Table 4

In Table 4 we use seven different input states in the polarimetric experiment simulations: three fundamental orthogonal pairs used to construct the sample’s imaging PM matrix and unpolarized input to compare the propagation of unpolarized light with both types of programs. We re-employ the definitions and conventions introduced in [1]; particular attention should be paid to the definition of the circular states.

Table 4. Jones and Stokes Vectors versus Polarization State

View Table | View all tables in this article

Appendix B: List of Input States Used in the Simulations: Table 5

Table 5. Specifications Relative to Both Programs

View Table | View all tables in this article

This research has been supported by the Swiss National Science Foundation, Grant Nos. PBBEP2_142142 and 205320_144432/1. Some of the authors would like to thank Joachim Wiest for his valuable comments and help with respect to the implementation of Program U.

References and Notes

1. H. G. Akarçay, A. Hohmann, A. Kienle, M. Frenz, and J. Rička, “Monte Carlo modeling of polarized light propagation. Part I. Stokes versus Jones,” Appl. Opt.53, 7576–7585 (2014).

2. F. Voit, A. Hohmann, J. Schäfer, and A. Kienle, “Multiple scattering of polarized light: comparison of Maxwell theory and radiative transfer theory,” J. Biomed. Opt 17, 045003 (2012). [CrossRef]  

3. J. Rička and M. Frenz, “Polarized light: electrodynamic fundamentals,” in Optical-Thermal Response of Laser-Irradiated Tissue, A. J. Welch and M. J. C. van Gemert, eds., 2nd ed. (Springer, 2011), Chap. 4.

4. J. Rička and M. Frenz, “From electrodynamics to Monte Carlo simulations,” in Optical-Thermal Response of Laser-Irradiated Tissue, A. J. Welch and M. J. C. van Gemert, eds., 2nd ed. (Springer, 2011), Chap. 7.

5. H. G. Akarçay and J. Rička, “Simulating light propagation: towards realistic tissue models,” Proc. SPIE 8088, 80880K (2011). [CrossRef]  

6. H. G. Akarçay, “Polarized light propagation in biological tissue: towards realistic modeling,” Ph.D. dissertation (University of Bern, 2011), http://www.iap.unibe.ch/publications/pub-detail.php?lang=en&id=3706.

7. H. G. Akarçay, J. Rička, and M. Frenz are preparing a manuscript to be called “jaMCp3: towards the realistic modeling of light propagation in biological tissues.”

8. D. Côté and 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). [CrossRef]  

9. R. T. Holm, “Convention confusions,” in Handbook of Optical Constants of Solids: Index, E. D. Palik, ed. (Academic, 1998), Vol. 3, Chap. 2.

10. R. H. Muller, “Definitions and conventions in ellipsometry,” Surf. Sci. 16, 14–33 (1969). [CrossRef]  

11. Our simulated dataset is available on the following website (within the “Polarized light propagation in biological tissue” project): www.iapbp.unibe.ch/content.php/home/projects/.

12. G. W. Rolfe, The Polariscope in the Chemical Laboratory: An Introduction to Polarimetry and Related Methods (Macmillan, 1919).

13. R. M. A. Azzam, “A perspective on ellipsometry,” Surf. Sci. 56, 6–18 (1976). [CrossRef]  

14. A. H. Hielscher, J. R. Mourant, and I. J. Bigio, “Influence of particle size and concentration on the diffuse backscattering of polarized light from tissue phantoms and biological cell suspensions,” Appl. Opt. 36, 125–135 (1997). [CrossRef]  

15. A. H. Hielscher, A. A. Eick, J. R. Mourant, D. Shen, J. P. Freyer, and I. J. Bigio, “Diffuse backscattering Mueller matrices of highly scattering media,” Opt. Express 1, 441–453 (1997). [CrossRef]  

16. R. S. Gurjar, V. Backman, L. T. Perelman, I. Georgakoudi, K. Badizadegan, I. Itzkan, R. R. Dasari, and M. S. Feld, “Imaging human epithelial properties with polarized light-scattering spectroscopy,” Nat. Med. 7, 1245–1248 (2001). [CrossRef]  

17. S. L. Jacques, J. C. Ramella-Roman, and K. Lee, “Imaging skin pathology with polarized light,” J. Biomed. Opt. 7, 329–340 (2002). [CrossRef]  

18. 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, 3399–3408 (1999). [CrossRef]  

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

20. B. D. Cameron, M. J. Raković, M. Mehrübeoğlu, G. W. Kattawar, S. Rastegar, L. V. Wang, and G. L. Coté, “Measurement and calculation of the two-dimensional backscattering Mueller matrix of a turbid medium,” Opt. Lett. 23, 485–487 (1998). [CrossRef]  

21. W. S. Bickel and W. M. Bailey, “Stokes vectors, Mueller matrices, and polarized light scattering,” Am. J. Phys. 53, 468–478 (1985). [CrossRef]  

22. V. Maxia, “Light polarization problems,” Appl. Opt. 15, 2576–2578 (1976). [CrossRef]  

23. As stated above, $n_s=4$ would actually be sufficient to construct a sample’s PM matrix. Here, we prefer to solve an overdetermined system primarily for didactic purposes, thereby remaining consistent with respect to our PM matrix definition given in [1].

24. C. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley, 1998).

25. R. Simon, “Nondepolarizing systems and degree of polarization,” Opt. Commun. 77, 349–354 (1990). [CrossRef]  

26. J. J. Gil and E. Bernabeu, “A depolarization criterion in Mueller matrices,” Opt. Acta 32, 259–261 (1985). [CrossRef]  

27. J. Ramella-Roman, S. Prahl, and S. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part I,” Opt. Express 13, 4420–4438 (2005). [CrossRef]  

28. S. L. Jacques, “Monte Carlo modeling of light transport in tissue (steady state and time of flight),” in Optical-Thermal Response of Laser-Irradiated Tissue, A. J. Welch and M. J. C. van Gemert, eds., 2nd ed. (Springer, 2011), Chap. 5.

29. E. Collett, “Mueller-Stokes matrix formulation of Fresnel’s equations,” Am. J. Phys. 39, 517–528 (1971). [CrossRef]  

30. E. Landi Degl’Innocenti, “The physics of polarization,” in Astrophysical Spectropolarimetry, J. Trujillo-Bueno, F. Moreno-Insertis, and F. Sánchez Martinez, eds. (Cambridge University, 2000).

31. J. M. Bennett, “Polarization,” in Handbook of Optics, M. Bass, ed. (McGraw-Hill, 1995), Vol. 1, Chap. 2.5.

32. D. H. Goldstein, “Polarized Light: Second Edition, Revised and Expanded” (CRC Press, 2003).

33. H. Goldstein, C. P. Poole, and J. L. Safko, Classical Mechanics, 3rd ed. (Addison Wesley, 2001), Chap. 4, Appendix A.

34. http://mathworld.wolfram.com/EulerAngles.html.

35. 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). [CrossRef]  

36. I. L. Maksimova, S. V. Romanov, and V. F. Izotova, “The effect of multiple scattering in disperse media on polarization characteristics of scattered light,” Opt. Spectra 92, 915–923 (2002). [CrossRef]  

37. A. Al-Qasimi, O. Korotkova, D. James, and E. Wolf, “Definitions of the degree of polarization of a light beam,” Opt. Lett. 32, 1015–1016 (2007). [CrossRef]  

38. F. C. MacKintosh, J. X. Zhu, D. J. Pine, and D. A. Weitz, “Polarization memory of multiply scattered light,” Phys. Rev. B 40, 9342–9345 (1989). [CrossRef]  

39. D. Bicout, C. Brosseau, A. S. Martinez, and J. M. Schmitt, “Depolarization of multiply scattered waves by spherical diffusers: influence of the size parameter,” Phys. Rev. E 49, 1767–1770 (1994). [CrossRef]  

40. G. G. Stokes, “On the composition and resolution of streams of polarized light from different sources,” Trans. Cambridge Philos. Soc. 9, 399–416 (1852).

41. F. Perrin, “Polarization of light scattered by isotropic opalescent media,” J. Chem. Phys. 10, 415–427 (1942). [CrossRef]  

42. H. C. van de Hulst, “Light Scattering by Small Particles” (Wiley, 1957).

43. C. Schwartz, “Probing random media with singular waves,” Ph.D. dissertation (University of Central Florida, 2006).

44. M. H. Smith, “Interpreting Mueller matrix images of tissues,” in The International Symposium on Biomedical Optics (International Society for Optics and Photonics, 2001).

45. N. Ghosh, I. A. Vitkin, and M. F. G. Wood, “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, 044036 (2008). [CrossRef]  

46. M. R. Antonelli, A. Pierangelo, T. Novikova, P. Validire, A. Benali, B. Gayet, and A. De Martino, “Mueller matrix imaging of human colon tissue for cancer diagnostics: how Monte Carlo modeling can help in the interpretation of experimental data,” Opt. Express 18, 10200–10208 (2010). [CrossRef]  

47. P. Sun, X. Cao, H. Sun, M. Sun, and M. He, “Spatial pattern characterization of linear polarization-sensitive backscattering Mueller matrix elements of human serum albumin sphere suspension,” J. Biol. Phys. 39, 501–514 (2013). [CrossRef]  

48. S.-M. F. Nee and T.-W. Nee, “Polarization of dipole scattering by isotropic medium,” Proc. SPIE 7065, 70650P (2008). [CrossRef]  

49. Perhaps, it would not be preposterous to claim that the “phenomenological” character attributed to the Stokes vector is an appropriate one, as linear states seem to be preponderant in nature, due not only to the passing of the sunlight through the atmosphere, but also to the diverse reflections that occur on various surfaces.

50. G. W. Kattawar, “A search for circular polarization in nature,” Opt. Photon. News 5(9), 42–43 (1994). [CrossRef]  

51. Random number generators in C++ intended for Monte Carlo applications and released under the Gnu general public license can be found at: http://www.agner.org/random/.

52. W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes in Pascal: The Art of Scientific Computing, 1st ed. (Cambridge University, 1991), Chap. 7.

53. D. E. Knuth, The Art of Computer Programming Vol. 2: Seminumerical Algorithms, 3rd ed. (Addison-Wesley, 1997), Chap. 3.