## 1.

## Introduction

Many biological tissues contain large amounts of fibrous components that exhibit birefringent properties; examples of such components are collagen fibers in tendons, intracellular microtubules, and microfilaments.^{1}2.^{–}^{3} In a recent study,^{4} we showed that scattering by well-aligned cylindrical scatterers displayed anisotropic effects, including a retardance that was similar to that of a piece of birefringent media whose extraordinary axis, and the fast axis, was along the axial direction of the aligned cylinders. We also examined the different contributions to the polarization features provided by the well-aligned cylindrical scatterers and the anisotropic interstitial medium.^{5}6.^{–}^{7} Several other groups have tried to quantify the direction of a birefringent media based on polarized light imaging and studying the three-dimensional (3-D) orientation of the optical axis of birefringent turbid media.^{8}^{,}^{9} It has been shown that for biological tissues containing several layers of different anisotropic properties, the axial variation in the birefringence can result in ambiguity when the measurement is averaged over the sampling volume.^{10}11.^{–}^{12}

In this study, we examined in detail how well-aligned cylindrical scatterers and birefringent media contributed to the value and direction of the total retardance of the media. Using Monte Carlo simulations,^{13} we simulated how the total retardance and the equivalent extraordinary axis of the bilayer sample changed with the change in the value of the retardance and directions of the extraordinary axis of the two anisotropic layers, respectively. The results showed that the retardance of the well-aligned cylindrical scatterers could be regarded as the retardance generated by a piece of birefringent media, whose equivalent extraordinary axis is aligned with the axial direction of the cylinders. We started with the superposition of two layers of uniaxial birefringent media whose axes were in the same plane, perpendicular to the propagation direction of the incident light. Furthermore, we showed that the characteristic features of the bilayer birefringent sample (BBS) were maintained even if one or both of the birefringent layers were replaced by well-aligned cylindrical scatterers.

## 2.

## Theory and Method

## 2.1.

### Mueller Matrix Decomposition Method

Lu and Chipman^{14} proposed a Mueller matrix containing three physical properties that can be extracted via a polar decomposition process. A Mueller matrix $M$ can be regarded as the product of three constituent “basic” matrices, representing depolarization (${M}_{\mathrm{\Delta}}$), retardance (${M}_{\mathrm{R}}$), and diattenuation (${M}_{\mathrm{D}}$), respectively

The depolarization coefficient $\mathrm{\Delta}$ can be calculated from the elements of ${M}_{\mathrm{\Delta}}$ and the linear retardance $\delta $ can be calculated from the elements of ${M}_{R}$

## (2)

$$\mathrm{\Delta}=1-\frac{|\mathrm{trace}({m}_{\mathrm{\Delta}})|}{3},\phantom{\rule[-0.0ex]{1em}{0.0ex}}0<\mathrm{\Delta}<1,$$^{15}

## (4)

$${r}_{i}=\frac{1}{2\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\delta}\sum _{j,k=1}^{3}{\u03f5}_{ijk}{({m}_{\mathrm{R}})}_{jk},$$## 2.2.

### Monte Carlo Simulations and Scattering Model

In earlier studies, we developed polarized photon scattering models and the corresponding Monte Carlo simulation program for examining the propagation of polarized photons in an anisotropic, turbid sample.^{13} The simulation program tracks the trajectory and polarization states of the photons; corresponding illuminations of at least four different polarization states of the forward scattering photons are recorded to evaluate the Mueller matrix of the sample

## (6)

$${S}_{\text{out}}=M\xb7{S}_{\text{in}}\phantom{\rule{0ex}{0ex}}{\left(\begin{array}{l}I\\ Q\\ U\\ V\end{array}\right)}_{\text{out}}=\left(\begin{array}{cccc}{m}_{11}& {m}_{12}& {m}_{13}& {m}_{14}\\ {m}_{21}& {m}_{22}& {m}_{23}& {m}_{24}\\ {m}_{31}& {m}_{32}& {m}_{33}& {m}_{34}\\ {m}_{41}& {m}_{42}& {m}_{43}& {m}_{44}\end{array}\right){\left(\begin{array}{l}I\\ Q\\ U\\ V\end{array}\right)}_{\text{in}}.$$^{4}5.6.

^{–}

^{7}

^{,}

^{13}The diameter and density of the spherical and cylindrical scatters, the degree of the alignment and spatial orientation of the cylindrical scatters, the optical properties of the scatterers, and the interstitial media were all variables used in the simulations to mimic the properties of turbid media such as tissues.

## 3.

## Results and Discussion

## 3.1.

### Bilayer Birefringent Sample

Typically, the Mueller matrix of a linear retarder can be written as in Ref. 16:

## (7)

$${M}_{\mathrm{R}}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& {\mathrm{cos}}^{2}\text{\hspace{0.17em}}2\theta +{\mathrm{sin}}^{2}\text{\hspace{0.17em}}2\theta \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\delta & \mathrm{sin}\text{\hspace{0.17em}}2\theta \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}\hspace{0.17em}}2\theta (1-\mathrm{cos}\text{\hspace{0.17em}}\delta )& -\mathrm{sin}\text{\hspace{0.17em}}2\theta \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\delta \\ 0& \mathrm{sin}\text{\hspace{0.17em}}2\theta \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}2\theta (1-\mathrm{cos}\text{\hspace{0.17em}}\delta )& {\mathrm{sin}}^{2}\text{\hspace{0.17em}}2\theta +{\mathrm{cos}}^{2}\text{\hspace{0.17em}}2\theta \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\delta & \mathrm{cos}\text{\hspace{0.17em}}2\theta \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\delta \\ 0& \mathrm{sin}\text{\hspace{0.17em}}2\theta \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\delta & -\mathrm{cos}\text{\hspace{0.17em}}2\theta \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\delta & \mathrm{cos}\text{\hspace{0.17em}}\delta \end{array}\right].$$We considered a sample consisting of two parallel birefringent layers with retardance values of ${\delta}_{1}$ and ${\delta}_{2}$ and extraordinary axes at ${\theta}_{1}$ and ${\theta}_{2}$, respectively. Figure 1 shows a schematic illustration of the two layers of birefringent media superposed in the same optical path; here, the planes including both the ordinary and extraordinary axes are parallel to each other and perpendicular to the light propagation direction.

The Mueller matrix of such a BBS is

## (8)

$${M}_{\mathrm{R}}(\delta ,\theta )={M}_{{\mathrm{R}}_{2}}({\delta}_{2},{\theta}_{2})\xb7{M}_{{\mathrm{R}}_{1}}({\delta}_{1},{\theta}_{1}).$$Figure 2 shows how the total retardance $\delta $ and the equivalent extraordinary axis orientation $\theta $ varied with changes in the retardance and orientation of the two birefringent layers. In the calculations, the retardance and the orientation of the first layer (${\delta}_{1},{\theta}_{1}$) remained constant, the retardance of the second layer (${\delta}_{2}$) changed from 0 to 1 rad, while the angle of intersection between the two extraordinary axes took different values in the range from 0 deg to 90 deg. In Figs. 2(a) and 2(b), ${\delta}_{1}=0.25\text{\hspace{0.17em}\hspace{0.17em}}\text{rad}$, and in Figs. 2(c) and 2(d), ${\delta}_{1}=0.5\text{\hspace{0.17em}\hspace{0.17em}}\text{rad}$. The quantitative variation of the retardance and equivalent extraordinary axis from the BBS was independent of the retardance produced by the first layer.

As shown in Figs. 2(a) and 2(c), if the extraordinary axes of the two anisotropic layers were parallel, the total retardance $\delta $ was equal to the sum of ${\delta}_{1}$ and ${\delta}_{2}$. The contributions from the two layers to the total retardance were additive when the angle of intersection between the two extraordinary axes was smaller than 45 deg, but cancelled each other out when the angle of intersection exceeded 45 deg. If the extraordinary axes of the two layers were perpendicular to each other, or if $\theta =90\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$, the anisotropy due to the two layers cancelled out exactly at ${\delta}_{2}={\delta}_{1}$, resulting in an isotropic bilayer sample. When ${\delta}_{2}$ was significantly larger than ${\delta}_{1}$, the slope of the total retardance curve approached 1, indicating that layer 2 provided the dominant contribution to the total retardance $\delta $ and to the anisotropic properties of the bilayer sample.

Figures 2(b) and 2(d) show the orientation of the equivalent extraordinary axis of the bilayer sample as a function of ${\delta}_{2}$. The equivalent extraordinary axis of the bilayer sample rotated gradually from the direction of the first layer axis (when ${\delta}_{2}$ was much smaller than ${\delta}_{1}$) to the direction of the second layer axis (when ${\delta}_{2}$ became much larger than ${\delta}_{1}$). When the extraordinary axes of the two layers were perpendicular to each other, the equivalent extraordinary axis of the bilayer sample jumped abruptly from the direction of layer 1 (for ${\delta}_{2}<{\delta}_{1}$) to the direction of layer 2 (for ${\delta}_{2}>{\delta}_{1}$). The orientation of the extraordinary axis was undefined at ${\delta}_{2}={\delta}_{1}$, because the anisotropy of the two layers cancelled each other out completely when the angle of intersection was 90 deg. In other words, under such a condition, the bilayer medium became isotropical.

Following the analysis in the previous section, we further applied Monte Carlo simulations and Mueller matrix decomposition to two types of bilayer anisotropic samples. Figure 3(a) shows a sample with a layer of birefringent media and a layer of well-aligned cylindrical scatterers [Birefringence and cylindrical scatterer sample (BCS)], and Fig. 3(b) shows a sample with two layers of cylindrical scatterers (CCS).

It was shown in our previous study that well-aligned cylindrical scatterers contribute to the retardance.^{4} The extraordinary axis of the scattering-induced birefringence is the extraordinary axis and is along the direction of the well-aligned cylinders. The magnitude of the birefringence is related to the density, size, and order of alignment of the cylindrical scatterers. Next, we investigated whether the characteristic features in Fig. 2 still applied to the total retardance and equivalent extraordinary axis of a bilayer anisotropic scattering sample.

## 3.2.

### Birefringence and Cylindrical Scatterer Sample

For the BCS, we varied the angles between the orientation of the well-aligned cylinders and the extraordinary axis of the birefringent layer from 0 deg to 90 deg. In the simulation, the diameter of the cylinders was $1\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, the scattering coefficient of the aligned cylinder layer was ${\mu}_{\mathrm{s}}=15\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$, and the thickness of the sample was 1 cm. In Fig. 4(a), the vertical axis represents the total retardance of the BCS and the horizontal axis is the retardance produced by the birefringent layer. The different shapes in Fig. 4(a) represent different angles between the aligned cylinders and the extraordinary axis of birefringent media, and the solid lines represent the calculated results for the BBS. Here, the nonzero starting value of 0.38 rad corresponded to the retardance produced by the well-aligned cylindrical scatterers. Figure 4 shows that there was good agreement between the simulated results for BCS and the calculated curves for BBS. The equivalent extraordinary axis of BCS also rotated gradually from the direction of aligned cylinders to the direction of the extraordinary axis of the birefringent layer as the retardance of the birefringent layer increased. The very similar anisotropic behaviors observed for the BCS and BBS confirmed that the well-aligned cylindrical scatterers generated birefringence, whose extraordinary axis was along the direction of the aligned cylinders.

## 3.3.

### Two Layers of Cylindrical Scatterers Sample

Since many real fibrous tissues have multilayer microstructures with layers aligned in different directions, we considered further samples consisting of two layers of cylindrical scatterers with different orientations and alignments, as shown in Fig. 3(b). It is known from a previous study^{4} that the birefringence induced by aligned cylindrical scatterers is affected by the diameter, the optical properties, and the alignment of the cylinders, as well as the scattering coefficient and thickness of the layer. In the simulations, we kept the parameters of the first layer constant and changed the scattering coefficient and the other parameters of the second layer to examine in detail how the scattering-induced anisotropy of the two layers affected the total retardance and the equivalent extraordinary axis of the sample. In Figs. 5(a) and 5(b), the cylinders were perfectly aligned without any fluctuations in the orientation. The horizontal axis represents the retardance of the second layer of cylindrical scatterers obtained by decomposing the simulated Mueller matrices. For Figs. 5(c) and 5(d), the orientations of the cylindrical scatterers followed a Gaussian distribution, whose full width at half maximum (FWHM) was $\eta =20\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$. In the simulations, the diameter of the cylinder was $1\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, the first layer cylinder scattering coefficient was ${\mu}_{\mathrm{s}}=15\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$, the second layer cylinder scattering coefficient changed from ${\mu}_{\mathrm{s}}=0\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$ to ${\mu}_{\mathrm{s}}=40\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$, and the thickness of the two layers was 1 cm, corresponding to a change in retardance from 0 to 1.14 rad. The different shapes show data for different intersection angles between the directions of the two layers of cylindrical scatterers. The solid lines represent the calculated results for BBS. We also investigated a mixed model with two types of cylindrical scatterers and different alignment orientations in the two layers. The simulation results in Fig. 5 were very similar to those shown Figs. 4(a) and 4(b). This good agreement confirmed that the anisotropic optical features from the coupling of the two layers of cylindrical scatterers with different orientations were approximately equivalent to the features observed for the combination of two pieces of birefringent media with different extraordinary axes.

We also studied the variation of the depolarization in the bilayer model; the first layer cylinder scattering coefficient was ${\mu}_{\mathrm{s}}=15\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$, and the second layer cylinder scattering coefficients were ${\mu}_{\mathrm{s}}=20\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$ and ${\mu}_{\mathrm{s}}=40\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$. In Fig. 6, the horizontal axis is the angle between the two layers of cylinders, and the vertical axis is the variation of the retardance and depolarization in this model. Similar to our previous study,^{4} an increased disorder in the orientation of the cylinders resulted in smaller retardance and larger depolarization, which could explain the decreasing retardance and increasing depolarization with increases in the angle between the cylinder scatterers.

## 3.4.

### Experiment Results and Analysis of Bilayer Glass Fibers

We made a bilayer fibrous phantom to validate the simulated and calculated results experimentally. In this experiment, we used two slabs of well-aligned glass fibers as the cylindrical scatterers. The two slabs were placed together with four different intersection angles; these angles were 0 deg, 30 deg, 45 deg, and 90 deg, as shown in Fig. 7. The cylindrical scatterers were made from nonbirefringent glass fibers with a refractive index of 1.547 and a radius of $5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$. The thickness of each glass fiber layer was $\sim 1.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$, the scattering coefficient was approximately ${\mu}_{\mathrm{s}}=50\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{cm}}^{-1}$, and the absorption coefficient ${\mu}_{\mathrm{a}}$ was $\sim 0$. The experimental setup for the Mueller matrix detection was a typical forward measurement configuration. The details of the experimental setup can be found in the literature.^{4} The mean experimentally measured retardance values corresponding to the intersection angles of 0 deg, 30 deg, 45 deg, and 90 deg are shown in Fig. 8; the results were consistent with the results of the theoretical calculations and the simulation results as shown in Fig. 6(a).

## 4.

## Conclusion

Our previous research showed that there are two sources of anisotropy in tissue: one is cylindrical scatterers and the other is birefringent interstitial media. In this study, we investigated the coupling between two-layer anisotropic samples with different extraordinary axis orientations. We applied Monte Carlo simulations and a Mueller matrix decomposition method for three types of bilayer anisotropic samples, and examined how the total retardance and equivalent extraordinary axis varied as a function of the rotation of one of the layers and the consequent variation in its extraordinary axis and retardance values. The calculations for the two anisotropic layers showed that the total retardance was not a linear superposition of the retardance of the two layers. When the birefringence of the first layer was fixed, and the birefringence of the second layer was increased from zero to well above that of the first layer, the equivalent extraordinary axis of the bilayer sample rotated gradually from the first layer to the second layer. During the same process, the total retardance always increased if the intersection angle between the extraordinary axes of the two layers was smaller than 45 deg, but first decreased and then increased when the intersection angle was larger than 45 deg. When the axes of the two layers were orthogonal, the total retardance became zero when the retardance of the two layers was equal. When one or both of the birefringent layers were replaced by well-aligned cylindrical scatters, Monte Carlo simulation results for the bilayer samples displayed the same characteristic features as those shown by the BBS in terms of both the total retardance and the equivalent extraordinary axis. These results confirmed that the well-aligned cylindrical scatterers generated retardance and behaved the same as a layer of birefringent media. A preliminary experiment using glass fibers as cylindrical scatterers was conducted and the results confirm the above conclusions. This study helped to reveal the connection between observed anisotropy and microstructures and provided more evidence for the feasibility of the quantitative characterization of complex tissues.

## Acknowledgments

This work has been supported by the National Natural Science Foundation of China, Grants Nos. 11174178, 11374179, 61205199, and 41475125.