1. INTRODUCTION

The polarization state of an optical wave undergoes a nonnegligible change as light propagates in a highly scattering medium, which can be modeled as a volume containing small randomly oriented birefringent compartments [1–3]. The intrinsic information about the scattered sample can be evaluated by measuring the depolarization degree, retardance, and other vectorial properties [4,5] of the input polarized light or by directly detecting the appearance of the new polarization component in scattered vectorial light. Modeling this birefringent scattering is essential and can provide a useful tool to quantitatively analyze the interaction between a complex birefringent sample and incident polarized wave, which is urgently needed in the fields that study the optical properties of semiconducting materials [6,7], design polarization-related integrated devices [8–11], fabricate micro/nanocomponents [12,13], and image biomedical samples for clinical applications [14,15]. The commonly used numerical finite-difference time-domain (FDTD) method [16,17] and the Monte Carlo method [18,19] can handle the propagation for electromagnetic problems of arbitrary computation sizes. However, the accuracy of the computation is sensitive to the sampling density of the computation volume, resulting in the requirement for high-speed access storage for large-scale simulation [16]. Moreover, the irreversible nature of these nonanalytical models makes solving the physically related inverse problems difficult, especially in the iterative reconstruction of 3D learning-based computational imaging [20] and end-to-end optimization of optical devices [8]. Therefore, development of an efficient and reversible computation model to deal with complex birefringent scattering is necessary.

The scattering process of a birefringent sample is more complicated than the common scalar scattering [21–26], and two factors must be considered: The first is the natural characteristics of the birefringent objects, which will exhibit different responses according to the external polarization state and will induce polarization changes for the incident polarized light field [27,28]; the second is depolarization of the far-field propagation of the polarized light field [29]. Several studies of the vectorial beam propagation method (VBPM) have been reported in the past few years. The VBPM is based on the vectorial Helmholtz equations, and the birefringent effects of the material and polarization properties of the electric field are all considered in this method. However, due to the complexity of solving the vectorial wave equation, the VBPM is not feasible for numerical computation. Some approximations are usually made in the VBPM, including a weak longitudinal polarization component [20], paraxiality [30], and small birefringence [31], and these approximations largely limit its accuracy and applicability. More importantly, for ill-posed inverse scattering problems, this inaccurate method will introduce unacceptable errors, making its use in the field of 3D computational imaging difficult [32–34].

Here, we present a multislice computational model for birefringent scattering. The proposed model considers the complete polarization components of the vectorial field and includes the full elements of the $3 \times 3$ scattering potential tensor without neglecting the longitudinal (i.e., z polarized) component. We exploit the dyadic Green’s function and define what we believe, to the best of our knowledge is a new polarization transfer function tensor (PTFT) to describe the polarization changes and depolarization during the scattering process. We verify the validity of our birefringent scattering model by simulating the vectorial light passing through the anisotropic objects and compare the results with those of the FDTD and VBPM method. Then, we use an optical-tweezer-assisted polarimeter to measure the full Mueller matrix of synthetic birefringent sample to confirm the theoretical proposal. In both simulations and experiments, the proposed computational model is proven to be highly accurate and efficient.

2. THEORY

A. Birefringent Scattering Model

The proposed model includes the three orthogonal polarization components by writing the electric field in the form of a column vector:

The polarization coupling problem during the scattering process is complicated and highly nonlinear. The first Born approximation [37] is used here, and it only considers the single scattering of the incident wave. This is valid when the scattered electric field is much smaller than the incident field. The total vectorial electric field $\vec u({\vec r})$ is composed of an unscattered incident vectorial field ${\vec u_{\text{in}}}({\vec r})$ and a single scattered vectorial field ${\vec u_s}({\vec r})$, where ${\vec u_s}({\vec r})$ can be linearly related to the defined scattering potential tensor ${\stackrel{\leftrightarrow} V} ({\vec r})$ by the dyadic Green’s function. The total field can be formulated by the 3D dyadic Green’s function-related first-order Born approximation, which is written as

B. Polarization Transfer Function Tensor

Since the birefringent scattering model [Eq. (3)] describes the sample-induced polarization changes by the convolution relationship between the 3D point spread function tensor and scattering potential tensor, we can give a comprehensive description of the polarization coupling in the frequency domain by using the Fourier representation of the 3D dyadic Green’s function. After imposing the constraints on the propagation vector and using Cauchy’s residual theorem for Fourier integration, we obtain the forward 3D dyadic Green’s function in terms of a continuous angular plane wave:

 

Fig. 1. Illustration of proposed multislice computational model for birefringent scattering. (a)–(b) Amplitude and phase distribution of the PTFT. (c) 3D birefringent samples are decomposed into multiple thin slices. The vectorial electric field is composed of an unscattered field and a single scattered field for every scattering slice, and its polarization will change during the scattering process.

Download Full Size | PDF

C. Multislice Birefringent Scattering Algorithm

In the scattering model [Eq. (3)], the scattering potential only interacts with the unscattered vectorial incident field, and the multiple scattering during the propagation is neglected; this approximation is only suitable for a thin sample. To overcome this restriction and make the model computationally accurate, we use the previously reported multislice method [23] to improve the applicability of the birefringent scattering model. We divide the 3D scattering potential into multiple 2D slices with a finite thickness ${{\Delta z}}$ along the axial direction, as shown in Fig. 1(c). For each slice, we apply the 3D dyadic Green’s function-related, first-order Born approximation, and we can then obtain the vectorial multislice model suitable for birefringent scattering. The forward propagation model can be written as

The vectorial free-space propagation to the ${(n + 1)_{\text{th}}}$ slice must first be calculated based on the given boundary field, which is the incident vectorial field $\vec u_{\text{in}}^n({x,y,n\Delta z})$ at slice ${n_{\text{th}}}$. The solution of this problem belongs to the classical Dirichlet problem [38], and it can be approached by applying the method of images [39] to the vectorial field. In this case, a modified 3D dyadic Green’s function has the form

The modified 3D dyadic Green’s function subtracts the image term ${R_1}$ to satisfy the boundary condition. By using Eq. (10) and considering the radiation condition in the far field, the vectorial free-space propagation process can be formulated as

Let ${\vec r} = ({x,y,({n + 1})\Delta z})$, ${{\vec r} ^\prime} = ({x^\prime ,y^\prime ,n\Delta z})$, and the normal vector $\vec n$ be on the boundary points in the ${-}{\rm{z}}$ direction. Substituting Eqs. (6), (7), and (10) into Eq. (11), the vectorial free-space propagation process can be written as

Meanwhile, based on the convolution theorem, and assuming that the scattering potential tensor does not vary axially within each slice, we can simplify the forward birefringent scattering in Eq. (9) as

The vectorial scattering operator ${\textbf{H}}$ describes the sample scattering-induced polarization change. By sequentially using Eqs. (11) and (14) on each slice, the 3D forward vectorial field can be efficiently evaluated. See Supplement 1, Section 3, for the derivation of the forward birefringent scattering process.

3. SIMULATION VERIFICATION

We confirmed the accuracy of the proposed computational model by simulating the vectorial light passing through the birefringent objects and comparing the scattering results with those of the FDTD and VBPM methods. In our simulation, a set of birefringent phantoms with known birefringent tensors were chosen to demonstrate the verification, and we used the results calculated by FDTD as the ground truth. As shown in Fig. 2(a), four 3 µm diameter birefringent beads were placed in a homogenous background medium with ${n_m} = 1.33$. We set the principal RI as ${n_{\textit{xx}}} = 1.37$, ${n_{\textit{yy}}} = 1.40$, and ${n_{\textit{zz}}} = 1.44$. The RI contrast with respect to the background medium is nearly 0.11, which is largely beyond the Born approximation [40]. To study the effect of the birefringence orientation on the incident polarized light, we made four birefringent beads have different birefringence orientations. The optical axes of the first three spheres were set parallel to the coordinate axes but mutually orthogonal; hence, the dielectric tensors of these three spheres only have the on-diagonal components. The fourth sphere was set to ${\theta _x} = {\theta _y} = {\theta _z} = \frac{\pi}{4}$, and it has off-diagonal dielectric components. Figure 2(b) shows the 3D organization of the four spheres, where the double yellow arrows indicate the birefringence orientations. The total computational volume was set to $11.7\;{\rm{\unicode{x00B5}{\rm m}}} \times 11.7\;{\rm{\unicode{x00B5}{\rm m}}} \times 4.55{\rm{\;\unicode{x00B5}{\rm m}}}$ with a sampling size of 65 nm in all three directions. In this simulation, we used an x-polarized 405 nm plane wave with an equal amplitude of 1 as the incident light.

 

Fig. 2. Simulation of the birefringent scattering by FDTD, VBPM, and the proposed model in normal illumination. (a) Parameter used to define the birefringent spheres. (b) 3D rendering of the birefringent spheres. The double yellow arrows indicate the birefringence orientations. (c)– (d) Amplitude and phase distribution of three orthogonal polarization components ${E_x}$, ${E_y}$, and ${E_z}$ after the light passes through the birefringent samples. The first row shows the ground truth calculated by FDTD, the second row shows the simulated results from VBPM, the third row is the results given by the proposed model, and the final row shows the discrepancies between our method and FDTD.

Download Full Size | PDF

Figures 2(c) and 2(d) show the amplitude and phase distribution of the three orthogonal polarization components ${E_x}$, ${E_y},$ and ${E_z}$ after the light passes through the birefringent samples. Not being able to include depolarization elements during the scattering, VBPM fails to give accurate results for #1–#3 spheres. The calculated scattering components ${E_y}$ and ${E_z}$ are zero, as shown in the second row in Fig. 2(c) and differs significantly from the FDTD results. The phase distribution calculated by VBPM further shows a large variation with the FDTD calculation as shown in Fig. 2(d). On the contrary, the proposed method agrees well with FDTD results in all three polarization components. The calculated vectorial fields follow the same distribution in both the amplitude and phase compared to the FDTD method, and only have the maximum 6% calculated errors. Importantly, running on a workstation equipped with $2 \times$ Inter Xeon E5-2697A V4, 2.6 GHz 128 GB, the computation time of our proposed model was only 0.77 s, which is significantly faster than that of the FDTD method (35 s) and VBPM (67.7 s).

Furthermore, we compared the calculated results between the three methods with the oblique illumination of ${\rm NA} = 0.8$. The paraxial approximation based VBPM gives more severe errors in oblique illumination. However, the proposed method shows little discrepancies in both amplitude and phase. Supplement 1, Section 4, presents the detailed derivation and error analyses for VBPM and the simulated results under highly oblique illumination.

We also gave an anisotropic phantom with radial birefringence distributions to show the accuracy of our method, and the results are shown in Supplement 1, Section 5. To test the validity of our method for the high birefringent RI contrast sample, more detailed analysis is discussed in Supplement 1, Section 6.

4. EXPERIMENTAL VERIFICATION

To verify the proposed method, we built an optical tweezer-assisted polarimeter to measure the full Mueller matrix. The experimental setup is shown in Fig. 3. For the imaging part, to reduce the coherent and environmental noise, we used a fiber-coupled LED (M405FP1, Thorlabs) as the imaging source. The x-polarized collimated imaging light first passed through a 4f system comprising an achromatic lens L1 (AC254-100-A, Thorlabs) and a water-immersion objective Obj1 ($60 \times$, ${\rm NA} = 1.0$, LUMPLFLN60XW). The half-wave plate HWP1 (WPH10M-405, Thorlabs) and the quarter-wave plate QWP1 (WPQ10M-405, Thorlabs) were used to generate four independent polarization states, including three linearly polarizations horizontal (H), vertical (V), and 45° (P), and one circularly polarization right circular (R). The polarized light passing through the birefringent sample was collected using an oil-immersion objective Obj2 ($100 \times$, ${\rm NA} = 1.4$, UPlanSApo, Olympus) and collimated by L2 (ACT508-180-A-ML, Thorlabs). A polarizer P (LPUV100, Thorlabs) and another quarter-wave plate QWP2 (WPQ10M-405, Thorlabs) were placed in front of an sCMOS camera (ORCA-Flash4.0 V3, Hamamatsu), and they were rotated to measure the Stoke vector by filtering the H, V, P, and R polarization components of the scattering light. The exposure time of the camera was set at 100 ms for all the measurements. The full Muller matrix of the scattering sample can be calculated by using all 16 intensity measurements. For the trapping part, a single longitudinal mode laser (${{\lambda}} = 473\;{\rm{nm}}$, 200 mW, CNI) was used as the trapping beam. After collimation by lens L3 (AC254-040-A, Thorlabs), the polarized trapping light passed through a half-wave plate HWP1 (WPH10-405, Thorlabs), which could control the polarization of the light. The beam was then reflected by a pair of mirrors to the back focal plane of the objective lens and formed optical trapping after being focused by Obj2. The imaging light and trapping light were separated by a long-pass dichroic mirror DM (DMLP425R, Thorlabs). Using the trapping beam, we can orient the optical axis (slow axis) of the vaterite particles to any angle about the beam axis and study the effect of the birefringent orientation for the vectorial scattering process.

 

Fig. 3. Experimental setup of the optical tweezer-assisted polarimeter to detect a full Mueller matrix.

Download Full Size | PDF

The birefringent samples used in our experiment were vaterite particles of 4–6 µm in diameter. The vaterite particles have a polycarbonate structure of calcium carbonate composed of nanocrystals that are 20–30 nm with positive uniaxial birefringence [41–43]. The alignment of the optical axes of the vaterite particles has a hyperbolic distribution throughout the volume, which makes the vaterite have a high birefringence coefficient of up to 0.1 [44,45]. The size of the vaterite particles is essential to detect the vectorial scattering light. Particles with small sizes cannot induce sufficient scattering and will bring very large discrepancies compared to the numerical results. Hence, we modified the synthesis of the vaterite particles based on a previously published protocol [41,42]. Aqueous solutions of ${\rm CaCl}_{2}$, ${\rm MgSO}_{4}$, and ${\rm K}_{2}$${\rm CO}_{3}$were prepared with a molarity of 0.1 M. First, 1.5 mL of ${\rm CaCl}_{2}$was mixed with 60 µL of ${\rm MgSO}_{4}$ in a 5 mL reaction vessel, followed by 90 µL of ${\rm K}_{2}$${\rm CO}_{3}$. The solution was agitated by violently pipetting the solution for 4 min. 100 µL of a surfactant solution (XYS-3500, Yancheng Yunfeng Chemical Co., Ltd.) was added to the solution to stabilize the reaction. Then, the synthetic particles in the first step were used as the seeds to promote the crystal growth. 5 mL of ${\rm CaCl}_{2}$ and 1.5 mL of ${\rm MgSO}_{4}$ were mixed in a 10 mL reaction flask, and then 20 µL of the seed solution and 1 mL of ${\rm K}_{2}{\rm CO}_3$ were added. The solution was gently mixed and kept undisturbed for 5 min. After that, 100 µL of a surfactant solution was added to stop the crystal growth. The concentration of the surfactant solution is essential to the growth of the vaterite particles. We dissolved 2.2 mg of the surfactant in 200 mL of distilled water and achieved the expected particle size.

The experimental and simulated results are shown in Fig. 4 (all elements have been normalized by ${m_{11}}$). The optical axis of the trapped vaterite has been indicated by the yellow double axis. The Mueller matrix the exhibits a mirror symmetry [46] that can be written as

 

Fig. 4. Simulated and experimental full Mueller matrix. (a) and (b) Slow axis of vaterite particle was rotated to the horizontal direction. (c) and (d) Slow axis of vaterite particle was rotated to the vertical direction. The colorbar for the diagonal elements ${m_{34}}$ and ${m_{43}}$ is [${-}1$, 1], and for other elements is [${-}0.1$, 0.1]. All the elements are normalized by ${m_{11}}$.

Download Full Size | PDF

As shown in Fig. 4, this symmetric distribution indicates that the vaterite particle is an anisotropic material [46]. The values of the diagonal Mueller matrix vary across the imaging, indicating that the depolarization happens during the scattering process [3]. The depolarization-induced reduction of ${m_{22}}$ and ${m_{23}}$ represents linear depolarization, and the reduction of ${m_{44}}$ shows circular depolarization. Moreover, the appearance of the off-diagonal elements ${m_{34}}$ and ${m_{43}}$ indicate the existence of birefringence in the scattering sample. When birefringent orientation of the vaterite particle is rotated by the trapping light, the elements ${m_{12}}$, ${m_{21}}$, ${m_{34}}$, and ${m_{43}}$ will change greatly, and these elements can be used to distinguish the vaterite particles with different birefringent orientations.

In the simulation, the phantoms were assumed to be the vaterite microspheres with uniformly orientated birefringent alignment. However, the experimentally synthesized vaterite particles were having polycarbonate structures of the calcium carbonate composed of many nanocrystals with positive uniaxial birefringence. The 3D dielectric distribution inside the synthesized particles were not perfectly homogenous, and the birefringent orientations of the particle were not uniformly aligned. The proposed model, however, still exhibits good agreement with the experimental results, especially the birefringent related elements ${m_{34}}$ and ${m_{43}}$. Moreover, Supplement 1, Section 7, gives the intensity distribution of the scattering light field, and the experimental and theoretical results of the cross-polarized field ${| {{E_y}} |^2}$ agree well in value and distribution, which further proves the accuracy of the proposed model.

5. DISCUSSION AND CONCLUSION

In conclusion, we propose a multislice computational model for birefringent scattering. Taking the complete vectorial polarization components into consideration, our proposed model is unambiguous by using the dyadic Green’s function. A new PTFT is defined to describe the polarization changes and depolarization during the scattering process. The vectorial properties of the electric field associated with the birefringent material are fully considered without neglecting any longitudinal polarization component in the model. The theoretical and experimental results clearly support the general applicability of this computational model and prove that can perform various scattering calculations with high accuracy and efficiency. Moreover, the proposed vectorial model can also be used to analyze the reflected vectorial scattering field. Supplement 1, Section 8, gives the analytical model for reflection in birefringent scattering.

The proposed method can be used to solve the 3D reconstruction of the scattering potential tensor. Currently, the main challenge of tensor imaging is the fact that ${E_z}$ cannot not be measured on the camera plane; hence, the $3 \times 3$ tensor becomes the $2 \times 2$ by neglecting the longitudinal component. Progress on tensor imaging was recently reported using the novel angular spectrum decomposition method [47]. This method, however, needs a polarization-sensitive holographic setup to measure the complex field and calculate the ${E_z}$ pixel by pixel, which is time-consuming and complicated for large-scale 3D tensor imaging. Alternatively, even if the measurement is incomplete (i.e., the experimental ${E_z}$ component is missing), the multiple-angle illumination can produce different transverse polarization components (as described by the polarization transform matrix ${\textbf{M}}$ in Supplement 1, Section 4) which can be used to solve the ill-posed inverse problem by iteratively minimizing the cost function between the theoretical results $\vec u({\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over V}})$ and experimental measurements $\vec y$ by following

The proposed method can be applied not only in the 3D reconstruction of computational imaging, but also in research fields related to vector light engineering [51]. The proposed method, for example, would be a useful tool to study the interaction of the special mode optical fields with trapped particles to explore special phenomenon and effects [52,53]. The dynamic characteristic optical force and toque on the particles is essential to design the optical tweezer [54]. Moreover, the proposed method also can be deployed to model and evaluate the polarization-related end-to-end design of photonic devices [55–57], providing a general analytical tool to tackle the challenges of the vectorial scattering problem.