Light scattering and the modeling thereof plays an important role for a wide range of scientific topics and engineering applications ranging from atmospheric research to medicine. It has been demonstrated that certain properties of a turbid medium can be determined from measurements of reflected polarized light.1, 2, 3, 4, 5, 6, 7, 8 This kind of diagnostics, however, is based on reverse engineering, which relies on accurate and efficient simulation tools. Two major theories exist to describe such phenomena. The more rigorous approach (analytical theory) is based on Maxwell’s electromagnetic equations, but due to its mathematical complexity and associated high computational cost, the simpler and computationally more efficient transport theory is preferred for many applications. This is based on a vector Boltzmann equation for the Stokes vector, which determines light intensity and polarization state as a function of position, propagation direction, and time. An established approach to solve this equation is the Monte Carlo method, where a large number of computational particles with individual properties are employed to represent the Stokes vector distribution. One possibility is to assign a Stokes vector, position, and propagation direction to each particle and use the so-called Stokes-Mueller formalism to evolve these properties.4, 5, 6, 7 In another equivalent approach, which is adopted in this paper, the evolution of complex electric field vectors is computed for each particle.9
Recently, there is a large interest in applying polarized light to examine properties of biological tissues.1, 2, 3, 10 Due to anisotropy of the tissue structure and due to the presence of chiral molecules (e.g., glucose and proteins), effects such as linear birefringence and optical activity have to be incorporated into the Monte Carlo methods. The usual assumption regarding linear birefringence is that a tissue is uniaxial, which means that there exists only one axis of anisotropy. The main axis along which the refractive index differs from those along the other main directions is called the extraordinary axis. Note that the ordinary axes along which the refractive index is uniform are perpendicular to the extraordinary axis. In most biological tissues, the birefringence value (difference between refractive indicies and along the extraordinary and ordinary axes, respectively) is small and, e.g., in Refs. 1, 2, where birefringence was modeled, the assumption that it has no significant influence on the scattering phase function was adopted. However, the effect of linear birefringence on the components of the electric field vector between scattering events is captured by the Jones N-matrix formalism as described in Refs. 2, 11 (retardation). In this paper, we want to examine the impact of on the backscattering Mueller matrices. Assuming spherical scattering particles embedded in a nonabsorbing, homogeneous, isotropic medium allows us to employ the Lorenz-Mie theory for the computation of the phase functions. While it is not straightforward to also take anisotropic background media into account, it is shown how the scattering phase function alters if the refractive index along the extraordinary instead of the ordinary axis is used for its computation. From theory, it also follows that the same value of has a stronger impact on the scattering phase function as the radius of the scattering particles grows.
In Sec. 2, the governing equations and some basic theory are presented together with the method for modeling linear birefringence. In Sec. 3, validation and computational studies are presented, and last, conclusions are given in Sec. 4.
In the transport theory for polarized light, one has to consider a vector Boltzmann equation, which readsis the Stokes vector, the position in physical space, the propagation direction, the wavelength, and the time. The vectors and form a so-called scattering plane, for which the single-scattering Mueller matrix is defined. Note that describes the interaction between the electromagnetic wave and an isolated particle; a detailed description of the single-scattering Mueller matrix can be found in Ref. 12. Extinction and scattering coefficients are denoted by and , respectively. A single scattering event requires rotation of the Stokes vector into the scattering plane, and its multiplication with the single-scattering Mueller matrix is required. The components of the Stokes vector are defined as and are the complex parts of the phasors , and . Here, and are the parallel and perpendicular components of the electric field vector with respect to the scattering plane. The phase difference between and is denoted by , and are the amplitudes of and , is the frequency, and is the phase of at ( is the speed of light).
The change of the electric field vector components and , which are in fact the Jones vector parameters, has to be computed at every scattering event; a complete description of the method can been found in Refs. 9, 13. Note, however, that the Stokes vector can be computed from and at any time for any scattering plane.
Monte Carlo Framework for Polarized Light
If the vector Boltzmann Eq. 1 is solved with a Monte Carlo method, a large number of particles has to be employed. To evolve a particle, first the random time until the next scattering event has to be sampled assuming exponential distribution with expectation . The new position then becomes , where and are particle position and propagation direction, respectively. The superscript denotes the state after collisions. Consistently, the particle clock time is updated as , and its weight is decreased by the absorbed energy . In order to compute the scattered propagation direction , which together with the incident propagation direction defines the scattering plane, the scattering angles and have to be sampled from a joint probability density function given asand are the elements of the single-scattering Mueller matrix. Then, and of the particle are changed as described in Ref. 9, whereas the new phasors of the particle, i.e., and , point in the direction of and , respectively. Note that not all particles experience the same number of collisions, and it is important to synchronize them. The local coordinate systems of all detected particles are first rotated into a so-called laboratory coordinate system, and then the Stokes vector is computed from the phasors and .
Last, in order to correctly capture effects at the interface between two media with different refractive indicies, Fresnel’s and Snell’s laws were implemented. An implementation of surface effects in the case of unpolarized light, where the scalar Boltzmann equation is solved, is described in Ref. 14. A nice description of the physics of light crossing an interface can be found in Refs. 15. Implementation in the Monte Carlo framework was already described in Refs. 16, 17.
Linear Birefringence in a Monte Carlo Framework
Media with anisotropic structure cause light to experience linear birefringence. Here, uniaxial media with a single axis of anisotropy are considered. The birefringence magnitude is defined by , where and are the refractive indicies along the extraordinary and ordinary axes, respectively. The Monte Carlo method was extended to include the effect of linear birefringence through the Jones N-matrix formalism. Here, we assume that linear birefringence has no effect on reflection and transmission of light at the medium interfaces and focus on its impact on the scattering phase function.
The refractive index is a function of the angle between the photon’s propagation direction and the extraordinary axis , and it is given by the expressioncaused by the difference in refractive indicies , i.e., is the traveled distance. As described in Ref. 2, the photon’s reference frame must be rotated first such that becomes parallel with the projection of the extraordinary axis onto the plane. Note that this rotation of the reference frame changes only the orientation angle of the polarization ellipse; the ellipticity stays the same. Now, the phasors and are expressed such that the M-matrix can be employed, i.e., the complex parts of the phasors change between two scattering events as
Jones N-Matrix Formalism
Details about this method can be found in Refs. 11, but for completeness here, the most important part is described. The effect of linear birefringence over an infinitely small optical path segment can be described by the so-called differential N-matrix and integration along the trajectory between two scattering events leads to the M-matrix. Since we trace the phasors and instead of the Stokes vector, the M-matrix can directly be applied. In our study, the depolarization effect occurs only due to multiple scattering, i.e., no depolarization occurs between scattering events.
The N-matrix for linear birefringence is given asis the retardation per unit distance. If a photon’s reference frame is chosen such that is parallel to , one can write and are defined as in Eq. 10 denotes the distance between two successive scattering events; i.e., the distance between two collisions. In Ref. 2, it was shown how two effects like linear birefringence and optical activity can be combined in one N-matrix. Note that various polarizing effects can be modeled with the N-matrix formalism. Since we consider only linear birefringence, the elements and are equal to zero. Once computed, the elements of the M-matrix are applied to the complex part of the phasors and as and have to be rotated first as described in Sec. 2.2 in order to apply Eq. 12. Although the use of the N-matrix formalism is unnecessarily complex for this study, since only one polarization effect is present, i.e., linear birefringence, it is used to make the developed model more flexible and extendable for future developments.
Before performing more relevant numerical investigations, models and implementations of the linear birefringence and surface effects (Fresnel and Snell laws) are validated. Then, scattering phase functions are shown for particles with different radii and refractive indices and for different refractive indices of the background medium, reflecting the difference between the refractive indices along extraordinary and ordinary axes. For many applications, it is of interest to estimate properties of scattering media by comparing measured and calculated Mueller matrices, which involves reverse engineering. Therefore, it is essential that the Monte Carlo method properly takes into account all relevant physical effects. So far, however, linear birefringence has typically been neglected. In the following, we want to demonstrate that this simplification may be questionable. For this paper, and are considered, and it is shown, based on the Lorenz-Mie theory, that the scattering phase functions based on and differ significantly, thus indicating that should not be ignored.
The Monte Carlo implementation described in the previous section was validated with experimental data.2 The experimental and numerical results published in Figs. 5 and 6 of Ref. 2 were kindly provided by the first author thereof. The simulation setup is shown here in Fig. 1 . A turbid medium is illuminated with circularly polarized laser light of wavelength propagating in the positive -direction, entering the turbid medium at the point (2, 0.5, 0). All geometrical dimensions are in cm, and the medium size is with the extraordinary axis parallel to the axis of the global coordinate system (see Fig. 1). The detectors have circular shape with an area of , and their centers are at (2, 0.5, 1) (detector 1) and at (2, 0, 0.5) (detector 2). The acceptance angle (angle between photon’s propagation direction and surface normal) is . Note that in Ref. 2, the same detectors were used in the experiments, but for the numerical studies, detectors with a slightly different geometry (rectangular area with a surface area of ) and an acceptance angle of were considered. The turbid medium, surrounded by air with refractive index of 1, is composed of spherical scattering particles embedded in a nonabsorbing host medium, which allows us to employ the Lorenz-Mie theory. The resolution to precompute and tabulate the scattering phase function is with . All spherical particles have the same radius and the same refractive index . The background medium has an ordinary refractive index of , and the birefringence value varies from 0 to . For the computations presented here, if not mentioned otherwise, the refractive index was used to calculate the scattering phase function; linear birefringence was considered only through retardation (as in Refs. 1, 2). Two test cases were considered: for the first, whose results are shown in Fig. 2 , a scattering coefficient of was chosen, and for the second, whose results are shown in Fig. 3 , was set to . The number of simulated particles was and for the media with lower and higher scattering coefficients, respectively. Again, note that in Ref. 2, the number of simulated particles was .
Figures 2 and 3 represent results “measured” with detector 1, while Figs. 2 and 3 show results for detector 2. In both figures, curves labeled with circles, squares, and asterisks represent changes of the V, U, and Q components, respectively, of the detected Stokes vectors as functions of . Solid and dashed lines show numerical and experimental results published in Ref. 2, while dashed-dotted lines represent results obtained with our method, i.e., with the code Scatter3D.13, 14, 18, 19 It can be observed that the results obtained with Scatter3D are in the same range as the numerical results published in Ref. 2. The small difference between the numerical simulations may be attributed mainly to the slightly different detector geometries.
Influence of Linear Birefringence on Phase Function
For the numerical studies presented in this and the following sections, the scattering phase functions were computed according to the Lorenz-Mie theory, whose input consists of the particle size parameter and the relative refractive index . While not included here, the influence of linear birefringence on the scattering phase function becomes apparent by comparing the scattering phase functions computed with and . Note, however, that this comparison delivers only an indication of the error due to neglecting this influence of linear birefringence, which can be further highlighted by comparing the corresponding backscattering Mueller matrices. Therefore, a domain of size ( is the optical mean free path length) is considered. Here, the global coordinate system is defined by the unit vectors , , and pointing in the directions of , the upper surface normal and , respectively. The local coordinate system of a particle is defined by the unit vectors , and . The incident laser beam intersects the medium surface at its center at an angle of and is characterized by the electric field vector , where and are uniformly distributed random numbers. The corresponding Stokes vector of the incident laser beam (defined in the local coordinate system) is , and the Stokes vector of the reflected particles is defined in the system with the unit vectors , , and . The number of simulated particles was , and all reflected particles are sampled on a structured grid at the surface of the turbid medium. Studies also reveal, as theoretically expected, that the influence of linear birefringence is more prominent for larger scattering particles and shorter wavelengths. Note that in all simulations, the influence of linear birefringence on reflection/transmission at the substrate surface was neglected. For retardation, however, it was considered according to Eqs. 4, 5.
For the first test case, phase functions for constant and , but two different values for , i.e., and , were computed. Note that this corresponds to . It can be observed in Fig. 5 that the difference between the two phase functions is negligible. The probability density function of the deflection angle [normalized product of scattering phase function and ] is plotted in Fig. 6 .
While the same , , and as for test case 1 were used in the second test case, was increased to 0.01, resulting in and . As can be observed in Figs. 7 and 8 , the two resulting phase functions differ significantly.
The plots in Fig. 9 show the 16 components of the backscattering Mueller matrix computed with on a patch on the surface centered around the incident point. Note that the values of each matrix element were normalized by the maximum value of the (1,1) component. The scattering medium is the one used in test case 2, and more details are provided in Ref. 13. As can be seen by comparing Fig. 9 with Fig. 10 , where no retardation is considered (isotropic medium) but with same , linear birefringence results in more diffuse patterns for all components except for (1,1), (1,2), (2,1), and (2,2). Note that the scattering coefficient was the same in all directions. A quantitative comparison with the backscattering Mueller matrix based on for the phase function computation is presented in Fig. 11 , where the normalized difference is shown for each of the 16 components. In the calculations of both backscattering Mueller matrices, was used for retardation.
Next, the difference between two backscattering Mueller matrices (one based on and one based on for the computation of the scattering phase function) for scattering media having the properties of test case 1 is shown in Fig. 12 . In both simulations, was used for retardation. As can be observed, the difference for the elements (1,2) and (2,1) is less than 10%.
In order to look more closely at the difference between the scattering phase functions based on and , the relative difference computed asand are the scattered light intensities at the angle for and , respectively, is plotted in Fig. 13 . The following parameters were used for these investigations: and . As expected, the average discrepancy results in a significant difference between the two corresponding backscattering Mueller matrices (see Fig. 11). On the other hand, for case 1, the average discrepancy and the difference between the backscattering Mueller matrices are much smaller (see Fig. 12). An interesting effect can be observed in Fig. 13, which shows that for the same but different particle radii (case 1) and (case 3), the average discrepancies are and , respectively.
Next, scattering phase functions are calculated for two test cases with but (test case 5) and (test case 6). The refractive indices of the host medium in ordinary and extraordinary directions are and , respectively. The dependence of on is presented in Fig. 14 , and it can be seen that the average discrepancy does not necessarily grow with increasing relative refractive index.
For light scattering in biological tissues, the influence of linear birefringence on the scattering phase function is usually neglected, which was the motivation for the performed study. It can be concluded that the influence of linear birefringence in scattering phase function computations depends on the particle size parameter and the relative refractive index . The difference in the scattered light intensity governed by phase functions based on and , as the Lorenz-Mie theory shows, is larger for larger scattering particles and also depends on the relative refractive index . Thus, not just the value of linear birefringence has to be considered in order to neglect its influence on scattering, but instead the comparison between two scattering phase functions has to be made. All this results in a variable so-called average discrepancy , which reflects both dependencies and quantifies the difference between two scattering phase functions computed with and . Note that the modeling of the effect of linear birefringence on a scattering phase function would require an extension of the Lorenz-Mie theory for an anisotropic refractive index of the background medium. The present study is just an indication of the problem and shows errors produced if one neglects the influence of linear birefringence on scattering. In this study, the proposed boundary value for the average discrepancy is 0.01. In all calculations of backscattering Mueller matrices, an influence of linear birefringence on reflection/transmission at the sample boundaries was neglected.
Birefringence is incorporated in our Monte Carlo framework through the Jones N-matrix formalism, and since we deal with the components of the electric field vector and , a M-matrix can be directly applied. The implementation was validated with results from an established reference before numerical studies were conducted.