Speckle-based determination of the polarisation state of single and multiple laser beams

Laser speckle is generated by the multiple interference of light through a disordered medium. Here we study the premise that the speckle pattern retains information about the polarisation state of the incident field. We analytically verify that a linear relation exists between the Stokes vector of the light and the resulting speckle pattern. As a result, the polarisation state of a beam can be measured from the speckle pattern using a transmission matrix approach. We perform a quantitative analysis of the accuracy of the transmission matrix method to measure randomly timevarying polarisation states. In experiment, we find that the Stokes parameters of light from a diode laser can be retrieved with an uncertainty of 0.05 using speckle images of 150×150 pixels and 17 training states. We show both analytically and in experiment that this approach may be extended to the case of more than one laser field, demonstrating the measurement of the Stokes parameters of two laser beams simultaneously from a single speckle pattern and achieving the same uncertainty of 0.05.


I. INTRODUCTION
When coherent light is diffused by a disordered medium, it produces a typical granular pattern called speckle. Despite their random and uncontrollable nature, speckle patterns encode information about both the diffuser and the light, and can therefore be used to perform a range of measurements [1]. Two approaches are possible: if one considers the incident light to be timeinvariant, the speckle pattern can be harnessed to probe properties of the diffuser. This is the dominant idea in speckle metrology. This has been applied to many types of measurements, such as displacement [2][3][4][5], vibration and sound [6,7], blood flow mapping in tissues [8], among many others. Another approach is to consider the diffuser to be constant in time, in which case the speckle pattern can be harnessed to probe properties of the incident light. This more recent concept has been applied to the measurement of wavelength variations and laser stabilisation [9][10][11][12][13][14][15][16], spectroscopy [17][18][19][20], and transverse mode characterisation of structured light [21,22].
Measurement of polarisation of a light field is a key requirement in a breadth of photonics applications, and studies related to polarisation measurement with speckle have also been carried out. A previous study derived an explicit expression for estimating the Stokes parameters of a laser beam, in terms of the cross-correlations between four particular speckle patterns corresponding to the four classical polarisation filters [23]. Later studies included a generalised approach using Jones-like transmission matrices [24], and Mueller-like transmission matrices allowing spatially resolved polarimetry [25], and spectropolarimetry [26]. * mf225@st-andrews.ac.uk In this paper we focus on quantitative analysis of the transmission matrix method performance when applied to polarisation measurement. The measurements are performed on randomly time-varying polarisation states of the incident field. Importantly we also extend the method to the case of multiple beams, where the polarisation state of multiple beams can be measured simultaneously from one single speckle pattern. We demonstrate this for the case of two light fields. We also provide a demonstration verifying the linearity between the Stokes vector of the input beam and the resulting speckle pattern, and extend this result to the case of multiple beams. Such an approach may have applications in optical telecommunications, optical manipulation of birefringent particles and polarisation microscopy.

II. BACKGROUND
Changing the polarisation state of a laser beam, for example by manually rotating a waveplate, induces a visible change in the speckle pattern produced after diffusion. In this section we derive an expression for the linearity that exists between the polarisation state and the speckle pattern, which is at the core of the transmission matrix method.
We consider the geometry displayed in figure 1 where an input beam of spatially constant polarisation state is diffused by a single reflective surface, which we model as an assembly of discrete elements. The light wave is described by its electric component E, which is a 1 × 3 complex vector. We assume that the electric field incident upon point j of the observation plane after being diffused by the ith element of the diffuser is of the form where E i is the electric field at point i, and α ij is a com- plex tensor containing all the details of the field's transformation from i to j, which includes the contribution of the diffusion itself, as well the transport from i to j [24]. This relation expresses the linear nature of the diffusion process, which is our assumption. The total field at j is the sum of the fields propagating from all the illuminated points of the diffuser, which gives As the input beam has a spatially constant polarisation state, the field incident upon i can be written as E i = eρ i e iφi+iϕ(t) , where ρ i is the amplitude of the electric field at point i, φ i is the spatial part of the phase at point i, ϕ(t) is the temporal part of the phase (containing the ωt term and any kind of temporal fluctuations, common to all i), and e is the normalised Jones vector. The resultant field at j is then where α j contains all the details of the field's transformation terminating at point j.
It is established that when the field undergoes a linear transformation, as the last relation indicates, the Stokes vector also undergoes a linear transformation [27]. We show later in this paper that we can extend this to the case of multiple beams. First, we express the coherency matrix, which is defined as C = E ⊗ E * , where the brackets denote time averaging, ⊗ the outer (or Kronecker) product, and * the complex conjugate. Using the fact that the outer product can be expressed as a usual matrix product when the vectors are properly shaped, the coherency matrix at point j is with the prime denoting transposition, and C 0 the coherency matrix of the input beam normalised to the intensity (as it is computed from the normalised Jones vector). Note that, by definition, the electric field in the expression of the coherency matrix is no longer a 1 × 3 vector, but a 1 × 2 vector, expressed in the plane perpendicular to the direction of propagation. As this direction is arbitrary (the observation plane can be anywhere), a change of basis in which the electric field is expressed is required. The change of coordinate, as well as the cropping of the electric field, can be performed by a matrix multiplication of E j , which we implicitly absorb in α j for clarity. Now that we know how the coherency matrix transforms, we wish to determine how the Stokes parameters subsequently transform. The Stokes parameters are related to the coherency matrix via the relation S m = T r(Cσ m ), where T r is the trace, S m is the mth Stokes parameter, and σ m is the mth Pauli matrix. This operation is analogous to a projection of the coherency matrix onto the three Pauli matrices (to which is added the unit matrix), as the Stokes parameters are defined as twice the coefficients of decomposition of the coherency matrix in this basis, expressed as C = n 1 2 S n σ n [27]. Applying these two relations at point j we have Finally, we find that the Stokes vector at point j is linearly related to the Stokes vector of the input beam S, through a 4 × 4 matrix M j , the elements of which are given by M j,nm = 1 2 T r(α j σ n α * j σ m ). In other words, M j is the Mueller matrix associated to the diffuser at point j. Note that S is the normalised Stokes vector, as it is computed from the normalised coherency matrix.
As a camera only records the intensity, given by the first Stokes parameter S 0 j , only the first column of M j is needed, and the intensity observed at point j is given by I j = SM j,1 , with M j,1 the first column of M j .

III. METHOD
The last relation found above can be extended to any set of points on the observation plane, and leads to the following central relation with I the 1 × L image of the speckle pattern reshaped into a row vector, and S the 1 × 4 Stokes vector of the input beam. M is a 4 × L matrix making the connection between the two, and is usually referred to as transmission (or transfer, or measurement) matrix [25,26,28]. The transmission matrix is unknown and depends on many parameters, such as the spectrum, beam profile, and angle of incidence of the input beam, its position on the diffuser, and the detailed structure of the diffuser. However, assuming all those conditions to be timeinvariant, we can determine the transmission matrix using a set of N (with any N ≥ 4) training polarisation states and their corresponding speckle patterns. Applying (6) to the training sets and performing a simple matrix inversion leads to with I 0 a N × L matrix containing the training speckle patterns stacked in rows, and S 0 a N × 4 matrix containing the corresponding Stokes parameters (S + 0 being the 4 × N pseudo-inverse of S 0 ). Here the pseudo-inverse is used, rather than the standard inverse, because S 0 is not square. Indeed, four training states would be enough to close the system, but in practice better results are obtained using more training states (see section 4 and 5), in which case S 0 is rectangular. With N > 4, expression (7) is an over-determined system, and solving by (8) corresponds to the minimisation of the Euclidean (or Frobenius, or L2) norm I 0 − S 0 M 2 , defined by Armed with M , we now have all we need in (6) to find the Stokes parameters of any new polarisation state, given its corresponding speckle pattern. By solving for the Stokes vector we find As a passing remark, one can see from this last relation that what is needed in practice is M + and not M . To avoid unnecessary interim calculations, one can directly compute M + from the training sets given by

IV. EXPERIMENTAL IMPLEMENTATION
Our setup is shown in figure 2. A laser beam of 780 nm wavelength (7 mW power and 1 MHz linewidth) is initially set to a highly stable linear polarisation state by means of a 40 dB optical isolator (Isowave I-80-T-5-H). It then passes through three successive waveplates (respectively quarter wave, half wave, and quarter wave), mounted on independent motorised rotating stages (Thorlabs KPRM1E/M), to allow generation of an arbitrary and dynamically controlled polarisation state. The beam is then split into two paths using a nonpolarising beam splitter. One path leads to a commercial polarimeter (Thorlabs PAX1000IR1/M), and one path leads to a rough surface where the light impinges at 45 • and is diffused to form a speckle pattern directly on the camera (Mikrotron MotionBLITZ EoSens mini2) without intermediate lenses. The rough surface is a 12.5 mm-diameter, 1 mm-thick disk of a PTFE-based material with high reflectivity and highly Lambertian reflectance in the 250 -2500 nm wavelength range (Thorlabs SM05CP2C). The distance between rough surface and camera is chosen to be 12 cm, so that individual grains in the speckle pattern cover approximately 15×15 pixels. We simultaneously record the polarisation measured by the commercial polarimeter and the speckle pattern. An example of an obtained speckle pattern is shown in figure 2. In order to test the method described above, we prepared the laser beam in a random and continuously time-varying state of polarisation by rotating the waveplates with incommensurate angular speeds (5.77 • s −1 , 10.77 • s −1 , and 20.77 • s −1 ). As the plates rotated, we simultaneously recorded the measurements of the commercial polarimeter and 150×150-pixels images of the speckle patterns, at regular intervals of 0.2 s. We picked 17 states at the beginning of the time series to make up our training sets S 0 and I 0 respectively in (7). Once M + was determined, we estimated the Stokes parameters of the subsequent states by applying equation (9) to the speckle patterns and compared to the measurements of the commercial polarimeter. We show the results in figure 3, and give a quantitative analysis of the uncertainty in terms of image size and number of training images in figure 4.
We define the uncertainty on the Stokes parameters retrieval as the standard deviation of the residuals, given by the difference between the measurements of the commercial polarimeter and the estimation from the speckle patterns. We find that the uncertainty rapidly reaches a value of 0.05 for about 15 training images and an image size of 100×100 pixels. For comparison, the resolution of the commercial polarimeter is 0.01.
It is worth pointing out that what is retrieved from the speckle patterns is the polarisation of the light that is incident upon the polarimeter, not necessarily that of the light incident upon the diffuser. These may differ due to the reflection in the beam splitter, and are linearly related by the Mueller matrix associated with the reflection. For any unknown beam, the polarisation retrieved from the speckle patterns is the one that would be measured by the polarimeter. If this is of any importance in a given application, for example if one needs to determine the polarisation of the light incident upon the diffuser, one would simply have to determine the Mueller matrix of the beam splitter. That task would be strictly analogous to what is done in the method section, where the unknown matrix M would be the Mueller matrix of the beam splitter, and the two linearly related vectors would be the polarisation states of the two paths.

V. MULTIPLEXING
As the method relies on the acquisition of speckle images, it can be extended to the simultaneous measurement of multiple beams. Indeed, if multiple beams are diffused on the same surface, all the corresponding speckle patterns superpose on the camera. If, in addition, the different beams originate from different sources, then the speckle patterns superpose without interference. It follows that the polarisation state of each beam is linearly encoded in the resulting speckle pattern, analogously to the single-beam case, as the polarisation state of each beam is already linearly encoded in its own speckle pattern. In this section we show that a relation equivalent to (6) holds in the case of two beams, allowing the same Trajectory of the polarisation state across the Poincaré sphere from t=160 s to t=190 s, measured by the commercial polarimeter (black) and retrieved from the speckle patterns (red). (c-e) The Stokes parameters S1 to S3 as a function of time, measured by the commercial polarimeter (black) and retrieved from the speckle patterns (red). (f) The error as the absolute residual, was averaged over the Stokes parameters. The estimation was performed using 150×150-pixels images and 17 training states. The uncertainty is given by the standard deviation of the residuals. It is shown as a function of the number of training images (four being the minimum required), for different image sizes ranging from 20×20 to 150×150 pixels. We see that the uncertainty reaches a minimum of 0.05 after about 15 training images and an image size of 100×100 pixels.
method to be applied, and we test it experimentally. The electric field is now a superposition of two fields The electric field at point j is then E j = e 1 α 1,j e iϕ1(t) + e 2 α 2,j e iϕ2(t) , the alpha tensor being different for each beam, as they may hit different diffusers with different phases and amplitudes. Injecting this in the expression of the coherency matrix at point j gives where the cross terms are zero as the beams are not coherent with each other, C 1 and C 2 are the normalised coherency matrices of each individual input beam. Now expressing the Stokes vector at point j we have where M 1,j and M 2,j are the Mueller matrices associated to point j for each beam, and S 1 and S 2 are the normalised Stokes vectors of each beam. Again, as only the intensity is observed on the camera, the intensity at point j is given by I j = S 1 M 1,j,1 + S 2 M 2,j,1 , where M 1,j,1 and M 2,j,1 are the first column of M 1,j and M 2,j . I j can actually be expressed in one single dot product as where the two Stokes vectors S 1 and S 2 are concatenated in one 1 × 8 vectorS, and M 1,j,1 and M 2,j,1 are concatenated in one 8 × 1 vectorM j . Generalising again to a set of points on the observation plane, we finally find the two-beams equivalent of relation (6): where I is the 1 × L speckle image,M is the two-beams 8 × L transmission matrix. As this relation is of the same form as relation (6), it implies that all the analysis performed in the method section can be applied in the exact same way. The only difference is that the inversion (9) gives both Stokes vectors in a single 1 × 8 vector. In principle this approach can be extended to more beams, a mathematical limit being one quarter of the number of pixels in the speckle image, and a physical limit being contrast reduction as the number of beams increases.
To test this relation, we added a second laser beam to the setup described in figure 2, of the same wavelength and power as the first one. This is also prepared in an initial high purity linear polarisation state by means of a 60 dB optical isolator (Toptica DSR780). We inserted a quarter waveplate in a fixed arbitrary orientation on the path of the second laser, so that its polarisation state was different than that of the first laser before entering the three rotating waveplates. This ensured that both beams hit the diffuser with different time-varying polarisation states. Also, we needed to know the polarisation state of each individual beam at any given time, which is not possible when they enter the commercial polarimeter simultaneously. Therefore we rotated the waveplates by small discrete increments, between which the waveplates were left stationary for 5 seconds. During those 5 seconds, three measurements were made: each beam was sequentially blocked while the polarisation state of the other one was measured by the commercial polarimeter, and an image of the speckle pattern was recorded when both beams hit the diffuser. This way of proceeding took an extended period of time, implying a timescale difference with our single-beam experiment and less data points. For consistency with the single-beam experiment, we also picked 17 training states at the beginning of the time series and applied the retrieval to the subsequent states, using 150×150-pixels images. We show the results in figure 5, and the uncertainty analysis in 6. Interestingly, we find that the uncertainty reaches a minimum after approximately the same number of training images than in the single-beam case, i.e. 15 training images. Stokes parameters S1 to S3 are given as a function of time, measured by the commercial polarimeter (black) and retrieved from the speckle patterns (red and blue, one for each beam). The estimation was performed using 150×150-pixels images and 17 training states.
The method could be applied to beams of different and time-varying powers, without needing to train for different powers. In equation (6) and (13)  by the power, then the transmission matrix found by (8) would be power-independent. It follows that the estimation (9) would give the non-normalised Stokes vector, for any input power. In the multiple-beams case, each training Stokes vector ofS in (13) would need to be individually multiplied by the power of its corresponding beam, and again the estimation would provide the nonnormalised Stokes vectors, for any input powers.

VI. ACQUISITION SPEED AND REGULARITY
Another advantage to our speckle-based polarisation measurement technique is that, when performed with a fast-framing camera, it allows a higher sampling rate compared to commercial polarimeters based on mechanically-rotating polarisers.
The Thorlabs PAX1000IR1/M we use as a benchmark can theoretically achieve a sampling rate of 400 Hz, but in practice we found it limited to 110 Hz, with an irregular sampling. In this section we explore the high speed capability by applying faster polarisation changes.
We performed the same single-beam experiment as in section 4 but replacing the three waveplates by an electrooptic modulator (EOM, Thorlabs EO-AM-NR-C4), allowing a very rapid, electrically tunable modulation of polarisation. We applied a periodic modulation by applying a sinusoidal voltage to the EOM, and chose its frequency so that the commercial polarimeter made at least 10 measurements per period of modulation. Below that number of points per period, the undersampled waveform resembled random noise. As the sampling rate of the commercial polarimeter is at maximum 110 Hz, we applied a modulation of 10 Hz, leading to 11 points per period. While the polarisation state of the beam was modulated, we simultaneously recorded the measurements of the commercial polarimeter and the speckle patterns, with acquisition rates of 110 Hz and 1000 Hz respectively. We compare the measurements in figure 7.
We further explored the speed capability by increasing the modulation frequency to 500 Hz, which is above the maximum acquisition rate of the commercial polarimeter but still clearly visible when retrieved from the speckle patterns using a camera frame rate of 5000 Hz. This was the maximum achievable frame rate in our setup before the intensity of the speckle pattern became too low. We show the corresponding measurements in figure 7.
Although this sampling rate is about 50 times higher than that of the commercial polarimeter, it is worth pointing out that it is still much lower than what can be achieved with polarimeters based on intensity measurements after four polarising elements, sometimes called four-detector polarimeters [29]. In that case the sampling rate is that of the photodiodes used, which can reach GHz. The Stokes parameters S1 to S3 as a function of time, measured by the commercial polarimeter (black) and retrieved from the speckle patterns (red), when a 10 Hz modulation is applied. The acquisition rates are respectively 150 Hz (on average) and 1000 Hz. Left column: The Stokes parameters as a function of time retrieved from the speckle patterns, when a 500 Hz modulation is applied. At this point the modulation is no longer visible on the commercial polarimeter.

VII. SUMMARY AND CONCLUSION
A one-to-one relation exists between the Stokes vector of a laser beam (describing its state of polarisation) and the speckle pattern it produces after diffusion. This relation takes a simple linear form, which we derived based on a linear diffusion model. We also derived an analogous relation in the case of two beams which are not coherent with each other. It involves a transmission matrix, which is unknown but can be determined through a training stage. Exploiting such a linear relation in a measurement purpose constitutes what is usually referred to as a transmission matrix method. In our case, the method essentially transfers the measurement process from the polarimeter to the camera, the knowledge being passed on during the training stage. Then the polarimeter is no longer needed and the polarisation can be retrieved from the speckle patterns alone, which unlocks several advantages.
The main advantage we investigated is the possibility of multiplexing, i.e. measuring the polarisation state of several beams simultaneously from one single image. Another advantage we explored is the possibility of higher acquisition rate and more regular sampling than commercial polarimeters. We achieved a sampling rate of 5000 Hz, which is about 50 times higher than commercial available polarimeters. We demonstrated that the Stokes parameters of one and two beams could be retrieved with an uncertainty of 0.05, using in both cases 17 training states and 150×150-pixels images, the typical resolution of a commercial polarimeter being 0.01.
Multiplexing is interesting in that it allows the embedding of information from several beams into one single image. This may find applications in optical imaging and manipulation of multiple birefringent particles. For example, in the field of levitated optomechanics, optical binding has been studied with two vaterite microparticles [30]. Our approach would allow the detailed analysis of the polarisation change of the light field from each particle, and hence the particle rotation and dynamics, in a facile, informative manner.