24 January 2017 Daytime sky polarization calibration limitations
Author Affiliations +
J. of Astronomical Telescopes, Instruments, and Systems, 3(1), 018001 (2017). doi:10.1117/1.JATIS.3.1.018001
The daytime sky has recently been demonstrated as a useful calibration tool for deriving polarization cross-talk properties of large astronomical telescopes. The Daniel K. Inouye Solar Telescope and other large telescopes under construction can benefit from precise polarimetric calibration of large mirrors. Several atmospheric phenomena and instrumental errors potentially limit the technique’s accuracy. At the 3.67-m AEOS telescope on Haleakala, we performed a large observing campaign with the HiVIS spectropolarimeter to identify limitations and develop algorithms for extracting consistent calibrations. Effective sampling of the telescope optical configurations and filtering of data for several derived parameters provide robustness to the derived Mueller matrix calibrations. Second-order scattering models of the sky show that this method is relatively insensitive to multiple-scattering in the sky, provided calibration observations are done in regions of high polarization degree. The technique is also insensitive to assumptions about telescope-induced polarization, provided the mirror coatings are highly reflective. Zemax-derived polarization models show agreement between the functional dependence of polarization predictions and the corresponding on-sky calibrations.
Harrington, Kuhn, and Ariste: Daytime sky polarization calibration limitations



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%12.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. 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. 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,2021.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.2324.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 DST4,2829.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.3233.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,3536.37 Complex modulation and calibration strategies are required for such a multi-instrument system.35,36,3839.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.4243.44.45



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: Si=[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: DoP=Q2+U2+V2I=q2+u2+v2. 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×4 set of transfer coefficients that describes how an optic changes the input Stokes vector (Siinput) to the output Stokes vector (Sioutput): Sioutput=MijSiinput. 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




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.5051. Anisotropic scattered sunlight from reflections off land or water can be highly polarized and temporally variable.6768. Aerosol particle optical properties and vertical distributions also vary.58,7374. The polarization can change across atmospheric absorption bands or can be influenced by other scattering mechanisms.8384.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.8889. Conditions at twilight with low solar elevations can present some spectral differences.7980.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.


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 (δmax) scales the polarization pattern across the sky.

The all-sky model requires knowing the solar location and the scale factor (δ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 (δmax) in this model occurs at a scattering angle (γ) of 90 deg. The Rayleigh sky model predicts the DoP (δ) at any telescope pointing (azimuth and elevation) as


and the spherical geometry is computed as cos(γ)=sin(θs)sin(θ)cos(ψ)+cos(θs)cos(θ), where γ is the angular distance between the telescope pointing and the sun, θs is the solar zenith angle, θ is the angular distance between the telescope pointing and the zenith, and ψ is the azimuthal angle between the solar direction and the telescope pointing. The geometry comes from the law of cosines with γ, θs, and θ as the angular distances and ψ representing the interior angle. The spherical triangle formed by the solar location, zenith, and telescope pointing can be seen in Fig. 1. A detailed example of this single-scattering model can be seen in Fig. 2, which is computed on January 27th 2010. This figure shows several model parameters either in altitude–azimuth projections or in orthographic projects. The single-scattering model has the highest DoP in a band of 90-deg scattering angle.

Fig. 1

The celestial triangle representing the geometry for the sky polarization computations at any telescope pointing. γ is the angular distance between the telescope pointing and the sun. θs is the solar angular distance from the zenith. θ is the angular distance between the telescope pointing and the zenith. φ is the angle between the zenith direction and the solar direction at the telescope pointing. The angle ψ 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 ψ, 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% (δ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 +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 +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.



Solving for Telescope Mueller Matrix Elements

We model the 3×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×3 cross-talk elements as representing the telescope Mueller matrix. We denote the three Euler angles as (α,β,γ) and use a short-hand notation where cos(γ) is shortened to cγ. We specify the rotation matrix (Rij) using the ZXZ convention for Euler angles as


With this definition for the rotation matrix, we solve the Euler angles assuming a linearly polarized daytime sky scaled to 100% DoP as calibration input. If we denote the measured Stokes parameters, Si, as (qm,um,vm) with i=1,2,3 and the input sky Stokes parameters, Rj, as (qr,ur,0), then the 3×3 quv Mueller matrix elements at each wavelength are


We have no V input from the daytime sky to constrain the VQ, VU, and VV terms. Nevertheless, two measurements at different input polarization angles are sufficient to fully specify the rotation matrix. Thus, we use the fact that the sky polarization changes orientation with time and take measurements at identical telescope pointings separated by enough time for the solar sky illumination to change. A set of observations with a changing input AoP yields an overconstrained solvable problem for all six linear polarization terms in the Mueller matrix.

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 (α, β, γ). This system of equations can be solved using a normal nonlinear least-squares minimization by searching the (α, β, γ) 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 (α, β, γ) space for minima in squared error. With the measured Stokes vector (Si), i=(1,2,3), the Rayleigh sky input vector (Rj), j=(1,2), and a rotation matrix (Rij), we define the error (ε) as ε2(α,β,γ)=i=13j=12[SiRiRij(α,β,γ)]2. For n measurements, this gives us 3×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.102103.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 (Rij) for i independent observations and j input Stokes parameters. The measured Stokes parameters (Si) become individual column vectors. The unknown Mueller matrix elements are arranged as a column vector by output Stokes parameter (Mj). If we write measured Stokes parameters as (qmi,umi,vmi) and the Rayleigh input Stokes parameters as (qri,uri), we can explicitly write a set of equations for two Mueller matrix elements


We have three such equations for each set of Mueller matrix elements sampled by sky measurements. We can express the residual error (εi) for each incident Stokes parameter (Si) with an implied sum over j as εi=SiRijMj. The normal solution of an overspecified system of equations is easily derived in a least-squares sense using matrix notation. With the total error E as the sum of all residuals for m independent observations, we get E=i=1mεi2. We solve the least-squares system for the unknown Mueller matrix element (Mj) by minimizing the error with respect to each equation. The partial derivative for εi with respect to Mj is just the sky input elements Rij. Taking the partial with respect to each input Stokes parameter, we get


We have inserted a dummy sum over the index k. Multiplying out the terms and rearranging gives us the normal equations: ikRijRikMk=iRijSi. This is written in matrix notation, which is the familiar solution of a system of equations via the normal method M=RTSRTR.

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 A with an implied sum over i observations for each term. As an example for a single element, if we compute the inverse of A and multiply out A1 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.


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.


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 ζ to denote the location of a point on a stereographic projection of the sky. In Cartesian geometry, ζ=x+iy. In polar coordinates, ζ=reiφ. 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 γ(ζ), they represent the stereographic projection for the sky polarization pattern as w(ζ)=|w|e2iγ(ζ). For the single-scattering case, this simple relation behaves as ζ2 and can be scaled to an amplitude of 1 and written in polar coordinates (r,φ) as w(ζ)ζ2=r2e2i(φπ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 ζ=±iA, which corresponds to a Cartesian y value of ±A. To make the singularities at the antisun location, the equation was generalized to w(ζ)(ζ2+A2)(ζ2+1A2).

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 δ=4ATAN(A)=27  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 white=1 and 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 δ=4ATAN(A)=12  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 α=0  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 δ=4ATAN(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.


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.106107. There have been many implementations of achromatic and polychromatic designs in both stellar and solar communities.40,41,112113. In the notation of these studies, the instrument modulates the incoming polarization information into a series of measured intensities (Ii) for i independent observations via the modulation matrix (Oij) for j input Stokes parameters (Sj): Ii=OijSj. 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: Si=RijMj.

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 (O) via the normal least squares formalism: S=OTIOTO. The demodulation matrix is typically defined as Dij=[OTO]1OT.

In our daytime sky technique, the Rayleigh sky input parameters become the modulation matrix (Oij=Rij) 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 σ and there are n total measurements, then the noise on each demodulated parameter (σi) becomes σi2=nσ2j=1nDij2. The efficiency of the observation becomes ei=(nj=1nDij2)12.

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.


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.

Fig. 6

The computed solar azimuth and elevation for all HiVIS daytime sky observations. May data had the sun near the zenith while winter observations (October and December) had the sun low and in the south. The observing allocations were dominated by sunrise to noon times, giving far more observations with the sun east and overhead.



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 (±q, ±u, ±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.

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 6300  Å 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 (γ) 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.


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.


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 (ε 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 (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, ε Aurigae that we observed in 2015 to illustrate a typical altitude–azimuth combination required for calibration.



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 (δmax) from the HiVIS measured DoP (δ) and the scattering angle (γ): δmax=δ1+cos2γsin2γ.

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

Fig. 10

(a) The measured sky DoP for HiVIS daytime sky observations as a function of angular distance from the sun (γ). Only points passing a 10% DoP filter and a 30% δmax data filter are included. Colored curves show the Rayleigh sky polarization δ as function of scattering angle for δmax amplitudes of 40%, 60%, and 80%. (b) The estimated maximum atmospheric DoP (δ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 δ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 δmax values computed from the HiVS measurements. There are a large number of points showing a high maximum degree of sky polarization (δ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 (δ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.


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 (δ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× 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 δmax points for agreement with Rayleigh model (40% typical δmax for a clear day).

  • Reject pointings with low qu input diversity (>20  deg for each pointing).

  • Reject observations where the calibrated Mueller matrix calibrations show high error outliers (Scal·Mangle<0.1).

  • Reject iteratively by calibrated residual angle between measurements and theory (Scal·R) 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%.

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 40  deg and azimuths 90  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.


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 sin(2Az+2El), sin(Az+El), sin(2El), sin(Az), 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 sin=S, cos=C and subscripts for A=azimuth, E=elevation. The best function we found to fit for each Mueller matrix estimate contains these six terms: S2AC2E,S2EC2A,C2AC2E,S2AS2E,S2E,C2E.

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.

Fig. 12

(a) The six-term trigonometric function coefficients fit to the rotation matrices derived from all six Mueller matrix estimate terms for spectral orders 5 at a wavelength of 6881 Å. Each color corresponds to a different Mueller matrix element. For instance, black shows the QQ term, which is dominated by coefficients 2 and 3 giving a functional form of C2AC2E+S2AS2E. Demodulation was done using the + inputs, and a minimum 15% measured DoP filter has been applied. The iterative consistency and Rayleigh minimum DoP maximum filters have also been applied. (b) The cumulative distribution function (CDF) for the errors between trig-based Mueller matrix estimates and the corresponding rotation matrix fits for spectral order 3. Each color shows one of the six Mueller matrix estimate residual CDFs. The 15% DoP and consistency-filters have been applied.


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.

Fig. 13

The difference between the trigonometric function fit rotation matrix elements and each of the corresponding Mueller matrix element estimates derived directly from the data. The gray-scale has been highly stretched to highlight differences at ±0.1 amplitudes. The differences were computed only where valid observations are recorded after all the filters are applied. Linear interpolation was performed to all other pointings to make smooth maps. The 15% minimum measured DoP filter and iterative consistency-filters have been applied. These represent the disagreement between the empirical Mueller matrix estimates and the trigonometric functions.



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.


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 ±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 ±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 ±1 amplitudes at the shorter wavelengths <7000  Å 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 ±1 amplitudes with strong functional dependence on azimuth and elevation at the shorter wavelengths <7000  Å, 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.


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/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×Ex+Ey×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×ExEy×Ey). Stokes U is computed from X and Y electric field amplitudes as well as phase variations: 2.×Ex×Ey×cos(δ). Stokes V is similarly computed with both XY field amplitudes and phases: 2.×Ex×Ey×sin(δ). The term δ 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 (Al2O3), fused silica (SiO2), 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 ±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 ±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 ±180, elevation–180–elevation). The intensity to quv terms have been linearly scaled to ±1.5%. The quv to quv terms have been scaled to ±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 Ag-SiO2. Black is MgF2 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.



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


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.


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.

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



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



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.


Limitation Summary

In this section, we presented several other instrument calibration limitations specific to HiVIS. The angular error histograms presented31 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 6  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 0.1 in a Stokes parameter. Liquid crystal temperature instability, reduced spectral sampling and optical misalignment could all contribute to errors in calibrated data.



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 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 (δ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  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 (Scal·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.


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 ζ to denote the complex location of a point on a stereographic projection of the sky. In Cartesian geometry, ζ  =x+iy. In polar coordinates, ζ=reiφ. 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 γ(ζ), they represent the stereographic projection for the sky polarization pattern as w(ζ)=|w|e2iγ(ζ). For the single-scattering case, this simple relation behaves as ζ2 can be scaled to an amplitude of 1 and written in polar coordinates (r,φ) as w(ζ)ζ2=r2e2i(φπ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 ζ=±iA, which corresponds to a Cartesian y value of ±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 α and use the stereographic projection where the sun is on the y axis at the location ζ=iys=i(1tanα2)(1+tanα2), the four zero polarization singularities are located



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


The denominator was chosen with the complex modulus terms to ensure the amplitude |w| is always 1. With this simple equation, you can relate the constant A to the angular separation of the polarization zero points as δ=4 ATAN(A).

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.


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 20  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× higher at these telescope pointings.


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.

Fig. 21

The measured data over the entire campaign for a telescope pointing of azimuth 60 deg and elevation 25 deg. Stokes quv parameters have been normalized by the measured DoP to create 100% polarized quv points. The measured DoP used to do the normalization is shown in the lower right panel. All points are shown as functions of the input polarization angle. May data are shown in black. October data are in orange. November data are in red. December data are in blue. Over the October to May timeframe, the sun moves across the sky and generates a wide range of input linear polarization angles (plotted as the x axis). There is a direct functional relationship between the input polarization angle and the output quv data. The smooth functional relationship (with some outliers to be rejected) illustrates the mapping between input qu angles and output quv angles on the Poincaré sphere. The demodulation and 100% DoP scaling has been applied. The <10% measured DoP filter has been applied to reject low DoP points.


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.


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 ξi(α,β,γ).

  • 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 (S·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: Qi=ACOS(Ri·Si).

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.


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.

Fig. 22

The measured sky quv spectra as functions of input AoP for different coudé windows (BK7 and Infrasil). Blue symbols show the Infrasil window; red symbols are for the BK7 window substrate. The telescope was at a pointing of azimuth 330-deg elevation 20 deg. Similar results are seen at other telescope pointings. The diamonds show measured Stokes q, the triangles show measured Stokes u, and the asterisks show measured Stokes v. All points have been scaled to 100% DoP and have been demodulated using the + state. Observations have been minimally filtered using a minimum 15% measured DoP threshold and a 30% minimum δmax filter to guarantee a reasonably polarized sky signal (rejecting clouds, obvious outliers, or highly contaminated data sets). The measured quv points are consistent between different window substrates, different times, and different input polarization angles. We consider the two windows to not limit the calibration precision for this study.


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×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×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 I2Q2U2V2>0. This set of rules means Mueller matrices are part of a group, formally called SO(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 (α,β,γ) 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 (di) for i=1,2,3: M=exp(αS1+βS2+γS3)=exp(ΣdiSi). 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 S1 is often denoted J3 and is


If we write a matrix using this notation for a rotation about the Z axis by and Euler angle α, we recover a rotation matrix, which is also the standard form for a circular retarder with a retardance α via exp(αJ3)



The standard infinitesimal generator for an X rotation is denoted J1 and if we rotate by an angle β, we get a Mueller matrix for a linear retarder via exp(βJ1)


We can use the same exponential notation for generating IQ and QI terms of the Mueller matrix using the generator K1 and a small diattenuation term ε. The matrix would be computed as exp(εK1) For this Mueller matrix, we take the approximation that sinh(ε)ε and cosh(ε)1, so the Mueller matrix simplifies.



The arbitrary form for the generators can be represented by the terms applicable to the first row and column of the Mueller matrix, Ki and the terms corresponding to rotations, Ji. The full set of generators can be written with a set of infinitessimal operators (often called boosts) in quv denoted as ζi and a set of rotations in quv with rotations denoted as θi


We can show how errors in the first row and column of the Mueller matrix propagate into cross-talk elements for quv to quv. Consider a Mueller matrix of the form Mij=exp(εK1+αS1+βS2+γS3)=exp(εK1+ΣdiSi). Using the infinitessimal generators outlined above, and the rotation matrix Rij as defined in Eq. (3) (in main text) we get a Mueller matrix:


Each element Rij is part of the rotation matrix, and ε is a small error caused by dichroism in the system. We will neglect the other dichroism terms (IU and IV) for simplicity. In the limit of small ε, we can derive a sensitivity of the Rij elements to ε. Note that for AEOS, the Rij terms are as large as 1 while the induced polarization and depolarization terms ε are <0.05.

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: exp(X+Y)=exp(X)exp(Y)exp(12[X,Y]). Additional terms grow complex quickly but have small amplitudes. The next correction to this equation involves nested commuting: 112[X,[X,Y]]+[Y,[Y,X]]. For our case, the noncommuting terms would be [X,Y]=[K1,Si]. We write our Mueller matrix where we denote the sum of noncommuting rotations exp(ΣdiSi) with Si as the ZXZ rotation group as a single matrix Rij of the three Euler angles di=(α,β,γ) for i=1,2,3.


This gives us three terms multiplied to create the Mueller matrix. Since ε is small, we approximate the first term exp(εK1) as the identity matrix (II) plus an additional term representing the diattenuation (εK1) to first order.


For the third term representing [X,Y], we need to use the commutation relation [Xi,Yj]=εijkXk where εijk is the 3-D Levi–Civita symbol. The commutation relation for the three Euler angles (i=1,2,3) gives [K1,Si]εijkKk. The Levi–Civita symbol will be +1 for the cyclic (1, 2, 3) ordering, adding the term βK3. The Levi–Civita symbol will be 1 for the anticyclic (1, 3, 2) ordering, adding the term γK2. If we collect terms and include the proper Euler angles in the sum over generators ΣdiSi, we find that the term [K1,ΣdiSi] includes only (βK3γK2). Therefore, the third term in Eq. (18) for the Mueller matrix becomes


Using a similar approximation for small ε we can write the [X,Y] term as



Combining the three terms for Eq. (18), we get


Multiplying out this equation and neglecting terms that are ε2, we collect four terms for the Mueller matrix with the original IQ perturbation ε, the Euler angles β, γ and the generators K1, K2, and K3


This expression gives the approximation for the Mueller matrix under our assumptions for small errors ε in the polarization to and from intensity cross-talk level terms of the Mueller matrix. The Mueller matrix can be represented as a rotation matrix (R) multiplied by a group of four correction terms. For the case we examined here, a small error in the QI and IQ terms gives us a first-order correction to the Mueller matrix. When expanding out the terms in Eq. (23), we get first-order corrections that have no impact inside the cross-talk terms R from errors in ε.

We can write the first Mueller matrix correction term K1R as



This term is scaled by ε 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 K2R


A similar Mueller matrix correction would be generated for the K3 term. From these corrections to the Mueller matrices in Eqs. (24) and (25), we see that there is no correction in the rotation matrix terms R due to this first-order approximation. In the limit of small ε, we see that neglecting the IQ and QI terms only impacts the first row of the Mueller matrix. From this sensitivity analysis, we can conclude that our method of approximating the 3×3 cross-talk elements of the Mueller matrix as a rotation matrix is relatively insensitive to errors in the first row and column of the Mueller matrix. Corrections to M to first-order in ε as above only affect the intensity to polarization terms (the dichroism in M). The rotation matrix fit terms are second-order in ε and, as such, have errors much smaller than other typical limiting noise sources. For AEOS and HiVIS, the cross-talk terms are of order 1. The induced polarization and depolarization terms are of order 5%. The second-order corrections from neglecting the first row and column in the rotation matrix fitting are thus of order ε0.052, which is 2.5×103. The method of fitting rotation matrices to the cross-talk elements of the Mueller matrix is robust against dichroism type errors in induced polarization and depolarization.


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.



J.-C. Hsu and M. Breger, “On standard polarized stars,” Astrophys. J. 262, 732–738 (1982).http://dx.doi.org/10.1086/160467Google Scholar


R. Gil-Hutton and P. Benavidez, “Southern stars that can be used as unpolarized standards,” Mon. Not. R. Astron. Soc. 345, 97–99 (2003).MNRAA40035-8711http://dx.doi.org/10.1046/j.1365-8711.2003.06957.xGoogle Scholar


L. Fossati et al., “Standard stars for linear polarization observed with FORS1,” in The Future of Photometric, Spectrophotometric and Polarimetric Standardization, Vol. 364, p. 503 (2007).Google Scholar


H. Socas-Navarro et al., “Characterization of telescope polarization properties across the visible and near-infrared spectrum. Case study: the Dunn solar telescope,” Astron. Astrophys. 531, A2 (2011).AAEJAF0004-6361http://dx.doi.org/10.1051/0004-6361/201015804Google Scholar


H. Socas-Navarro et al., “Spinor: visible and infrared spectro-polarimetry at the national solar observatory,” Sol. Phys. 235, 55–73 (2006).http://dx.doi.org/10.1007/s11207-006-0020-xGoogle Scholar


H. Socas-Navarro et al., “High precision polarimetry with the advanced technology solar telescope,” Proc. SPIE 5901, 590105 (2005).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.616000Google Scholar


H. Socas-Navarro, “Polarimetric calibration of large-aperture telescopes. II. Subaperture method,” Opt. Soc. Am. J. 22, 907 (2005).http://dx.doi.org/10.1364/JOSAA.22.000907Google Scholar


H. Socas-Navarro, “Polarimetric calibration of large-aperture telescopes. I. Beam-expansion method,” Opt. Soc. Am. J. 22, 539 (2005).http://dx.doi.org/10.1364/JOSAA.22.000539Google Scholar


J. R. Kuhn et al., “Removing instrumental polarization from infrared solar polarimetric observations,” Sol. Phys. 153, 143–155 (1994).http://dx.doi.org/10.1007/BF00712497Google Scholar


D. F. Elmore et al., “Utilization of redundant polarized solar spectra to infer the polarization properties of the new generation of large aperture solar telescopes,” Proc. SPIE 7735, 77354E (2010).http://dx.doi.org/10.1117/12.857061Google Scholar


J. Sánchez Almeida, V. Martinez Pillet and A. D. Wittmann, “The instrumental polarization of a Gregory-Coude telescope,” Sol. Phys. 134, 1–13 (1991).http://dx.doi.org/10.1007/BF00148738Google Scholar


E. Giro et al., “Polarization properties at the Nasmyth focus of the alt-azimuth TNG telescope,” Proc. SPIE 4843, 456 (2003).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.458607Google Scholar


F. Patat and M. Romaniello, “Error analysis for dual-beam optical linear polarimetry,” Publ. Astron. Soc. Pac. 118, 146–161 (2006).PASPAU0004-6280http://dx.doi.org/10.1086/pasp.2006.118.issue-839Google Scholar


J. Tinbergen, “Accurate optical polarimetry on the Nasmyth platform,” Publ. Astron. Soc. Pac. 119, 1371–1384 (2007).PASPAU0004-6280http://dx.doi.org/10.1086/509224Google Scholar


F. Joos et al., “Reduction of polarimetric data using Mueller calculus applied to Nasmyth instruments,” Proc. SPIE 7016, 70161I (2008).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.788915Google Scholar


G. van Harten, F. Snik and C. U. Keller, “Polarization properties of real aluminum mirrors, I. Influence of the aluminum oxide layer,” Publ. Astron. Soc. Pac. 121, 377–383 (2009).PASPAU0004-6280http://dx.doi.org/10.1086/597155Google Scholar


R. Roelfsema et al., “The ZIMPOL high-contrast imaging polarimeter for SPHERE: design, manufacturing, and testing,” Proc. SPIE 7735, 77354B (2010).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.857045Google Scholar


K. Ichimoto et al., “Polarization calibration of the solar optical telescope onboard Hinode,” Sol. Phys. 249, 233–261 (2008).http://dx.doi.org/10.1007/s11207-008-9169-9Google Scholar


B. W. Lites et al., “The Hinode spectro-polarimeter,” Sol. Phys. 283, 579–599 (2013).http://dx.doi.org/10.1007/s11207-012-0206-3Google Scholar


P. G. Judge et al., “Evaluation of seeing-induced cross talk in tip-tilt-corrected solar polarimetry,” Appl. Opt. 43, 3817 (2004).APOPAI0003-6935http://dx.doi.org/10.1364/AO.43.003817Google Scholar


J. Sánchez Almeida, “Instrumental polarization in the focal plane of telescopes. 2: effects induced by seeing,” Astron. Astrophys. 292, 713–721 (1994).AAEJAF0004-6361Google Scholar


J. Sánchez Almeida and V. Martinez Pillet, “Instrumental polarization in the focal plane of telescopes,” Astron. Astrophys. 260, 543–555 (1992).AAEJAF0004-6361Google Scholar


M. Stangalini et al., “The effects of AO systems on polarized light,” Proc. SPIE 9148, 91486P (2014).http://dx.doi.org/10.1117/12.2056470Google Scholar


D. Harrington et al., “InnoPOL: an EMCCD imaging polarimeter and 85-element curvature AO system on the 3.6-m AEOS telescope for cost effective polarimetric speckle suppression,” Proc. SPIE 9147, 91477C (2014).http://dx.doi.org/10.1117/12.2056667Google Scholar


M. van Noort, “Spatially coupled inversion of spectro-polarimetric image data,” Astron. Astrophys. 548, A5–A14 (2012).AAEJAF0004-6361http://dx.doi.org/10.1051/0004-6361/201220220Google Scholar


M. De Juan Ovelar et al., “Modeling the instrumental polarization of the VLT and E-ELT telescopes with the M&m’s code,” Proc. SPIE 8449, 844912 (2012).http://dx.doi.org/10.1117/12.926588Google Scholar


D. M. Harrington, J. R. Kuhn and K. Whitman, “The new HiVIS spectropolarimeter and spectropolarimetric calibration of the AEOS telescope,” Publ. Astron. Soc. Pac. 118, 845–859 (2006).PASPAU0004-6280http://dx.doi.org/10.1086/pasp.2006.118.issue-844Google Scholar


G. Barrick, T. Benedict and D. Sabin, “Correcting polarization crosstalk in the ESPaDOnS spectro-polarimeter,” Proc. SPIE 7735, 77354C (2010).http://dx.doi.org/10.1117/12.856084Google Scholar


G. Barrick and T. Benedict, “Experimental results from using two laminated film polarizers to make absolute measurements of polarization crosstalk in an optic,” Proc. SPIE 7735, 773548 (2010).http://dx.doi.org/10.1117/12.856085Google Scholar


D. M. Harrington et al., “Correcting systematic polarization effects in Keck LRISp spectropolarimetry to 0.05 percent,” Publ. Astron. Soc. Pac. 127, 757–775 (2015).PASPAU0004-6280http://dx.doi.org/10.1086/682323Google Scholar


D. Harrington, J. R. Kuhn and R. Nevin, “Calibrating and stabilizing spectropolarimeters with charge shuffling and daytime sky measurements,” Astron. Astrophys. 578, A126 (2015).AAEJAF0004-6361http://dx.doi.org/10.1051/0004-6361/201322791Google Scholar


J. P. McMullin et al., “Construction status of the Daniel K. Inouye solar telescope,” Proc. SPIE 9145, 914525 (2014).http://dx.doi.org/10.1117/12.2055483Google Scholar


S. L. Keil et al., “ATST: the largest polarimeter,” in Solar Polarization 6. Proc. of a Conf. Held in Maui, Vol. 437, p. 319 (2011).Google Scholar


T. R. Rimmele et al., “Instrumentation for the advanced technology solar telescope,” Proc. SPIE 5492, 944 (2004).http://dx.doi.org/10.1117/12.551853Google Scholar


D. F. Elmore, S. R. Sueoka and R. Casini, “Performance of polarization modulation and calibration optics for the Daniel K. Inouye solar telescope,” Proc. SPIE 9147, 91470F (2014).http://dx.doi.org/10.1117/12.2054610Google Scholar


D. F. Elmore et al., “The Daniel K. Inouye solar telescope first light instruments and critical science plan,” Proc. SPIE 9147, 914707 (2014).http://dx.doi.org/10.1117/12.2057038Google Scholar


W. Schmidt et al., “A two-dimensional spectropolarimeter as a first-light instrument for the Daniel K. Inouye solar telescope,” Proc. SPIE 9147, 91470E (2014).http://dx.doi.org/10.1117/12.2056322Google Scholar


S. R. Sueoka, R. A. Chipman and D. F. Elmore, “Characterization of DKIST retarder components with polarization ray tracing,” Proc. SPIE 9293, 929308 (2014).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.2071279Google Scholar


W. H. Schubert, E. Petrak and T. G. Baur, “Measurement of polarization assemblies for the Daniel K. Inouye solar telescope,” Proc. SPIE 9369, 93690N (2015).http://dx.doi.org/10.1117/12.2077749Google Scholar


A. G. de Wijn, “Preliminary design of the visible spectro-polarimeter for the advanced technology solar telescope,” Proc. SPIE 8446, 84466X (2012).http://dx.doi.org/10.1117/12.926497Google Scholar


A. G. de Wijn et al., “The polychromatic polarization modulator,” Proc. SPIE 7735, 77354A (2010).http://dx.doi.org/10.1117/12.857745Google Scholar


J. Sánchez-Capuchino et al., “Current concept for the 4m European Solar Telescope (EST) optical design,” Proc. SPIE 7652, 76520S (2010).http://dx.doi.org/10.1117/12.871604Google Scholar


F. C. M. Bettonvil et al., “The polarization optics for the European Solar Telescope,” in Solar Polarization 6. Proc. of a Conf. Held in Maui, Vol. 437, p. 329 (2011).Google Scholar


F. C. M. Bettonvil et al., “The polarization optics for the European Solar Telescope (EST),” Proc. SPIE 7735, 77356I (2010).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.857817Google Scholar


M. Collados et al., “European Solar Telescope: project status,” Proc. SPIE 7733, 77330H (2010).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.856994Google Scholar


E. Collett, Polarized Light. Fundamentals and Applications, Vol. 1, CRC Press (1992).Google Scholar


D. Clarke, Stellar Polarimetry, CLARKE:STELLAR POLARIS O-BK, John Wiley & Sons, Weinheim, Germany (2009).Google Scholar


K. L. Coulson, “Characteristics of skylight at the zenith during twilight as indicators of atmospheric turbidity. 1: degree of polarization,” Appl. Opt. 19, 3469 (1980).APOPAI0003-6935http://dx.doi.org/10.1364/AO.19.003469Google Scholar


K. L. Coulson, Polarization and Intensity of Light in the Atmosphere, A Deepak Publication (1988).Google Scholar


G. Horvath et al., “First observation of the fourth neutral polarization point in the atmosphere,” Opt. Soc. Am. J. 19, 2085 (2002).http://dx.doi.org/10.1364/JOSAA.19.002085Google Scholar


G. Horváth, I. Pomozi and J. Gál, “Neutral points of skylight polarization observed during the total eclipse on 11 August 1999,” Appl. Opt. 42, 465 (2003).APOPAI0003-6935http://dx.doi.org/10.1364/AO.42.000465Google Scholar


H. Horvath, “Basic optics, aerosol optics, and the role of scattering for sky radiance,” J. Quant. Spectrosc. Radiat. Transfer 139, 3–12 (2014).JQSRAE0022-4073http://dx.doi.org/10.1016/j.jqsrt.2013.08.009Google Scholar


R. L. Lee, “Digital imaging of clear-sky polarization,” Appl. Opt. 37, 1465 (1998).APOPAI0003-6935http://dx.doi.org/10.1364/AO.37.001465Google Scholar


Y. Liu and K. Voss, “Polarized radiance distribution measurement of skylight. II. Experiment and data,” Appl. Opt. 36, 8753 (1997).APOPAI0003-6935http://dx.doi.org/10.1364/AO.36.008753Google Scholar


B. Suhai and G. Horváth, “How well does the Rayleigh model describe the E-vector distribution of skylight in clear and cloudy conditions? A full-sky polarimetric study,” Opt. Soc. Am. J. 21, 1669 (2004).http://dx.doi.org/10.1364/JOSAA.21.001669Google Scholar


J. Gál et al., “Polarization of the moonlit clear night sky measured by full-sky imaging polarimetry at full Moon: comparison of the polarization of moonlit and sunlit skies,” J. Geophys. Res. 106, 22647–22653 (2001).http://dx.doi.org/10.1029/2000JD000085Google Scholar


J. Gál et al., “Polarization patterns of the summer sky and its neutral points measured by full-sky imaging polarimetry in Finnish Lapland north of the Arctic Circle,” in Royal Society of London Proc. Series A, p. 1385, The Royal Society (2001).Google Scholar


A. Vermeulen, C. Devaux and M. Herman, “Retrieval of the scattering and microphysical properties of aerosols from ground-based optical measurements including polarization. I. Method,” Appl. Opt. 39, 6207 (2000).APOPAI0003-6935http://dx.doi.org/10.1364/AO.39.006207Google Scholar


L. Li et al., “A method to calculate Stokes parameters and angle of polarization of skylight from polarized CIMEL sun/sky radiometers,” J. Quant. Spectrosc. Radiat. Transfer 149, 334–346 (2014).JQSRAE0022-4073http://dx.doi.org/10.1016/j.jqsrt.2014.09.003Google Scholar


Z. Li et al., “Improvements for ground-based remote sensing of atmospheric aerosol properties by additional polarimetric measurements,” J. Quant. Spectrosc. Radiat. Transfer 110, 1954–1961 (2009).JQSRAE0022-4073http://dx.doi.org/10.1016/j.jqsrt.2009.04.009Google Scholar


S. Liang and P. Lewis, “A parametric radiative transfer model for sky radiance distribution,” J. Quant. Spectrosc. Radiat. Transfer 55(2), 181–189 (1996).JQSRAE0022-4073http://dx.doi.org/10.1016/0022-4073(95)00155-7Google Scholar


I. Pomozi, G. Horváth and R. Wehner, “How the clear-sky angle of polarization pattern continues underneath clouds: full-sky measurements and implications for animal orientation,” J. Exp. Biol. 204, 2933–2942 (2001).Google Scholar


T. W. Cronin, E. J. Warrant and B. Greiner, “Celestial polarization patterns during twilight,” Appl. Opt. 45, 5582 (2006).APOPAI0003-6935http://dx.doi.org/10.1364/AO.45.005582Google Scholar


T. W. Cronin, E. J. Warrant and B. Greiner, “Polarization patterns of the twilight sky,” Proc. SPIE 5888, 58880R (2005).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.613053Google Scholar


R. Hegedüs, S. Åkesson and G. Horváth, “Anomalous celestial polarization caused by forest fire smoke: why do some insects become visually disoriented under smoky skies?” Appl. Opt. 46, 2717 (2007).APOPAI0003-6935http://dx.doi.org/10.1364/AO.46.002717Google Scholar


M. V. Berry, M. R. Dennis and R. L. J. Lee, “Polarization singularities in the clear sky,” New J. Phys. 6, 162–162 (2004).http://dx.doi.org/10.1088/1367-2630/6/1/162Google Scholar


P. Litvinov et al., “Reflection models for soil and vegetation surfaces from multiple-viewing angle photopolarimetric measurements,” J. Quant. Spectrosc. Radiat. Transfer 111(4), 529–539 (2010).JQSRAE0022-4073http://dx.doi.org/10.1016/j.jqsrt.2009.11.001Google Scholar


J. Peltoniemi et al., “Polarised bidirectional reflectance factor measurements from soil, stones, and snow,” J. Quant. Spectrosc. Radiat. Transfer 110(17), 1940–1953 (2009).JQSRAE0022-4073http://dx.doi.org/10.1016/j.jqsrt.2009.04.008Google Scholar


X. He et al., “A vector radiative transfer model of coupled ocean-atmosphere system using matrix-operator method for rough sea-surface,” J. Quant. Spectrosc. Radiat. Transfer 111(10), 1426–1448 (2010).JQSRAE0022-4073http://dx.doi.org/10.1016/j.jqsrt.2010.02.014Google Scholar


S. Salinas and S. Liew, “Light reflection from a rough liquid surface including wind-wave effects in a scattering atmosphere,” J. Quant. Spectrosc. Radiat. Transfer 105(3), 414–424 (2007).JQSRAE0022-4073http://dx.doi.org/10.1016/j.jqsrt.2007.01.051Google Scholar


Y. Ota et al., “Matrix formulations of radiative transfer including the polarization effect in a coupled atmosphere-ocean system,” J. Quant. Spectrosc. Radiat. Transfer 111, 878–894 (2010).JQSRAE0022-4073http://dx.doi.org/10.1016/j.jqsrt.2009.11.021Google Scholar


V. Kisselev and B. Bulgarelli, “Reflection of light from a rough water surface in numerical methods for solving the radiative transfer equation,” J. Quant. Spectrosc. Radiat. Transfer 85, 419–435 (2004).JQSRAE0022-4073http://dx.doi.org/10.1016/S0022-4073(03)00236-XGoogle Scholar


B. Wu and Y. Jin, “Twilight polarization and optical depth of stratospheric aerosols over Beijing after the Pinatubo volcanic eruption,” Appl. Opt. 36, 7009 (1997).APOPAI0003-6935http://dx.doi.org/10.1364/AO.36.007009Google Scholar


G. van Harten et al., “Atmospheric aerosol characterization with a ground-based SPEX spectropolarimetric instrument,” Atmos. Meas. Tech. 7, 4341–4351 (2014).http://dx.doi.org/10.5194/amt-7-4341-2014Google Scholar


A. K. Shukurov and K. A. Shukurov, “Field studies of the correlation between the atmospheric aerosol content and the light polarization at the zenith of the daytime sky,” Izvestiya 42, 68–73 (2006).http://dx.doi.org/10.1134/S0001433806010063Google Scholar


O. S. Ougolnikov and I. A. Maslov, “Wide-angle polarimetry of the night sky: measurements of atmospheric glow and zodiacal light,” Cosmic Res. 43, 17–24 (2005).CSCRA70010-9525http://dx.doi.org/10.1007/s10604-005-0015-7Google Scholar


O. S. Ougolnikov and I. A. Maslov, “Multicolor polarimetry of the twilight sky: the role of multiple light scattering as a function of wavelength,” Cosmic Res. 40, 224–232 (2002).CSCRA70010-9525http://dx.doi.org/10.1023/A:1015920819478Google Scholar


O. Ougolnikov and I. Maslov, “Polarization studies of contribution of aerosol scattering to the glow of twilight sky,” Cosmic Res. 43(6), 404–412 (2005).CSCRA70010-9525http://dx.doi.org/10.1007/s10604-005-0064-yGoogle Scholar


O. S. Ugolnikov, “Twilight sky photometry and polarimetry: the problem of multiple scattering at the twilight time,” Kosm. Issled. 37, 168–175 (1999).KSIUAWGoogle Scholar


O. S. Ugolnikov and I. A. Maslov, “Scattering function of tropospheric aerosol according to the data of polarimetry of the twilight and night sky background,” Cosmic Res. 49, 187–193 (2011).CSCRA70010-9525http://dx.doi.org/10.1134/S0010952511030105Google Scholar


O. S. Ugolnikov and I. A. Maslov, “Optical properties of the undisturbed mesosphere from wide-angle twilight sky polarimetry,” Cosmic Res. 51, 235–240 (2013).CSCRA70010-9525http://dx.doi.org/10.1134/S0010952513040084Google Scholar


O. S. Ugolnikov and I. A. Maslov, “Studies of the stratosphere aerosol layer based on polarization measurements of the twilight sky,” Cosmic Res. 47, 198–207 (2009).CSCRA70010-9525http://dx.doi.org/10.1134/S0010952509030022Google Scholar


E. Boesche et al., “Polarization of skylight in the O2A band: effects of aerosol properties,” Appl. Opt. 47, 3467 (2008).APOPAI0003-6935http://dx.doi.org/10.1364/AO.47.003467Google Scholar


E. Boesche et al., “Effect of aerosol microphysical properties on polarization of skylight: sensitivity study and measurements,” Appl. Opt. 45, 8790 (2006).APOPAI0003-6935http://dx.doi.org/10.1364/AO.45.008790Google Scholar


J. Zeng, Q. Han and J. Wang, “High-spectral resolution simulation of polarization of skylight: sensitivity to aerosol vertical profile,” Geophys. Res. Lett. 35, L20801 (2008).GPRLAJ0094-8276http://dx.doi.org/10.1029/2008GL035645Google Scholar


I. Aben, D. M. Stam and F. Helderman, “The ring effect in skylight polarization,” Geophys. Res. Lett. 28, 519–522 (2001).GPRLAJ0094-8276http://dx.doi.org/10.1029/2000GL011901Google Scholar


I. Aben et al., “Spectral fine-structure in the polarization of skylight,” Geophys. Res. Lett. 26, 591–594 (1999).GPRLAJ0094-8276http://dx.doi.org/10.1029/1999GL900025Google Scholar


N. Pust and J. A. Shaw, “Imaging spectropolarimetry of cloudy skies,” Proc. SPIE 6240, 624006 (2006).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.670612Google Scholar


N. J. Pust and J. A. Shaw, “How good is a single-scattering model of visible-NIR atmospheric skylight polarization?” Proc. SPIE 7461, 74610B (2009).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.828343Google Scholar


N. J. Pust and J. A. Shaw, “Digital all-sky polarization imaging of partly cloudy skies,” Appl. Opt. 47, H190–H198 (2008).APOPAI0003-6935http://dx.doi.org/10.1364/AO.47.00H190Google Scholar


N. J. Pust and J. A. Shaw, “All-sky polarization imaging,” Proc. SPIE 6682, 668204 (2007).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.736330Google Scholar


N. J. Pust and J. A. Shaw, “Dual-field imaging polarimeter using liquid crystal variable retarders,” Appl. Opt. 45, 5470 (2006).APOPAI0003-6935http://dx.doi.org/10.1364/AO.45.005470Google Scholar


N. J. Pust and J. A. Shaw, “Dual-field imaging polarimeter for studying the effect of clouds on sky and target polarization,” Proc. SPIE 5888, 588812 (2005).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.618773Google Scholar


J. A. Shaw et al., “Continuous outdoor operation of an all-sky polarization imager,” Proc. SPIE 7672, 76720A (2010).http://dx.doi.org/10.1117/12.851374Google Scholar


R. Swindle, “Towards improved diagnostics in terrestrial and solar spectropolarimetry,” MS Thesis [PhD thesis], University of Hawai’i at Manoa (2014).Google Scholar


R. Swindle and J. R. Kuhn, “Haleakalā sky polarization: full-sky observations and modeling,” Publ. Astron. Soc. Pac. 127(956), 1061–1076 (2015).PASPAU0004-6280http://dx.doi.org/10.1086/683177Google Scholar


O. S. Ugolnikov, O. V. Postylyakov and I. A. Maslov, “Effects of multiple scattering and atmospheric aerosol on the polarization of the twilight sky,” J. Quant. Spectrosc. Radiat. Transfer 88, 233–241 (2004).JQSRAE0022-4073http://dx.doi.org/10.1016/j.jqsrt.2003.12.033Google Scholar


M. P. Fetrow et al., “Results of a new polarization simulation,” Proc. SPIE 4481, 149 (2002).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.452902Google Scholar


A. Berk et al., “MODTRAN5: 2006 update,” Proc. SPIE 6233, 62331F (2006).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.665077Google Scholar


D. M. Harrington, J. R. Kuhn and S. Hall, “Deriving telescope Mueller matrices using daytime sky polarization observations,” Publ. Astron. Soc. Pac. 123, 799–811 (2011).PASPAU0004-6280http://dx.doi.org/10.1086/660894Google Scholar


D. M. Harrington et al., “Achromatizing a liquid-crystal spectropolarimeter: retardance vs. Stokes-based calibration of HiVIS,” Publ. Astron. Soc. Pac. 122, 420–438 (2010).PASPAU0004-6280http://dx.doi.org/10.1086/651621Google Scholar


C. R. Givens and A. B. Kostinski, “A simple necessary and sufficient condition on physically realizable Mueller matrices,” J. Mod. Opt. 40, 471–481 (1993).http://dx.doi.org/10.1080/09500349314550471Google Scholar


A. B. Kostinski, C. R. Givens and J. M. Kwiatkowski, “Constraints on Mueller matrices of polarization optics,” Appl. Opt. 32, 1646 (1993).APOPAI0003-6935http://dx.doi.org/10.1364/AO.32.001646Google Scholar


Y. Takakura and M.-P. Stoll, “Passivity test of Mueller matrices in the presence of additive Gaussian noise,” Appl. Opt. 48, 1073 (2009).APOPAI0003-6935http://dx.doi.org/10.1364/AO.48.001073Google Scholar


A. Adhikari, K. Dev and A. Asundi, “Subwavelength metrological characterization by Mueller matrix polarimeter and finite difference time domain method,” Opt. Lasers Eng. 86, 242–247 (2016).http://dx.doi.org/10.1016/j.optlaseng.2016.06.014Google Scholar


E. Compain, S. Poirier and B. Drevillon, “General and self-consistent method for the calibration of polarization modulators, polarimeters, and Mueller-matrix ellipsometers,” Appl. Opt. 38, 3490 (1999).APOPAI0003-6935http://dx.doi.org/10.1364/AO.38.003490Google Scholar


J. C. del Toro Iniesta, Introduction to Spectropolarimetry Cambridge University Press, Cambridge (2003).Google Scholar


J. C. del Toro Iniesta and M. Collados, “Optimum modulation and demodulation matrices for solar polarimetry,” Appl. Opt. 39, 1637 (2000).APOPAI0003-6935http://dx.doi.org/10.1364/AO.39.001637Google Scholar


A. de Martino et al., “Optimized Mueller polarimeter with liquid crystals,” Opt. Lett. 28, 616 (2003).http://dx.doi.org/10.1364/OL.28.000616Google Scholar


K. Nagaraju et al., “An efficient modulation scheme for dual beam polarimetry,” Bull. Astron. Soc. India 35, 307 (2007).Google Scholar


S. Tomczyk et al., “Wavelength-diverse polarization modulators for Stokes polarimetry,” Appl. Opt. 49, 3580–3586 (2010).APOPAI0003-6935http://dx.doi.org/10.1364/AO.49.003580Google Scholar


P. G. Nelson et al., “The visible spectro-polarimeter (ViSP) for the advanced technology solar telescope,” Proc. SPIE 7735, 77358C (2010).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.857610Google Scholar


F. Snik et al., “Design of a full-Stokes polarimeter for VLT/X-shooter,” Proc. SPIE 8446, 844625 (2012).http://dx.doi.org/10.1117/12.926163Google Scholar


A. G. de Wijn et al., “The polychromatic polarization modulator,” Ground-based and Airborne Instrumentation for Astronomy III Ed., pp. 7735, 143 (2010).http://dx.doi.org/10.1117/12.857745Google Scholar


A. G. D. Wijn et al., “Wavelength-diverse polarization modulators for Stokes polarimetry,” Solar Polarization 6. Proc. of a Conf. Held in Maui, Vol. 437, p. 413 (2011).Google Scholar


D. Gisler, A. Feller and A. M. Gandorfer, “Achromatic liquid crystal polarisation modulator,” Proc. SPIE 4843, 45 (2003).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.458835Google Scholar


Y. Hanaoka, “Ferroelectric liquid crystal polarimeter for high-cadence Hα imaging polarimetry,” Sol. Phys. 222, 265–278 (2004).http://dx.doi.org/10.1023/B:SOLA.0000043581.05390.2fGoogle Scholar


C. Xu et al., “Polarimeter with two ferroelectric liquid-crystal modulators attached to the Yunnan solar tower,” Appl. Opt. 45, 8428 (2006).APOPAI0003-6935http://dx.doi.org/10.1364/AO.45.008428Google Scholar


A. R. Dahlberg, N. J. Pust and J. A. Shaw, “Effects of surface reflectance on skylight polarization measurements at the Mauna Loa observatory,” Opt. Express 19, 16008–16021 (2011).http://dx.doi.org/10.1364/OE.19.016008Google Scholar


A. R. Dahlberg, N. J. Pust and J. A. Shaw, “All-sky imaging polarimeter measurements of visible and NIR skylight at Mauna Loa, Hawaii,” Proc. SPIE 7461, 746107 (2009).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.826537Google Scholar


D. M. Harrington and J. R. Kuhn, “Spectropolarimetric observations of Herbig Ae/Be stars. I. HiVIS spectropolarimetric calibration and reduction techniques,” Publ. Astron. Soc. Pac. 120, 89–117 (2008).PASPAU0004-6280http://dx.doi.org/10.1086/528881Google Scholar

Biographies for the authors are not available.

David M. Harrington, Jeffrey R. Kuhn, Arturo López Ariste, "Daytime sky polarization calibration limitations," Journal of Astronomical Telescopes, Instruments, and Systems 3(1), 018001 (24 January 2017). http://dx.doi.org/10.1117/1.JATIS.3.1.018001
Submission: Received 16 November 2016; Accepted 30 November 2016







Atmospheric modeling

Back to Top