9424

Research Article

Vol. 59, No. 30 / 20 October 2020 / Applied Optics

Dynamic calibration for enhancing the stability ofa channeled spectropolarimeter
Benjamin D. Chrysler,1, *Kazuhiko Oka,2 Yukitoshi Otani,33

AND Nathan Hagen

1

University of Arizona, James C. Wyant College of Optical Sciences, 1630 E. University Blvd., Tucson, Arizona 85721, USA
Faculty of Science and Technology, Hirosaki University, 3 Bunkyo-cho, Hirosaki 036-8561, Japan3
Utsunomiya University, Department of Optical Engineering, Center for Optical Research and Engineering (CORE), 7-1-2 Yoto,
Utsunomiya, Tochigi 321-0954, Japan*Corresponding author: bchrysler@email.arizona.edu
2Received 2 September 2020; revised 20 September 2020; accepted 22 September 2020; posted 23 September 2020 (Doc. ID 409054);
published 15 October 2020

Channeled spectropolarimeters are optical instruments that measure the spectral dependence of the polarization oflight without any mechanically moving parts. An important factor in achieving stable and accurate measurements
is the calibration process, especially in dynamic environments where temperature fluctuations or other factors
affect the retardance of the components in the polarimeter. In previous research, a self-calibration algorithm that
accounts for these variations was developed, without any additional reference measurements. In this paper, we identify an ambiguity in the self-calibration phase, which limits the allowed temperature changes to surprisingly small
ranges. We show how to adaptively estimate and correct for the phase ambiguity using a polynomial curve-fitting
algorithm, extending the temperature range to virtually all practical scenarios. Lastly, we demonstrate the abilityof the modified self-calibration algorithm to provide stable reconstruction of the Stokes vector for a temperature
range >40◦ C, using an experimental channeled spectropolarimeter. © 2020 Optical Society of America
https://doi.org/10.1364/AO.409054

1. INTRODUCTION
Polarization is an important property of light, and it has a variety
of applications [1]. The polarization state of incoherent light*is described by the Stokes vector S (λ) as a function of wavelength λ, and it provides important information that cannot
be obtained using conventional intensity measurements orimages. Stokes polarimeters are instruments that measure thepolarization state of incoherent light, and they are useful inastronomy [2], remote sensing [3], biology [4,5], and industrial
metrology [6].
Several types of polarimeters have been developed, withpolarizing filter array cameras [7] and rotating element polarizers [8] being the most common. Polarizing filter array cameras
are advantageous for their low cost and the ability to capture thepolarization information of a scene. Unfortunately, the cameras
do not typically capture information about circular polarizationor spectral variation and thus do not measure potentially important information. In contrast, rotating element polarimeters
often measure spectral variation and capture the full Stokes
vector. While they have found application in industry, they are
relatively slow and bulky, and the use of mechanically moving
parts increases the chance of breaking down over time.
Channeled spectropolarimeters (CHSP) are another class ofpolarimeters using mechanically fixed thick retarders and linear
1559-128X/20/309424-10 Journal © 2020 Optical Society of America

polarizers to modulate the Stokes vector in the wavenumber
domain [9]. The Stokes vector components are reconstructedby taking the Fourier transform of the wavenumber spectrum,
applying a set of notch filters in the optical path difference(OPD) domain, and taking the inverse Fourier transform.
While currently less developed than the rotating-element polarizer, the CHSP has the potential advantages of being faster, lessexpensive, more compact, and more robust to vibration.
One of the challenges in obtaining accurate measurementswith CHSPs is the sensitivity to fluctuations in environmental parameters such as temperature, which can vary between−40◦ C and +100◦ C for some out-of-laboratory applications.
One approach to solving this problem is to reduce the temperature sensitivity of the thick retarders by designing athermal
composite retarders [10]. However, these specialized retarders
are more expensive and cannot completely eliminate the problem of thermal phase fluctuations. Other approaches have been
to dynamically compensate the changing phase delay of theretarder using a “self-calibration” approach, which compensatesfluctuations without any additional calibration measurements[11–13]. A more recent approach uses phase-shifting interferometry to correct for phase fluctuations, but this method
uses mechanically moving parts and multiple measurements
per reconstruction cycle [14], which detracts from the CHSP
advantages of being compact, durable, and fast.

Research Article

Vol. 59, No. 30 / 20 October 2020 / Applied Optics

While the self-calibration algorithm is an effective calibration
technique for CHSPs in dynamic environments, improvements
still need to be made to enhance the accuracy and stability of thealgorithm. In this paper, we first review the process for implementing self-calibration in a CHSP system, including somesteps taken in practice to provide accurate Stokes vector reconstructions [12,13] not explicitly stated when the algorithm wasinitially developed [11], such as the incorporation of occasionalreference calibration measurements and phase-unwrapping
algorithms. Next, we identify a source of instability resulting
from a mathematical ambiguity in the bias value of the reconstructed phase. Large errors in the reconstructed Stokes vector
are shown to occur when the temperature exceeds a thresholddetermined by the temperature sensitivity of the birefringent
material and the order of the retarder. We show the importanceof correcting the phase ambiguity since many advantageous
CHSP configurations are especially sensitive to this problem.
For example, quartz retarders have high temperature sensitivity[15] but are otherwise a practical material for use in CHSPs
since many quarter- and half-wave plates sold on the market are
sufficiently high order for use in CHSPs [16] and are relatively
inexpensive. Another important example is the Mueller-Matrix
CHSP, which uses very thick retarders to attain sufficiently highspectral resolution [17], but accordingly suffers from greater
instability. We propose a modified version of the self-calibration
algorithm that adaptively corrects for the phase ambiguity, providing stable reconstruction for most CHSP configurations inpractical scenarios. This approach works by comparing the slopeof the self-calibration phase during the reference calibrationand in real-time, by using curve-fitting polynomial approximations and determining deviations greater than those that canbe detected using self-calibration directly. Lastly, we show howtemperature differences between the CHSP retarders affect thereconstructed phase and how to correct for these errors if thetemperature of both retarders is known.
2. CHANNELED SPECTROPOLARIMETRYA CHSP consists of two thick retarders, a linear polarizer, anda spectrometer, as illustrated in Fig. 1. The first thick retarder
has a thickness d1 , and the orientation of the fast axis of thecrystal determines the horizontal axis of the coordinate system.
The second thick retarder has a thickness d2 , and the fast axis isoriented at 45◦ . The transmission axis of the linear polarizer isaligned with the horizontal axis. The polarization of the input
beam is described by the Stokes vector:*TS (σ ) = S0 (σ ) S1 (σ ) S2 (σ ) S3 (σ ) ,(1)

where σ is the wavenumber of light defined as the inverse of
the wavelength λ, σ = 1/λ. The first component of the Stokes
vector S0 (σ ) is the intensity of the light, while S1 (σ ), S2 (σ ),and S3 (σ ) describe the statistical tendency of the beam tohave a horizontal/vertical linear, +45◦/− 45◦ linear, or righthanded/left-handed circular polarization. The spectrum P (σ )measured by the CHSP spectrometer [an example spectrum isshown in Fig. 2(a)] can be determined analytically using Mueller
calculus, in which Mueller matrices describe the transformationof the Stokes vector of light propagating through polarization

9425

Fig. 1. Configuration of the optical components in a channeledspectropolarimeter. The fast axis of the second retarder is oriented at
45◦ with respect to the first retarder, and the transmission axis of thelinear polarizer is horizontal. The components modulate the Stokes
vector components of the input beam onto the spectrum measured by
the spectrometer.

components. The results of applying these calculations to the
CHSP components give [9]11
P (σ ) = S0 (σ ) + S1 (σ ) · cos(φ2 (σ ))221+ S2 (σ ) · sin(φ2 (σ )) · sin(φ1 (σ ))2
1− S3 (σ ) · sin(φ2 (σ )) · cos(φ1 (σ )),2(2)where φ1 (σ ) and φ2 (σ ) are the phase delays between the fast andslow axes of the first and second retarders caused by the birefringence B:
φ(σ ) = 2πB(σ )σ d .

(3)

By properly selecting the thicknesses of the two thick retarders,the Stokes vector components are separated into frequency
channels C i (h) in the Fourier domain:C i (h) = F{P (σ )} · Wi (h),(4)

where i is the channel number, h is the OPD, and Wi (h) is arectangular windowing function as depicted in Fig. 2(b). The
channels are represented in wavenumber space by taking theinverse Fourier transform of the channels in the OPD domain:
Fi (σ ) = F −1 {C i (h)}.

(5)

When the thicknesses of the first and second retarder are configured in a 1:2 ratio, the channel content relates to the Stokes vector components:F0 (σ ) = 21 S0 (σ ) F (σ ) = 1 S ∗ (σ )exp( j [φ (σ ) − φ (σ )])
1
2
1
8 23
,
(6)
1

F
(σ
)
=
S
(σ
)exp(
j
[φ
(σ
)])
2
1
2

4


F3 (σ ) = − 18 S23 (σ )exp( j [φ2 (σ ) + φ1 (σ )])where S23 (σ ) = S2 (σ ) − j S3 (σ ). The amplitude and phaseterms in Eq. (6) can be combined and written as the calibration
coefficients K i (σ ) for each channel [12]:Research Article

Vol. 59, No. 30 / 20 October 2020 / Applied Optics

9426

Fig. 2. (a) Modulated spectrum P (λ) measured by the spectrometer and (b) channel structure obtained by taking the Fourier transform of P (σ ).
A total of seven channels are obtained using a retarder thickness ratio of 1:2. Individual channels are filtered using a rectangular notch filter Wi (h).

K (σ ) = 2 0
K 1 (σ ) = 8 · exp(− j [φ2 (σ ) − φ1 (σ )])K i (σ ) =. (7) K 2 (σ ) = 4 · exp(− j [φ2 (σ )])K 3 (σ ) = −8 · exp(− j [φ2 (σ ) + φ1 (σ )])Multiplying the calibration coefficients with the channel
contents then yields the Stokes vector components:
F0 (σ ) · K 0 (σ )*F2 (σ ) · K 2 (σ ) .
(8)S (σ ) = 
Re{F3 (σ ) · K 3 (σ )} 
−Im{F3 (σ ) · K 3 (σ )}In this paper, the channel i = 3 is used for reconstructing bothS2 (σ ) and S3 (σ ) to reduce channel crosstalk from the central
channel and the effects of dichroism [18], although it is possible
to reconstruct them using channel i = 1 as well:
F0 (σ ) · K 0 (σ )* F2 (σ ) · K 2 (σ ) .
(9)S (σ ) = 
Re{F1 (σ ) · K 1 (σ )} 
Im{F1 (σ ) · K 1 (σ )}The compensation factors K i (σ ) must be accurately calibrated for a precise reconstruction of the Stokes vector. The
compensation factor can be viewed as the product of an initialcompensation factor K r ,i (σ ), determined during the initialreference calibration of the system. and a dynamic compensation factor K d ,i (σ ), determined in real-time during the Stokes
vector measurement. It accounts for changes in the phase of theretarders due to fluctuations in the temperature, beam angle, orother factors:
K i (σ ) = K r ,i (σ ) · K d ,i (σ ).

(10)


T
(22.5◦ )*1 1
S (σ )= S0 (σ ) · 1 √ √ 0 .
2 2

The initial compensation factors are calibrated based on thismeasurement:K r ,0 (σ ) = 2
(22.5◦ )
1 K r ,0 (σ )·{F0 (σ )}◦) K r ,1 (σ ) = √2
(22.5
{F1 (σ )}◦(12)
K r ,0 (σ )·{F0 (σ )}(22.5 ) .
1

√
K
(σ
)
=
◦)

r
,2
(22.5

2
{F
(σ
)}
2


(22.5◦ ) K r ,3 (σ ) = √1 K r ,0 (σ )·{F0 (σ )}◦(22.5 )
2
{F3 (σ )}While this method is simple since it uses a single measurement,
higher accuracy can be obtained using two reference measurements taken with the polarizer oriented at 0◦ and at 45◦ ,respectively:*{ S (σ )}(22.5 ) :
◦

*
(0◦ )
S (σ )= S0 (σ ) · [1 1 0 0]T ,

(13)

*
(45◦ )S (σ )= S0 (σ ) · [1 0 1 0]T .

(14)

The 0◦ measurement is used to calibrate the i = 2 channel, andthe 45◦ measurement is used to calibrate the i = 1, 3 channels:K r ,0 (σ ) = 2



◦


K (σ )·{F0 (σ )}(45 ) K r ,1 (σ ) = r ,0◦{F1 (σ )}(45 )◦(15)
K r ,0 (σ )·{F0 (σ )}(0 ) .
K r ,2 (σ ) =
◦)

(0

{F2 (σ )}(45◦ )

 K r ,3 (σ ) = K r ,0 (σ )·{F0 (σ )}◦(45 ){F3 (σ )}3. REFERENCE CALIBRATION
While the self-calibration algorithm is sometimes understoodas eliminating the need for reference measurements, in practice, it is usually implemented in combination with occasionalreference measurements [12,13]. The reference measurements
may be taken when initializing the system before use, or when
vibrations or other disturbances misalign system components.
The initial compensation factors K r ,i (σ ) are obtained usingreference measurements of a linear polarizer with a known stateof polarization. In previous literature, the calibration polarizer
was oriented at 22.5◦ [18] with the resulting polarization state(11)

The two-point reference calibration is more tolerant of calibration polarizer misalignment errors when compared to the
one-point reference calibration. This can be understood byconsidering the sensitivity of the Stokes vector components tothe angle of the calibration polarizer. A polarizer oriented at anangle θ with error  generates the following Stokes vector:
*

(θ+)


T
= S0 (σ ) · 1 cos(2 · (θ + ε)) sin(2 · (θ + ε)) 0 .
(16)
The Taylor series approximation of the alignment error centeredat θ is{ S (σ )}Research Article

Vol. 59, No. 30 / 20 October 2020 / Applied Optics

9427

the channels F1 (σ ), F2 (σ ), and F3 (σ ), as shown in Eq. (6),such that the argument of the complex function is a functionof the phase φ2 (σ ), and the magnitude is non-zero for any
polarized beam [11]:1  2
S + S22 + S32 e j 2φ2 (σ ) .
16 1
(20)
Taking the complex argument of Eq. (20) yields a function
φ0 (σ ) = 2φ2 (σ ) for ideal systems with band-limited signals andinfinite spectral resolution. In practice, only approximate estimates ψ0 (σ ) of the actual phase functions are obtained, whereψ0 (σ ) ≈ φ0 (σ ). The change in phase 1ψ 0 (σ ) between thecurrent and calibration measurement is determined by dividing
Eq. (20) by a similarly defined function using the reference
coefficients [12]:F22 (σ ) − 4 · F1 (σ ) · F3 (σ ) =Fig. 3. Error in the reconstructed Stokes vector components as afunction of calibration polarizer misalignment error using a “one-pointreference” method, consisting of a single calibration measurement of apolarizer oriented at approximately 22.5◦ , and a “two-point reference”
method, consisting of a set of two calibration measurements at 0◦ andat 45◦ , respectively.

1ψ0 (σ ) = UNWRAP
( (
× arg

1
(θ+)
*
cos(2θ ) − 2 sin(2θ ) − 4 2 cos(2θ )
.
S (σ )≈ S0 (σ ) sin(2θ ) − 2 cos(2θ ) − 4 2 cos(2θ )
0(17)
The relative error between the measured value of the Stokes vec**
tor { S (σ )}(m) and the actual value of the Stokes vector { S (σ )}(t)depends on the misalignment error of the calibration polarizer(1 pr )and can be calculated for the one-point method *andS(2 pr )
by combining Eqs. (8), (12),the two-point method *S(15), and (17):




(m)
*

S (σ )(1 pr )
*
=1−
(t)
S
*
S (σ )

22= 0 2 + 4 2 + 4 2 + 4


2 T

,

(m)
*

T

S (σ )(2 pr )
= 0 4 2 4 2 4 2 ,
*
=1−
(t)
S
*
S (σ )

F2 2 (σ ) − 4 · F1 (σ ) · F3 (σ )∗ ,4 · K r ,2 (σ )2 − K r ,1 (σ ) · K r ,3 (σ )(21)

where the arg function returns the angle of the complex function
between the bounds [x − π, x + π ), where x is the bias of
the function bounds and is commonly x = 0 or x = π . The
UNWRAP function adds or subtracts integer multiples of 2πfrom each data point to generate an unwrapped phase function
without phase differences of greater than π between data points.
While the phase-unwrapping step is not explicitly stated inliterature [11–13], it is essential for avoiding phase ambiguity
errors that cause discontinuities in the reconstructed phasesince the unwrapped phase function needs to be divided toobtain the retarder phase functions. The change in phase of thesecond retarder 1ψ 2 (σ ) is determined directly by dividing themeasured phase function by 2:
1ψ2 (σ ) =

(18)

(19)

where it is assumed that in the two-point method, both calibration measurements have the same misalignment error . The
error of the Stokes vector reconstructed using the two calibration methods is compared in Fig. 3, showing that the two-pointreference method has a lower sensitivity to misalignment errors
compared to the one-point method, resulting from the difference in the sensitivity of S1 (σ ) and S2 (σ ) to misalignment
errors for angles around 22.5◦ compared to 0◦ or 45◦ .

))1ψ 0 (σ )
.
2

(22)

The change in phase of the first retarder is indirectly determinedby dividing the change in phase of the second retarder by theratio of retarder thicknesses ρ = d2 /d1 :
1ψ1 (σ ) =1ψ 2 (σ )
.
ρ

(23)

The changes in the phase of the retarders are then incorporated
in the dynamic calibration coefficients, which are used in reconstructing the Stokes vector in Eqs. (8)–(10):



1

 K d ,1 (σ ) = exp − j1ψ 2 (σ ) 1 − ρ
K d ,2 (σ ) = exp − j 1ψ 2 (σ ).
(24)
 K (σ ) = exp − j 1ψ (σ ) 1 + 1 
d ,3
2
ρ

4. DYNAMIC CALIBRATION
A. Review of the Self-Calibration Algorithm

B. Phase Bias Ambiguity

The “self-calibration” technique dynamically determines
the compensation factors K d ,i (σ ) to account for changes
in the phase of the retarders. The method works by combining

An important consideration not previously discussed in
the literature is the existence of an ambiguity in the offset
value of the measured phase functions, existing even if the9428

Research Article

Vol. 59, No. 30 / 20 October 2020 / Applied Optics

phase-unwrapping algorithm has been implemented. Whencalculating the phase functions 1ψ 0 (σ ), the arg function
returns values within the range of x − π ≤ 1ψ 0 (σ ) < x + π .
However, the actual value of the phase functions do not necessarily lie within that range and can vary by integer multiplesof 2π . The phase-unwrapping algorithm eliminates therelative phase ambiguity between adjacent data points byadding integer multiples of 2π to ensure phase differences
between adjacent data points do not exceed π . However, theunwrapped function is relative to the phase value at the reference wavenumber σ0 , which is the wavenumber for which the
UNWRAP function directly returns the value of its argumentx − π ≤ 1ψ 0 (σ0 ) < x + π and for which the relative valuesof surrounding points are subsequently added or subtracted toreconstruct the unwrapped function. The actual value of thephase at the reference wavenumber itself may vary by integer
multiples of 2π from the value returned by the arg function.
Therefore, an “adjusted” phase function 1θ 0 (σ ) varies withrespect to the self-calibration phase function 1ψ 0 (σ ):1θ0 (σ ) = 1ψ0 (σ ) + m · 2π ,(25)

where m is the bias integer and determines the function bias
level. The adjusted phase difference of the two retarders is thencalculated by dividing by 2 and by ρ:
1θ2 (σ ) = 1ψ2 (σ ) + m · π ,

(26)

m·π
.
ρ

(27)

1θ1 (σ ) = 1ψ1 (σ ) +The adjusted set of dynamic calibration coefficients can then be
written in terms of the measured phase functions and the bias
integer:
 



1
1
exp
−
j
mπ
1
−
K
(σ
)
=
exp
−
j
1ψ
(σ
)
1
−
d
,1

2
ρ
ρ

K d ,2 (σ ) = exp[− j 1ψ 2 (σ )]exp[−jmπ]
.


 


1
1
K d ,3 (σ ) = exp − j 1ψ 2 (σ ) 1 + ρ exp − j mπ 1 + ρ(28)While the first term in each of the adjusted calibration coefficients in Eq. (28) is the same as the calibration coefficients in
Eq. (24), the additional exponential factor significantly changes
the overall value of the phase for some values of m. For example, if the value of m = 1 and the retarder ratio ρ = 2, then thechannel phase in Eq. (28) will differ by π2 , π , and 3π
compared2
to that in Eq. (24), resulting in errors in the reconstructed Stokes
vector. The erroneous values of the reconstructed Stokes vector*0

S (σ ) when using the conventional self-calibration coefficients[Eq. (24)] can be calculated in terms of the actual value of the*

Stokes vector S (σ ) and the value of m using Eqs. (8), (24), and
(28):


S0
S1 · cos( j mπ )


*0


S (σ ) = S2 · cos  j 1 + 1 mπ  + S3 · sin j 1 + 1 mπ .
ρ
ρ



 

 
1
1
S3 · cos j 1 + ρ mπ − S2 · sin j 1 + ρ mπ
(29)

C. Stability of the Self-Calibration Algorithm

The phase ambiguity in the self-calibration algorithm has
potential for causing significant errors in the reconstructed
Stokes vector when m 6 = 0, as indicated by Eq. (29). This condition occurs when the phase retardance deviates too muchfrom the initial reference conditions and exceeds a thresholdcondition for stability, potentially caused by changes in theenvironmental temperature, input angle, or other disturbance.
The range in which the phase can deviate without causing
error depends on the reference wavenumber σ0 . While mostphase unwrapping algorithms use the first data point in theseries as the reference, the choice is arbitrary, and in most cases,
it is best to use a reference wavenumber in a part of the spectrumwith the highest signal-to-noise ratio. If the reference wavenumber is selected in a noisy region of the spectrum, then any errors
occurring in the phase-unwrapping algorithm may affect thereconstructed phase, even in regions where there is a high signalto-noise ratio. For example, if the modulated spectrum obtained
from the spectrometer contains 1000 data points, and a noise
value greater than π exists between data points 5 and 6, thenthe phase-unwrapping algorithm will erroneously increase or
decrease by an integer value of 2π . If the first data point is used as
the reference, the remaining 994 data points after the noisy data
points will all be offset by integer multiples of 2π and exhibit theerrors of Eq. (29). However, if any data point after data point 6
is used as a reference, then a total of 995 data points will be error
free, with only the first several points having error. While moresophisticated phase-unwrapping algorithms that filter noise
may alleviate some of this risk, it is still a good practice.
To analyze the extent to which the retarders can deviate
and not cause errors, the bounds of the functions in Eqs. (25)and (26) are evaluated at the reference wavenumber. Since
the arg function returns an absolute value within the bounds[x − π, x + π ), and the UNWRAP function keeps the valueof its argument unchanged at the reference wavenumber, thebounds arex − π ≤ 1ψ 0 (σ0 ) < x + π .

(30)

The bounds of 1ψ 2 (σ0 ) are determined by dividing by 2:
x +π
x −π
≤ 1ψ 2 (σ0 ) <
.
2
2

(31)

Therefore, the range in which the functions 1ψ 0 (σ0 ) and1ψ 2 (σ0 ) can vary without error is 2π and π , respectively, withthe bounds being offset by the bias value x . The phase bounds in
Eq. (31) can be used to derive equations predicting the bounds
in which specific environmental parameters can vary withoutcausing error. The stable temperature bounds are determined byfinding the relationship between changes in environmental temperature and phase retardance. The material specific coefficientsγ given by Etzel et al. [15] define a proportionality between thechange in temperature 1T and normalized change in phase:1φ= γ 1T,φ(32)

where 1φ is the change in phase retardance and φ is the phaseretardance, which is often specified by the manufacturer inResearch Article

Vol. 59, No. 30 / 20 October 2020 / Applied Optics

Fig. 4. Best-fit approximation of the coefficient of normalizedchange in birefringence as a function of temperature γ . The redasterisks correspond to the data points measured by NREL [15].

terms of the number of waves of retardance w at a given wavelength λ. A quarter-wave plate designed for the He–Ne laser
wavelength, for example, has w = 1/4 and λ = 632.8. Since thereference wavelength is not necessarily the same as the reference
wavenumber σ0 , the phase retardance should first be determinedat the reference wavenumber:φ(λ) = 2π wλσ0 .

(33)

By combining Eqs. (30)–(32), the bounds over which thetemperature T can deviate from the reference calibration
temperature T0 is determined:
T0 +

x −π
x +π
≤ T < T0 +.
4π wλσ0 γ
4π wλσ0 γ(34)

The following example illustrates the effect the stability bounds
in Eqs. (31) and (34) have on the stability of the Stokes vector
reconstruction. We simulate a CHSP system using Mueller
calculus with ideal matrices for retarders and linear polarizers [1]. The birefringence of the quartz retarders is modeledusing Sellmeier coefficients for the ordinary and extraordinaryrefractive indices [19]. The effect of temperature changes on thebirefringence is modeled using the temperature coefficients γ .
Since the parameter γ is a function of wavelength λ, but it has
only been measured at four wavelengths [15], a linear best-fit
approximation in the wavenumber domain is used to derivea function for γ (λ) and is plotted in Fig. 4. The function isthen used to adjust the birefringence of the CHSP retarders,
depending on the difference between the simulated temperatureand reference calibration temperature.

9429

The simulated system consists of 60- and 120-wave retarders, specified at λ = 633 nm with a reference wavenumber
value of σ0 = 1.36 µm−1 . Simulated reference measurements
at room temperature are taken using the two-point reference measurement technique described in Section 3. The
numpy.angle function from the NumPy package in Python was
implemented, which has a bias value of x = 0 and returns thecomplex angle between [−π, +π ). Based on these parameters, the stable temperature range is calculated to be between−1.2◦ C ≤ T < +43.2◦ C.
Simulated measurements of a quartz half-wave retarderoriented at 22.5◦ with respect to a linear polarizer are taken at
42◦ C and 45◦ C, showing the accuracy of the self-calibration
method when the phase and temperature lie within the boundsof Eqs. (31) and (34) and the errors that occur outside thesebounds. The results for the simulated measurement at 42◦ C
are shown in Fig. 5. In Fig. 5(a) the self-calibration phase functions 1ψ 0 (σ ), 1ψ 1 (σ ), and 1ψ 2 (σ ) deviate by <0.0029 π ,0.0007 π , and 0.0015 π , respectively from the actual phasefunctions 1φ0 (σ ), 1φ1 (σ ), and 1φ2 (σ ). At this temperature,
the value of 1ψ 0 (σ0 ) lies just within the stable bounds marked
by the grey box. The reconstructed Stokes vector components
S1 , S2 , and S3 are plotted in Fig. 5(b) and deviate by <0.006,
0.011, and 0.007, respectively, from the actual values. After
increasing the temperature to 45◦ , the self-calibration phasefunctions 1ψ 0 (σ ), 1ψ 1 (σ ), and 1ψ 2 (σ ) deviate by 2π, π/2,
and π from the actual phase functions 1φ0 (σ ), 1φ1 (σ ), and
1φ2 (σ ), as shown in Fig. 6(a), while the adjusted phase functions (m = −1) 1θ 0 (σ ), 1θ 1 (σ ), and 1θ 2 (σ ) deviate by<0.0036 π, 0.0009 π, 0.0018 π , from the actual phase functions. At this temperature, 1ψ 0 (σ0 ) lies just outside the stablephase range marked by the grey box. The actual and reconstructed values of the Stokes vector components are plotted in
Figs. 6(b)–6(d). The components S1 , S2 , and S3 reconstructedusing self-calibration have errors as large as 0.44, 1.37, and1.37, respectively, when compared to the actual values, while
the adjusted (m = −1) reconstructed components deviate lessthan 0.0059, 0.0110, and 0.0084, respectively. This example
shows that the large errors described in Eq. (29) occur when thephase increases or decreases beyond the stability thresholds ofthe system given in Eqs. (31) and (34).

Fig. 5. Simulation results of a quartz half-wave plate oriented at 22.5◦ with respect to an ideal linear polarizer, for a CHSP configuration using
60- and 120-wave quartz retarders at a temperature of 42◦ C. (a) The actual phase functions 1φ0 (σ ), 1φ1 (σ ), and 1φ2 (σ ) (solid lines) deviate byless than 0.0029 π, 0.0015 π, and 0.0007 π , respectively, from the self-calibration phase 1ψ0 (σ ), 1ψ1 (σ ), and 1ψ2 (σ ) (dashed lines), making any
differences difficult to see at this scale. The gray area indicates the phase range in which 1ψ0 (σ0 ) must lie to provide stable results. In this example,
1ψ0 (σ0 ) lies just within the stable range. (b) The normalized Stokes vector components S1 , S2 , and S3 (solid lines) deviate less than 0.006, 0.011, and
0.007, respectively, from the actual values, and any differences are difficult to see at this scale.

9430

Research Article

Vol. 59, No. 30 / 20 October 2020 / Applied Optics

Fig. 6. CHSP simulation results of a quartz half-wave plate rotated at 22.5◦ with respect to an ideal linear polarizer, for a CHSP configuration
using 60- and 120-wave quartz retarders at a temperature of 45◦ C. (a) The actual phase functions 1φ0 (σ ), 1φ1 (σ ), and 1φ2 (σ ) (solid lines) deviate
by 2π, π , and π/2, respectively, from the self-calibration phase functions 1ψ0 (σ ), 1ψ1 (σ ), and 1ψ2 (σ ) (dashed lines). The self-calibration phasefunctions adjusted for m = −1 (dotted lines) vary by less than 0.0036 π, 0.0018 π, and 0.0009 π , respectively, from the actual phase functions, making any differences difficult to detect at this scale. In this example, 1ψ0 (σ0 ) lies just outside the gray area that indicates the stable phase range, resulting in erroneous reconstructed phase. (b)–(d) The normalized Stokes vector components (solid lines) deviate as much as 0.44, 1.37, and 1.37 withrespect to the self-calibration phase (dashed lines), while the adjusted self-calibration (dotted lines) deviate by less than 0.0059, 0.0110, and 0.0084,
respectively, from the actual phase functions.

5. ENHANCED-STABILITY DYNAMIC
CALIBRATIONThe Stokes vector reconstruction errors given in Eq. (29) can
be compensated if the value of m is known and used in the calibration coefficients in Eq. (28). In this section, we describe amethod for extending the stability of the CHSP by dynamicallyobtaining the value of m. The method utilizes a polynomial
curve-fitting model of 1ψ 0 (σ ) to evaluate the absolute deviation of the phase retardance from the initial reference condition.
The polynomial curve-fitting function A(σ ) is given by theformulaA(σ ) =KX

Ak · σ k ,

(35)

k=0

where K is the order of the polynomial fit, and A k are thecoefficients that give the best fit to the function 1ψ 0 (σ ). The
coefficients A k can be determined using functions in standard
libraries such as numpy.polyfit() in Python. Next, an assumption
is made that the phase functions should approximately be 0
for wavenumbers close to the origin, i.e., 1θ 0 (0) ≈ 0. This
assumption is based on physical models of phase retardancesuch as Eq. (3), which states the retardance is proportional to
the wavenumber and is therefore 0 at the origin. However, theself-calibration phase 1ψ 0 (σ ) that the curve-fitting model isbased on is offset from the adjusted phase function 1θ 0 (σ )by a constant, depending on the value of m. The polynomial
coefficients A k for the curve-fitting model A(σ ) are therefore
accurate for k > 0, but the constant term A 0 is incorrect in
the cases that m 6 = 0. Given that the A k coefficients are always

accurate for k > 0 and the physical constraint that 1θ 0 (0) ≈ 0,
an accurate model is obtained by removing the constant term
A 0 . The adjusted model for A(σ ) is compared to 1ψ 0 (σ )to find the value of m that best adjusts the difference between
the measured phase function 1ψ 0 (σ ) and the modeled phasefunction:
( K
)
X


k
A k · σ0 − 1ψ 0 (σ0 ) + n · 2π
m = argminn, (36)
k=1

where argminn is a function that returns the value of the integer n that minimizes the value of the argument and can be
implemented using functions in standard libraries such as
numpy.argmin() in Python. In most cases, a linear polynomial(i.e., K = 1) is sufficient for accurately obtaining the value of m,eliminating the errors described in Eq. (29), and extending thetemperature and phase stability beyond the thresholds shown in
Eqs. (31) and (34).
6. EXPERIMENT
A demonstration CHSP was assembled in a 1:2 retarder configuration (Fig. 7) using a set of three 60-wave quartz retarders
and a Glan–Thompson linear polarizer. The first retarder wasoriented at 0◦ , and the other two oriented at 45◦ to effectively
form a 120-wave retarder. The linear polarizer was oriented at
0◦ , and a fiber coupled Ocean Optics USB4000+ spectrometer[20] was used to measure the modulated spectrum. A polarization state generator was constructed using a fiber-coupled
broadband halogen lamp and a linear polarizer. The CHSP
was placed in a heat chamber consisting of a thick protective

Research Article

Vol. 59, No. 30 / 20 October 2020 / Applied Optics

9431

Fig. 7. Demonstration CHSP consists of three 60-wave retarders
assembled in a 1:2 configuration. A polarization state generator on
the left consists of a fiber coupled broadband lamp and polarization
components such as retarders and linear polarizers. The CHSP isplaced inside a heat controlled chamber consisting of a thick covering,a burner stove, and a fan.

covering, a burner stove to generate heat, and fan to circulate air,as well as an alcohol thermometer to record the temperature.
The CHSP was calibrated at 21◦ C using the two-point reference method described in Section 3. A linear polarizer oriented
at 18◦ was placed in front of the broadband lamp, and spectra
were recorded at 14 different temperatures as the burner heatedthe chamber from 21◦ C to 62◦ C in one hour. The spectra were
processed using both the conventional and modified versionsof the self-calibration algorithm assuming a linear fit (K = 1).
For both methods, the reconstructed normalized Stokes vector
components were averaged across the spectrum and plotted asa function of temperature alongside the theoretical values ofeach component in the dotted lines [Figs. 8(a) and 8(b)]. The
results for the conventional self-calibration method in Fig. 8(a)show that for some temperature conditions the self-calibration
yields correct results, while for other temperature conditions,the results have substantial error. In Fig. 8(b), the adjusted
self-calibration method is used, and the large errors of Fig. 8(a)do not occur. In Fig. 8(c), the reconstructed phase functions
1ψ 0 (σ0 ) and 1θ 0 (σ0 ) are evaluated at the reference wavenumber and plotted as a function of temperature, alongside the valueof m. The dotted grey lines show the bounds in which the theory
from Section 4 predicts the value of m to change. In this experiment, the value of m started at m = 0 and decreased to m = −1
between the measurements taken at 36◦ C and 39◦ C. The valueof m further decreased to m = −2 between the measurements at
58◦ C and 62◦ C.
Each of the reconstructed Stokes vector components for the
adjusted self-calibration method are plotted in Fig. 8(d) andshow good agreement for S1 and S2 and fair agreement for S3 .
We account for errors in S3 by generalizing the phase ratio ρin Eq. (23) for situations in which the first and second retarder
are held at different temperatures T1 and T2 , respectively. The
new phase ratio ρ1 is determined by using the temperature
coefficients γ defined in Eq. (32):
ρ1 =

1T20
1φ2 1T20 · φ2 · γ==ρ,
1φ1 1T10 · φ1 · γ
1T10

(37)

where 1T10 = T1 − T0 and 1T20 = T2 − T0 are the temperature differences between the first and second retarders and theFig. 8. (a) Spectrally averaged Stokes vector reconstruction of alinear polarizer as a function of temperature when using the conventional self-calibration algorithm. Large errors begin to occur when thetemperature exceeds 37◦ C due to the phase ambiguity. (b) Spectrally
averaged Stokes vector reconstruction of a linear polarizer as a functionof temperature when using the adjusted self-calibration algorithm,
resulting in no abrupt errors. (c) The phase functions 1ψ0 (σ0 ) and1θ0 (σ0 ) evaluated at the reference function alongside the value of m.
(d) The Stokes vector components as a function of wavelength for the
adjusted algorithm are plotted and overlaid.

calibration temperature. In the case where retarders are in thermal contact, we assume 1T10 = 1T20 and ρ1 = ρ. However,in the demonstration system, the retarders are not in thermal
contact, and temperature differences between the retarders
likely exist.
The primary effect of having an inaccurate value of ρ1 is
to change the phase of S23 (σ ), resulting in some S2 (σ ) lightbeing erroneously measured as S3 (σ ) and vice versa. Although
we did not independently measure the temperature of the tworetarders, we show that the errors in S3 (σ ) can be corrected byadjusting ρ1 , the resulting temperature difference between thefirst and second retarders 1T21 = T2 − T1 can be estimated,
and a simulated version of the system with the estimated temperature differences show similar errors compared to the ones
measured experimentally. For each temperature point, ρ1 wasadjusted until S3 (σ ) ≈ 0 since the Stokes generator used a linear
polarizer and should not produce circular polarization. An estimated temperature difference 1T21 was then calculated using
Eq. (37) by assuming the value of ρ = 2 and T0 = 21◦ C and thatthe temperature measured by the thermometer is equal to T2 .
The value of ρ1 required to correct the S3 (σ ) measurement isplotted in Fig. 9(a), the corresponding temperature difference
1T21 is plotted in Fig. 9(b), and the corrected Stokes vector
reconstruction is shown in Fig. 9(c). The results show thatthe temperature difference is the greatest in the early stages ofheating and that it then reduces as the entire chamber heats up.
The system is simulated using parameters from the demonstration example, assuming the temperature difference 1T21

9432

Research Article

Vol. 59, No. 30 / 20 October 2020 / Applied Optics

Fig. 9. (a) Value of the generalized phase ratio parameter ρ1required to correct the experimental data from Fig. 8 so that S3 (σ ) ≈ 0.
We assume the temperature of the second retarder T2 is equal to thetemperature measured in the chamber. (b) Temperature difference 1T21 as determined by combining the data from (a) with
Eq. (36), showing a maximum temperature difference of around 5◦ C.
(c) Experimental Stokes vector data when corrected using ρ1 from(a). (d) The spectrally averaged Stokes vector components from theexperimental data (dotted lines) is compared to data simulated using
the temperature difference from (b) (solid lines), showing that thetemperature differential 1T21 is expected to cause similar errors insimulation as observed experimentally in Fig. 8.

from Fig. 9(b). The Stokes vector was reconstructed assuming
ρ1 ≈ ρ ≈ 2 for both the simulated and experimental systems
and then compared by averaging the data across the spectrum foreach temperature point, as shown in Fig. 9(d). The agreement
between the simulated and experimental results shows that the
primary source of errors in the adjusted self-calibration resultsof the experimental system were due to temperature differencesof a few degrees centigrade between the retarders. While previous literature [11] eliminated the problem of temperature
difference between the retarders by placing them in thermal
contact, our analysis shows that the problem can also be solved
by measuring the temperature of each retarder and correcting
for the temperature differential by using the generalized phaseratio parameter ρ1 .
7. DISCUSSION
The self-calibration instability described in this paper has notbeen, to the best of our knowledge, noted in literature during
prior demonstrations of the algorithm. In Taniguchi et al. [11]the Stokes vector was successfully reconstructed between 5◦ C
and 45◦ C using the self-calibration algorithm. In that experiment, calcite retarders with approximately 35- and 70-wave
retardances were used. Based on the birefringence temperature
coefficients γ for calcite given by Kimura et al. [21] and Eq. (34),this configuration should have a maximum temperature rangeof around 130◦ C. While this is sufficient for stable operation
across the 40◦ C range and might explain why these errors were
not noted, other configurations of practical interest using different birefringent materials or higher-order retarders are notstable over this temperature range and would experience large
errors in the same environmental conditions. In a simulated
example using CdS retarders, Altaqui et al. [14] reconstructed
the Stokes vector using self-calibration without experiencing
threshold stability errors when taking measurements between
5◦ C and 25◦ C. However, in this simulation the second retarder
was held constant at room temperature, and only the temperature of the first retarder was decreased. Due to the nature of theself-calibration algorithm, only the retardance of the second
retarder is directly reconstructed and changes in the first retarder
are indirectly determined. In this case, the changing retardanceof the first retarder does not affect the ability of self-calibration
to reconstruct the retardance of the second retarder, and thethreshold boundary errors do not occur.
As described in Eq. (34), the stable temperature range
depends on the temperature sensitivity of the birefringent
material γ and the order w of the thickest retarder. Both of theseparameters have potential to vary significantly for different
configurations. The maximum temperature range is listed forseveral different configurations in Table 1. For each configuration, it is assumed the retardance is specified at λ = 633 nm, andthe reference wavenumber is σ0 = 1.72 µm−1 . Configurations
using quartz or CdS have the lowest stable operation range,
while the polymer materials from Meadowlark Optics [22]have a range 10× greater than quartz. Based on the results inTable 1, it can be seen that the need to implement the algorithm
in Section 5 depends on the birefringent material, retarder
thickness, and expected temperature operating range.
The algorithm described in Section 5 extends the stabletemperature thresholds for most configurations and practicalapplications. In simulation, the algorithm provided stable andaccurate Stokes vector reconstructions over a 800◦ C range whenusing a linear fit (K = 1), which is a 26× enhancement, andwhen using a quadratic fit (K = 2), the system was stable over
a 1000◦ C range, which is a 32× enhancement. Based on theseresults, we expect that the curve-fitting compensation algorithm
will be suitable for attaining stable results for temperature ranges
greater than 100◦ C for nearly all configurations listed in Table 1.

Table 1. Maximum Stable Temperature Range ina
Degrees Celsius for Different CHSP Configurations

Quartz [15]Calcite [22]MgF2 [15]CdS [23]
Polymer [22]

50λ

100λ

250λ

500λ

75.2
180.1
167.0
24.2
874.7

37.6
90.0
83.5
12.1
437.4

15.0
36.0
33.4
4.8
174.9

7.5
18.0
16.7
2.4
87.5

a
The rows indicate the birefringent material, and the columns indicate theorder of the second retarder.

Research Article
If the algorithm in Section 5 is not implemented, it is important that the value of the bias offset x is 0 or else the temperaturebounds will be offset according to Eq. (34). In this case, theinitial reference temperature may be very close to the upperor lower temperature bound, resulting in limited stable temperature increases or decreases regardless of the birefringent
material being used. For example, a CHSP calibrated at 21◦ C
using a 50λ MgF2 retarder and reconstructed using an algorithm that returns the complex phase within a [0, 2π ) range(x = π) would have a stable temperature operation range of−146◦ C ≤ T < +21◦ C. In this example, the system does
not tolerate temperature increases at all, even though it has alarge overall temperature range. Most complex angle functions
in standard libraries, such as numpy.angle in Python, assumex = 0, but this should be verified, especially when using custom
libraries.
8. CONCLUSION
In this paper, we analyzed the channeled spectropolarimeter
self-calibration algorithm and identified a source of instability
resulting from an ambiguity in the bias value of the reconstructed phase function. The ambiguity leads to large errors in
the reconstructed phase function when the temperature exceeds
a threshold determined by the birefringent material and orderof the retarders. In the case of a system using 60-wave:120-wave
quartz retarders, the temperature range is limited to around
30–40◦ C. The errors can be eliminated by increasing or decreasing the bias value of the reconstructed phase function by integer
multiples of 2π . The exact integer multiple is determinedby comparing the slope of the phase function of the current
measurement to the initial reference measurement. In an experimental system, we demonstrate the existence of large errors inthe self-calibration phase and reconstructed Stokes vector when
the temperature exceeds the temperature thresholds, and correct
it by implementing a slope-comparison algorithm. We identified additional errors as the result of temperature differences
between the retarders and proposed a method for compensating
this type of error if the temperature difference between theretarders is known.
Funding. Japan Society for the Promotion of Science(GR-18001); National Science Foundation (DGE-1746060);
NSF Graduate Research Opportunities Worldwide program.
Acknowledgment. Any opinions, findings, and conclusions or recommendations expressed in this material are thoseof the authors and do not necessarily reflect the views of the
National Science Foundation.
Disclosures. The authors declare no conflicts of interest.

Vol. 59, No. 30 / 20 October 2020 / Applied Optics

9433

REFERENCES
1. D. H. Goldstein, Polarized Light (CRC press, 2017).
2. R. R. Antonucci and J. S. Miller, “Spectropolarimetry and the natureof NGC 1068,” Astrophys. J. 297, 621–632 (1985).
3. F. Waquet, B. Cairns, K. Knobelspiesse, J. Chowdhary, L. D. Travis,
B. Schmid, and M. I. Mishchenko, “Polarimetric remote sensing ofaerosols over land,” J. Geophys. Res. Atmos. 114, D01206 (2009).
4. N. Ghosh, M. F. Wood, S. Li, R. D. Weisel, B. C. Wilson, R. Li, and I. A.
Vitkin, “Mueller matrix decomposition for polarized light assessmentof biological tissues,” J. Biophoton. 2, 145–156 (2009).
5. E. Du, H. He, N. Zeng, M. Sun, Y. Guo, J. Wu, S. Liu, and H. Ma,“Mueller matrix polarimetry for differentiating characteristic featuresof cancerous tissues,” J. Biomed. Opt. 19, 076013 (2014).
6. S. Liu, X. Chen, and C. Zhang, “Development of a broadband Mueller
matrix ellipsometer as a powerful tool for nanostructure metrology,”Thin Solid Films 584, 176–185 (2015).
7. V. Gruev, R. Perkins, and T. York, “CCD polarization imaging sensor with aluminum nanowire optical filters,” Opt. Express 18,
19087–19094 (2010).
8. D. Goldstein, “Mueller matrix dual-rotating retarder polarimeter,”Appl. Opt. 31, 6676–6683 (1991).
9. K. Oka and T. Kato, “Spectroscopic polarimetry with a channeledspectrum,” Opt. Lett. 24, 1475–1477 (1999).
10. J. Craven-Jones, B. M. Way, M. W. Kudenov, and J. A. Mercier,“Athermalized channeled spectropolarimetry using a biaxial
potassium titanyl phosphate crystal,” Opt. Lett. 38, 1657–1659(2013).
11. A. Taniguchi, K. Oka, H. Okabe, and M. Hayakawa, “Stabilization ofa channeled spectropolarimeter by self-calibration,” Opt. Lett. 31,
3279–3281 (2006).
12. H. Okabe, M. Hayakawa, J. Matoba, H. Naito, and K. Oka, “Errorreduced channeled spectroscopic ellipsometer with palm-size
sensing head,” Rev. Sci. Instrum. 80, 083104 (2009).
13. X. Ju, C. Yan, J. Zhang, B. Yang, and W. Xing, “Adaptive correction ofretardations with immunity to alignment errors for a channeled spectropolarimeter,” Appl. Opt. 57, 8134–8142 (2018).
14. A. Altaqui and M. W. Kudenov, “Phase-shifting interferometry-based
Fourier transform channeled spectropolarimeter,” Appl. Opt. 58,
1830–1840 (2019).
15. S. M. Etzel, A. H. Rose, and C. M. Wang, “Dispersion of the temperature dependence of the retardance in SiO2 and MgF2 ,” Appl. Opt. 39,
5796–5800 (2000).
16. N. Hagen, P. Buranasiri, and Y. Otani, “Compensating for nonlinear dispersion in channeled spectropolarimetry,” Appl. Opt. 59,
5032–5040 (2020).
17. N. Hagen, K. Oka, and E. L. Dereniak, “Snapshot Mueller matrix
spectropolarimeter,” Opt. Lett. 32, 2100–2102 (2007).
18. J. Craven and M. W. Kudenov, “False signature reduction inchanneled spectropolarimetry,” Opt. Eng. 49, 053602 (2010).
19. G. Ghosh, “Dispersion-equation coefficients for the refractive index
and birefringence of calcite and quartz crystals,” Opt. Commun. 163,
95–102 (1999).
20. “Ocean Insight,” 2020, https://www.oceaninsight.com/products/
spectrometers/usb-series/.
21. T. Kimura and M. Saruwatari, “Temperature compensation of birefringent optical filters,” Proc. IEEE 59, 1273–1274 (1971).
22. “Meadowlark Optics,” 2020, https://www.meadowlark.com/store/
PDFs/Retarders.pdf.
23. M. P. Lisitsa, L. F. Gudymenko, V. N. Malinko, and S. F. Terekhova,“Dispersion of the refractive indices and birefringence of CdSx Se1-x
single crystals,” Phys. Status Solidi B 31, 389–399 (1969).