Eigenvalue calibration methods for polarimetry

Complex polarisation sensitive systems such as imaging Mueller matrix polarimeters are commonly calibrated using the eigenvalue calibration method. In this paper we present an extensive review of the method and an existing variant. We also introduce two more variants of the method to calibrate imaging polarimeters that use high numerical aperture optics. The calibration methods are tested using a Mueller matrix confocal microscope of high numerical aperture, and the effect of the pinhole size on the polarisation is also assessed experimentally. [DOI: http://dx.doi.org/10.2971/jeos.2012.12004]


INTRODUCTION AND OVERVIEW
Several calibration methods have been developed for polarimetry [5,6,8,7], all with the drawback of requiring a priori information such as the polarisation states produced or the angular orientation of the polarisation optics in the system.This drawback was overcome with the introduction of the eigenvalue calibration method (ECM) [1,4].The ECM has shown recently an increase in popularity and it is expected to become the preferred method to calibrate polarimeters; however, due to the development of more complex experimental layouts, it has been necessary to modify the ECM.A variant of the ECM was introduced [9] to calibrate systems configured in doublepass, yet restricted to low numerical aperture (NA) imaging systems because of the assumption it requires -it assumes reflection from a perfect mirror at normal incidence; assumption that can only be meet using low NA optics and by removing the objective lens of the imaging system.In high NA imaging systems and particularly in confocal microscopes, the longitudinal component of the light contributes significantly to the imaging process, and this contribution is disregarded when removing the objective lens.Two more variants are thus introduced in this paper to calibrate systems using high NA optics.These optical systems are routinely used in high-resolution microcopy and imaging, which provides the motivation for our work.One of the variants is based in a two-step procedure, whilst the other is on a three-step procedure.The double-pass, two-step, and three-step variants are explained next after reintroducing the ECM.

Eigenvalue calibration method
The ECM [1] is an ingenious method based in control theory and uses linear algebra extensively.The objective of the method is to find the matrices W, T and M i in the linear system depicted in Figure 1 when only D is known.To achieve this, M i is changed until the system is overdetermined and the solution becomes unique.When the linear system is a Mueller matrix polarimeter as in this case, M i represents a polarisation element such as a polariser or a waveplate and is referred to as calibration sample.To reach a solution in this kind of systems, three assumptions are necessary: the noise in the system is negligible, the calibration samples have a given Mueller matrices, and any polarisers used as calibration samples must be perfect polarisers.
To explain the method in more detail, consider the mathematical expression of the Mueller matrix polarimeter depicted in Figure 1 to be given by where M i is the Mueller matrix of the ith polarisation element used as calibration sample, D i is in general a q × n matrix with irradiance readings as entries, T is the q × 4 matrix known as instrument matrix (polarisation state analyser: PSA), and W is a 4 × n matrix formed by n generated Stokes vectors (polarisation state generator: PSG); however, for simplicity and without the loss of generality, assume that q = 4 and n = 4.
If Eq. ( 1) is left multiplied by the inverse of the measurement of air, i.e.
one obtains where C i and M i are similar matrices.This means that C i and M i have the same characteristic polynomial, and hence the same eigenvalues [12].Given that M i can be written in general as [10] where and ∆ i are the ellipsometric parameters [10], and that the eigenvalues of M i are it is possible to partially reconstruct M i from the eigenvalues of C i .The reconstruction is partial because the angular orientation of the polarisation sample θ is found later by means of an optimisation.Note that Eq. ( 2) is only valid if the Mueller matrix of air is assumed to be the identity matrix; hence the assumption of the noise being negligible.If the assumption is correct, the problem is simplified to find W from Eq. ( 3), and then use Eq. ( 2) to find T.
To find W, one needs to manipulate Eq. ( 3) and solve the resulting equation.Multiplying Eq. ( 3) on both sides by W and rearranging one obtains which is a special case of the Sylvester equation (MW − WC + K = 0 with K = 0).The Sylvester equation is used in the fields of control theory [13], neuronal networks [15], and for wavefront reconstruction in optics [16].A more complete literature review of this equation can be found in [14].Common methods to solve Eq. ( 5) are Schur decomposition with backward substitution [17], diagonalisation [18], and column stacking [20].The ECM uses the last two methods in the following way.
In order to solve Eq. ( 5) using the diagonalisation and the column stacking method, it is necessary to first introduce the vec operator [20].This is an operator that stacks the columns of a n × m matrix to form a 1 × n * m vector.A matrix multiplication can be rewritten using this operator as where J, K, and F are arbitrary matrices, T represents transpose, and ⊗ denotes Kronecker product [19,21,22].It is then possible to use Eq. ( 6) to rewrite Eq. ( 5) as [23] ( where At this point W can be located anywhere in the 16dimensional eigenspace of H T i H i .The null space method [24] is used to reduce the number of solutions and find W.In this method n measurements are taken with different polarisation elements as calibration samples.An equivalent equation to Eq. ( 9) is written for each measurement and then linearly combined to form where As also shown in the Appendix, if the n calibration samples are chosen appropriately vec(W) is the eigenvector of L corresponding to the null eigenvalue of L. The matrix L however, must have only one null eigenvalue to ensure that the solution is unique.If for example two eigenvalues of L are zero, vec(W) lies in a 16−dimensional hyperplane defined by the two eigenvectors corresponding to the null eigenvalues.If three eigenvalues are zero, vec(W) thus lies in a 16−dimensional hypervolume and so forth.
To ensure that W is the unique solution of Eq. ( 10), the n calibration samples are choosen such that the dimension of the null space of L is reduced to one 1 .In theory a perfect polariser and an ellipsometric surface can be used if the orientation of the polariser is different from zero or π/2, and the ellipsometric parameters of the surface are such that Φ = 0 or π/4, and ∆ = 0 or π [1].Once it is certain that only one eigenvalue of L is zero, the angular position of the calibration samples θ i can be found by minimising which is the ratio of the smallest and second smallest eigenvalue.Note that the minimisation yields only the relative angular orientation of the polarisation elements.Since the eigenvalues of a matrix are invariant to rotation, it is possible to write M i = R(θ )M i R(−θ ) for the n calibration elements and still obtain the same eigenvalues with a different W.An absolute angular orientation can be obtained if the frame of reference is set such that it coincides to the reference axis of the first polariser placed after the source.Once W is found, T can finally be obtained from Note that Eq. ( 12) is a metric that quantifies the errors in the calibration originated from imperfections in the polarisation elements, laser noise, and electronic noise.A statistical comparison of different measurements of D air can quantify the contribution of the laser and electronic noise to .
The ECM described so far can be used to calibrate a system set up in transmission.As explained in the Following section, calibration of a system in reflection or in double-pass requires the calibration to be complemented with more measurements.

VARIANTS OF THE ECM
To introduce the variants of the ECM consider the experimental setup depicted by the block diagram shown in Figure 2. In this figure the bold letters represent the Mueller matrix of each unit: N represents a confocal microscope formed by imaging optics F (including the pinhole) and a flat mirror E as sample; B is the Mueller matrix of the beamsplitter, and Q represents the setup excluding the PSA and the PSG.Three calibration procedures are considered to calibrate such a system.These are a two-step, a three-step, and the double pass calibration [9].The two-step calibration process yields W, T, and Q.The double pass procedure yields W, T, B, and F if E is assumed to be a perfect mirror.The three-step procedure is a combination of the two-step, the double pass calibration, and a fitting algorithm.It yields W, T, B, and N, and also F and E if the latter is assumed to be a perfect mirror.Regardless of the method used, the position in the setup where the polarisation is effectively being measured is given by the location of the calibration samples.The coloured arrows in Figure 2 denote such positions for the different methods.
Note that for high NA imaging systems, the system must be calibrated at the image plane.This is because the distribution of the field at the pinhole depends on the spatial variations of phase, amplitude, and polarisation occurring at the back focal plane of the high NA lens [28].At this plane the longitudinal component of the scattered field is intermixed with the transverse components.Calibrating at this plane thus disregards the contribution of the longitudinal component to the formation of the image because the calibration samples can not discriminate the transverse from longitudinal contributions.The image plane is located at the pinhole and hence inaccessible in practice, but if the pinhole is sufficiently small the calibration samples can be placed between the imaging optics and the beamsplitter (blue and green dot); this is because there is no longer contribution from the longitudinal component at this place [28].

Two-step calibration
The two-step calibration yields W, T, and the products TQ and QW from which Q can be easily found.This calibration is suitable for setups in which Q shows only retardance, and when one is interested in the relative polarisation properties of the sample.As long as Q does not exhibit diattenuation, polarisance or depolarisation, the information present in the polarisation states measured is mapped, i.e. rotated around the Poincare sphere, by the retardance in Q and thus conserved.
The appropriateness of this method should thus be assessed by performing a Lu-Chipman decomposition [11] on Q after calibration.
As its name suggests, this calibration is a two step process.
First W and the product TQ are obtained by calibrating with the polarisation elements M i placed between the PSG and the beamsplitter.Note that at this stage it is not possible to obtain T alone since D air = TQW.In the second step the polarisation elements are placed between the PSA and the beamsplitter.This step yields T and the product QW.The Mueller matrix of the system Q can then be calculated from In practice however, the equality does not hold due to electronic noise.The norm of the difference between Q 1 and Q 2 can be used as a relative error estimator as seen next.It is important to mention that the calibration sample that defines the reference axis of the polarisation must be the same in both steps.
The efficiency of the two-step calibration procedure was assessed using a flat mirror as object to be imaged.A high NA objective lens (NA = 0.95) was used, and the chosen polarisation elements for the calibration were [27] a polariser at = 0 • , one at ≈ 90 • , and a quarter waveplate at ≈ 30 • .The pinhole used was filtering 2/3 of the incident Airy disc.The measured Q showed a retardance of λ/24 (possibly introduced by the beamsplitter alone), polarisance and diattenuation of 0.015, and 0.008 of depolarisation.The relative error in Q was 8.5×10 −4 , which corresponds to the order of magnitude of the electronic noise.

Double pass calibration
The double pass calibration is a variant of the ECM introduced by Lara and Dainty [9] for systems where the light passes twice through a portion of the setup.To use this variant of the ECM it is necessary to assume that the light passing back through the polarisation elements is reflected by a perfect mirror at normal incidence, and that the behaviour of the polarisation elements is independent of the direction of propagation of the light.The practical validity of these assumptions can be assessed using an algebraic procedure explained at the end of this section.
In systems such as the one depicted in Figure 2 light passes twice through the same polarisation element M i , and in con-sequence Eq. ( 1) has no longer same form.Adapting Figure 2 to Figure 3 following the direction of propagation, it can be seen that Eq. ( 1) is replaced by where B f and B b denote the Mueller matrices of the forward and backward propagating path of the beamsplitter, respectively.The equation equivalent to Eq. ( 3) is where Note from Eq. ( 15) that C and K are similar matrices, however, K and M b i are not.The latter means that the eigenvalues of K are no longer independent of the angular orientation of M i .Unless M b i and N commute, it is not possible to find M i and hence H can not be constructed.To find W using Eq. ( 15) one needs to assume that N = E, where E is the Mueller matrix of a perfect mirror at normal incidence.One also needs to assume that the polarisation elements used for calibration are independent of the propagation direction.The latter means that the Mueller matrices of each polarisation element can be written as and Substituting Eqs. ( 16) and ( 17) into Eq.( 15) yields if E is chosen to commute to the left.If E is chosen to commute to the right the substitution yields In both cases it is possible to reconstruct H and solve the problem using the null space method described in section (1.1).The choice of the direction of commutation depends on the angular orientation of the calibration samples input to the calibration algorithm.If the samples are placed at θ and the data is processed with the angular orientation of the samples at θ, E commutes to the left naturally.For E to commute to the right the calibration samples are still placed at θ, but the data is processed as if the element was placed at −θ.In the first case the calibration yields B f W and TB b i .In the second case the calibration yields TB b i E and B f W. If the assumptions were valid, it follows that Two scenarios were used to test the double pass calibration method.One scenario consisted in the calibration of the system using a flat mirror and removing both the pinhole and the objective lens.In this case the assumptions mentioned above were valid.The relative error between the matrices obtained from the left and right commutation was 1.2 × 10 −4 .The second scenario consisted in calibrating the system with a flat mirror and with the objective lens in place and using a pinhole of 2/3 of the Airy disc.In this case the relative error was 0.03.Even though the first scenario yielded satisfactory results, one needs to guarantee that the path of the light is the same with and without the pinhole and the objective lens.The results from the second scenario were unsatisfactory, however, these results can be improved with the three-step calibration process as discussed next.

Three-step calibration
The three-step calibration method can be used if the polarisation properties of the polarisation elements M i depend on the direction of propagation of the light or if the mirror is not perfect.This calibration method is a combination of the double pass calibration, the two-step method and a fitting algorithm.The aim is to find T, W, the polarisation elements M f i and M b i used as calibration samples, B f , B b , and N.
The calibration procedure is as follows.First a set of measurements are taken placing the n calibration samples in three different places of the setup: between the PSG and the beamsplitter, between the beamsplitter and the PSA, and between the beamsplitter and the pinhole.Since the first two places correspond to the positions where the measurements of the twostep calibration are taken, it is possible to obtain the Mueller matrices of the calibration samples in both directions M f i and M b i , T, W, and Q.Note that Q = B b NB f and can be written using Eq. ( 6) as where B = (B f ) T ⊗ B b .The measurements, taken between the beamsplitter and the pinhole correspond to the the double pass method, yield After multiplying by the inverse of T and W Eq. ( 21) can be simplified to Equation ( 6) is again used together with a property of the Kronecker product2 to write Eq. (22) as The linear combination of the n measurements can thus be given by where , and the 16 × 16 identity matrix from the measurement of air is represented by I.Note the similarity between Eqs. (20) and (24).The remaining matrices B f , B b , and N can be found from these two equations using a Levenberg-Marquardt algorithm [26].With this algorithm it was possible to solve Eqs. ( 20) and ( 24) using 4 × 4 random matrices as entries with an accuracy comparable with the electronic noise of the experiment.
The three-step calibration method was tested using four polarisation elements as calibration samples.These where a polariser at = 0 • , one at ≈ 90 • , a quarter waveplate at ≈ 30 • , and a λ/8 waveplate at ≈ 30 • .The relative error between the measured Q and the retrieved product B b NB f was 6.5 × 10 −6 .It was found however that the result depends on the starting points of the algorithm.Rough a priori information of B f and B b is thus essential.This information can be obtained by performing Mueller matrix polarimetry on the beamsplitter in reflection and transmission.Also note that the Mueller matrix measured (N) comprises both the sample and the imaging system.Disentanglement between these two should not be sought unless a PSA capable of measuring the x, y, and z components of the field and a 3D notation is used instead.This is because the scanning confocal imaging system projects the 3D interaction of the polarisation into 2D, hence validating the use of the 2D notation.
As an empirical validation of the three-step calibration, the next section shows the dependance of the polarisation in the system due to the size of the pinhole.

EFFECT OF THE PINHOLE SIZE ON THE POLARISATION
The effect of the pinhole size on the polarisation was measured to corroborate the calibration of the system when using the three-step method.With a flat mirror in the object plane of the objective lens, the Mueller matrix N was measured for different sizes of the pinhole.Figure 4 shows the results from the polar decomposition [11] of the different N measured.The top left and right plots show the magnitude of the depolarisation and polarisance, respectively.The bottom left and right plots show the magnitude of the retardance and diattenuation, respectively.
Apart from the depolarisation plot, the rest of the plots in Figure 4 show that the measurement tends towards a perfect mirror as the pinhole is reduced.One would expect that a perfect mirrors yields a diattenuation, depolarisation and polarisance equal to zero, and a retardance equal to 180 degrees; however, the measured polarisation properties are different from the expected due to the spatial variations of the polarisation introduced by the high NA objective lens [3], which are minimised by the pinhole [28,2], and due to electronic noise.The trend in the depolarisation plot is different from the rest due to a decrease in the signal-to-noise ratio o

CONCLUSION
It was seen that the choice of the method depends on the polarisation properties of the system.The two-step method is ideal for systems that only show retardance and if the user is only interested on the relative polarisation properties of the sample.This method can be used for systems configured in reflection or in transmission.The double pass variant can be used for systems where the light passes twice through the calibration samples M i , i.e. systems in reflection configuration at normal incidence.To use this variant, it is necessary to assume that the light passing back through M i was reflected by a perfect mirror at normal incidence and that the beam path remains unchanged with or without the objective lens.The three-step calibration method can be used if the behaviour of the polarisation elements M i depends on the direction of propagation of the light or if the mirror is not perfect.This calibration method is a combination of the double pass calibration, the two-step method and a Levenberg-Marquardt fitting algorithm.It was also mentioned that, despite of the chosen calibration method, it is imperative to calibrate the system after the pinhole when using high NA optics.
Lemma The matrix W written as a 1 × 16 vector is the eigenvector that corresponds to the null eigenvalue of H [30].
Proof Let X = {x T 1 , ..., x T m } be the eigenvectors of M with corresponding eigenvalues µ 1 , ..., µ m , and Y = {y It then follows from equation ( 25) that W is now in the eigensystem of L, and that the dimension of the null space of L can be reduced to one as long as W and M i does not commute.

FIG. 1
FIG. 1 Block diagram of a Mueller matrix polarimeter.

FIG. 4
FIG. 4 Parameters from the polar decomposition of the confocal signal produced by a mirror for different pinhole diameters.

T 1 ,YΛY − 1 = W − 1 XΛX − 1 W
..., y T m } be the eigenvectors of C with corresponding eigenvalues κ 1 , ..., κ m .Since H = M ⊗ I − I ⊗ C T , it then follows that[22] (M ⊗ I − I ⊗ C T )(x ⊗ y) = (Mx ⊗ y) − (x ⊗ C T y) = (µx ⊗ y) − (x ⊗ κy) = (µ − κ)(x ⊗ y).Therefore the non trivial solution is given when (µ − κ) is zero.LemmaThe matrix W can be written as(x ⊗ y) if C is Hermitian.Proof Note that M and C are similar matrices and assume that at least C is Hermitian.Let M = XΛX −1 and C = YΛY −1 be the spectral decomposition of M and C, respectively.The matrix Λ is diagonal and contains the eigenvalues of M and C. Rewriting C = W −1 MW as equation (26) in terms of Kronecker products as vec(W) = (X ⊗ Y −1T )vec(I), and since the eigenvectors of a Hermitian matrix are orthogonal, i.e.Y −1 = Y T , then vec(W) = (X ⊗ Y)vec(I), = (x ⊗ y).Lemma A linear system of equations denoted by L can be constructed to contain W as the unique solution.ProofNote that the 16 eigenvalues of H are given by the possible combinations of µ 1..4 − κ 1..4 .The number of times that µ − κ = 0 appears can thus be given by[1] α is the multiplicity of each eigenvalue, and ℵ represents the dimension of the null space of H. Consider now n different measurements denoted byC M i = W −1 M i W,for i = 1, ..., n.It is then possible to produce n different H M i and write is the 4 × 4 identity matrix and H i is a 16 × 16 matrix.