Figure 1 shows the Mueller matrix testing setup. A light beam emitted by an 83.06-mW laser at a wavelength of 532 nm passes through an attenuator (ATT) to adjust the illumination power. A rotatable linear polarizer (LP1, LPVISC050-MP2, Thorlabs, US) and a quarter wave-plate (W1, AQWP10M-580, Thorlabs, US) composed of the polarization state generator (PSG) are placed in front of the ATT to control the incident polarization state as one of horizontal (H), vertical (V), 45 ∘ (P), 135 ∘ (M) linear polarization states, and right-hand (R) and left-hand (L) circular polarization states. After beam expansion by the expander (BE, BE02-05-A, Thorlabs, US), the photons forward-scattered from the particles in the environment modulation system (EMS) pass through a polarization state receiver (PSR). The PSR consists of a quarter-wave plate (W2, AQWP10M-580, Thorlabs, US) and a linear polarizer (P2, LPVISC050-MP2, Thorlabs, US). Finally, the patterns can be recorded by using a CCD camera (MER-502-79U3C, Daheng Optics, China) with a pixel size of 3.45 μ m ×3.45 μ m and 2448 × 2048 pixels. In the environment modulation system (EMS), the oil smoke environment was generated by glycerol ultrasonic atomization with a refractive index of 1.47 (sample data provided by Antari Company, Taiwan). Before relative humidity testing, we measured the smoke particle size by using Malvern Spraytec (STP5311, made by Malvern Instruments Limited, Britain), shown in Section 2.2. The monitoring method of the relative humidity and its effect on the optical thickness is described in Section 2.3.

During testing, to ensure the consistency of the concentration of dry particles at each relative humidity level, we control the amount of dry particles in the EMS by using the particle filling time. In this paper, we choose a dry particle filling time of 5 s, and then the volume distribution of the particle size can be measured directly by Spraytec, which is shown in Figure 2A. To ensure suitability of the testing data for theory simulation analysis, we converted the cumulative volume distribution to the number distribution in our previous study [16]. Figure 2B shows the cumulative number distribution of oil smoke, and the oil smoke particle has a mean radius value of 0.77 μ m .

For monitoring of the relative humidity in the EMS, this paper uses ultrasonic atomization of water (ZS-10Z, Hangzhou Zhengdao Electrical Equipment Co., Ltd, China) for humidity adjustment and the relative humidity calibrator (RHCL-2, OMEGA Engineering Inc, China) for real-time measurement in the EMS. The water fog filling will stop once the relative humidity meets the testing requirement. During the experiments, the relative humidity levels of 0%, 10%, 40%, and 95% were chosen as the low, medium, and high typical relative humidity. We characterize the medium concentration by the optical thickness according to Beer’s law to learn the concentration variation at different levels of relative humidity. The optical thickness measurement schematic diagram is shown in Figure 3, in which the light power meter (LPM, PM100D, made by Thorlabs, US) is applied for outgoing light intensity testing. The optical thickness τ can be obtained from the incident light intensity I 1 and emergent light intensity I 2 , which can be expressed as [17] I 2 = I 1 ⁡ exp − τ = I 1 ⁡ exp − ρ k e l , (1) where l is the thickness of the medium in the EMS (units of cm), ρ is the medium concentration (units of #/cm³), and k e = k s + k a is the extinction cross-section (units of cm²).

Figure 4 shows the stability of the emergent light intensity as a function of sampling time for one of the groups at 40% relative humidity. The result shows that the smoke environment is stable 3–10 min after the smoke filling, so we can conduct testing during this time. When the relative humidity levels are 0%, 10%, 40%, and 95%, the corresponding optical thicknesses are 2.02, 2.08, 2.45, and 3.26, respectively. Increase in the optical thickness is due to the larger relative humidity resulting in the increase of the extinction cross-section k e in Eq. 1.

Due to the hygroscopicity of the oil smoke particles, the radius and refractive index of the oil smoke particles after water absorption can be obtained from the Hanel’s empirical growth equation [18]. When the relative humidity H r varies from 0 to 1, the particle size R and refractive index (real part M r e and imaginary part M i e ) can be calculated by the empirical growth equation and the volume weighted average method, respectively, which can be expressed as R R 0 = 1 − H r − 1 u , (2) M r e = M r w + M r o − M r w R / R o − 3 , (3) M i e = M i w M r w 2 + M i o M r o 2 + 2 − M i w M r w 2 R R o − 3 M r e 2 + 2 , (4) where M r w and M i w are the real and the imaginary parts of the refractive index of pure water particles, respectively; M r o and M i o denote the real and imaginary parts of the refractive index of the dry oil smoke particles, respectively; μ is the quality growth factor, which takes the value of 4 [19]; R 0 is the radius of the dry oil smoke particle. Figure 5 shows the variation of particle radius and the refractive index with the relative humidity. The imaginary part of the particle refractive index after moisture absorption is tiny, so it can be ignored.

A Monte Carlo method based on Ramella-Roman et al. [3] is used to perform the calculations for the multiple scattering of oil smoke particles. It involves the change of the transmission energy and the direction of the photons, which is finally characterized by polarization.

The Stokes vector uses four parameters I , Q , U , and V to describe the polarized light [20], in which I , Q , U , and V correspond to the total intensity of the light, linear polarization in the horizontal or vertical direction, linear polarization in + 45 o or − 45 o direction, and circular polarization, respectively. The Stokes vector of the light can be expressed as [21] I = E 1 2 + E 2 2 Q = E 1 2 − E 2 2 U = 2 E 1 E 2 ⁡ cos ⁡ θ V = 2 E 1 E 2 ⁡ sin ⁡ θ , (5) where E 1 and E 2 denote parallel and perpendicular components of plane wave amplitudes, respectively, and θ is the phase difference between E 1 and E 2 .

When the polarized light passes through the medium, the Mueller matrix M θ , ϕ of the medium connects the incident Stokes vector S in with the outgoing vector S out of the transmission light, which is S out = M θ , ϕ ⋅ S in , (6) I o u t Q o u t U o u t V o u t = 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 I i n Q i n U i n V i n , (7) where M θ , ϕ = R − ϕ S θ R − ϕ , R − ϕ is the rotation matrix for the Stokes vector of light transformation between the reference plane and the scattering plane.

For the spherical, isotropic, and homogeneous particles, the Mueller matrix is S θ = S 11 S 12 0 0 S 12 S 11 0 0 0 0 S 33 S 34 0 0 − S 34 S 33 , (8) and the rotation matrix can be written as R − ϕ = 1 0 0 0 0 cos 2 ϕ − sin 2 ϕ 0 0 sin 2 ϕ cos 2 ϕ 0 0 0 0 1 , (9) so M θ , ϕ can be expressed as [22] M θ = 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 = 1 0 0 0 0 cos 2 ϕ − sin 2 ϕ 0 0 sin 2 ϕ cos 2 ϕ 0 0 0 0 1 S 11 S 12 0 0 S 12 S 11 0 0 0 0 S 33 S 34 0 0 − S 34 S 33 1 0 0 0 0 cos 2 ϕ − sin 2 ϕ 0 0 sin 2 ϕ cos 2 ϕ 0 0 0 0 1 = S 11 S 12 0 0 S 12 ⁡ cos 2 ϕ S 11 ⁡ cos 2 ϕ − S 33 ⁡ sin 2 ϕ − S 34 ⁡ sin 2 ϕ S 12 ⁡ sin 2 ϕ S 11 ⁡ sin 2 ϕ S 33 ⁡ cos 2 ϕ S 34 ⁡ cos 2 ϕ 0 0 − S 34 S 33 1 0 0 0 0 cos 2 ϕ − sin 2 ϕ 0 0 sin 2 ϕ cos 2 ϕ 0 0 0 0 1 = S 11 S 12 ⁡ cos 2 ϕ − S 12 ⁡ sin 2 ϕ 0 S 12 ⁡ cos 2 ϕ S 11 cos 2 2 ϕ − S 33 sin 2 2 ϕ − S 11 ⁡ cos 2 ϕ sin 2 ϕ − S 33 ⁡ sin 2 ϕ cos 2 ϕ − S 34 ⁡ sin 2 ϕ S 12 ⁡ sin 2 ϕ S 11 ⁡ cos 2 ϕ sin 2 ϕ + S 33 ⁡ sin 2 ϕ cos 2 ϕ − S 11 cos 2 2 ϕ + S 33 sin 2 2 ϕ S 34 ⁡ cos 2 ϕ 0 − S 34 ⁡ sin 2 ϕ − S 34 ⁡ cos 2 ϕ S 33 . (10)

Photon energy loss can occur in the process of the transmission; once the photon energy weight decreases below a certain threshold or the photon flies off the boundary, the photon transmission is terminated. After the photon has been scattered n times, the energy weight W n is W n = W n − 1 ⋅ k s k s + k a , (11) where k s and k a are the scattering and absorption cross-sections of the polydispersed system, respectively, which can be obtained from k s = π ∫ r 1 r 2 r 2 Q s c n r d r , (12) k a = π ∫ r 1 r 2 r 2 Q a b n r d r , (13) in which r 1 and r 2 are the lower and upper limits of the particle radius, respectively; Q s c and Q a b are the scattering and absorption factors, respectively; and n r is the number distribution of the particle size.

The Mueller matrix can be acquired in various ways; in this paper, we use a combination of 36 groups to adjust the polarization state generator (PSG) and the polarization state receiver (PSR) parts to H , V , P , M , L , and R , respectively. The detailed combination of the 36 groups is shown in Table 1.

According to the aforementioned combinations, the Mueller matrix can be derived as [23] M = 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 = 1 2 S HH + S HV + S VH + S VV S HH + S HV − S VH − S VV S PH + S PV − S MP − S MM S RH + S RV − S LH − S LV S HH + S VH − S HV − S VV S HH + S VV − S HV − S VH S PH + S MV − S MH − S PV S RH + S LV − S RV − S LH S HP + S VP − S HM − S VM S HP + S VM − S VP − S HM S PP + S MM − S PM − S MP S RP + S LM − S LP − S RM S HR + S VR − S RL − S LL S HR + S VL − S VR − S HL S PR + S ML − S MR − S PL S RR + S LL − S RL − S LR , (14) where S indicates the light intensity and the subscripts H , V , P , M , L , and R are the corresponding polarization states described in Table 1.

With the polar decomposition method, the Mueller matrix M can be decomposed into the product of the depolarization matrix M Δ , retardance matrix M R , and diattenuation matrix M D [24], which is M = M Δ M R M D , (15) in which M R and M D represent the transformation process of polarizations from one type of fully polarized light to another one, as well as M Δ indicates the conversion of fully polarized light to partially polarized light [25,26]. M ′ contains only the retarder and depolarizer, which can be written as M ′ = M M D − 1 = M Δ M R = 1 0 → T P → Δ m Δ 1 0 → T 0 → T m R = 1 0 → T P → Δ m Δ m R = 1 0 → T P → Δ m ′ , (16) where P → Δ is the polarizance of the depolarizer and m ′ = m Δ m R , m Δ , and m R are the 3 × 3 submatrix of the depolarization matrix and retardance matrix.

We set λ 1 , λ 2 , and λ 3 as the eigenvalues of m Δ , which can be calculated as m Δ = ± m ′ m ′ T + λ 1 λ 2 + λ 2 λ 3 + λ 3 λ 1 I − 1 × λ 1 + λ 2 + λ 3 m ′ m ′ T + λ 1 λ 2 λ 3 I . (17)

The total depolarizer coefficient can be written as Δ = 1 − t r m Δ 3 = 1 − t r M Δ − 1 3 , 0 ≤ Δ ≤ 1 , (18) where M Δ can be obtained from Eqs 16, 17.

Our goal in this work is to investigate the effect of relative humidity on the polarization Mueller matrix. Since the empty EMS without particles has little impact on the Mueller matrix of the oil smoke environment, we chose dry particles and low, medium, and high typical humidity for the study. The Mueller matrix patterns of oil smoke particles at different relative humidity are simulated by the Monte Carlo simulation method in Section 3.2 with the initial parameters of the wavelength of 532 nm and photon number of 10⁷. The results are shown in Figure 6, in which Figures 6A, B, C, D show the Mueller matrix simulation with each element pattern pixels of 100 × 100 at the relative humidity of 0%, 10%, 40%, and 95% respectively. We can see that the 16 elements of the Mueller matrix exhibit symmetry, which conforms to the cosine distribution according to Eq. 10. With the increasing relative humidity, the scattering pattern size of the Mueller matrix is found to increase. This is due to different relative humidity leading to variations in particle size and refractive index, which leads to a rise in the scattering and absorption coefficients of the medium and a more significant deflection of the photons scattered in the medium. Moreover, the number of photons in the pattern that deviate from the incident center point increases after transmission, which increases pattern sizes.

To further quantify the effect of different relative humidity on each element of the Mueller matrix, we first normalize the pixel intensity values of each element at different relative humidity by m₁₁ to eliminate the effects of light intensity and keep the intensity value of each pixel in the range of [−1 1] [27], and then according to the appearance frequency of the intensity values, the interval segmentation statistics of intensity values are performed for each normalized matrix element m₁₁–m₄₄ to get the number of intensity values in a small interval of different central values. Finally, the normalized intensity frequency distribution of each element is established, shown in Figure 7. In Figure 7, the horizontal coordinate is the central value of the intervals and the vertical coordinate is the number of intensity values occurring in this division interval. As can be seen from the figures, the three diagonal elements m₂₂, m₃₃, and m₄₄ change more significantly as the relative humidity increases. The depolarization phenomenon is more significant as observed from the Mueller matrix polarization decomposition in Section 3.4. This phenomenon is also consistent with the conclusion that the scattering depolarization of spherical, isotropic, and homogeneous particles is more significant than the retardance and diattenuation [28].

The depolarization coefficient pattern of Figure 6 under different relative humidity is solved by Eq. 15. Figure 8 shows the Mueller matrix depolarization coefficient patterns' three-dimensional simulation results for 0%, 10%, 40%, and 95% relative humidity. The figures show that the minimum value of the oil smoke depolarization coefficient pattern increases gradually with increasing relative humidity, which are 0.426, 0.511, 0.565, and 0.816, respectively. This indicates that the depolarization phenomenon becomes apparent as the relative humidity of the oil smoke increases.

The PSG and PSA parts are adjusted in the measurement according to the data in Table 1, and 16 sets of Mueller matrix element patterns can be obtained by calculating Eq. 14. To correspond to the image size of the Mueller matrix elements in the simulations, the raw patterns captured by the CCD with the size of 2448 × 2048 are cropped to the 500 × 500 pixel images with the center position as the original point. Based on the simulation analysis and our previous research for spherical, isotropic, and homogeneous particles [10], in this paper, we focus on the pattern variations of the three matrix elements m₂₂, m₃₃, and m₄₄ related to depolarization. Figures 9–12 show the testing results under the relative humidity of 0%, 10%, 40%, and 95%, respectively. From the results, we can see that with the increasing relative humidity, the patterns calculated by Eq. 14 keep getting more prominent, which is consistent with the trend of simulation results. Figure 13 shows the schematic diagram of the quantitative analysis process, the edge points of the m₂₂, m₃₃, and m₄₄ patterns in the horizontal leftward and rightward and vertical upward and downward directions are used as reference points, and the absolute values of the pixel intensities in the four directions are summed for the m₂₂, m₃₃, and m₄₄ patterns under the relative humidity of 0%, 10%, 40%, and 95%. According to calculation, at the four relative humidity levels, the absolute sum of the pixel intensities of m₂₂ in four directions are 1.09, 1.12, 2.46, and 2.81, for m₃₃ are 0.86, 0.92, 0.99, and 1.34, and for m₄₄ are 0.44, 0.61, 1.11, and 1.71, respectively. It is further verified that as the relative humidity increases, the photon motion becomes more deflected with the rise of the medium scattering coefficient, leading to the more significant deviation of the pattern depolarization terms m₂₂, m₃₃, and m₄₄ from the central point of propagation.

To verify the accuracy of the simulations, we compared the simulation result of Mueller matrix elements with the experimental images, and the mean square error (MSE) was calculated [29]. Each Mueller matrix element’s simulated and measured patterns under the same relative humidity are subjected to a point-taking operation. During the calculation, the center pixel is used as the center of the circle, and the points are taken at 45° intervals in the 360° direction of the circle. To make the values consistent, the edge points of the pattern are taken in these eight directions to complete the values compatible with the simulation and experiment. Finally, the normalized intensity values of the taken points are calculated for the MSE, which can be expressed as M S E = 1 n ∑ 1 n R − R m 2 , (19) where R and R m are the simulated and measured normalized intensity values, respectively, and n = 8 is the number of samples at different directions. Table 2 shows the MSE calculation of simulation and measurement for matrix elements m₂₂, m₃₃, and m₄₄ at the relative humidity of 0%, 10%, 40%, and 95%.

It is apparent that the average value of the MSE with different relative humidity levels is 0.00538. This is due to the fact that the relative humidity controlled in the experiment cannot reach the ideal state, resulting in the simulation and experimental patterns not being fully consistent. However, the simulation and the observed images have the exact change trend, which can indicate the correctness of the simulation.