## 2 Principle
The polarization transfer characteristics of the waveplate can be described by the Muller matrix (\(M_{\rm S}\)) as follows [1,12]. \[M_{\rm S}= \begin{bmatrix}1&0&0&0&0\\ 0&\cos^{2}(2\theta)+\sin^{2}(2\theta)\cdot\cos(\phi)&\sin(4\theta)\cdot\sin^{2 }(\phi/2)&\sin(2\theta)\sin(\phi)\\ 0&\sin(4\theta)\cdot\sin^{2}(\phi/2)&\sin^{2}(2\theta)+\cos^{2}(2\theta)\cdot \cos(\phi)&-\cos(2\theta)\sin(\phi)\\ 0&-\sin(2\theta)\sin(\phi)&\cos(2\theta)\sin(\phi)&\cos(\phi)\end{bmatrix}  \tag{3}\]  
where \(\phi\) is the phase retardance. It can be further expressed as \(\phi=\frac{2\pi}{\lambda}(n_{e}-n_{o})d\), where \(\lambda\) is the wavelength of incident light, \(d\) is the thickness of the waveplate, and \(n_{o}\) and \(n_{e}\) are the refractive indices of ordinary and extraordinary light, respectively. \(\theta\) is the fast-axis azimuth. It can be observed from the above formula that the polarization characteristics of the waveplate are fully determined by the two parameters of phase retardance \(\phi\) and fast-axis azimuth \(\theta\). They must be strictly measured and calibrated, especially the parameters over the entire clear aperture. Furthermore, the transmission axis of the analyzer is set at 45\({}^{\circ}\), and its Muller matrix (\(M_{\rm A}\)) is expressed as follows. \[M_{\rm A}=\frac{1}{2} \begin{bmatrix}1&0&1&0\\ 0&0&0&0\\ 1&0&1&0\\ 0&0&0&0\end{bmatrix}  \tag{4}\]  
Here, we insert into Equation (1) the Muller matrices of the PEMs, waveplate, and analyzer expressed in Equations (2)-(4), together with the Stokes vector of incident light. We consider that the light intensity perceived by the detector is the first component of the Stokes vector [13,14]. It can, thus, be solved to obtain the output light intensity detected by the detector to the measurement system as follows. \[ \begin{array}{l}I=\frac{I_{0}}{2}(1+\cos\delta_{1}\cos\delta_{2}\sin(4\theta )\sin^{2}(\frac{\phi}{2})+\sin\delta_{1}\sin\delta_{2}\cos\phi+\\ +\cos\delta_{1}\sin\delta_{2}\cos(2\theta)\sin\phi+\sin\delta_{1}\cos\delta_{2 }\sin(2\theta)\sin\phi)\end{array}  \tag{5}\]  
In the formula above, the two PEM modulation items, \(\sin\delta_{i}=\sin(\delta_{i0}\sin\omega_{i}t)\) and \(\cos\delta_{i}=\cos(\delta_{i0}\sin\omega_{i}t)\), are expanded as \(\sin\delta_{i}=2\sum\limits_{2k-1}J_{2k-1}(\delta_{i0})\sin((2k-1)\omega_{i}t)\) and \(\cos\delta_{i}=J_{0}(\delta_{i0})+\sum\limits_{2k}2J_{2k}(\delta_{i0})\cos(2 k\omega_{i}t)\) using the first kind of Bessel functions; \(k=1,2,\cdots\) are positive integers; \(J_{0}\), \(J_{2k-1}\), and \(J_{2k}\) represent the Bessel series of the 0 th order, the \((2k-1)\)th order, and the \(2k\)th order, respectively [11,15,16]. Additionally, PEM1 and PEM2 correspond to \(i=1\) and \(2\), respectively. Taking \(k=1\) as the low-order Bessel series, Equation (5) can be further rewritten as follows. \[ \begin{array}{l}I=\frac{\omega_{1}}{2}(1+(J_{0}(\delta_{i0})+2I_{2}(\delta_{ i0})\cos(2\omega_{1}t))\cdot(J_{0}(\delta_{20})+2J_{2}(\delta_{20})\cdot \cos(2\omega_{2}t))\sin(4\theta)\sin^{2}(\frac{\phi}{2})\\ +4J_{1}(\delta_{i0})I_{1}(\delta_{20})\sin(\omega_{1}t)\sin(\omega_{2}t)\cos \phi+(J_{0}(\delta_{i0})+2I_{2}(\delta_{i0})\cos(2\omega_{1}t))\cdot(2J_{1}( \delta_{20})\sin(\omega_{2}t))\cdot(\cos(2\theta)\sin(\omega_{2}t))\cdot\\ \cos(2\theta)\sin\phi+(J_{0}(\delta_{20})+2J_{2}(\delta_{20})\cos(2\omega_{2}t ))\cdot(2J_{1}(\delta_{10})\sin(\omega_{1}t))\sin(2\theta)\sin\phi)\end{array}  \tag{6}\]  
From the analysis of Equation (6), it can be observed that the relevant terms for the waveplate phase retardance and fast-axis azimuth are loaded in different frequency signals. The relevant term \(\sin(4\theta)\sin^{2}(\phi/2)\) is loaded into \(2\omega_{1}\), \(2\omega_{2}\), \(2\omega_{1}+2\omega_{2}\), and \(2\omega_{2}-2\omega_{1}\)signals, and the relevant term \(\cos\phi\) is loaded into \(\omega_{2}-\omega_{1}\) and \(\omega_{1}+\omega_{2}\) signals.