## 1.

## Introduction

Polarization calibration of large telescopes and modern instruments is often limited by the availability of suitable sources for calibration. Several calibration techniques exist using stars or the sun, internal optical systems, or *a priori* knowledge of the expected signals, but each technique has limitations. For altitude–azimuth telescopes, coudé or Nasmyth instruments, or telescopes with off-axis primaries, the polarization calibration usually requires bright, highly polarized sources available over a wide range of wavelengths, altitude–azimuth pointings. In night-time astronomy, polarized standard stars are commonly used, but they provide very limited altitude–azimuth coverage are faint, and have low polarization amplitudes (typically $<5\%$^{1}2.^{–}^{3}). Unpolarized standard stars also exist, but they are also faint and provide limited altitude–azimuth coverage. Solar telescopes can use solardisk-center as a bright, zero-polarization target, provided there is no magnetic field activity. Solar observations often lack bright, significantly polarized targets of known properties. Smaller telescopes can use fixed polarizing filters placed over the telescope aperture to provide known input states that are detected and yield terms of the Mueller matrix as for the Dunn Solar Telescope.^{4}5.6.7.^{–}^{8} Symmetries of spectropolarimetric signatures from the Zeeman effect have been used in solar physics to determine terms in the Mueller matrix.^{9}^{,}^{10} Many studies have either measured and calibrated telescopes, measured mirror properties, or attempted to design instruments with minimal polarimetric defects (cf. Refs. 1112.13.14.15.16.–17). Space-based polarimetric instruments such as Hinode also undergo detailed polarization calibration and characterization.^{18}^{,}^{19}

Many telescopes use calibration optics such as large polarizers, polarizer mosaic masks, polarization state generators, or optical injection systems at locations in the beam after the primary or secondary mirror. For some systems, a major limitation is the ability to calibrate the primary mirror and optics upstream of the calibration system. These systems include telescopes with large primary mirrors, systems without accessible intermediate foci, or many-mirror systems without convenient locations for calibration optics. System calibrations are subject to model degeneracies, coherent polarization effects in the point spread function, and other complex issues such as fringes or seeing-induced artifacts.^{4}^{,}^{5}^{,}^{20}21.^{–}^{22} Modern instrumentation is often behind adaptive optics systems requiring detailed consideration of active performance on polarization artifacts in addition to deconvolution techniques and error budgeting.^{23}24.25.^{–}^{26} Modeling telescope polarization is typically done either with simple single-ray traces using assumed mirror refractive indices or with ray-tracing programs such as Zemax.^{11}^{,}^{24}^{,}^{27}

Every major observatory addresses a diversity of scientific cases. Often, cross-talk from the optics limits the polarization calibration to levels of 0.1% to $>1\%$ in polarization orientation (e.g., ESPaDOnS at CFHT, LRISp at Keck, and SPINOR at DST^{4}^{,}^{28}29.^{–}^{30}). Artifacts from the instruments limit the absolute degree of polarization (DoP) measurements from backgrounds or zero-point offsets. Hinode and the Daniel K Inouye Solar Telescope (DKIST) project outline attempts to create error budgets, calling for correction of these artifacts to small fractions of a percent. The calibration techniques presented here aim to calibrate the cross-talk elements of the Mueller matrix to levels of roughly 1% of the element amplitudes, consistent with internal instrument errors. We also show that the limitation of the method is not the model for the polarization patterns of the sky but the other instrumental and observational issues.

The High-resolution Visible and Infrared Spectrograph (HiVIS) is a coudé instrument for the 3.67-m AEOS telescope on Haleakala, Hawaii. The visible arm of HiVIS has a spectropolarimeter, which we recently upgraded to include charge shuffling synchronized with polarization modulation using tunable nematic liquid crystals.^{31} In Ref. 31, hereafter called H15, we outline the coudé path of the AEOS telescope and details of the HiVIS polarimeter.

The DKIST is a next-generation solar telescope with a 4-m-diameter off-axis primary mirror and a many-mirror folded coudé path.^{32}33.^{–}^{34} This altitude–azimuth system uses seven mirrors to feed light to the coudé lab.^{32}^{,}^{35}^{,}^{36} Its stated scientific goals require very stringent polarization calibration. Operations involve four polarimetric instruments spanning the 380- to 5000-nm wavelength range with changing configuration and simultaneous operation of three polarimetric instruments covering 380 to 1800 nm.^{6}^{,}^{35}36.^{–}^{37} Complex modulation and calibration strategies are required for such a multi-instrument system.^{35}^{,}^{36}^{,}^{38}39.40.^{–}^{41} With a large off-axis primary mirror, calibration of DKIST instruments requires external (solar, sky, and stellar) sources. The planned 4-m European Solar Telescope, though on-axis, will also require similar calibration considerations.^{42}43.44.^{–}^{45}

## 1.1.

### Polarization

The following discussion of polarization formalism closely follows Refs. 46 and 47. In the Stokes formalism, the polarization state of light is denoted as a four-vector: ${\mathbf{S}}_{i}={[I,Q,U,V]}^{T}$. In this formalism, $I$ represents the total intensity, $Q$ and $U$ are the linearly polarized intensities along polarization position angles 0 deg and 45 deg in the plane perpendicular to the light beam, respectively, and $V$ is the right-handed circularly polarized intensity. The intensity-normalized Stokes parameters are usually denoted as ${[1,q,u,v]}^{T}={[I,Q,U,V]}^{T}/I$. The DoP is the fraction of polarized light in the beam: $\mathrm{DoP}=\frac{\sqrt{{Q}^{2}+{U}^{2}+{V}^{2}}}{I}=\sqrt{{q}^{2}+{u}^{2}+{v}^{2}}$. For this work, we adopt a term angle of polarization (AoP) from the references on daytime sky polarimetry, which defines the angle of linear polarization (AoP) as $ATAN(Q/U)/2$. The Mueller matrix is a $4\times 4$ set of transfer coefficients that describes how an optic changes the input Stokes vector (${\mathbf{S}}_{{i}_{\text{input}}}$) to the output Stokes vector (${\mathbf{S}}_{{i}_{\text{output}}}$): ${\mathbf{S}}_{{i}_{\text{output}}}={\mathbf{M}}_{ij}{\mathbf{S}}_{{i}_{\text{input}}}$. If the Mueller matrix for a system is known, then one inverts the matrix to recover the input Stokes vector. One can represent the individual Mueller matrix terms as describing how one incident polarization state transfers to another. In this paper, we will use the notation

## 1.2.

### Daytime Sky as a Calibration Target

The daytime sky is a bright, highly linearly polarized source that illuminates the telescope optics similar to distant targets (sun, stars, satellites, and planets) starting with the primary mirror. A single-scattering Rayleigh calculation is often adequate to describe the sky polarization to varying precision levels and is introduced in great detail in several text books (e.g., Refs. 48 and 49).

There are many atmospheric and geometric considerations that change the skylight polarization pattern. The linear polarization amplitude and angle can depend on the solar elevation, atmospheric aerosol content, aerosol vertical distribution, aerosol scattering phase function, wavelength of the observation, and secondary sources of illumination such as reflections off oceans, clouds, or multiple scattering.^{50}51.52.53.54.55.56.57.58.59.60.61.62.63.64.65.^{–}^{66} Anisotropic scattered sunlight from reflections off land or water can be highly polarized and temporally variable.^{67}68.69.70.71.^{–}^{72} Aerosol particle optical properties and vertical distributions also vary.^{58}^{,}^{73}74.75.76.77.78.79.80.81.^{–}^{82} The polarization can change across atmospheric absorption bands or can be influenced by other scattering mechanisms.^{83}84.85.86.^{–}^{87} Deviations from a single-scattering Rayleigh model grow as the aerosol, cloud, ground, or sea-surface scattering sources affect the telescope line-of-sight. Clear, cloudless, low-aerosol conditions should yield high linear polarization amplitudes and small deviations in the polarization direction from a Rayleigh model. Observations generally support this conclusion.^{88}89.90.91.92.93.94.95.^{–}^{96} Conditions at twilight with low solar elevations can present some spectral differences.^{79}80.81.^{–}^{82}^{,}^{97}

An all-sky imaging polarimeter deployed on Haleakala also shows that a single-scattering sky model is a reasonable approximation for DKIST and AEOS observatories.^{95}^{,}^{96} The preliminary results from this instrument showed that the AoP agreed with single-scattering models to better than 1 deg in regions of the sky more than 20% polarized. We will show later how to filter data sets based on several measures of the daytime sky properties to ensure that second-order effects are minimized. The daytime sky DoP was much more variable, but as shown in later sections, the DoP variability has minimal impact on our calibration method. More detailed models include multiple scattering and aerosol scattering and are also available using industry standard atmospheric radiative transfer software such as MODTRAN.^{98}^{,}^{99} However, the recent studies on Haleakala applied measurements and modeling techniques to the DKIST site and found that the AoP was very well described by the single-scattering model for regions of the sky with DoP greater than 15%.^{95}^{,}^{96} The behavior of the DoP was much more complex and did not consistently match the single-scattering approximation. The technique we developed uses only the AoP.

## 1.3.

### Single-Scattering Sky Polarization Model

Sky polarization modeling is well represented by simple single-scattering models with a few free parameters. The simplest Rayleigh sky model includes single scattering with polarization perpendicular to the scattering plane. A single scale factor for the maximum degree of linear polarization (${\delta}_{\mathrm{max}}$) scales the polarization pattern across the sky.

The all-sky model requires knowing the solar location and the scale factor (${\delta}_{\mathrm{max}}$) to compute the DoP and the AoP projected onto the sky. The geometry of the Rayleigh sky model is shown in Fig. 1. The geometrical parameters are the observer’s location (latitude, longitude, and elevation) and the time. The solar location and relevant angles from the telescope pointing are computed from the spherical geometry in Fig. 1. The maximum DoP (${\delta}_{\mathrm{max}}$) in this model occurs at a scattering angle ($\gamma $) of 90 deg. The Rayleigh sky model predicts the DoP ($\delta $) at any telescope pointing (azimuth and elevation) as

## (2)

$$\delta =\frac{{\delta}_{\mathrm{max}}\text{\hspace{0.17em}}{\mathrm{sin}}^{2}\gamma}{1+{\mathrm{cos}}^{2}\gamma},$$## 1.4.

### Solving for Telescope Mueller Matrix Elements

We model the $3\times 3$ cross-talk elements ($quv$ to $quv$ terms) of the Mueller matrix as a rotation matrix.^{100} This method makes the assumption that a telescope with weakly polarizing optics can have a Mueller matrix that is well represented by a rotation matrix. We find cross-talk of 100%, but the induced polarization and depolarization terms are less than 5%. A rotation matrix has been a good fit to our past data, is predicted by our Zemax modeling, and is easily described with three Euler angles to produce the nine terms of the cross-talk matrix.^{31}^{,}^{101} We also perform a sensitivity analysis in later sections to show that this approximation is reasonable. We find in the appendices that we can neglect the first row and column of the Mueller matrix as the correction to the inner $quv$ to $quv$ terms is second-order in these neglected terms.

For our procedure, all Stokes vectors are scaled to unit length (projected onto the Poincaré sphere) by dividing the Stokes vector by the measured DoP. This removes the residual effects from changes in the sky DoP, telescope-induced polarization, and depolarization. Since we ignore the induced polarization and depolarization, we consider only the $3\times 3$ cross-talk elements as representing the telescope Mueller matrix. We denote the three Euler angles as ($\alpha ,\beta ,\gamma $) and use a short-hand notation where $\mathrm{cos}(\gamma )$ is shortened to ${c}_{\gamma}$. We specify the rotation matrix (${\mathbb{R}}_{ij}$) using the $ZXZ$ convention for Euler angles as

## (3)

$${\mathbb{R}}_{ij}=\left(\begin{array}{ccc}{c}_{\gamma}& {s}_{\gamma}& 0\\ -{s}_{\gamma}& {c}_{\gamma}& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& {c}_{\beta}& {s}_{\beta}\\ 0& -{s}_{\beta}& {c}_{\beta}\end{array}\right)\left(\begin{array}{ccc}{c}_{\alpha}& {s}_{\alpha}& 0\\ -{s}_{\alpha}& {c}_{\alpha}& 0\\ 0& 0& 1\end{array}\right)=\left(\begin{array}{ccc}{c}_{\alpha}{c}_{\gamma}-{s}_{\alpha}{c}_{\beta}{s}_{\gamma}& {s}_{\alpha}{c}_{\gamma}+{c}_{\alpha}{c}_{\beta}{s}_{\gamma}& {s}_{\beta}{s}_{\gamma}\\ -{c}_{\alpha}{s}_{\gamma}-{s}_{\alpha}{c}_{\beta}{c}_{\gamma}& -{s}_{\alpha}{s}_{\gamma}+{c}_{\alpha}{c}_{\beta}{c}_{\gamma}& {s}_{\beta}{c}_{\gamma}\\ {s}_{\alpha}{s}_{\beta}& -{c}_{\alpha}{s}_{\beta}& {c}_{\beta}\end{array}\right).$$## (4)

$${\mathbf{S}}_{i}=\left(\begin{array}{c}{q}_{m}\\ {u}_{m}\\ {v}_{m}\end{array}\right)={\mathbf{M}}_{ij}{\mathbf{R}}_{j}=\left(\begin{array}{ccc}QQ& UQ& VQ\\ QU& UU& VU\\ QV& UV& VV\end{array}\right)\left(\begin{array}{c}{q}_{r}\\ {u}_{r}\\ 0\end{array}\right).$$When using this rotation matrix approximation for the telescope Mueller matrix, the Rayleigh sky input Stokes parameters multiply each term of the rotation matrix to give a system of equations for the three Euler angles ($\alpha $, $\beta $, $\gamma $). This system of equations can be solved using a normal nonlinear least-squares minimization by searching the ($\alpha $, $\beta $, $\gamma $) space for minima in squared error. This direct solution of this set of equations using standard minimization routines is subject to several ambiguities that affect convergence using standard minimization routines. The details of our methods for deriving Euler angles and an example of how one could plan sky calibration observations are outlined in the Appendix of Ref. 100.

Equating Mueller matrix elements to rotation matrix elements, we can write the system of equations for the three Euler angles. This system of equations can be solved using a normal nonlinear least-squares minimization by searching the ($\alpha $, $\beta $, $\gamma $) space for minima in squared error. With the measured Stokes vector (${\mathbf{S}}_{i}$), $i=(\mathrm{1,2},3)$, the Rayleigh sky input vector (${\mathbf{R}}_{j}$), $j=(\mathrm{1,2})$, and a rotation matrix (${\mathbb{R}}_{ij}$), we define the error ($\epsilon $) as ${\epsilon}^{2}(\alpha ,\beta ,\gamma )=\sum _{i=1}^{3}\sum _{j=1}^{2}{[{\mathbf{S}}_{i}-{\mathbf{R}}_{i}{\mathbb{R}}_{ij}(\alpha ,\beta ,\gamma )]}^{2}$. For $n$ measurements, this gives us $3\times n$ terms. This solution is easily solvable in principle but has ambiguities. An alternative method for the direct least-squares solution for Euler angles is done in two steps. First, we directly solve a system of equations for the Mueller matrix elements that are not subject to rotational ambiguity. With the estimated Mueller matrix elements in hand, we can then perform a rotation matrix fit to the derived Mueller matrix element estimates. This two-step process allows us to use accurate starting values to speed up the minimization process and to resolve Euler angle ambiguities. When deriving the Mueller matrix elements of the telescope, one must take care that the actual derived matrices are physical. For instance, there are various matrix properties and quantities one can derive to test the physicality of the matrix.^{102}103.104.^{–}^{105} Noise and systematic errors might give overpolarizing or unphysical Mueller matrices. By fitting a rotation matrix, we avoid unphysical matrices.

The normal solution for Mueller matrix elements can be computed via the normal least-squares method. We can rearrange the time-varying Rayleigh sky inputs to (${\mathbf{R}}_{ij}$) for $i$ independent observations and $j$ input Stokes parameters. The measured Stokes parameters (${\mathbf{S}}_{i}$) become individual column vectors. The unknown Mueller matrix elements are arranged as a column vector by output Stokes parameter (${\mathbf{M}}_{j}$). If we write measured Stokes parameters as (${q}_{{m}_{i}},{u}_{{m}_{i}},{v}_{{m}_{i}}$) and the Rayleigh input Stokes parameters as (${q}_{{r}_{i}},{u}_{{r}_{i}}$), we can explicitly write a set of equations for two Mueller matrix elements

## (5)

$${\mathbf{S}}_{i}=\left(\begin{array}{c}{q}_{{m}_{1}}\\ {q}_{{m}_{2}}\\ {q}_{{m}_{3}}\end{array}\right)={\mathbf{R}}_{ij}{\mathbf{M}}_{j}=\left(\begin{array}{cc}{q}_{{r}_{1}}& {u}_{{r}_{1}}\\ {q}_{{r}_{2}}& {u}_{{r}_{2}}\\ {q}_{{r}_{3}}& {u}_{{r}_{3}}\end{array}\right)\left(\begin{array}{c}QQ\\ UQ\end{array}\right).$$## (6)

$$\frac{\partial E}{\partial {\mathbf{M}}_{j}}=2\sum _{i}{\epsilon}_{i}\frac{\partial {\epsilon}_{i}}{\partial {\mathbf{M}}_{j}}=-2\sum _{i}{\mathbf{R}}_{ij}({\mathbf{S}}_{i}-\sum _{k}{\mathbf{R}}_{ik}{\mathbf{M}}_{k})=0.$$This solution is stable provided a diverse range of input states are observed to give a well-conditioned inversion. The noise properties and inversion characteristics of this equation can be calculated in advance of observations and optimized. We can write the matrix $\mathbf{A}$ with an implied sum over $i$ observations for each term. As an example for a single element, if we compute the inverse of $\mathbf{A}$ and multiply out ${\mathbf{A}}^{-1}$ for the $QQ$ term, we can write

## (7)

$$\mathbf{A}={\mathbf{R}}^{T}\mathbf{R}=\left(\begin{array}{cc}{q}_{{r}_{i}}{q}_{{r}_{i}}& {q}_{{r}_{i}}{u}_{{r}_{i}}\\ {q}_{{r}_{i}}{u}_{{r}_{i}}& {u}_{{r}_{i}}{u}_{{r}_{i}}\end{array}\right),{A}_{11}=QQ=\frac{({q}_{{r}_{i}}{q}_{{m}_{i}})({u}_{{r}_{i}}{u}_{{r}_{i}})-({u}_{{r}_{i}}{q}_{{m}_{i}})({q}_{{r}_{i}}{u}_{{r}_{i}})}{({q}_{{r}_{i}}{q}_{{r}_{i}})({u}_{{r}_{i}}{u}_{{r}_{i}})-({q}_{{r}_{i}}{u}_{{r}_{i}})({q}_{{r}_{i}}{u}_{{r}_{i}})}.$$The solution to the equations for the three sets of Mueller matrix elements is outlined in the Appendix of Ref. 100. In this manner, we can easily implement the usual matrix formalism with a time-series of daytime sky observations to measure six Mueller matrix elements.

## 2.

## Single-Scattering Model Limitations

The assumption of a single-scattering model for computing the daytime sky polarization is incorrect under some circumstances. Multiple scattering, contributions from multiple light sources (upwelling, cloud reflections, and ocean reflections) all complicate the computation of the DoP and associated linear polarization angle. In this section, we outline a second-order scattering model and show how this model can be used to choose calibration observations to avoid such issues. By planning observations in regions of the sky where multiple scattering issues are minimized, this calibration technique can be efficiently used with a simple single-scattering model.

## 2.1.

### Multiple-Scattering Models

We show here that the common two-component multiple-scattering model imparts minimal changes to the AoP in wide regions of the sky. Several additions to the single-scattering model are possible but behave similarly. Along any line of sight in the sky, there are contributions from the single-scattered sunlight along with multiply scattered light off a range of airborne- and ground-based sources as well as extinction. Contributions from Mie scattering of water droplets, ice crystals, or large aerosols modify the models in complex ways. As an example of the variations between the single-scattered Rayleigh model and a simple multiple-scattering model, we follow the mathematical formalism of Ref. 66 to derive general properties of the polarization imprinted from the most common multiple scattering source.

In their notation, they use $\zeta $ to denote the location of a point on a stereographic projection of the sky. In Cartesian geometry, $\zeta =x+iy$. In polar coordinates, $\zeta =r\hspace{0.17em}{e}^{i\phi}$. In Ref. 66, they used the term $w$ to represent the polarization pattern across the sky. By breaking the exponential equation into an amplitude term $|w|$ and a complex orientation term $\gamma (\zeta )$, they represent the stereographic projection for the sky polarization pattern as $w(\zeta )=|w|{e}^{2i\gamma (\zeta )}$. For the single-scattering case, this simple relation behaves as ${\zeta}^{2}$ and can be scaled to an amplitude of 1 and written in polar coordinates ($r,\phi $) as $w(\zeta )\sim {\zeta}^{2}={r}^{2}{e}^{2i(\phi -\frac{\pi}{2})}$.

To add multiple scattering to this equation, we must consider the shift of the zero polarization points away from the solar and antisolar locations. These zero points are Brewster and Babinet points near the sun as well as the Arago and second Brewster point near the antisolar location. Several empirical results show that the singularities are found above and below the sun along the solar meridian. This generally follows from the empirical result that double scattering is the dominant contribution to multiple scattering in the typical locations surveyed. This double-scattering contribution is generally polarized in the vertical direction as it represents the light scattered into the line of sight from the integrated skylight incident on all points along the line of sight. When the sun is low in the horizon, the low DoP regions of the sky are also low on the horizon. This double-scattering contribution is of the same amplitude as the single-scattered light when the single-scattered light is weak and horizontally polarized, which occurs above and below the sun at low solar elevations during sunrise and sunset.

The simplest perturbation to the model is to add a constant that represents a small additional polarization of assumed constant orientation, denoted as $A$. Following Ref. 66, the zero polarization singularities fall at the locations of $\zeta =\pm iA$, which corresponds to a Cartesian $y$ value of $\pm A$. To make the singularities at the antisun location, the equation was generalized to $w(\zeta )\sim ({\zeta}^{2}+{A}^{2})({\zeta}^{2}+\frac{1}{{A}^{2}})$.

A simple example of this two-term scattering model is shown in Fig. 3. The stereographic projection convention has been used. In this case, we put the sun on the horizontal axis to match the North = up convention of Fig. 2. However, in this formalism, the Stokes $qu$ parameters are not referenced to the altitude–azimuth frame and there is no singularity at the zenith. An angular splitting of 27 deg was chosen, and the sun is at an elevation of 89 deg. This solar elevation puts the sun in the center of the image with the horizon projected on the edge of the circle.

The calibration method we have pursued is based on the assumption that the AoP of the sky polarization pattern is known as a modeled input parameter with a high degree of accuracy. Variations between the single Rayleigh scattering model and the real input Stokes vector can cause errors in our calibration methodology.

Figure 4 shows the AoP variations between the simple single-scattering model and the multiple-scattering model, considering the double scattering term in a stereographic projection for a range of multiple-scattering models. In the regions of highest DoP at scattering angles of 90 deg, the difference between this second-order model and the simple Rayleigh model is less than 0.001 deg, as can be seen in Fig. 5. The band of high DoP following the 90-deg scattering angle arc shows similar agreement in AoP. Regions near the neutral points show strong angular variation. This is in agreement with the all-sky imaging polarimetry on Haleakala.^{95}^{,}^{96}

In this section, we outlined a second-order scattering model that included two components contributing to the polarization pattern of the sky. We showed that by choosing regions of the sky with high DoP, one can avoid several contaminations of the AoP to a small fraction of a degree as shown in Fig. 5. Choose points near the 90 deg scattering plane and away from the horizon with high airmass to avoid multiple scattering contamination when using this calibration technique.

## 2.2.

### Planning Sky Observations for Diversity and Efficiency

This technique requires a diversity of input polarization angles to minimize noise propagation when deriving telescope Mueller matrices. There is an analogy between the time-dependent Rayleigh sky input polarization and the retardances chosen to create an efficient modulation scheme for polarization measurements. By choosing telescope pointings and observing times such that the solution for Mueller matrix elements is well conditioned (efficient modulation by the daytime sky), a good calibration can be derived. Polarimeters typically produce intensity modulations by changing the incident polarization state with retardance amplitude and orientation changes. This retardance modulation translated into varying intensities using an analyzer such as a polarizer, polarizing beam splitter, or crystal blocks such as Wollaston prisms or Savart plates. These modulation schemes can vary widely for various optimizations and schemes to maximize or balance polarimetric efficiency over user-chosen Stokes parameters, wavelengths, and instrumentation systems.^{106}107.108.109.110.^{–}^{111} There have been many implementations of achromatic and polychromatic designs in both stellar and solar communities.^{40}^{,}^{41}^{,}^{112}113.114.115.116.117.^{–}^{118} In the notation of these studies, the instrument modulates the incoming polarization information into a series of measured intensities (${\mathbf{I}}_{i}$) for $i$ independent observations via the modulation matrix (${\mathbf{O}}_{ij}$) for $j$ input Stokes parameters (${\mathbf{S}}_{j}$): ${\mathbf{I}}_{i}={\mathbf{O}}_{ij}{\mathbf{S}}_{j}$. This is analogous to our situation in which we changed the matrix indices to be $i$ independent Stokes parameter measurements for $j$ different sky input Stokes parameters: ${\mathbf{S}}_{i}={\mathbf{R}}_{ij}{\mathbf{M}}_{j}$.

In most night-time polarimeters, instruments choose a modulation matrix that separates and measures individual parameters of the Stokes vector, typically called a Stokes definition modulation sequence

## (8)

$${\mathbf{O}}_{ij}=\left(\begin{array}{cccc}1& +1& 0& 0\\ 1& -1& 0& 0\\ 1& 0& +1& 0\\ 1& 0& -1& 0\\ 1& 0& 0& +1\\ 1& 0& 0& -1\end{array}\right).$$Other instruments choose a wide range of modulation schemes to balance the efficiencies over a number of exposures. One recovers the input Stokes vector from a series of intensity measurements by inverting the modulation matrix ($\mathbf{O}$) via the normal least squares formalism: $\mathbf{S}=\frac{{\mathbf{O}}^{T}\mathbf{I}}{{\mathbf{O}}^{T}\mathbf{O}}$. The demodulation matrix is typically defined as ${\mathbf{D}}_{ij}={[{\mathbf{O}}^{T}\mathbf{O}]}^{-1}{\mathbf{O}}^{T}$.

In our daytime sky technique, the Rayleigh sky input parameters become the modulation matrix (${\mathbf{O}}_{ij}={\mathbf{R}}_{ij}$) and the formalism for noise propagation developed in many studies such as Refs. 107 and 108 apply. If each measurement has the same statistical noise level $\sigma $ and there are $n$ total measurements, then the noise on each demodulated parameter (${\sigma}_{i}$) becomes ${\sigma}_{i}^{2}=n{\sigma}^{2}\sum _{j=1}^{n}{\mathbf{D}}_{ij}^{2}$. The efficiency of the observation becomes ${\mathbf{e}}_{i}={(n\sum _{j=1}^{n}{\mathbf{D}}_{ij}^{2})}^{-\frac{1}{2}}$.

One must take care with this technique to build up observations over a wide range of solar locations, so the inversion is well conditioned, as outlined in the Appendices of Ref. 100. The path of the sun throughout the day will create regions of little input sky Stokes vector rotation causing a poorly constrained inversion with high condition number. For instance, at our location in the tropics, the sun rises and sets without changing azimuth until it rises quite high in the sky. We are constrained to observing in early morning and late evening with the dome walls raised since we may not expose the telescope to the sun. This causes input vectors at East-West pointings to be mostly $q$ oriented with little rotation over many hours. Observations at other times of the year or at higher solar elevations are required to have a well-conditioned inversion. One can easily build up the expected sky input polarizations at a given observing site with the Rayleigh sky polarization equations. Then, determining the modulation matrix and noise propagation for a planned observing sequence to ensure a well-measured telescope matrix with good signal-to-noise (SNR) is straightforward.

## 3.

## HiVIS Daytime Sky Observing Campaign

From October 2014 to May 2015, we collected a large data set of daytime sky observations with HiVIS using the new liquid crystal charge-shuffling configuration.^{31} We obtained over 1700 measurements in our standard setup with 17 spectral orders and 4000 pixels per order. The daytime sky was observed in a grid of telescope pointings (azimuth elevation combinations).

The first subset of telescope pointings was chosen starting North-South with a 60-deg spacing for azimuths of [060, 120, 180, 240, 300, 360] and elevations of [10, 25, 50, 75]. The second subset of telescope pointings was chosen starting East-West with azimuths of [030, 090, 150, 210, 270, 330] and elevations of [20, 35, 60, 89]. See Ref. 31 for a schematic and optical layout of HiVIS. The solar azimuth and elevation for all observations are shown in Fig. 6. The sun was low in the south for October to December 2014, while nearly passing through the Zenith in May of 2015. These telescope pointings were used during daytime sky observations over several days: October 19, 24, 25, 29, 30; December 01, 11, 14, 15; and May 9, 10, 11, 16, 17, 18, for a total of 15 days spread over 7 months.

## 3.1.

### HiVIS Data Extraction

As part of routine calibration, modulation matrix elements were derived using our polarization calibration unit.^{31} This unit is a wire grid polarizer and a Bolder Vision Optik achromatic quarter-wave plate on computer-controlled rotation-translation stages. This polarization state generator unit is mounted immediately in front of the HiVIS slit and dichroics slit window. An alignment procedure was done during initial installation to find the stepper motor rotation positions where the polarizing axis of the polarizer and the fast axis of the quarter-wave retarder are aligned with the Savart plate at nominal wavelengths. By using a standard sequence of polarizer and quarter-wave plate retarder orientations, six pure Stokes inputs ($\pm q$, $\pm u$, $\pm v$) are used to derive redundant calibration sets.

There is cross-talk in the quarter-wave retarder that can be compensated for by additional fitting techniques, but this effect is also removed by using the daytime sky calibrations.^{31} By demodulating the polarization state generator calibration data at the slit, we decouple the spectrograph polarization response from the telescope. The average system modulation matrix as the average of all October 2013, December 2013, and May 2014 modulation matrices is shown in Fig. 7. There is little variation in the derived modulation matrix within the main observing periods of October, December, or May. For clarity, only the median modulation matrix value for each spectral order is shown. Calibrations are derived by doing spectral averaging (binning) to 50 spectral pixels per order. The variation in individual modulation matrix elements is small ($<0.05$). To remove any effect by varying system modulation matrix elements, we used calibration observations taken for each major observing season. Typically within each run, full calibration sequences were taken daily with little change shown over timescales of days to 2 weeks.

These modulation matrices are used to demodulate the dual-beam charge shuffled measurements (3 exposures, 12 intensity spectra) into individual $quv$ measurements. The algorithm for computing the telescope cross-talk elements assumes that the measurements are projected onto the Poincaré sphere. Each individual demodulated spectrum is divided by the measured DoP to create scaled Stokes vectors with 100% DoP.

## 4.

## Filtering the Data for High Accuracy Calibrations

There are several sources of error present when using daytime sky measurements for computing telescope Mueller matrices. In this section, we outline techniques to reject observations.

One limitation is that single-scattering Rayleigh sky model is only an approximation. In areas of the sky with low DoP, the computed AoP can vary substantially. This fact immediately suggests removing data points with low measured DoP as well as avoiding using the low DoP region of the sky for this technique. Cirrus clouds have been shown to rotate the AoP and also cause strong departures from the Rayleigh single-scattering model. Cirrus clouds are known to decrease the measured DoP in addition to rotating the polarization by a large angle. On Haleakala, occasional small patches of low-laying cumulus can blow over the telescope aperture a few hundred feet above the ground. If a patch of cloud depolarizes a single exposure of a data set, strong deviations from the Rayleigh sky model can be seen. Figure 8 shows the measured DoP and our scattering angle coverage for this observing campaign.

We experienced and must compensate for several types of errors.

• Rayleigh sky model is inaccurate in low DoP regions.

• Cirrus clouds rotate AoP strongly.

• Optical window uncertainty (BK7/Infrasil encoder failure and motor replacement).

• Operators manually point telescope to wrong pointing (no computer feedback).

• Cumulus blow-by in single exposures (often on Haleakala marginal inversion layer).

The measured DoP for each data set is shown in Fig. 12. On three separate days there were thick cirrus clouds that impacted the measured DoP seen in measurements 200 to 300, 400 to 500, and around 800. Light cirrus clouds were present on May 17 and 18. Low-lying cumulus clouds blowing over from the Haleakala crater were possible in certain October and December days. There are several ways to identify and filter this large data set to select quality data.

## 4.1.

### Filter: Measured DoP Threshold

A simple data filter that improves calibration quality is discarding observations showing low measured DoP. Low detected polarization is often an indicator of either bad atmospheric conditions or issues with the data. This calibration technique requires knowledge of the AoP with reasonable precision to keep noise amplification low.

At low DoP values, the AoP uncertainties grow substantially. Figure 9 shows the number of daytime sky observations we have in the data set after several filtering processes are applied. The grid of (azimuth and elevation) observation points was linearly interpolated to a continuous map over all observable (azimuth and elevation) optical geometries. The small black triangles show the position of a stellar target ($\epsilon $ Aurigae) we use for calibration purposes. These triangles show a typical azimuth, elevation track for a target marking each individual data set where independent calibrations are required.

## 4.2.

### Filter: Rayleigh DoP Agreement Threshold

Several geometrical calculations are required for assessing and filtering data using the Rayleigh sky model. Figure 12 shows the angular separation between each measurement pointing and the computed solar location. Note that there are three exposures per complete full Stokes measurement set, so there are only $4890/3=1630$ unique data sets. From the pointing and solar geometry, we derive the input Rayleigh-sky stokes parameters. One way of checking the agreement of the HiVIS data is to compute the Rayleigh sky parameters from the HiVIS measurements at all pointings. We can rearrange the Rayleigh sky polarization equation to give the calculated maximum DoP (${\delta}_{\mathrm{max}}$) from the HiVIS measured DoP ($\delta $) and the scattering angle ($\gamma $): ${\delta}_{\mathrm{max}}=\delta \frac{1+{\mathrm{cos}}^{2}\text{\hspace{0.17em}}\gamma}{{\mathrm{sin}}^{2}\gamma}$.

From this equation we can use the data to calculate a measure of atmospheric conditions (${\delta}_{\mathrm{max}}$), and we can create a data filter to reject HiVIS observations on hazy days with low ${\delta}_{\mathrm{max}}$. Figure 10 shows the computed ${\delta}_{\mathrm{max}}$ as a function of scattering angle derived from the data set. The different color curves show 40%, 60%, and 80% ${\delta}_{\mathrm{max}}$ scalings. The ${\delta}_{\mathrm{max}}$ functions are reasonably constrained by all-sky polarimeter measurements and MODTRAN models.^{95}^{,}^{96} A simple function for the maximum sky DoP ${\delta}_{\mathrm{max}}$ on a clear day is used following Mauna Loa measurements: ${\delta}_{\mathrm{max}}=80\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}-20\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}\times \frac{90\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}}{\text{solar altitude}}$.^{119}^{,}^{120}

By taking this simple relation, a set of data filters can be created. Data points with low ${\delta}_{\mathrm{max}}$ predictions can be rejected as likely influenced by clouds, haze, multiple scattering, and other effects. The typical Rayleigh sky dependence on solar elevation is scaled down by 30% and calculated for every data set to show a minimum acceptable DoP for each measurement. Figure 10 shows the ${\delta}_{\mathrm{max}}$ values computed from the HiVS measurements. There are a large number of points showing a high maximum degree of sky polarization (${\delta}_{\mathrm{max}}$) as expected for a high, dry observing site such as Haleakala.^{119}^{,}^{120} There are clusters of points at low DoP values that correspond to days with cirrus clouds. These points are rejected by data filters.

The measured DoP roughly follows the expected Rayleigh patterns. The polarization is higher at scattering angles approaching 90 deg, and the predicted maximum sky DoP (${\delta}_{\mathrm{max}}$) matches Mauna Loa measurements on cloud-free days.^{119}^{,}^{120} This rule is only approximate as daytime sky polarization is modified near sunrise and sunset as well as by varying solar elevation.

## 4.3.

### Data Filtering Summary

We use several methods for ensuring data integrity and solution consistency. First, we require an SNR threshold for every spectrum at a nominal wavelength. Second, we require a minimum number of observations at each (azimuth and elevation) combination on the telescope pointing grid. Third, we reject observations with a low measured DoP after demodulation. Fourth, we reject observations where the computed sky polarization as estimated by the projected maximum DoP (${\delta}_{\mathrm{max}}$) suggests haze, cloud, or other data contamination. Fifth, we reject (azimuth and elevation) grid points with low AoP diversity to ensure a well-conditioned solution (fit) to each Mueller matrix element estimate. Sixth, an iteration is done to ensure that the remaining observations give consistent Mueller matrix estimates. The angular distance between observations and Rayleigh model is preserved for a system that is not depolarizing. Thus, data sets showing inconsistent angles between the bulk of the observations are rejected. Seventh, an iterative process is followed to ensure the Rotation matrix fits give consistently calibrated measurements. Observations with a residual angle between the calibrated observations and the Rayleigh sky model above a threshold are rejected.

• Reject observations by measured SNR ($>500$ for all $quv$, spectral order 3 after binned $80\times $ to sampling of 50 spectral pixels per order).

• Require several observations to compute the six Mueller matrix element estimates ($>4$ points optimal).

• Reject observations by measured DoP ($>15\%$ after demodulation).

• Reject low ${\delta}_{\mathrm{max}}$ points for agreement with Rayleigh model (40% typical ${\delta}_{\mathrm{max}}$ for a clear day).

• Reject pointings with low $qu$ input diversity ($>20\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$ for each pointing).

• Reject observations where the calibrated Mueller matrix calibrations show high error outliers (${S}_{\mathrm{cal}}\xb7\text{Mangle}<0.1$).

• Reject iteratively by calibrated residual angle between measurements and theory (${S}_{\mathrm{cal}}\xb7R$) until convergence below a threshold angle (e.g., 25 deg) for consistency.

As examples of some of these filters, Fig. 11 shows three different sets of Mueller matrix element estimates. We show the element estimates on the (azimuth and elevation) grid for (1) data from all seasons filtered as in the above list, (2) data from only October but with a minimum of 3 points per (azimuth and elevation) grid point, and (3) data from May with a minimum measured DoP of 30%.

As the filters reject observations based on season, DoP, diversity, or other minimum thresholds, some (azimuth and elevation) grid points become excluded. Inspection of Fig. 11 also shows that there are some points where there is seasonal disagreement. The May and October data sets disagree at elevations of 89 deg.

We find that calibrations are consistent when the data filters are set to the parameters above. Knowing in advance the impact of several atmospheric factors, we outlined here can help when planning a calibration observing campaign. To know how many telescope pointings must be observed, we now assess the functional dependence and errors when interpolating the telescope model to intermediate pointings.

## 5.

## Mueller Matrices With Azimuth, Elevation, and Wavelength

Once the telescope Mueller matrix elements have been estimated on a grid of azimuth and elevation points, the full telescope Mueller matrix must be interpolated to every possible azimuth and elevation combination for every target that requires calibration. As seen in Ref. 31, the Mueller matrix is smooth trigonometric functions of azimuth and elevation. This is caused by the fold mirror axes crossing in the $f/200$ coudé path.

The Mueller matrix elements can be modeled as sin and cos functions with a range of possible forms. As a simple test, we included fits to functions of azimuth and elevation with the forms $\mathrm{sin}(2Az+2El)$, $\mathrm{sin}(Az+El)$, $\mathrm{sin}(2El)$, $\mathrm{sin}(Az)$, $\mathrm{sin}(El)$, constant.

We will show in later sections that these azimuth–elevation trigonometric functions are a natural consequence of the optical design and are the only terms required to fit the Zemax-predicted Mueller matrix elements from the design.

These equations can be expanded using trigonometric identities to include 13 possible terms. Since the domain of the fit is restricted to azimuths of 0 deg to 360 deg and elevations of only 0 to 90, care must be taken regarding the uniqueness of the fit parameters given various combinations of the functions. After testing several functional forms, we found that only the sin($2Az+2El$) and sin($2Az$) terms had significant amplitude. These functional forms are fit to the Mueller matrix element estimates as well as the rotation matrix fits to the estimates. Typically, only very small differences between the individual Mueller matrix estimates and their corresponding trigonometric function fits are seen. We use a shorthand notation $\mathrm{sin}=S$, $\mathrm{cos}=C$ and subscripts for $A=\mathrm{azimuth}$, $E=\text{elevation}$. The best function we found to fit for each Mueller matrix estimate contains these six terms: ${S}_{2A}{C}_{2E},{S}_{2E}{C}_{2A},{C}_{2A}{C}_{2E},{S}_{2A}{S}_{2E},{S}_{2E},{C}_{2E}$.

As an example of the coefficients found when fitting each of the six Mueller matrix element estimates, Fig. 12 shows the terms for spectral order five at a wavelength of 6881 Å. Each particular Mueller matrix element seems to be dominated by only 1 or 2 terms.

The small differences between the trigonometric function fits to the rotation matrices and the original Mueller matrix element estimates show how well the interpolation method reproduces all Mueller matrix element estimates over the azimuth–elevation grid. Figure 13 shows the difference between the trigonometric fits and the original Mueller matrix elements. The differences were computed only for (azimuth and elevation) points where all filters were applied and a valid result was obtained. Some slight variation at elevations of 0 and 90 are seen in a few Mueller matrix elements.

## 5.1.

### Interpolation Scheme Errors: Rotation Refits to Trig Functions

The trigonometric function fits to the rotation matrix element azimuth–elevation dependences cause some interpolation errors. The rotation matrices are refit to ensure that the Mueller matrices are strictly rotation matrices. This adds one more step in processing to ensure that any data calibrated at an arbitrary azimuth–elevation is not corrupted by the interpolation process. In addition, interpolation from our chosen grid to the actual telescope pointing of any desired target adds uncertainty. This interpolation introduces another source of error. The trigonometric fitting functions create errors in the Mueller matrix estimates at all intermediate pointings because they are derived from interpolated rotation matrices. We chose a finely sampled azimuth–elevation grid spacing of 1 deg for this campaign.

To quantify this fitting error, differences between interpolated Mueller matrix elements and the corresponding refit rotation matrix elements were derived for telescope pointings in between the nominal observed azimuth–elevation grid at the maximum angular distance. The cumulative distribution function for these Rotation matrices minus the Mueller matrix estimate residuals is shown in Fig. 12.

## 5.2.

### Wavelength Dependence

The Mueller matrices for HiVIS are smooth functions of wavelength. Figure 14 shows the azimuth–elevation dependence for HiVIS in four spectral orders numbered [0,5,10,15] corresponding to wavelengths of [6260, 6880, 7650, and 8600 Å], respectively.

For wavelengths short of 7500 Å, the linear to circular and circular to linear cross-talk terms are quite large. The $VV$ terms show elevation dependence and are much less than 1. For the last spectral order at 8600 Å, the $VV$ term is nearly 1 and shows negligible dependence on elevation. The $VQ$ and $VU$ terms show nearly $\pm 1$ amplitudes with strong functional dependence on azimuth and elevation at the shorter wavelengths $<7000\text{\hspace{0.17em}\hspace{0.17em}}\AA $, but are nearly 0 at 8600 Å. As wavelengths increase, the polarization response of HiVIS goes from severe linear-to-circular cross-talk to very benign cross-talk approaching the nominal geometrical $qu$ variation expected for an altitude–azimuth referenced coordinate frame.

## 5.3.

### Zemax System Modeling of Azimuth-Elevation Dependence

The optical ray tracing program, Zemax, has the ability to perform fully polarized ray propagation using a Jones formalism. With this program, we can verify that the trigonometric functions in azimuth and elevation fully capture the behavior of the telescope Mueller matrix. We have used this program in the past to create polarization models of HiVIS.^{24}^{,}^{27} With the polarized ray trace function called POLTRACE, a Zemax user can propagate rays from any pupil coordinate ($Px,Py$) to any field coordinate ($Hx,Hy$) in the Zemax file. By selecting a set of fully polarized inputs in the Jones formalism that also correspond to the Stokes vectors ($quv$), one can determine the polarization response of an optical design. With the Zemax programming language, we trace many rays propagated from a grid of coordinates ($Px,Py$) across the pupil through the optical system to the corresponding focal plane. Typical sampling of 10% where the pupil coordinates ($Px,Py$) are scanned in steps of 0.1 achieves 0.0001 level or better agreement to Mueller matrix calculations using more fine pupil sampling. This match between Mueller matrix terms depends on the details of the optical system including optical power, tilted optics, and, in general, the symmetries in the polarization properties of the exit pupil. However, a 0.1 step in $Px,Py$ seems to be a good compromise between model run speed and calculation sensitivity. Tests run at step sizes of 0.01, 0.025, and 0.05 do not vary by more than the fifth decimal place under typical, mostly symmetric, nonvignetted system configurations.

The Zemax electric field calculations in the Jones formalism are turned into Stokes vector formalism for each of the pure $quv$ input states. The POLTRACE function outputs electric field vector amplitudes ($Ex,Ey,Ez$) and phases following the Jones formalism for every ray traced. For simplicity in large $f/\text{number}$ beams, we project this three-dimensional (3-D) field onto a two-dimensional surface ignoring the $z$ components along the direction of propagation. As the HiVIS polarimeter operates at $f/40$, this is a reasonable assumption for this analysis.

The computed Stokes intensity is the square of the electric field amplitudes ($Ex\times Ex+Ey\times Ey$). This incoherent average is also a reasonable approximation for seeing-limited systems or systems not fully sampling the polarized diffraction limited point spread function. Stokes $Q$ goes as the $X$ and $Y$ intensity difference: ($Ex\times Ex-Ey\times Ey$). Stokes $U$ is computed from $X$ and $Y$ electric field amplitudes as well as phase variations: $2.\times Ex\times Ey\times \mathrm{cos}(\delta )$. Stokes $V$ is similarly computed with both $XY$ field amplitudes and phases: $2.\times Ex\times Ey\times \mathrm{sin}(\delta )$. The term $\delta $ represents the phase difference.

A key parameter for determining the polarization response of an optical system is the coating formulation. The output polarization models are very sensitive to the coating model thicknesses. We performed some experiments with predicting AEOS and HiVIS polarization response using various coating formulations. As we do not have access to the coating or coating formulas for the AEOS mirrors, we show a range of representative functions derived from common coatings. Figure 15 shows some of the formulas commonly used in enhanced protected silver mirrors. The phase retardation and diattenuation for these coatings is shown in Fig. 15 is for a 45-deg reflection using enhanced protected silver-coated mirrors.

The retardance of these one- and two-layer coatings matches vendor-provided curves. For a typical two-layer coating, the retardance has two wavelengths where the nominal 180 deg phase crossing occurs. The exact thickness and materials in the coating determine which wavelengths, but overall the crossings typically occur in the blue-green and red-near-infrared regions of the spectrum. Usually the retardance is 10 deg to 30 deg above 180 deg for the intermediate bandpass and below 180 deg outside these wavelengths. Single-layer protective coatings such as the fused silica over silver have a retardation always below 180 deg. This is seen as the red curve in Fig. 15.

As expected, all the enhanced protected silver formulations with ZnS as an over-coating show diattenuation values below 1% for wavelengths longer than 550 nm. This includes all of the double-layer coatings. The formulas with a fused silica over-coating show higher induced polarization levels.

Most multilayer coatings have polarization responses that are strong functions of the angle of incidence. At near-normal angles of incidence, the coatings shown in Fig. 15 will all display the nominal 180 deg phase and minimal diattenuation. As the angle of incidence is increased past 45 deg, the retardation will be over 40 deg above nominal. As Zemax propagates light ray-by-ray through an optical system across both pupil and field, the angle of incidence is accounted for in the ray-trace. Any variations in angle of incidence across a beam footprint from asymmetries, off-axis rays, decentered optics, and/or vignetting will be propagated to the focal plane.

To demonstrate the typical functional dependence of polarization, we computed the system Mueller matrix as functions of azimuth and elevation for all spectral orders HiVIS samples. Figure 16 shows a typical output Mueller matrix. The intensity to polarization and the polarization to intensity terms were linearly scaled to $\pm 1.5\%$. Note that this predicted induced polarization and depolarization is quite small, further supporting the assumption that the telescope is only weakly diattenuating. The polarization to polarization terms ($quv$ to $quv$) have been scaled to $\pm 1$. The functional dependence is exactly as found using our trigonometric functions above. The $qu$ to $quv$ terms have a strong azimuthal dependence. The $qu$ to $v$ terms are dominated by elevation dependence.

To demonstrate how the coating formula changes polarization predictions, we computed Mueller matrix predictions for all coating formulas in Fig. 15. The full range of optical motion possible in the HiVIS system was used to show the dependence of coating retardation on the Mueller matrix element zero-crossings as well as the relationship between coating formula predictions. Figure 17 shows the functional dependence of the predicted Mueller matrices at chosen azimuth–elevation locations. Some terms are dominated by azimuthal dependence and are plotted versus azimuth while other terms are strongly dependent on elevation only and are plotted with elevation. Different coating formulas have very different telescope pointings where minimal or maximal values occur. The sign and even functional form of some Mueller matrix elements can change with coating formula.

## 5.4.

### Mueller Matrix Functional Dependence Summary

The system Mueller matrices are smooth functions of telescope pointing and wavelength. Using relatively simple trigonometric functional dependencies, calibration of any data set at an arbitrary telescope pointing is possible. The errors inherent in projecting from a coarsely sampled azimuth–elevation grid to any arbitrary location are less than or comparable to other errors presented in this paper. More accurate calibration can be obtained by calibrating the telescope at the actual pointings for a priority target.

Zemax modeling has been performed to verify the amplitude and basic functional dependence of the system Mueller matrix with azimuth and elevation given common coating formulas. Good agreement is seen between Zemax predictions and measured amplitudes. The induced polarization and depolarization terms are consistent with results previously presented for HiVIS.^{24}^{,}^{27}^{,}^{31}^{,}^{100}^{,}^{101}^{,}^{121}

## 6.

## Other Instrument Limitations

The HiVIS system also has other limitations to calibration precision arising from instrumental issues. As an example, during this campaign we discovered an optical misalignment that leads to unstable continuum polarization. As part of the InnoPOL campaign, we designed, built, and installed a diffraction limited $f/200$ laser simulator system.^{24} It was discovered with this laser simulator that the HiVIS fore-optics were delivering the beam to the spectrograph such that the pupil image was being vignetted by the echelle grating. Small changes in the illumination caused by guiding errors and atmospheric seeing caused a variable vignetting that influenced the continuum polarization. This vignetting is likely the source of many continuum polarization instabilities reported in Ref. 31. If the continuum polarization is unstable for point sources, stellar continuum estimates are subject to an additional source of error. This error is reduced for continuum sources. Vignetting and illumination of the edges of the optics caused some of the scattered light issues reported in Ref. 31, degrading the polarization calibration of highly polarized sources. This was shown as the asymmetry of continuum polarization between charge shuffled beams shown in Ref. 31.

## 6.1.

### Intensity to quv Cross-Talk

With such a large data set, we were able to investigate the median intensity to $quv$ cross-talk following a simple cross correlation.^{30} The $quv$ spectra are known to contain very large continuum variations with measured DoP above 85%. However, with so many spectra, we can compute the average continuum-subtracted daytime sky polarization spectrum to very high shot noise statistical limits. This allows us as a very sensitive test of the intensity to polarization cross-talk from artifacts in the data reduction pipeline.

These intensity spectra were normalized by a continuum fitting process. The spectra were median-filtered in wavelength. Note that HiVIS has two amplifiers used in reading out each CCD and there are two CCDs in the mosaic focal plane. The median smoothed intensity spectra for each amplifier was fit with an individual polynomial. The intensity spectra were then divided by these polynomial fits to create continuum-normalized spectra. There was some additional small level fluctuation with wavelength that was subsequently removed with a 60-pixel boxcar smooth fit.

The corresponding intensity and average continuum subtracted median $v=V/I$ spectra for selected spectral orders are shown in Fig. 18. The intensity to polarization cross-talk is immediately apparent at levels of 0.1% to 0.3% for these continuum measurements. As shown in Ref. 31, there is substantial blending between the charge shuffled and modulated exposures when looking at continuum sources with a wider slit length. We recently installed several new dichroic masks in addition to slits of different lengths and widths to test and overcome these limitations. The high-sensitivity spatial profiles presented in the Appendix of Ref. 31 show that roughly 15 spatial pixels of separation was required in the old optical configuration to obtain 0.1% intensity contrast. An investigation after the new optical alignment may change these numbers substantially.

## 6.2.

### Calibration at Reduced Spectral Sampling

For this paper, we performed the demodulation and calibration analysis at very low spectral sampling. The HiVIS data were spectrally averaged by $4000\times $ to a single spectral measurement per order. This achieved high SNR for all spectral orders but it does neglect real spectral variation with wavelength. There are known variations with wavelengths for the demodulation matrices as well as the measured sky polarization across spectral orders that have been neglected for this study.

Figure 19 shows the median modulation matrix derived for the $UU$ term for spectral order 3 at 6620 Å. The variation across the spectral order is roughly 0.02 in amplitude. Since the calibration process uses only a single point per spectral order, there are potential errors at all wavelengths arising from using the average value when the modulation varies from 0.045 to 0.065 across the order. Similar amplitude and slowly varying spectral dependence is seen in all terms of the modulation matrix. As the other errors presented in this paper are of order 0.01 to 0.05, this spectral binning is a contributor to the calibration uncertainties. Performing calibration at increased spectral sampling should remove this error source.

## 6.3.

### Modulator Stability

The temporal drifts or other instabilities of an instrument can be a major limitation to any calibration effort. During this campaign, we did not change any HiVIS configurations to minimize system drifts.

The polarization calibration procedures were followed and standard modulation matrices were recorded. Longer term variations in the system modulation matrix recorded over the campaign showed variations in terms of up to 0.05, but this is compensated for by calibrations specific to each night of data collection. This figure shows the difference between the average modulation matrix and the modulation matrix derived for the three main observing sessions (October, December, and May). The variations are roughly 0.06 peak-to-peak for modulation matrices with amplitudes ranging from $-1$ to $+1$.

Although we did perform calibrations for every observing run, these liquid crystals are temperature-dependent and are known to drift. Our internal calibrations show the temperature sensitivity.^{101} Unfortunately, the temperature control of the AEOS coudé room resulted in temperature changes of up to 5°C in some cases. These drifts certainly contributed to calibration errors similar to those presented in Ref. 31.

## 6.4.

### Limitation Summary

In this section, we presented several other instrument calibration limitations specific to HiVIS. The angular error histograms presented^{31} showed the angle between theoretical sky polarization and the calibrated observations. We can compute the angular error between calibrated daytime sky spectra and the single-scattering model for our data set in a similar manner. Figure 20 shows such angular differences with wavelength for a data set at azimuth 180 deg and elevation 50 deg.

Angular variations of $\sim 6\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$ on the Poincaré sphere correspond to Mueller matrix element and Stokes vector values of 0.1. There are several errors presented in this section that can create calibration uncertainties of $\sim 0.1$ in a Stokes parameter. Liquid crystal temperature instability, reduced spectral sampling and optical misalignment could all contribute to errors in calibrated data.

## 7.

## Summary

We have presented a 6-month long campaign of 1600 polarized daytime sky spectra to test the algorithms to compute telescope Mueller matrix calibrations. There are several considerations for planning a calibration observing campaign including the AoP diversity and the functional dependence expected for the telescope Mueller matrix with azimuth and elevation. Plan observations to have a large input polarization angular diversity to provide a well-conditioned matrix inversion (least squares solution) and to ensure good estimates for the Mueller matrix elements. We presented several data-based filters necessary to avoid potential problems with data collection, sky model uncertainties, sky polarization variations, and instrumental limitations.

A single-scattering model for the linear polarization angle of the daytime sky can be used if calibration data are taken away from regions of naturally low sky polarization. This avoids contamination from multiple scattering as shown in our analysis of a second-order scattering model. We showed a two-term scattering model for the sky polarization that consists of a single-scattering model plus a constant polarized background from typical multiple-scattering sources. The regions of sky 90 deg from the sun are highly polarized and the AoP is much less sensitive to multiple scattering. Predictions with large multiple-scattering contributions showed that the AoP agreed to better than 0.01 deg linear polarization rotation in the highly polarized region of the sky near the 90-deg scattering plane.

We showed in the appendices that treating a weakly polarizing telescope Mueller matrix as a rotation matrix is insensitive to induced polarization and depolarization terms. By neglecting induced polarization and depolarization, we fit the $quv$ cross-talk elements as a rotation matrix. We showed that this assumption is second-order in small terms from the first row and column of the Mueller matrix. For our system with many enhanced protected silver-coated mirrors, the induced polarization is predicted by Zemax to less than 2% while cross-talk is 100%. Our observations of unpolarized standard stars support this low number, within the limitations of the optical misalignments reported here.

Several data rejection filters can be applied to ensure quality calibrations. We reject observations with low measured SNR. For this work, we used $\mathrm{SNR}>500$ for all $quv$ in spectral order 3 at 6620 Å after spectral binning by a factor of 80× to sampling of 50 spectral pixels per order. We require at least $n>3$ observations to compute the six Mueller matrix element estimates. We reject observations if the measured DoP is less than 15%. For every observation, a calculation of the Rayleigh sky maximum DoP (${\delta}_{\mathrm{max}}$ ) using the data is required to be above 40% to guarantee a clear sky. Telescope pointings where the data set provides a limited range in the input AoP are also rejected ($<20\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$) because the low angular diversity will amplify errors. Iterative filters based on convergence criteria were shown to reject statistical outliers. After initial calibrations are done to the data set, we reject individual outliers where there is a large residual angle on the Poincaré sphere between the calibrated measurements and the single-scattering theory (${S}_{\mathrm{cal}}\xb7\mathbf{R}$) until the group of points converge below a threshold angle (e.g., 25 deg) or the data become too sparse to compute a Mueller matrix.

We showed how observations taken on a sparse grid azimuth and elevation telescope pointing combinations can be interpolated onto a continuous function set. By using Zemax optical models of AEOS and HiVIS with representative enhanced protected silver mirror coatings, we can show the expected functional dependence of the Mueller matrix with azimuth and elevation. The Zemax models show that simple functional dependence is expected, and the assumption of weakly polarizing optics is both predicted in Zemax and observed by HiVIS using daytime sky calibrations. The interpolation from a sparse grid of (azimuth and elevation) measurements adds some errors, but with careful planning of calibration observations, the errors from interpolation can be minimized or removed by calibrating along the same pointings as the observations.

## Appendices

## Appendix A:

### Second-Order Scattering Model for the Sky

We summarize here the mathematics used to expand a scattering model with an additional constant polarization term in addition to the single-scattering Rayleigh term. In the notation of Ref. 66, they used $\zeta $ to denote the complex location of a point on a stereographic projection of the sky. In Cartesian geometry, $\zeta \text{\hspace{0.17em}\hspace{0.17em}}=x+iy$. In polar coordinates, $\zeta =\mathrm{r}\hspace{0.17em}{e}^{i\phi}$. In Ref. 66, they used the term $w$ to represent the polarization pattern across the sky. By breaking the exponential equation into an amplitude term $|w|$ and a complex orientation term $\gamma (\zeta )$, they represent the stereographic projection for the sky polarization pattern as $w(\zeta )=|w|{e}^{2i\gamma (\zeta )}$. For the single-scattering case, this simple relation behaves as ${\zeta}^{2}$ can be scaled to an amplitude of 1 and written in polar coordinates ($r,\phi $) as $w(\zeta )\sim {\zeta}^{2}={r}^{2}{e}^{2i(\phi -\frac{\pi}{2})}$.

To add multiple scattering to this equation, we must consider the shift of the zero polarization points away from the solar and antisolar locations. These zero points are Brewster and Babinet points near the sun as well as the Arago and second Brewster point near the antisolar location. Several empirical results show that the singularities are found above and below the sun along the solar meridian. This generally follows from the empirical result that double scattering is the dominant contribution to multiple scattering in the typical locations surveyed. This double-scattering contribution is generally polarized in the vertical direction as it represents the light scattered into the line of sight from the integrated skylight incident on all points along the line of sight. When the sun is low in the horizon, the low DoP regions of the sky are also low on the horizon. This double-scattering contribution is of the same amplitude as the single-scattered light when the single-scattered light is weak and horizontally polarized, which occurs above and below the sun at low solar elevations during sunrise and sunset.

The simplest perturbation to the model is to add a constant representing a small additional polarization of assumed constant orientation denoted $A$. Following Ref. 66, the zero polarization singularities fall at the locations of $\zeta =\pm iA$, which corresponds to a Cartesian $y$ value of $\pm A$. To make the singularities at the antisun location, the equation was generalized to

The Reference 66 notation showed the four zero polarization singularities as simple functions of the solar position and the constant $A$. If you denote the solar elevation as $\alpha $ and use the stereographic projection where the sun is on the $y$ axis at the location $\zeta =i{y}_{s}=i(1-\mathrm{tan}\frac{\alpha}{2})(1+\mathrm{tan}\frac{\alpha}{2})$, the four zero polarization singularities are located## (10)

$$(1){\zeta}_{+}=i\frac{{y}_{s}+A}{1-A{y}_{s}}\phantom{\rule[-0.0ex]{1em}{0.0ex}}(2){\zeta}_{-}=i\frac{{y}_{s}-A}{1+A{y}_{s}}\phantom{\rule[-0.0ex]{1em}{0.0ex}}(3)\frac{-1}{{\zeta}_{+}^{*}}(4)\frac{-1}{{\zeta}_{-}^{*}}.$$To make the polarization function symmetric across the sky (antipodally invariant) and to scale the DoP to 100%, the polarization equation can be normalized and written in terms of these singularity locations as

## (11)

$$w(\zeta )=-4\frac{(\zeta -{\zeta}_{+})(\zeta -{\zeta}_{-})(\zeta +\frac{1}{{\zeta}_{+}^{*}})(\zeta +\frac{1}{{\zeta}_{-}^{*}})}{{(1+{r}^{2})}^{2}|{\zeta}_{+}+\frac{1}{{\zeta}_{+}^{*}}||{\zeta}_{-}+\frac{1}{{\zeta}_{-}^{*}}|}.$$## Appendix B:

### Details of Data Set Filters

We outline in this section details of some of the filters applied to the HiVIS data set as we enforce consistency and quality across the daytime sky observing campaign.

## B.1.

#### Filter: Input Angular Diversity Threshold

This computational algorithm relies on the assumption that we can solve for a set of rotation matrix angles to match the input data with the Rayleigh sky model through a least squares process. If the input data lacks enough angular diversity, the least squares solution becomes badly conditioned, leading to very large noise propagation errors. We require that all (azimuth and elevation) grid points in the computation have sufficient range of input polarization angles (large angular diversity). We find that a minimum of $\sim 20\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$ angular diversity is required to give a well-conditioned solution.

The high-elevation telescope pointings near the zenith only have 30 deg of input diversity in our survey largely because the telescope cannot observe the Zenith while the sun is high above the horizon. Given the nearly East-West oriented rise and set azimuths for the sun during the summer, there is not much diversity in the input polarization angles along East-West telescope azimuths. For telescope pointings at low elevations in the north, the telescope sees a large input AoP diversity for all seasons. Given the winter data set with the sun rotating through all azimuths lower in the south, the input diversity is at least $3\times $ higher at these telescope pointings.

## B.2.

#### Filter: Consistency and Convergence of Matrix Element Estimates

By comparing the agreement between individual calibrated measurements and the least-squares solution for the calibration itself, we can identify outliers and reject them from the data set. We call this filter consistency because it checks the agreement between each individual measurement and the resulting average over all measurements. Statistical outliers can be rejected given a threshold. One would expect smooth variation in the rotation of the polarization caused by the telescope-induced cross-talk. The measured $quv$ parameters should be a smooth function of input polarization angle. As an example, the measured (demodulated, 100% scaled) Stokes parameters ($quv$) are shown in Fig. 21 for a telescope pointing of azimuth of 60 deg and elevation of 25 deg. There is strong rotation of the measured cross-talk as the input polarization angle varies with solar location throughout the year.

The angular rotation in $quv$ space between the initial input Stokes vector and subsequent measured Stokes vectors is a useful test of the Mueller-matrix-as-a-rotation-matrix assumption. Assuming the cross-talk is confined to a rotation matrix uninfluenced by the induced polarization or depolarization as shown by our sensitivity analysis, the angular rotation between subsequent measurements should also match the angular rotation of the modeled Rayleigh sky $qu$ inputs. In most cases, the modeled angular rotation in $quv$ space matches the measured angular rotation with errors of less than 10 deg.

## B.3.

#### Filter: Iterative Filtering and Rotation Matrix Solution Consistency

As another independent assessment of the data, we can ensure that the individual data points used to create Mueller matrix estimates give consistent results and are not statistical outliers compared to the average. We can apply a data filter by requiring that the derived rotation matrix calibrates each individual measurement to within some threshold tolerance of the average. To check the consistency of the telescope Mueller matrix solution, we create an iterative process. We apply the above processes of fitting the telescope Mueller matrix elements and calibrating all individual polarization measurements. Once we derive an initial set of calibrations, we apply the calibration to all individual measurements. With calibrated measurements, we can check the differences between each calibrated individual measurement and the associated Rayleigh sky model. If these differences are above a threshold, we can identify and reject the data as an outlier. An iterative process was written to follow these steps:

• 1. Compute Mueller matrix estimates from all measurements.

• 2. Fit for Euler angles ${\xi}_{i}(\alpha ,\beta ,\gamma )$.

• 3. Compute rotation matrix elements for telescope Mueller matrix approximation.

• 4. Calibrate all data points used to compute Mueller matrix estimates.

• 5. Compute angular difference ($\mathbf{S}\xb7\mathbf{R}$) between all calibrated data and predicted Rayleigh sky.

• 6. Check if all data points are below a threshold value for angular differences (Convergence?).

• 7. If no points are rejected, stop.

• 8. If points are above the threshold, reject the data point with highest angular difference (nominally 25 deg).

• Repeat all steps with newly filtered data set until convergence criteria are obtained.

Outliers are rejected and the process is repeated until convergence criteria are met. Note that in this formalism, the Rayleigh sky model and the measured demodulated, scaled Stokes parameters both have 100% DoP and are by definition vectors with a length of 1. We can calculate the angle between the calibrated measurements and the model sky polarization at each measurement ($i$) with a simple dot product as: ${\mathrm{Q}}_{i}=\mathrm{ACOS}({\mathbf{R}}_{i}\xb7{\mathbf{S}}_{i})$.

The first steps in computing Mueller matrix elements were shown in Fig. 11. These Mueller matrix estimates are then fit to a rotation matrix via the least squares process. These rotation matrices are used to calibrate all individual $quv$ measurements. The rotation error among these calibrated measurements and the Rayleigh sky is computed. As shown in Fig. 14 of Ref. 31, the fitting of the Mueller matrix elements to a rotation matrix results in some differences. The six individual Mueller matrix element estimates here are typically within 0.1 of the final rotation matrix fit values. Depending on the number of points, filters, input polarization angular diversity, and other factors, the cumulative error distribution functions show that 80% of the Mueller matrix estimates are between 0.03 and 0.1 of the final fit value.

## B.4.

#### Uncertainty in Optical Window

We have two different window options (BK7 and Infrasil). As outlined in the schematic of Ref. 31, the telescope has a wheel for different window substrates that separate the coudé lab from the telescope optical feed. The windows are mounted in between the sixth and seventh mirrors in the optical train. The windows are permanently mounted in the wheel and are in the vertical orientation as the incoming $f/200$ beam travels downward. These two windows were used on alternating days. The AEOS operators had a mechanical issue discovered only after our October observing run. The issue was with the window wheel and caused uncertainty about which window was used on which day in October. A new motor and computer controlled wheel rotation were installed between the December and May runs. We intentionally tested both windows. With operator confusion and uncertainty for the observing run, we will treat the window parameters as unknown in our October measurements. Certainly the window was fixed in one location on a particular day, but it could have changed configurations between observing days.

Several observing runs were done to assess the impact of changing windows from the BK7 window to the Infrasil window. Data were collected on several days, often interchanging windows. Figure 22 shows the measured $quv$ spectra taken in May where blue represents Infrasil and red represents BK7. Some data points show deviations from the smooth curves, but no systematic difference between window substrates is seen. Outliers from the smooth curves are likely caused by a variety of observing-related issues and can be easily removed via data filters.

## Appendix C:

### Sensitivity to Assuming Polarization Preservation: A Perturbation Analysis

Another assumption of this method is neglecting consideration of the first row and first column of the Mueller matrix. These intensities to/from $quv$ terms can introduce errors. In this section, we show a simple sensitivity analysis and demonstrate that neglecting these terms is a good assumption. The measured induced polarization and depolarization are low, giving us some confidence in this assumption. For high reflectivity mirrors with low diattenuation, as our mirrors are, the predicted Mueller matrix first row and column are typically low. Cross-talk values are 100% while induced polarization is less than 1% to 3%. Our method of fitting the $3\times 3$ cross-talk elements of a Mueller matrix with a rotation matrix is relatively robust against errors in the first row and column of the Mueller matrix. The analysis below will show that the $3\times 3$ cross-talk elements of the Mueller matrix are minimally impacted to second-order when including nonzero elements in the first row and column of the Mueller matrix.

A Mueller matrix is not a random combination of 16 numbers arranged in a square. A Mueller matrix is a group of numbers that must transform a Stokes vector into another Stokes vector through multiplication according to a set of rules that preserves certain properties of the vector. This means that any Mueller matrix is a transformation matrix that must belong to a group of transformations and must behave according to a set of rules. For instance, a familiar constraint would be that the Stokes vector cannot have greater than 100% polarization, expressed as ${I}^{2}-{Q}^{2}-{U}^{2}-{V}^{2}>0$. This set of rules means Mueller matrices are part of a group, formally called $\mathrm{SO}{(\mathrm{1,3})}^{+}$ as part of the orthochronous Lorentz group. This group describes a set of functions that transform Stokes vectors in a given four-dimensional space. Conveniently, it also provides many conditions that a Mueller matrix must obey to be physically realizable. It also provides specific frameworks (technically, a Lie algebra) for the convenient mathematical analysis of the Mueller matrix since you can generate any matrix as the exponential of a few simple constants multiplied by what are typically called generators. These generators are simple matrices that can be multiplied by a constant (for instance, a rotation angle) to define a physically realizable Mueller matrix. With these generators, we can do error propagation and we can test for the sensitivity of any approximation to small perturbations. In our method, we write the Mueller matrix in terms of the Euler angles ($\alpha ,\beta ,\gamma $) and the generators used to create the rotation matrix. In our $ZXZ$ convention for Euler angles, we rotate first about the $Z$ axis, then about the $X$ axis, then about the $Z$ axis in the rotated coordinates. To compute the Mueller matrix via exponentials, we use the generators for rotations about the $Z$ axis and about the $X$ axis.

We use the standard exponential notation to describe a Mueller matrix as an exponential of the three Euler angles (${d}_{i}$) for $i=\mathrm{1,2},3$: $M=\mathrm{exp}(\alpha {S}_{1}+\beta {S}_{2}+\gamma {S}_{3})=\mathrm{exp}(\mathrm{\Sigma}{d}_{i}{S}_{i})$. The $S$ matrices describe the standard infinitesimal generators for rotation for $1=Z$, $2=X$, $3=Z$ with the $ZXZ$ convention. The $Z$ rotation for ${S}_{1}$ is often denoted ${J}_{3}$ and is

## (12)

$${\mathbf{S}}_{i}={J}_{3}=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& -1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\end{array}\right).$$## (13)

$$M=\mathrm{exp}\left(\alpha \right[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& -1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\end{array}\left]\right)=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& {c}_{\alpha}& -{s}_{\alpha}& 0\\ 0& {s}_{\alpha}& {c}_{\alpha}& 0\\ 0& 0& 0& 1\end{array}\right).$$The standard infinitesimal generator for an $X$ rotation is denoted ${J}_{1}$ and if we rotate by an angle $\beta $, we get a Mueller matrix for a linear retarder via $\mathrm{exp}(\beta {J}_{1})$

## (14)

$${\mathbf{M}}_{ij}=\mathrm{exp}\left(\beta \right[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& -1\\ 0& 0& 1& 0\end{array}\left]\right)=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& {c}_{\beta}& -{s}_{\beta}\\ 0& 0& {s}_{\beta}& {c}_{\beta}\end{array}\right).$$## (15)

$${\mathbf{M}}_{ij}=\mathrm{exp}\left(\epsilon \right[\begin{array}{cccc}0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\left]\right)=\left(\begin{array}{cccc}\mathrm{cosh}(\epsilon )& \mathrm{sinh}(\epsilon )& 0& 0\\ \mathrm{sinh}(\epsilon )& \mathrm{cosh}(\epsilon )& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)=\left(\begin{array}{cccc}1& \epsilon & 0& 0\\ \epsilon & 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right).$$The arbitrary form for the generators can be represented by the terms applicable to the first row and column of the Mueller matrix, ${K}_{i}$ and the terms corresponding to rotations, ${J}_{i}$. The full set of generators can be written with a set of infinitessimal operators (often called boosts) in $quv$ denoted as ${\zeta}_{i}$ and a set of rotations in $quv$ with rotations denoted as ${\theta}_{i}$

## (16)

$${M}_{ij}=\mathrm{exp}({\zeta}_{x}{K}_{1}+{\zeta}_{y}{K}_{2}+{\zeta}_{z}{K}_{3}+{\theta}_{x}{J}_{1}+{\theta}_{y}{J}_{2}+{\theta}_{z}{J}_{3})=\mathrm{exp}\left(\right[\begin{array}{cccc}0& {\zeta}_{x}& {\zeta}_{y}& {\zeta}_{z}\\ {\zeta}_{x}& 0& -{\theta}_{z}& {\theta}_{y}\\ {\zeta}_{y}& {\theta}_{z}& 0& -{\theta}_{x}\\ {\zeta}_{z}& -{\theta}_{y}& {\theta}_{x}& 0\end{array}\left]\right).$$## (17)

$${\mathbf{M}}_{ij}=\left(\begin{array}{cccc}1& \epsilon & 0& 0\\ \epsilon & {R}_{11}& {R}_{21}& {R}_{31}\\ 0& {R}_{12}& {R}_{22}& {R}_{32}\\ 0& {R}_{13}& {R}_{23}& {R}_{33}\end{array}\right).$$Note that the Euler angles do not commute, and we have an equation for the Mueller matrix as an exponential of a sum of terms representing the rotation matrices. We can use the Zassenhaus formula, which expands noncommuting exponential functions of sums to an infinite series of terms, similar to other familiar expansions. The formula below includes the first of the additional terms and represents the expansion: $\mathrm{exp}(X+Y)=\mathrm{exp}(X)\mathrm{exp}(Y)\mathrm{exp}(\frac{1}{2}[X,Y])$. Additional terms grow complex quickly but have small amplitudes. The next correction to this equation involves nested commuting: $\frac{1}{12}[X,[X,Y]]+[Y,[Y,X]]$. For our case, the noncommuting terms would be $[X,Y]=[{K}_{1},{S}_{i}]$. We write our Mueller matrix where we denote the sum of noncommuting rotations $\mathrm{exp}(\mathrm{\Sigma}{d}_{i}{S}_{i})$ with ${S}_{i}$ as the $ZXZ$ rotation group as a single matrix ${\mathbb{R}}_{ij}$ of the three Euler angles ${d}_{i}=(\alpha ,\beta ,\gamma )$ for $i=\mathrm{1,2},3$.

## (18)

$${M}_{ij}=\mathrm{exp}(\int {K}_{1})\mathbb{R}{(\alpha ,\beta ,\gamma )}_{ij}\mathrm{exp}[\frac{1}{2}(\int {K}_{1},\mathrm{\Sigma}{d}_{i}{S}_{i})].$$## (19)

$${\mathbf{M}}_{ij}=\mathrm{exp}(\u03f5{K}_{1})=\mathbb{II}+\epsilon {K}_{1}=\left[\begin{array}{cccc}1& \epsilon & 0& 0\\ \epsilon & 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right].$$## (20)

$${\mathbf{M}}_{ij}=\mathrm{exp}[\frac{1}{2}(\epsilon {K}_{1},\mathrm{\Sigma}{d}_{i}{S}_{i})]=\mathrm{exp}[\frac{\epsilon}{2}(\beta {K}_{3}-\gamma {K}_{2})].$$## (21)

$${\mathbf{M}}_{ij}=\mathrm{exp}[\frac{1}{2}\epsilon (\beta {K}_{3}-\gamma {K}_{2})]\sim \mathbb{II}+\frac{\beta \epsilon}{2}{K}_{3}-\frac{\gamma \epsilon}{2}{K}_{2}.$$Combining the three terms for Eq. (18), we get

## (22)

$${\mathbf{M}}_{ij}=(\mathbb{II}+\epsilon {K}_{1})\mathbb{R}(\mathbb{II}+\frac{\beta \epsilon}{2}{K}_{3}-\frac{\gamma \epsilon}{2}{K}_{2}).$$## (23)

$${\mathbf{M}}_{ij}=\mathbb{R}(\mathbb{II}+\epsilon {K}_{1}+\frac{\beta \epsilon}{2}{K}_{3}-\frac{\gamma \epsilon}{2}{K}_{2}).$$We can write the first Mueller matrix correction term ${K}_{1}\mathbb{R}$ as

## (24)

$$\left(\begin{array}{cccc}0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& {R}_{11}& {R}_{21}& {R}_{31}\\ 0& {R}_{12}& {R}_{22}& {R}_{32}\\ 0& {R}_{13}& {R}_{23}& {R}_{33}\end{array}\right)=\left(\begin{array}{cccc}0& {R}_{11}& {R}_{12}& {R}_{13}\\ 1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right).$$This term is scaled by $\epsilon $ and does not include any terms in the $quv$ to $quv$ portion of the Mueller matrix. Similar corrections to the Mueller matrix are seen for the term ${K}_{2}\mathbb{R}$

## (25)

$$\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& {R}_{11}& {R}_{21}& {R}_{31}\\ 0& {R}_{12}& {R}_{22}& {R}_{32}\\ 0& {R}_{13}& {R}_{23}& {R}_{33}\end{array}\right)=\left(\begin{array}{cccc}0& {R}_{12}& {R}_{22}& {R}_{32}\\ 0& 0& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right).$$## Acknowledgments

Dr. Harrington acknowledges support from the InnoPol (Grant No. SAW-2011-KIS-7) from Leibniz Association, Germany, and from the European Research Council Advanced Grant HotMol (ERC-2011-AdG 291659) during the period 2013 to 2015 for this work. Dr. Kuhn acknowledges the NSF-AST DKIST/CryoNIRSP program. This program was partially supported by the Air Force Research Laboratory (AFRL) through partial salary support for Dr. Harrington until October 2015. This work made use of the Dave Fanning and Markwardt IDL libraries.