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–66 Anisotropic scattered sunlight from reflections off land or water can be highly polarized and temporally variable.^67–72 Aerosol particle optical properties and vertical distributions also vary.^58,73–82 The polarization can change across atmospheric absorption bands or can be influenced by other scattering mechanisms.^83–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–96 Conditions at twilight with low solar elevations can present some spectral differences.^79–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 }_{\max }$) scales the polarization pattern across the sky.

The all-sky model requires knowing the solar location and the scale factor (${\delta }_{\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 }_{\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

Fig. 1

The celestial triangle representing the geometry for the sky polarization computations at any telescope pointing. $\gamma$ is the angular distance between the telescope pointing and the sun. ${\theta }_s$ is the solar angular distance from the zenith. $\theta$ is the angular distance between the telescope pointing and the zenith. $\phi$ is the angle between the zenith direction and the solar direction at the telescope pointing. The angle $\psi$ is the difference in azimuth angle between the telescope pointing and the solar direction. The input $qu$ components are derived as sin and cos of $\psi$, respectively. The law of cosines is used to solve for any angles needed to compute DoP and position AoP.

Fig. 2

Various Rayleigh sky model parameters computed in a range of projections for mid-morning on January 27th, 2010, on Haleakala when the sun was at an elevation of 50.5 deg and an azimuth of 172 deg. The single-scattering model was scaled to a maximum DoP of 100% (${\delta }_{\max }$). (a) The model DoP with white at 100% and black as 0% for all altitudes and azimuths plotted on a Cartesian rectangular grid. (b) The same DoP model data but in an orthographic projection with North up and East left. (a, b) show equivalent DoP data just with differing projections. (d) The scattering angle in an orthographic projects. This scattering angle shows the input linear polarization angle in the reference frame of the telescope, which always has the +elevation axis point. (e, f) $q$ and $u$ in an orthographic projection where white is $+1$ and black is $-1$. The coordinate system for $qu$ was chosen to be $+1$ in the $+\text{altitude}$ direction of an altitude–azimuth system. This system has a singularity at the zenith where the telescope optics can degenerately point to the zenith with any azimuth. The $+\text{altitude}=+q$ system is referenced to the optical train through the orientation of the primary mirror mount against the sky. (g, h, i) Properties of the sky polarization model on the horizon. (g) The DoP with peaks in the East and West. The scattering angle in (h) shows sign changes at the solar azimuth of 172 deg. (i) shows $q$ as the solid black line and $u$ as the dashed line. Stokes $u$ changes sign at the solar azimuth of 172 deg.

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.¹⁰⁰ 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 $\cos (\gamma )$ is shortened to $c_{\gamma }$. We specify the rotation matrix (${\mathbb{R}}_{ij}$) using the $ZXZ$ convention for Euler angles as

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–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

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

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.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.

Fig. 3

The multiple-scattering model with a splitting constant of $\delta =4\mathrm{ATAN}(A)=27\mathrm{deg}$. All projections are stereographic with North up and East left. The sun is at the Zenith. (a) Stokes $q$ and (b) Stokes $u$ in the altitude–azimuth frame. The gray-scale corresponds to $\text{white}=1$ and $\text{black}=-1$. (c) The DoP with black as 0 and white as 100%. The AoP computed as 0.5ATAN($q,u$) is shown in (d). The two polarization zero points are seen as the split singularities near the zenith in the center of the AoP image. The AoP is linearly scaled from black to white from 0 deg to 180 deg.

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

Fig. 4

The angular differences (in deg) between the single-scattering Rayleigh model and the multiple-scattering model outlined above with the double-scattering term. We use a log color scale and chose a separation of $\delta =4\mathrm{ATAN}(A)=12\mathrm{deg}$ and the sun is at 10 deg elevation. A region within 20 deg of the sun has been masked and shows up as a black circle. The scale bar on the right shows the color scheme with white as 0.1 deg AoP angular difference and black as 0.001 deg AoP difference. The line of symmetry between the sun and the antisolar location is a region of minimal difference as is the 90-deg scattering plane shown as the curved black arc in this stereographic projection.

Fig. 5

The AoP variation between the single-scattering Rayleigh model and the multiple-scattering model here considering double scattering. For this figure, the sun was placed on the horizon at an elevation of $\alpha =0\mathrm{deg}$. Each curve shows a trace from horizon to the zenith (elevation 90 deg along the 90 deg scattering plane for maximum DoP). The different colors correspond to different splitting angles $\delta =4\mathrm{ATAN}(A)$ of 3 deg up to 27 deg in steps of 3 deg. As the double-scattered term grows stronger and the splitting angle increases, the AoP variations between single Rayleigh scattering and this multiple-scattering model increases from 0.005 deg up to and approaching 0.4 deg. However, at the 90-deg scattering location, the angular differences between single- and double-scattering models drop significantly.

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–111 There have been many implementations of achromatic and polychromatic designs in both stellar and solar communities.^{40,41,112–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

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.1.

HiVIS Data Extraction

As part of routine calibration, modulation matrix elements were derived using our polarization calibration unit.³¹ 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.³¹ 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.

Fig. 7

The $quv$ to $quv$ modulation matrix elements for the standard Stokes definition sequence liquid crystal voltages derived using the full-Stokes injection unit (a polarization state generator) in front of the HiVIS spectrograph slit. The unit is a wire grid polarizer and quarter-wave retarder on independently controlled rotation stages for creating known $quv$ inputs. The matrices shown here use only one (+) of the two polarization calibration unit input Stokes parameter sets. Wavelengths span the $\sim 6300\AA$ to 8800 Å range. For clarity, only the median value for each spectral order is shown (4000 spectral pixels per order). The liquid crystals were roughly tuned for a standard Stokes definition modulation set around 7000 Å. The diagonal elements are roughly 0.9 at these wavelengths. The nondiagonal elements are all nonzero, and some have amplitudes above 0.7 within the observed wavelength range.

Fig. 8

The computed angular separation ($\gamma$) between the telescope pointing and the sun for all 4890 exposures (1630 full-Stokes polarization measurements) shown in (a). The measured daytime sky DoP for each exposure set is shown in (b). Clouds, pointing, atmospheric properties, and time are all variables. Note that these measured DoP values are used to scale each measured Stokes vector to 100% DoP for use in our calibration algorithm.

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.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.

Fig. 9

The color coding shows the number of daytime sky polarization observations at each telescope pointing available for the telescope Mueller matrix calculations. The filtering by demodulated DoP ($\mathrm{DoP}>15\%$) is shown, but no iterative filter (convergence) has been applied. Interpolation (linear) between neighbors on the az–el grid has been applied for clarity. The small black triangles show the altitude–azimuth track for a star, $\epsilon$ Aurigae that we observed in 2015 to illustrate a typical altitude–azimuth combination required for calibration.

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 }_{\max }$) from the HiVIS measured DoP ($\delta$) and the scattering angle ($\gamma$): ${\delta }_{\max }=\delta \frac{1+{\cos }^2\gamma }{{\sin }^2\gamma }$.

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

Fig. 10

(a) The measured sky DoP for HiVIS daytime sky observations as a function of angular distance from the sun ($\gamma$). Only points passing a 10% DoP filter and a 30% ${\delta }_{\max }$ data filter are included. Colored curves show the Rayleigh sky polarization $\delta$ as function of scattering angle for ${\delta }_{\max }$ amplitudes of 40%, 60%, and 80%. (b) The estimated maximum atmospheric DoP (${\delta }_{\max }$) computed directly from the demodulated HiVIS exposures. The summer observing had a maximum solar elevation of 87 deg while the winter observing season had a maximum solar elevation of around 55 deg. As we observed in the afternoons more often (time allocation constraints) and had more time with the sun well above the horizon, there are less data points at angular distances larger than 90 deg.

By taking this simple relation, a set of data filters can be created. Data points with low ${\delta }_{\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 }_{\max }$ values computed from the HiVS measurements. There are a large number of points showing a high maximum degree of sky polarization (${\delta }_{\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 }_{\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 }_{\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.

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%.

Fig. 11

The six Mueller matrix element estimates computed at all azimuths and elevations with sufficient numbers of valid measurements. The grid of points corresponds to the (azimuth and elevation) sampling used for this study. Only data with valid observations are shown, giving rise to the grid pattern. Each of the four panels shows estimates for $QQ$ and $UQ$ on top, $QU$ and $UU$ in the middle, $QV$ and $UV$ on the bottom. (a) All data from October to May with a minimum 15% DoP measured as well as a minimum of five observations per telescope pointing requirement. (b) The October-only set with a 15% minimum DoP filter and a minimum of three observations per pointing. (c) The December data with a 20% minimum DoP filter and a minimum of three observations per pointing. (d) The May-only set with a 30% minimum DoP filter and a minimum of five observations per pointing. A close inspection of some points near elevations of 90 deg shows some disagreement between the October and both the December and May sets. In all the single-month October, December, and May data sets, some telescope pointings at elevations of $\sim 40\mathrm{deg}$ and azimuths $\sim 90\mathrm{deg}$ are not computed due to a lack of sufficient number of data sets. A calibration campaign must balance the requirements for efficiency, speed, sufficient pointing coverage, and having a large enough data set to reliably identify and reject contaminated outliers.

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.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.

Fig. 14

The best fit Mueller matrix derived for four different spectral orders as functions of azimuth and elevation. All panels are on a linear gray-scale of $\pm 1$. The final rotation matrices shown here are fits to the trigonometric function-based Mueller matrix maps on a fine altitude–azimuth grid. The 15% minimum measured DoP filter and iterative consistency-filters have been applied. The four spectral orders are [0,5,10,15], corresponding to wavelengths of 6260 Å, 6880 Å, 7650 Å, and 8600 Å, respectively. For the longest wavelength in the lower right corner, the $VV$ term is essentially $\pm 1$ with minimal dependence on elevation. The $QU$ terms show the expected geometrical projection from an altitude–azimuth-based reference frame to the fixed slit-based reference frame. The $VQ$ and $VU$ terms on the right side of each panel show nearly $\pm 1$ amplitudes at the shorter wavelengths $<7000\AA$ but are nearly 0 at 8600 Å.

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\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 \cos (\delta )$. Stokes $V$ is similarly computed with both $XY$ field amplitudes and phases: $2.\times Ex\times Ey\times \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.

Fig. 15

The top panel retardance in degrees phase for a 45-deg fold mirror coated with various flavors of enhanced protected silver formulas. The solid black line shows a typical enhanced protected silver specification. The other lines show the retardation caused by various coating formulations. Common materials and coatings include zinc sulfide (ZnS), sapphire (${\mathrm{Al}}_2{\mathrm{O}}_3$), fused silica (${\mathrm{SiO}}_2$), and over silver (Ag). The retardation of multilayer coatings often crosses the nominal 180 deg phase at two wavelengths in the visible. The red curve shows a single fused silica coating, and this curve never reaches 180 deg.

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.

Fig. 16

The Zemax simulated Mueller matrix at 626 nm for 360-deg azimuth variation and a fully articulated 0 deg to 180 deg elevation range using a nominal enhanced protected silver coating on all surfaces except the primary mirror (aluminum + aluminum oxide). The full elevation range was modeled to demonstrate the elevation dependence of some Mueller matrix elements. The mirror geometry provides two degenerate optical configurations when pointing to a particular point on the sky with different Mueller matrix calculations (under the substitution azimuth–azimuth $\pm 180$, elevation–180–elevation). The intensity to $quv$ terms have been linearly scaled to $\pm 1.5\%$. The $quv$ to $quv$ terms have been scaled to $\pm 1$. The II term has not been displayed because it is set to 1 in this model calculation.

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.

Fig. 17

Comparison of coatings for all Mueller matrix elements at 600 nm. Each panel shows either azimuth or elevation dependence at particular pointings as appropriate for the element. Telescope pointing locations were chosen to demonstrate how each coating formula changes the amplitude and relative zero-crossing for the different Mueller matrix elements. Solid lines are one limiting case in either azimuth or elevation with dashed being the opposite extreme to demonstrate maximal polarization effects (crossed mirrors and maximal differences). The polarization to intensity terms are dominated by elevation dependence. The $qu$ to $qu$ terms are strongly azimuth dependent. Red shows the fused-silica over silver single layer formula $\mathrm{Ag}\text{-}{\mathrm{SiO}}_2$. Black is ${\mathrm{MgF}}_2$ over silver. Dark blue is fused silica over sapphire over silver. Light blue is zinc sulfide over sapphire over silver. As we do not know the actual coating formula for the various AEOS mirrors, these functions represent the range of possibilities for an optical system using mirrors with common formulas.

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.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.³⁰ 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.

Fig. 18

The average intensity and $v$ polarization spectra for May 16th and 17th. (a) The continuum normalized intensity spectrum filtered by a 60-pixel boxcar smooth to remove low amplitude ripples. (b) Stokes $v$. select spectral orders are shown with a wide range in number of spectral lines and associated line depths. CCD edges were trimmed (from 2048 pixels down to 2000 pixels per CCD). The gap in the middle of each spectrum represents the CCD mosaic gap in addition to the 24 pixels trimmed from the edge of each device. Clear $v$ line polarization spectra are seen in all lines for all displayed wavelengths with amplitudes of 0.1% to 0.5%.

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.

Fig. 19

The modulation matrix element for the $UU$ term for spectral order 3 at 6620 Å. The chromatic variation across the 4000 spectral pixels accounts for variation of 0.2, which is comparable to other error sources presented in this paper.

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.¹⁰¹ 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³¹ 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.

Fig. 20

The angular residual variation between calibrated HiVIS measurements and the computed Rayleigh sky model on the Poincaré sphere. Calibrations and observations were taken at a telescope azimuth of 180 and elevation of 50. The minimum measured DoP filter of 30% was applied.

Angular variations of $\sim 6\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.