2.1.

Monte Carlo Model

A Monte Carlo program for simulation of polarized light propagation in scattering media was developed based on the radiative transfer theory.^{10,13,15–18,20,21,28} The Monte Carlo method simulates a multitude of random light propagation paths through the medium omitting the wave character of light, nonetheless considering polarization effects. The scatterers were assumed to be randomly distributed within a simulation volume of $V=10\times 10\times 10\mu {\mathrm{m}}^3$ and to have an angularly-resolved scattering probability (phase function) similar to the phase function of a single sphere.

This phase function was calculated using Mie theory (single homogeneous sphere in a homogeneous medium, see Bohren and Huffman³² and Mie³³) and used in case of a scattering event. For polarized light, this phase function depends on both the scattering angle as well as on the polarization state of the incident light. The optical properties of polystyrene microspheres (radius $r=1\mu \mathrm{m}$) in water were used as input parameters for the spherical scatterers in the Mie program: the refractive index of the spheres was $n=1.59$ in an infinite medium of $n_m=1.33$. There was no difference between the refractive index of the medium within the simulation volume and the surrounding medium. The vacuum wavelength of the light was set to $\lambda =600\mathrm{nm}$. This results in a size parameter $x=2\pi r{n}_m/\lambda \approx 13.928$.

The scattering efficiency $Q_{\mathrm{sca}}=C_{\mathrm{sca}}/(\pi r^2)$ of a single sphere resulting from the Mie program can be applied to calculate the scattering coefficient ${\mu }_s$ as in Eq. (2) used for the Monte Carlo simulations

A commonly used technique to gain the Müller matrix elements of a scatterer system via the Monte Carlo method is the subsequent execution of four simulations applying an incident Stokes vector ${\mathbf{S}}_i^{st}$ of ${(1,1,0,0)}^T$, ${(1,-1,0,0)}^T$, ${(1,0,1,0)}^T$ and ${(1,0,0,1)}^T$ which is labeled by H, V, P and R, respectively.^6,20 The exit (scattering) angle of each photon which left the simulation volume was registered for each simulation run together with its corresponding Stokes vector. The results are four scattered Stokes vectors ${\mathbf{S}}_f^{\mathrm{st}}$: ${(I_H,Q_H,U_H,V_H)}^T$, $(I_V,Q_V,U_V,V_V{)}^T$, ${(I_P,Q_P,U_P,V_P)}^T$ and ${(I_R,Q_R,U_R,V_R)}^T$. These four Stokes vectors, resulting from the average of many photon paths, can be used to calculate the Müller matrix of the system:

2.2.

Maxwell Theory

A Maxwell solution for the scattering problem is given by the Generalized Multiparticle Mie approach, an extension of Mie theory, which is able to handle multiple spherical scatterers and, as an exact solution, also takes into account dependent scattering effects. A Fortran code is provided by Xu et al.³⁴ and used as a reference Maxwell solver. Therein the scattering of an incident plane wave by a distinct arrangement of spheres is calculated by employing the wave equation as opposed to the statistical Monte Carlo model dealing with photon paths.

The total scattered field results from coherent superposition of the scattered field contributions ${\mathbf{E}}_{\mathrm{sca}}^j$ of each sphere $j$:

Since multiparticle Mie theory requires definitely positioned scatterers of finite extent for solution of the underlying boundary conditions, a collection of spheres is randomly distributed in a given simulation volume in order to achieve defined volume concentrations. To approximate homogeneous turbid media models, several random configurations of finite spheres (realizations) are generated where the results are averaged over these realizations. The spheres are randomly distributed under the restriction that no overlapping is allowed. For generation of these distributions a Metropolis shuffling algorithm^36,37 was programmed where each sphere in an initial periodic arrangement (simple cubic lattice) is randomly shuffled many times consi-dering periodic boundary conditions with respect to the cubic volume. This algorithm converges to an equilibrium distribution when performing a high number of displacement steps.

In order to improve convergence of the algorithm, the random step size is adjusted dynamically. An initial displacement step size is randomly chosen within a size interval in the order of the sphere radius. After shuffling each sphere for a predefined number of times, an acceptance rate of allowed (non-overlapping) with respect to tried sphere displacements is evaluated. In case the new position would lead to sphere overlap the tried displacement is rejected and the old position is kept.³⁷ A small acceptance rate indicates a poor probability to find allowed positions, then the step size is decreased for the next shuffling pass. On the other hand, a high acceptance rate indicates too small deviation from the initial arrangement, in that case the step size is increased for the next shuffling pass. Generation of sufficiently randomized sphere distributions is achieved if the acceptance rate stays around 50%.

For the Maxwell simulations, the same parameters as for the Monte Carlo simulations were used, i.e., the different volume concentrations are modeled by $N=12,24,48$ and 96 spheres ($r=1\mu \mathrm{m}$), respectively, in a simulation volume of $10\times 10\times 10\mu {\mathrm{m}}^3$. This volume size allows simulations to be performed in a reasonable amount of time even at higher volume concentrations. Apart from diffraction around forward scattering angles, the characteristics of an infinite scattering medium are reproduced in this subset simulation volume. For the four concentrations 10 distinct random realizations were generated where each sphere is shuffled ca. 10000 times. From these realizations the resulting Müller matrix elements are averaged to reduce interference effects. As above, the incident vacuum wavelength is set to $\lambda =600\mathrm{nm}$, the relative refractive index to $m=1.59/1.33$ resulting in a size parameter $x\approx 13.928$.

2.3.

Properties of the Müller Matrix

Both introduced methods are based on scatterers of spherical symmetry. This leads to a reduction in the number of independent Müller matrix elements. Moreover, random orientation of these scatterers has an additional impact on the matrix properties. For example, arbitrarily shaped scatterers of random orientation lead to 10 independent elements. Employing scatterers which have a plane of symmetry further reduces the number of independent elements to 6 (see van de Hulst³⁸). In case of randomly arranged spheres (labeled $\mathrm{rs}$), the Müller matrix can be denoted as

3.1.

Effects on Multiple and Dependent Scattering

This study was aimed at pointing out agreements and differences between the two light scattering simulation approaches for polarized light. It can be observed that especially for small concentrations the two methods are in good agreement for all Müller matrix elements. This was already described for unpolarized light²⁷ and was now observed for polarized light within this study.²⁸ With increasing scatterer concentrations slightly increasing differences can be observed for all Müller matrix elements, mainly for higher scattering angles. These differences are attributed to dependent scattering effects since multiple scattering is considered by both solution methods. In Fig. 6 the angular errors are depicted for the same Müller matrix elements as above and three different volume concentrations. For $M_{11}$ the relative error and for all normalized elements $i,j=1,\dots ,4$ the difference $M_{ij}^{GMM}-M_{ij}^{\text{MonteCarlo}}$ between the both methods is shown. Only for large investigated volume concentrations above 20 Vol.-% dependent scattering effects become predominant, especially in the backscattering hemisphere. It is not possible to reproduce dependent scattering by standard solutions of the radiative transfer theory due to near field interactions. This has to be kept in mind whenever results of Monte Carlo simulations of polarized light propagation are regarded.

Fig. 6

Concentration dependent deviations between Maxwell and Monte Carlo results for all relevant Müller matrix elements. In case of $M_{11}$ the relative error is calculated, otherwise the difference between normalized elements. For volume concentrations of $f_V=10.05\mathrm{Vol}$.-% (red), $f_V=20.11\mathrm{Vol}$.-% (blue) and $f_V=40.21\mathrm{Vol}$.-% (black) the resulting deviations are depicted versus the scattering angle $\theta$ in the range from 0 to 180 deg. Apart from forward- and backscattering peaks ($\theta \approx 20-170\text{deg}$) deviations do not exceed 30%.

All Müller matrix elements show oscillations in the whole angular range, especially for small concentrations, while the oscillations are decreasing with increasing scatterer concentrations. These oscillations for small concentrations result from the solution of scattering by a single sphere (Mie scattering). With increasing concentrations, multiple scattering becomes predominant which leads to a reduction of the oscillation amplitudes.

3.2.

Results for the Unpolarized Intensity

The $M_{11}$ element is a normalized measure for the scattered intensity for unpolarized light. The differences in the $M_{11}$ element between the two methods for small angles ($\theta \lesssim 10\text{deg}$) are due to forward interferences resulting from the limited simulation volume (Maxwell theory).²⁷ The angular width of these forward interferences would become infinitesimally small if a laterally infinite instead of a finite cubic simulation volume would be regarded. Another effect is the coherent backscattering peak⁴³ at $\theta \asymp 180\text{deg}$ which can be reproduced by Maxwell theory in contrast to the radiative transfer theory. It can be observed for all concentrations and for scattering angles $\theta$ near $180\text{deg}$ and becomes more pronounced for high scatterer concentrations. Forward interferences and coherent backscattering cannot be reproduced by the Monte Carlo method as the radiative transfer theory neglects the wave character of light.

Beside that, it can be observed that the angular distribution of the $M_{11}$ element flattens for higher sphere numbers, i.e., the scattered intensity becomes more isotropic with increasing concentrations. This effect appears for both methods. Further, a flattening of the oscillations attributed to the single Mie scattering function is found for higher concentrations, also for the other Müller matrix elements. Increasing concentrations lead to irregular deviations of the $M_{11}$ element between both methods being most apparent for the highest concentration (96 spheres). The Monte Carlo solution stays rather smooth while the exact Maxwell results show more oscillations, even if no speckles are considered.

3.3.

Polarization Specific Müller Matrix Elements

The absolute value of all normalized Müller matrix elements (except $M_{11}$) generally decreases for higher angles as the concentration increases. This effect, investigated with Monte Carlo simulations, was attributed to multiple scattering by Tuchin et al.⁶ The Maxwell solutions presented here show a similar behavior even if increasing deviations to the Monte Carlo solutions for higher concentrations are taken into account. The decrease of these Müller matrix elements indicates the depolarizing character of multiple scattering. It is important to mention that the described decrease is mostly present, though for some angles there are also exceptions of this behavior.

For angles near forward direction, the $M_{22}$, $M_{33}$ and $M_{44}$ elements take values near 1, while the off-diagonal elements take values near 0 for all concentrations. This means that the polarization properties are better maintained in the forward direction. For both simulation methods this effect is dominant for small concentrations.

The $M_{44}$ element contains information about the circular polarization. It also decreases when scatterer concentrations increase for forward scattering angles $\theta \lesssim 50\mathrm{deg}$, but the decrease is smaller than for the $M_{22}$ and $M_{33}$ Müller matrix elements. For higher scattering angles, no general decrease with higher concentrations of the $M_{44}$ element can be observed from the Monte Carlo simulations. It could be confirmed that this effect can mainly be attributed to multiple scattering as was assumed by Maksimova et al.,¹⁸ who considered constant scatterer concentrations (size parameters in the Mie regime), but increasing scattering volumes. This may indicate that circular polarization is better maintained than linear polarization for multiple scattering. However, for very high scatterer concentrations as is the case for 96 spheres dispersed in the medium it can be observed (see Fig. 5) that the $M_{44}$ element calculated by Maxwell theory is below the $M_{44}$ element calculated by the Monte Carlo method while the $M_{22}$ element shows an opposite behavior. This shows that the better maintenance of circular polarization for high scatterer concentrations is not as manifest as could be expected only due to multiple scattering. In this case dependent scattering leads to a further reduction of the $M_{44}$ element. It has to be mentioned that these results hold for the parameters used in this study with scatterers in the Mie regime and a relatively high refractive index mismatch between scatterer and surrounding medium. For scattering media in the Rayleigh regime or in the Mie regime in case of scatterers with low relative refractive index it is possible that linear polarization is better maintained than circular polarization.⁵