Simulation of polarized optical speckle fields: effects of the observation scale on polarimetry

: In this paper, we propose the simulation of polarized speckle ﬁelds using the Stokes formalism, which allows the description of partially polarized electromagnetic waves. We deﬁne a unique parameter which determines the partial decorrelation of the involved ﬁelds, allowing to simulate the polarized speckles produced by all types of scatterers, from simple to multiple scatterers. We validate this model by comparison with experimental measurements. We use that simulation model to study the impact of the imaging device parameters on polarimetric measurements: ﬁrst we emphasize a limit of resolution on retardance measurements, then we study the spatial depolarization, which appears when an observer is measuring any space-variant polarization map.


Introduction
Simulating the effects of speckle fields present an important interest in various domain of physics and applications, in addition to polarimetric imaging: metrology [1], telecommunications with polarization coding [2,3], biomedical imaging [4,5], remote sensing [6].
Usually, speckle fields are simulated through their scalar intensity, assuming that the fields are uniformly polarized. However in most of the case, there is multiple scattering leading to complex spatial polarization variation. These variations, which lead to a partial decorrelation of the fields [7], have to be taken into account if one wants to simulate the scalar intensity of a speckle field. Several approach have been used to simulate polarized speckle fields, for example based on the Jones formalism [8], or on multiple random wave superposition [9,10].
In that paper, we propose a phenomenological way to simulate polarized speckle fields, based on the plane wave description [11] of polarimetric states introduced by Gabriel Stokes (1852). That description is still valid for partially coherent light and allows us to set the type of scattering continuously from a single to a multiple scatterer. We compare the simulation results with experimental results for surface and multiple scattering [12,13,14,15]. We then study the impact of the imaging point spread function on polarimetry. Finally we propose a study of the spatial depolarization [16,17,18,19], which appears because of the spatial sampling by an observer of a continuous space-variant polarization map.

Independent speckle field simulation
In that paper, we focus our attention on subjective speckle fields i.e. speckle fields imaged by a lens. Moreover, we consider an imaging system with a circular aperture, and only in case of a monochromatic Gaussian and fully polarized illumination. Finally, we study fully developed speckle, produced by a large number of random phasors which phases are uniformly distributed between [0, 2π].
The simulation of that type of speckle is resumed in Fig. 1. We assume that the imaging device is in the far field region relatively to the sample, thus the Fraunhoffer diffraction equation can be applied and the field incoming on the pupil plane is the Fourier transform (FT) of the illumination beam phase shifted by the sample phasor map (we consider a sample with homogeneous reflectance). After the truncation by the pupil transmission function, we observe at the focal plane of the imaging system the Fourier transform of the pupil field. We refer to the coordinates in the object plane with the subscripts x o , y o ; x p , y p in the pupil plane and x i , y i in the image plane. The complex amplitude in the object plane E E E o o o is: With A 0 the Gaussian amplitude of the beam, and ϕ the random phasor map, uniformly distributed on [0, 2π]. We assume that the distances between each transverse plane are large enough to be in a far field approximation, i.e. the Fraunhoffer approximation is valid. Thus, the complex amplitude in the pupil plane E E E p p p (x p , y p ) is: With pup the transmittance function of the circular pupil, assumed to be 1 on the pupil surface, and 0 elsewhere. Finally, the complex amplitude imaged in the image plane We remind that the Point Spread Function (PSF) of the simulated imaging device is linked to the ratio of the wavelength over the pupil diameter. The diameter of the pupil, and so the PSF size, will determine the size of the imaged speckle grains, for a given wavelength.

Polarized speckle fields simulation
We can describe any State Of Polarization (SOP) by a sum of two plane waves perpendicularly polarized, and a depolarized intensity. Thus, if we consider a polarized wave propagating along the z axis, its SOP can be described by the following Stokes vector [11]: With I the total intensity (polarized or not), Q the amount of intensity polarized along the horizontal (resp. vertical) if Q > 0 (resp. Q < 0), U the amount of intensity polarized in the ±45 • direction (depending on the sign of U) and V the amount of intensity polarized along a right handed (resp. left handed) circular SOP if V > 0 (resp. V < 0); E E E i i ix x x the complex amplitude in the image plane of the first plane wave linearly polarized along the x axis, E E E i i iy y y the complex amplitude in the image plane of the second plane wave linearly polarized along the y axis, φ their relative phase shift, and C 2 a constant corresponding to the depolarized intensity [15]. The Degree Of Polarization (DOP) of this SOP is: Substituting Eq. (4) into Eq. (5) gives: Thus, in order to simulate a polarized speckle field, one have to determine the parameters E E E i i ix x x , E E E i i iy y y , C and φ , we will see in the following how to determine them.
By looking at the I parameter in Eq. (4), one can see that the polarized intensity is obtained by the incoherent sum of two plane waves intensities |E E E i i ix x x | 2 and |E E E i i iy y y | 2 polarized in orthogonal directions. Thus, in order to simulate a polarized speckle field, we can sum two independent speckle fields, respectively uniformly polarized in orthogonal directions.
One can notice that the convention for the plane wave decomposition is defined by Eq. (4): in the convention chosen here, the x and y directions are respectively corresponding to the ±45 • directions of polarization. We choose that convention in accordance with the simulated illumination SOP (chosen to be vertical), in order to produce two orthogonal plane waves with equivalent mean intensity levels (i.e. due to the projection of a vertical SOP on the ±45 • directions). In that case, the SOP variations are only produced by modifications of the relative phase shifts φ , and do not depend of a variation of the mean level of the intensities |E E E i i ix x x | 2 and |E E E i i iy y y | 2 .
Moreover, experimental measurements have shown that the more the scattered light comes from multiple scattering, the more the SOP are spreading on the Poincaré sphere surface, centered around the illumination SOP [12] [13]. Thus, the effect of the regime of diffusion (from surface to volume) can be modelized through the standard deviation of the relative phase shifts φ , around a mean phase value corresponding to the illumination SOP: π for a vertical incident SOP in the convention chosen here.
We notice that one could choose to modelize any other illumination SOP, in that case in order to be in the same configuration than the model presented here, he would have to use the plane waves convention in which the mean intensities |E E E i i ix x x | 2 and |E E E i i iy y y | 2 are equivalent, in order to ensure that the SOP variations are only due to relative phase shifts variations.
The variations of φ produce, after the fields propagation, a (spatial) partial decorrelation between the two speckle fields intensities |E E E i i ix x x | 2 and |E E E i i iy y y | 2 . The spatial correlation of these intensities is depending on the standard deviation of the relative phase shift variations.
In order to introduce a phase correlation between the two speckle fields, we first generate the phasor map of one of the speckle fields in the object plane, e.g. ϕ x (x o , y o ), with a uniform random distribution between 0 and 2π. The phasor map of the second field in the object plane E E E o o oy y y is obtained with: With G the normal distribution with mean ρ and standard deviation σ . One can notice that the mean value ρ lead to a uniform relative phase shift (a piston), e.g. produced by a uniform retarder. Thus, the degree of correlation between the two phasor maps is only depending on the standard deviation value σ .
At this point, we have simulated E E E i i ix x x (x i , y i ) and E E E i i iy y y (x i , y i ) obtained respectively by propagating the phasor maps ϕ x (x o , y o ) and ϕ y (x o , y o ), as described in section 2. Moreover, we compute φ , the relative phase shift in the image plane with: In the following of that paper, we neglect eventual experimental noises. As we modelize a monochromatic illumination, elastic scattering and temporally stable fields, there is no effects of spectral depolarization. We first neglect the spatial depolarization, i.e. C 2 = 0 and DOP=1 (see Eq. (6)). Thus, the [I Q U V] parameters of the polarization map can be computed with Eq. (4) in function of the modelized E E E i i ix x x , E E E i i iy y y and φ parameters. We will see in section 6 how to estimate the value C 2 when the spatial depolarization is no longer negligible.
One can notice that the adjustable parameters of the simulation are the characteristics of the imaging device i.e. the diameter of the entrance pupil, which will determine the speckle size, and the correlation between the two speckles depending on the standard deviation σ in Eq. (7), which determine the mode of diffusion. As explained above, the particular cases are for σ = 0 which represent a perfect surface scatterer, and σ >> 2π which represent a perfect multiple scatterer. Between that interval, all types of scatterers are described.

Validation of the model
We study experimental results obtained using SOPAFP method [13,14]. These results were produced by a black lambertian scatterer under a 532 nm Single Longitudinal Mode (SLM) illumination, vertically polarized. Using SOPAFP method, we measure the 4 parameters of the field scattered by the sample, using the convention defined by Eq. 4 for the plane wave decomposition (i.e. x and y directions are corresponding respectively to the ±45 • directions). The measured mean relative phase shift and its standard deviation value are mean(φ ) = 2.90 rad and std(φ ) = 0.73 rad. As the measured field is obtained by the surface scattering of a vertical incident SOP, one can notice that we should obtain a mean relative phase shift of π: we observe a global phase shift of π − 2.90 = 0.24 rad. This piston value is coming from the incidence and detection angles relatively to the sample interface, which have to be taken into account in the simulation. We determine the input parameters of the model in order to obtain simulation results comparable with the measurement: • We adjust the diameter of the simulated pupil in order to obtain roughly the same PSF size that our imaging system, and so the same speckle size in the image plane.
• We set the mean phase value ρ in Eq. (7) to 2.90 rad, in order to take into account the effect of the incidence and detection angles in the model.
• We adjust the input parameter σ in Eq. (7) for the simulation of the phase decorrelation in the object plane in order to obtain the same standard deviation value than experimentally in the image plane (std(φ ) = 0.73 rad). This value is roughly obtained with σ = π/6 rad.
The simulation and measurement results are represented in Fig. 2. We choose to represent the polarization maps with a RGB representation, defined by R = |Q|, G = |U|, B = |V |. One can see that the spatial variations of polarization and the SOP probability in the Poincaré sphere (Figs. 2(a)-2(b)), are equivalent. Moreover, the simulated (Fig. 2(c)) and measured ( Fig. 2(d)) phase spatial variations are equivalent at the speckle grain scale. Finally, one can see that the probability density function of the simulated phase in the image plane (Fig. 2(c), red) is close to the measured one (Fig. 2(c), blue), with a high cross correlation coefficient R 2 = 0.998. Thus, the simulation leads to results in accordance with experimentation, in terms of SOP statistics and spatial variations.
We simulate polarized speckle fields produced by surface and multiple scattering. The results are presented in Fig. 3 and 4. In case of surface scattering, the decorrelation between the orthogonal plane waves is simulated in the object plane through a normal distribution (Eq. (7)) with a σ = π 12 standard deviation. This deviation value is chosen arbitrarily. In case of multiple scattering, simulated with σ = 10 9 , the normal distribution G in Eq. (7) is quasi uniform on [0; 2π]. Thus there is no correlation between ϕ x and ϕ y and the relative phase shift in the image plane φ is uniformly distributed between 0 and 2π. One can see that we obtain simulation results in accordance with previous theoretical or experimental speckle fields studies (e.g. in [9,10,13,14,18]), in terms of polarization singularities and SOP statistics.

Impact of the PSF size on polarimetry
The aim of polarimetry is to measure the retardance and diattenuation produced by a sample. The spatial limit of resolution for diattenuation measurements is the Rayleigh criterion, as the diattenuation is depending on the spatial variations of the two plane waves intensities (|E E E i i ix x x | 2 and |E E E i i iy y y | 2 ). However, the retardance is defined as a spatial variation of the relative phase shift φ between these two plane waves.
Here, we study the effect of the PSF size on retardance measurements. We can simulate the effect of retardance by adding the corresponding phase variation in the object plane. Equation (7) becomes: With ψ the sample retardance map. We simulate the scattering by a given surface scatterer, with σ = π/12. We use two maps of retardance, spatially varying at different frequencies, but imaged with the same parameters and so the same PSF size. The spatial frequency of the retardance is 10 times higher in Figs. 5(d)-5(f) than in Figs. 5(a)-5(c). One can notice that the SOP variations in Fig. 5(c) are quasi only due to retardance variations ( Fig. 5(a)). However, the SOP are completely random in Fig. 5(f) when the spatial frequency of the retardance reaches a given value (Fig. 5(d)), similarly than in case of multiple scattering (see Fig. 4(c)). That effect can be explained by the fact that if the retardance ψ is spatially varying from 0 to 2π on a lower scale than the PSF size, it appears a loss of the phase correlation between the two orthogonal plane waves. In that case, the surface scatterer seems to have the same polarimetric behavior than a multiple scatterer. Thus, we point out a polarimetric limit of resolution, that can lead to false interpretations between the retardance variations and the mode of diffusion.

Spatial depolarization estimation
Another polarimetric parameter is the DOP of the scattered fields, which is determined here by the value C 2 in Eq. (6). As we study the case of monochromatic illumination and temporally stable fields, the only contributions to the C 2 value are the experimental noises, that are neglected here, and the spatial depolarization [16,17,18,19]. The spatial depolarization appears due to the sampling on the pixel surface of a continuously varying polarization map. Indeed, as the measurable quantity is the intensity integrated by a pixel, we can only measure the sum of all the polarimetric states incoming on its surface.
We simulate that effect, by generating a highly resolved speckle map produced by a perfect multiple scatterer, with a large PSF relatively to the sampling size. We then simulate its integration by a CCD, by computing the sum of the polarimetric states on larger pixels, which lowers the DOP. The spatially depolarized intensity (I sd p ) is obtained with the following equation:   I sd p = I − Q 2 +U 2 +V 2 (10) With [I Q U V ] the sum of [I Q U V] over the sampling surface. In the following we give relative sampling sizes, 1 × 1 being the dimension of the initial, highly resolved map. We simulate the effect of its integration on 2 × 2, 10 × 10 and 70 × 70 pixels, the results are presented in table 1 and Fig. 6. Table 1. Contrast value of the speckle scalar intensity I , defined by std(I )/mean(I ), and depolarization ratio, defined by mean(I sd p )/mean(I ), in function of the spatial sampling of a given field produced by multiple scattering. We can see in table 1 that the mean spatial depolarization produced by the integration on 10 × 10 pixels is 30 times higher than with 2 × 2 pixels. In case of classical imaging, one will lower the speckle noise by using a large pupil aperture in order to sum various speckle grains on the same pixel. One can notice that this summation leads to a huge spatial depolarization: for pixels at roughly the same size than the speckle grain (70 × 70), 70% of the total intensity is depolarized due to the field summation, in case of a multiple scatterer.
One can notice in Figs. 6(d)-6(f) that the structure of the spatially depolarized intensity is not depending on the pixel size. Indeed, the pixel size will only impact the absolute value of that intensity and the highest spatial frequency observable. Thus, the depolarization map is only depending on the spatial variations of the SOP map. We define an empirical equation in order to estimate the structure of that spatially depolarized intensity I sd p : With i = 2 : 4 the index value of the normalized Stokes vector component S i , g xi (resp. g yi ) the gradient of the i th Stokes component along the x (resp. y) dimension and I p the polarized part of the speckle scalar intensity (I p = |E E E i i ix x x | 2 + |E E E i i iy y y | 2 ). We used that equation to estimate the depolarized intensity from a polarization map sampled with 1 × 1 pixels. We compare that calculated depolarized intensity with the simulated one, obtained similarly than Figs. 6(d)-6(f), i.e. by the sum of the 1 × 1 polarization map on 45 × 45 pixels. The results are presented in Fig. 7, one can see that the spatial variations of the calculated and simulated maps are similar.
As a consequence, in order to simulate any measurement of monochromatic space-variant polarization map, one should first determine the polarization map and then compute the spatial depolarization. That intensity is the value C 2 which has to be summed with the intensities of the two orthogonal plane waves in Eq. (4). Thus, the DOP of the field simulated will be the highest DOP that one will obtain with his imaging system, i.e. for a given ratio of wavelength over pupil diameter and pixel size, if one neglects the experimental noises.

Conclusion
We propose a new mean of simulating polarized speckle fields, with a criterion allowing to evolve continuously from a simple scatterer to a multiple scatterer. We validate this algorithm by comparison of the simulation results with experimental measurements performed on the field scattered by a black lambertian scatterer. Then, we use that algorithm to simulate the effects of the PSF size on polarimetric measurements.
We show that if one performs a polarimetric analysis on a surface scatterer in which relative phase shifts (e.g. retardance) are varying with an higher spatial frequency than the PSF width of the imaging system, then the scatterer appears as a multiple scatterer. Reciprocally, one can see a multiple scatterer as a sample which phase variations are spatially varying more quickly than the PSF. Thus, as the Rayleigh criterion in scalar optics, one can define a polarimetric limit of resolution, which is increasing in frequency with the decreasing of the PSF size. However we saw that for a given pixel size, when the PSF size is decreasing the spatially depolarized intensity is increasing.
Finally, we show that the structure of the spatial depolarization is not depending on the pixel size, but only on the spatial variations of polarization. The pixel size only impacts its absolute value and its spatial sampling. If one knows the polarization map, he can estimate the intensity depolarized by spatial depolarization with an empirical equation.