© ESO, 2014

1. Introduction

Astronomical polarimetry and spectropolarimetry, which have historically been under-utilized, are now routinely being used in many experiments to search for and test new physics. Besides the diverse application of polarimetric techniques in solar and stellar astronomy, two areas of particular note are polarimetric observations of the cosmic microwave background (CMB) and spectropolarimetric observations of supernovae (SNe).

Detection of inflation signatures on the CMB using polarimetric imaging were recently announced by the BICEP2 team (Ade et al. 2014b,a). This result would confirm one of astrophysics most cherished but experimentally elusive theories. After a review of the public BICEP2 data, some have questioned the interpretation of the BICEP2 data and suggested that the foreground dust contamination had been underestimated or that the signal could have been produced by other mechanisms such as “radio loops” (Liu et al. 2014). The BICEP2 team have stood by their results although they have moderated their language in the peer-reviewed version of their paper (Ade et al. 2014b). If valid, the BICEP2 discovery continues the trend set by all-sky mapping satellites like COBE, WMAP, and Planck and various ground and balloon surveys of the CMB being perhaps the single most fruitful source to test cosmological theory and sometimes even fundamental physics. Further results related to CMB polarization are expected in the near future from BICEP3 (Kuo & BICEP3 and POLAR1 Collaborations 2013), the Keck Array (Staniszewski et al. 2012), Polarbear (The POLARBEAR Collaboration 2014), and the Planck team (Planck Collaboration XIX 2014). It is clear that the CMB will continue to be a vital source of cutting edge science; but, as the controversy around the BICEP2 announcement indicates, polarization work is tricky and extreme care must be made to handle properly statistical and systematic error in the data.

Other research groups have been using spectropolarimetric observations to investigate the systematic behavior and statistical variation of thermonuclear (Type Ia) and core-collapse SNe (all other types) arising from asymmetries. (The Wang & Wheeler 2008 Annual Review is a good introduction.) Spectropolarimetry is the only known way to probe asymmetries of unresolved sources such as extra-galactic supernovae so it will continue to be an effective research tool. The detection of asymmetries has helped place much needed constraints on SNe explosion physics/progenitor scenarios (recent examples include Maund et al. 2010b; Tanaka et al. 2012; Patat et al. 2012; Zelaya et al. 2013; Maund et al. 2013) and environments (Mauerhan et al. 2014). Asymmetries of Type Ia SNe are of particular interest as high-redshift events were famously used to discover cosmic acceleration (Riess et al. 1998; Perlmutter et al. 1999), a result for which leads of two large teams shared the 2011 Nobel Prize. This acceleration implies the existence of so-called “dark energy”, whose physical nature is completely puzzling and is perhaps the biggest unsolved problem in astrophysics. As such, any possible correlations between metrics of asymmetry in Type Ia SNe and other observables (Leonard et al. 2005; Wang et al. 2007; Maund et al. 2010a) are of high importance because correlations may be used to identify contaminating Type Ia sub-classes (Wang et al. 2013) or apply statistical corrections to the maximum brightness, both of which could help refine the measurement of the cosmological constant (Maeda et al. 2011). Spectropolarimetry of SNe is a blooming field but, just as with the CMB, handling the statistical and systematic error of polarimetric data for SNe is challenging.

As both polarimetric imaging and spectropolarimetry require dividing incoming light into bins, researchers in these fields often find themselves in photon-starved situations. Indeed, in both subjects discussed above, photon-limited data at low signal-to-noise must be confronted. (Low signal-to-noise is often the hallmark of forefront science by its very nature.) This problem is compounded in polarimetric imaging when a narrow-band filter is needed. Working at low signal-to-noise will continue to be an issue for polarimetrists.

Polarization data is often measured in Stokes parameters but expressed in human-digestible quantities like degrees of polarization and angles of polarization. Despite this, there has been relatively little attention given to the surprisingly complicated statistics involved when making transformations to these quantities at low signal-to-noise. Most of the published literature of this fundamental, practical topic has focused on the analysis of linear polarization alone that involves marginalizations over a two-dimensional gaussian distribution over the Q-U plane cross section of the Poincaré sphere (the “Poincaré disk” is a useful expression although it is already in use in other mathematical contexts). Marginalization over the linear polarization angle results in the Rice distribution (Rice 1945) while marginalization over the degree of linear polarization results in another distribution given by Vinokur (1965). The authors refer to this marginalized approach as the “one-dimensional problem”. A small number of papers such as Wang et al. (1997), Vaillancourt (2006), and Quinn (2012) have used a non-marginalized distribution over the Poincaré disk (the “two-dimensional problem”), with the latter being the first to put the subject on a rigorous mathematical footing. To date, there has been little to no attention paid to the three-dimensional problem defined by a three-dimensional distribution over the Poincaré sphere. Progress on higher dimensional generalizations of the two-dimensional problem have been made, for instance, by members of the Planck Collaboration (Planck Collaboration XIX 2014; Montier et al. 2014a,b) by including intensity in the analysis (see also Viola et al. 2014) but researchers have thus far generally ignored circular polarization when measuring linear polarization in astronomical sources (and vice versa); however, the measurement of one is not actually independent of the other. This effect is generally negligible except when one is working at very small signal-to-noise ratios or with objects with polarization near 100% (see Quinn 2012 for details on how large polarization introduces complication). The aim of this paper is to elucidate some of this interdependence and to see how it relates to the more common ways of treating linear polarization measurements at low signal-to-noise via the Rice distribution or through Bayesian techniques as outlined in Vaillancourt (2006) and Quinn (2012). The results should be of interest to all polarimetrists and those those studying the CMB or SNe in particular.

2. The sampling distribution

A very practical way of encoding the polarization state of electromagnetic radiation is the Stokes parameter formalism invented by George Stokes in 1852. Stokes parameters consist of four quantities: I (intensity), Q and U (related to linear polarization), and V (related to circular polarization)¹. They are convenient quantities to measure but it is easier for people to think about polarization in terms of percentages and angles. This may be done by introducing a spherical coordinate system (p,θ,φ) on points of the Poincaré sphere (a unit ball in Cartesian coordinates centered on the origin of $\frac{\mathit{Q}}{{\mathit{I}}_{\mathrm{0}}}$, $\frac{\mathit{U}}{{\mathit{I}}_{\mathrm{0}}}$, and $\frac{\mathit{V}}{{\mathit{I}}_{\mathrm{0}}}$ axes), which represents all possible polarization states. Accessible modern introductions to Stokes parameters and their theory are given in del Toro Iniesta (2003) and Landi Degl’Innocenti (2002).

If an astronomical source’s polarization state is measured with gaussian error (Σ_I,Σ_Q,Σ_U,Σ_V) on each Stokes parameter (I,Q,U,V), how do the corresponding spherical coordinates (p,θ,φ) relate to the “true” polarization state (p₀,θ₀,φ₀) determined uniquely from (I₀,Q₀,U₀,V₀)? (Unsubscripted variables will usually refer to measured values while “0”-subscripted values will refer to “true” values. Later, an overbar on some variables and distributions will indicate signal-to-noise related quantities.) As it turns out, a spherical coordinate system introduces non-trivial biases into the measured quantities and the measurement itself must be “corrected” to remove them. This paper extends Quinn (2012), which studied these effects in polar coordinates before discussing the Bayesian approach, to spherical coordinates. The reader is referred to that paper for a more extensive introduction on Stokes parameters.

Under the assumptions that the source intensity is large compared to the measurement error, i.e., I₀ ≫ Σ_I which implies I ≈ I₀ (Quinn 2012)², and the error on a Stokes measurement is described by a three-dimensional uncorrelated ellipsoidal Gaussian distribution (with semi-major axes Σ_Q, Σ_U, Σ_V), the normalized distribution, F_C, of the measured Stokes parameters (Q, U, V) given the “true” Stokes parameters (Q₀, U₀, V₀) is $\begin{array}{ccc}{\mathit{F}}_{\mathit{C}}\mathrm{\left(}\mathit{Q,U,V}\mathrm{\right|}{\mathit{Q}}_{\mathrm{0}}\mathit{,}{\mathit{U}}_{\mathrm{0}}\mathit{,}{\mathit{V}}_{\mathrm{0}}\mathit{,}{\mathrm{\Sigma }}_{\mathit{Q}}\mathit{,}{\mathrm{\Sigma }}_{\mathit{U}}\mathit{,}{\mathrm{\Sigma }}_{\mathit{V}}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathrm{\Sigma }}_{\mathit{Q}}{\mathrm{\Sigma }}_{\mathit{U}}{\mathrm{\Sigma }}_{\mathit{V}}}{\mathrm{e}}^{\mathrm{-}\left(\frac{\mathrm{(}\mathit{Q}\mathrm{-}{\mathit{Q}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathrm{\Sigma }}_{\mathit{Q}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{(}\mathit{U}\mathrm{-}{\mathit{U}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathrm{\Sigma }}_{\mathit{U}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{(}\mathit{V}\mathrm{-}{\mathit{V}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathrm{\Sigma }}_{\mathit{V}}^{\mathrm{2}}}\right)}\mathit{.}& & \end{array}$(1)These parameters are not usually worked with directly. Instead “normalized Stokes parameters” are defined as q ≡ Q/I₀, u ≡ U/I₀, v ≡ V/I₀ (with also q₀ ≡ Q₀/I₀, u₀ ≡ U₀/I₀, and v₀ ≡ V₀/I₀), and σ_q ≡ Σ_Q/I₀, σ_u ≡ Σ_U/I₀, and σ_v ≡ Σ_V/I₀, yielding a new normalized distribution $\begin{array}{ccc}{\mathit{f}}_{\mathit{C}}\mathrm{\left(}\mathit{q,u,v}\mathrm{\right|}{\mathit{q}}_{\mathrm{0}}\mathit{,}{\mathit{u}}_{\mathrm{0}}\mathit{,}{\mathit{v}}_{\mathrm{0}}\mathit{,}{\mathit{\sigma }}_{\mathit{q}}\mathit{,}{\mathit{\sigma }}_{\mathit{u}}\mathit{,}{\mathit{\sigma }}_{\mathit{v}}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{\sigma }}_{\mathit{q}}{\mathit{\sigma }}_{\mathit{u}}{\mathit{\sigma }}_{\mathit{v}}}{\mathrm{e}}^{\mathrm{-}\left(\frac{\mathrm{(}\mathit{q}\mathrm{-}{\mathit{q}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{q}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{(}\mathit{u}\mathrm{-}{\mathit{u}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{u}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{(}\mathit{v}\mathrm{-}{\mathit{v}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{v}}^{\mathrm{2}}}\right)}\mathit{.}& & \end{array}$(2)Lastly, one most often works with “signal-to-noise” versions of the Stokes parameters defined as $\overline{\mathit{q}}\mathrm{\equiv }\mathit{Q}\mathit{/}{\mathrm{\Sigma }}_{\mathit{Q}}\mathrm{=}\mathit{q}\mathit{/}{\mathit{\sigma }}_{\mathit{q}}$, $\overline{\mathit{u}}\mathrm{\equiv }\mathit{U}\mathit{/}{\mathrm{\Sigma }}_{\mathit{U}}\mathrm{=}\mathit{u}\mathit{/}{\mathit{\sigma }}_{\mathit{u}}$, and $\overline{\mathit{v}}\mathrm{\equiv }\mathit{V}\mathit{/}{\mathrm{\Sigma }}_{\mathit{V}}\mathrm{=}\mathit{v}\mathit{/}{\mathit{\sigma }}_{\mathit{v}}$ which gives another distribution ${\overline{\mathit{f}}}_{\mathit{C}}\mathrm{\left(}\overline{\mathit{q}}\mathit{,}\overline{\mathit{u}}\mathit{,}\overline{\mathit{v}}\mathrm{\right|}{\overline{\mathit{q}}}_{\mathrm{0}}\mathit{,}{\overline{\mathit{u}}}_{\mathrm{0}}\mathit{,}{\overline{\mathit{v}}}_{\mathrm{0}}\mathit{,}{\mathit{\sigma }}_{\mathit{q}}\mathit{,}{\mathit{\sigma }}_{\mathit{u}}\mathit{,}{\mathit{\sigma }}_{\mathit{v}}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{(}\overline{\mathit{q}}\mathrm{-}{\overline{\mathit{q}}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\overline{\mathit{u}}\mathrm{-}{\overline{\mathit{u}}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\overline{\mathit{v}}\mathrm{-}{\overline{\mathit{v}}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}}}\mathit{,}$(3)which is normalized, ${{}^{\mathrm{\int }}}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}{{}^{\mathrm{\int }}}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}{{}^{\mathrm{\int }}}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}{\overline{\mathit{f}}}_{\mathit{C}}\hspace{0.17em}\mathrm{d}\overline{\mathit{q}}\hspace{0.17em}\mathrm{d}\overline{\mathit{u}}\hspace{0.17em}\mathrm{d}\overline{\mathit{v}}\mathrm{=}\mathrm{1}$. This distribution no longer explicitly depends on the σ’s and naturally handles ellipsoidal distributions. The σ’s do, however, play a minor role: they limit the possible values of q₀, u₀, and v₀ so it is important not to forget their significance (Quinn 2012).

Fig. 1 Illustration of the spherical coordinate system used in the paper. The variable θ is used for the azimuthal angle with a range of (−π,π ], and the variable ϕ is measured from the positive v-axis with a range of [ 0,π ].

We wish to now transform our equations to spherical coordinates as the Poincaré sphere suggests. The transformation equations³ are $\begin{array}{ccc}& & \mathit{p}\mathrm{=}\sqrt{{\mathit{q}}^{\mathrm{2}}\mathrm{+}{\mathit{u}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}^{\mathrm{2}}}\\ & & \mathit{\theta }\mathrm{=}\arctan \mathrm{\left(}\mathit{u,q}\mathrm{\right)}\\ & & \mathit{\varphi }\mathrm{=}\arctan \mathrm{\left(}\sqrt{{\mathit{q}}^{\mathrm{2}}\mathrm{+}{\mathit{u}}^{\mathrm{2}}}\mathit{,v}\mathrm{\right)}\end{array}$(4)which has inverse $\begin{array}{ccc}& & \mathit{q}\mathrm{=}\mathit{p}\sin \mathrm{\left(}\mathit{\varphi }\mathrm{\right)}\cos \mathrm{\left(}\mathit{\theta }\mathrm{\right)}\\ & & \mathit{u}\mathrm{=}\mathit{p}\sin \mathrm{\left(}\mathit{\varphi }\mathrm{\right)}\sin \mathrm{\left(}\mathit{\theta }\mathrm{\right)}\\ & & \mathit{v}\mathrm{=}\mathit{p}\cos \mathrm{\left(}\mathit{\varphi }\mathrm{\right)}\mathit{,}\end{array}$(5)where θ is the angle in the q-u plane and ϕ is the angle from the positive v-axis (see Fig. 1). The “true values”, p₀, θ₀, and ϕ₀, are similarly related to q₀, u₀, and v₀.

It is crucial not to confuse the degree of polarization, p, with the degree of linear polarization in the following. The degree of linear polarization is a separate quantity equal to $\sqrt{{\mathit{q}}^{\mathrm{2}}\mathrm{+}{\mathit{u}}^{\mathrm{2}}}$ or psin(ϕ) and, although it is certainly important, is not discussed in this paper.

We now restrict ourselves to the case where all three standard deviations are equal (Σ_Q = Σ_U = Σ_V ≡ Σ) so that σ_q = σ_u = σ_v ≡ σ. The distribution in spherical coordinates is then $\begin{array}{ccc}{\mathit{f}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{p,\theta ,\varphi }\mathrm{\right|}{\mathit{p}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathit{,\sigma }\mathrm{)}\mathrm{=}\frac{{\mathit{p}}^{\mathrm{2}}\sin \mathit{\varphi }}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{\sigma }}^{\mathrm{3}}}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{p}}^{\mathrm{2}}\mathrm{+}{\mathit{p}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}\mathrm{2}\mathit{p}{\mathit{p}}_{\mathrm{0}}\left(\sin \mathit{\varphi }\sin {\mathit{\varphi }}_{\mathrm{0}}\cos \mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{\theta }}_{\mathrm{0}}\mathrm{)}\mathrm{+}\cos \mathit{\varphi }\cos {\mathit{\varphi }}_{\mathrm{0}}\right)}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\mathit{.}& & \end{array}$(6)The p²sinϕ in the numerator of the factor out front is due to the Jacobian of the transformation, i.e., dq du dv = p²sinϕ dp dθ dϕ⁴. It is rather promising that the analytic form of Eq. (6) is similar to and not much more complicated than the two-dimensional polar version. Already one can see that the functional form of the equation is not independent of the true value of the circular polarization.

The barred version is $\begin{array}{ccc}{\overline{\mathit{f}}}^{\mathrm{\prime }}\mathrm{\left(}\overline{\mathit{p}}\mathit{,\theta ,\varphi }\mathrm{\right|}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathit{,\sigma }\mathrm{)}\mathrm{=}\frac{{\overline{\mathit{p}}}^{\mathrm{2}}\sin \mathit{\varphi }}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}\frac{{\overline{\mathit{p}}}^{\mathrm{2}}\mathrm{+}\overline{\mathit{p}}\begin{array}{c}\mathrm{2}\\ \mathrm{0}\end{array}\mathrm{-}\mathrm{2}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\left(\sin \mathit{\varphi }\sin {\mathit{\varphi }}_{\mathrm{0}}\cos \mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{\theta }}_{\mathrm{0}}\mathrm{)}\mathrm{+}\cos \mathit{\varphi }\cos {\mathit{\varphi }}_{\mathrm{0}}\right)}{\mathrm{2}}}\mathit{,}& & \end{array}$(7)where $\overline{\mathit{p}}\mathrm{\equiv }\sqrt{{\mathit{Q}}^{\mathrm{2}}\mathrm{+}{\mathit{U}}^{\mathrm{2}}\mathrm{+}{\mathit{V}}^{\mathrm{2}}}\mathit{/}\mathrm{\Sigma }\mathrm{=}\mathit{p}\mathit{/}\mathit{\sigma }$ and ${\overline{\mathit{p}}}_{\mathrm{0}}\mathrm{\equiv }\sqrt{{\mathit{Q}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{+}{\mathit{U}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{+}{\mathit{V}}_{\mathrm{0}}^{\mathrm{2}}}\mathit{/}\mathrm{\Sigma }\mathrm{=}{\mathit{p}}_{\mathrm{0}}\mathit{/}\mathit{\sigma }$. This is still normalized, ${{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{\pi }}{{}^{\mathrm{\int }}}_{\mathrm{-}\mathit{\pi }}^{\mathit{\pi }}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathrm{1}\mathit{/}\mathit{\sigma }}{\overline{\mathit{f}}}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{d}\overline{\mathit{p}}\hspace{0.17em}\mathrm{d}\mathit{\theta }\hspace{0.17em}\mathrm{d}\mathit{\varphi }\mathrm{=}\mathrm{1}$.

3. Classical estimators

Now that ${\overline{\mathit{f}}}^{\mathrm{\prime }}\mathrm{\left(}\overline{\mathit{p}}\mathit{,\theta ,\varphi }\mathrm{\right|}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathit{,\sigma }\mathrm{)}$ (the sampling distribution in signal-to-noise variables) is at hand, two important classical estimators may be calculated. These are the so-called “Most Likely” and “Most Probable” estimators. (Note that the prime is not indicating a derivative in ${\overline{\mathit{f}}}^{\mathrm{\prime }}$ but serving as a warning that we are using the q-u plane angle, θ, instead of the sky angle for linear polarization.)

3.1. The “Most Likely” solution

The “Most Likely” (ML) solution is obtained by maximizing ${\overline{\mathit{f}}}^{\mathrm{\prime }}$ with respect to ${\overline{\mathit{p}}}_{\mathrm{0}}$, θ₀, and ϕ₀. Physically, this corresponds to finding the value of $\mathrm{\left(}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}$ which makes the measured value of $\mathrm{\left(}\overline{\mathit{p}}\mathit{,\theta ,\varphi }\mathrm{\right)}$ the most statistically likely, hence the name. A general solution may be found by solving the system $\frac{\mathit{\partial }{\overline{\mathit{f}}}^{\mathrm{\prime }}}{\mathit{\partial }{\overline{\mathit{p}}}_{\mathrm{0}}}\mathrm{=}\mathrm{0}$, $\frac{\mathit{\partial }{\overline{\mathit{f}}}^{\mathrm{\prime }}}{\mathit{\partial }{\mathit{\theta }}_{\mathrm{0}}}\mathrm{=}\mathrm{0}$, and $\frac{\mathit{\partial }{\overline{\mathit{f}}}^{\mathrm{\prime }}}{\mathit{\partial }{\mathit{\varphi }}_{\mathrm{0}}}\mathrm{=}\mathrm{0}$ for ${\overline{\mathit{p}}}_{\mathrm{0}}$, θ₀, and ϕ₀, which are then treated as estimators. Technically, a second derivative test must also be done to check that the point is actually a maximum. For functions of more than two variables, this test requires checking that all the eigenvalues of the function’s Hessian matrix are positive. (All negative would be a minima, while mixed values would be inconclusive.) The equations involved in the test are very large and cumbersome. From graphical investigation or intuition, it is obvious, however, that the solution to be found is a maximum.

From the azimuthal equation, $\frac{\mathit{\partial }{\overline{\mathit{f}}}^{\mathrm{\prime }}}{\mathit{\partial }{\mathit{\theta }}_{\mathrm{0}}}\mathrm{=}\mathrm{0}$, it is quickly deduced that θ₀ = θ everywhere except for values that correspond to an origin (where ${\overline{\mathit{p}}}_{\mathrm{0}}$ or $\overline{\mathit{p}}$ equals zero) or a circular polarization axis (where ϕ₀ or ϕ equals zero or π) of the true and measured Poincaré spheres. Further tedious but straight-forward algebra allows the full solution to be found, which is just

$\begin{array}{ccc}& & {\overline{\mathit{p}}}_{\mathrm{0}\mathit{,}\mathrm{ML}}\mathrm{=}\overline{\mathit{p}}\\ & & {\mathit{\theta }}_{\mathrm{0}\mathit{,}\mathrm{ML}}\mathrm{=}\mathit{\theta }\\ & & \end{array}$(8)

This solution involves no “bias corrections” at all. What you measure is your best estimate for the “true” polarization.

This solution is consistent with the polar coordinate case (that is, the two-dimensional case corresponding to ϕ → π/ 2, ϕ₀ → π/ 2, and σ_v → 0) presented in Quinn (2012). It is only when the Rice distribution (a one-dimensional distribution arising from marginalization) is used that a non-trivial solution is found for the ML method (Simmons & Stewart 1985). As marginalization intentionally removes information present in the original distribution, the logical basis for the use of non-trivial ML estimators found from a marginalized distribution to perform a “bias correction” is questionable. This is underscored by the fact that the p₀, θ₀, and ϕ₀ are not stochastic variables but input parameters in ${\overline{\mathit{f}}}^{\mathrm{\prime }}$ and should not be treated as such as is done in the ML approach. The logically correct way to “invert” the distribution is to use Bayes’ Theorem.

3.2. The “Most Probable” solution

The “Most Probable” (MP) solution is found by maximizing ${\overline{\mathit{f}}}^{\mathrm{\prime }}$ with respect to $\overline{\mathit{p}}$, θ, and ϕ. This corresponds to finding the $\mathrm{\left(}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}$ point that produces a distribution of measured points with a maximum at the actual measured point $\mathrm{\left(}\overline{\mathit{p}}\mathit{,\theta ,\varphi }\mathrm{\right)}$. A general solution may be found by solving the system $\frac{\mathit{\partial }{\overline{\mathit{f}}}^{\mathrm{\prime }}}{\mathit{\partial }\overline{\mathit{p}}}\mathrm{=}\mathrm{0}$, $\frac{\mathit{\partial }{\overline{\mathit{f}}}^{\mathrm{\prime }}}{\mathit{\partial \theta }}\mathrm{=}\mathrm{0}$, and $\frac{\mathit{\partial }{\overline{\mathit{f}}}^{\mathrm{\prime }}}{\mathit{\partial \varphi }}\mathrm{=}\mathrm{0}$ for ${\overline{\mathit{p}}}_{\mathrm{0}}$, θ₀, and ϕ₀. (As before, a second derivative test is technically necessary to determine if the solution is actually a maximum as opposed to a minimum or some point of mixed inflection but this test is similarly impractical to perform because the equations end up being quite large. Graphically it can be shown that the solutions are maximums.)

Again, one quickly finds θ₀ = θ from the azimuthal equation, $\frac{\mathit{\partial }{\overline{\mathit{f}}}^{\mathrm{\prime }}}{\mathit{\partial \theta }}\mathrm{=}\mathrm{0}$. Using this condition, the ϕ and $\overline{\mathit{p}}$ equations yield ${\overline{\mathit{p}}}^{\mathrm{2}}\mathrm{=}\mathrm{2}\mathrm{+}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\cos \mathrm{(}\mathit{\varphi }\mathrm{-}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{)}$ and $\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\cos \mathrm{\left(}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}\sin \mathrm{(}\mathit{\varphi }{\mathrm{)}}^{\mathrm{2}}\mathrm{=}\cos \mathrm{(}\mathit{\varphi }\mathrm{\left)}\mathrm{\right(}\mathrm{1}\mathrm{+}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\sin \mathrm{\left(}\mathit{\varphi }\mathrm{\right)}\sin \mathrm{\left(}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}\mathrm{)}$, respectively. This two-equation system is somewhat difficult to solve and the solution is best accomplished with a computer algebra system. The full analytic solution is

$\begin{array}{ccc}& & {\overline{\mathit{p}}}_{\mathrm{0}\mathit{,}\mathrm{MP}}\mathrm{=}\frac{\sqrt{\cot \mathrm{(}\mathit{\varphi }{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\overline{\mathit{p}}}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathrm{)}}^{\mathrm{2}}}}{\overline{\mathit{p}}}\\ & & {\mathit{\theta }}_{\mathrm{0}\mathit{,}\mathrm{MP}}\mathrm{=}\mathit{\theta }\\ & & \end{array}$(9)

which is only valid for points satisfying the condition $\overline{\mathit{p}}\mathit{>}\sqrt{\mathrm{1}\mathrm{+}\csc \mathrm{(}\mathit{\varphi }{\mathrm{)}}^{\mathrm{2}}}\mathit{,}$(10)where csc(x) is the cosecant function. Thus there are two coupled “bias corrections” that must be made to the data: one for $\overline{\mathit{p}}$ and one for ϕ. A graphical representation of the preceding solution is given in Fig. 2. The left-most extension of the region consisting of points satisfying Eq. (10) (shown in blue in the figure) approaches the point $\mathrm{\left(}\overline{\mathit{p}}\mathit{,\varphi }\mathrm{\right)}\mathrm{=}\mathrm{(}\sqrt{\mathrm{2}}\mathit{,\pi }\mathit{/}\mathrm{2}\mathrm{)}$. In the special case of ϕ = π/ 2 and $\overline{\mathit{p}}\mathit{>}\sqrt{\mathrm{2}}$, the correction is just ϕ_0,MP = π/ 2 and ${\overline{\mathit{p}}}_{\mathrm{0}\mathit{,}\mathrm{MP}}\mathrm{=}\overline{\mathit{p}}\mathrm{-}\mathrm{2}\mathit{/}\overline{\mathit{p}}$.

Fig. 2 Visualization of the $\mathrm{\left(}\overline{\mathit{p}}\mathit{,\varphi }\mathrm{\right)}$ bias correction field of Eq. (9). The tail of each arrow is located at a measured $\mathrm{\left(}\overline{\mathit{p}}\mathit{,\varphi }\mathrm{\right)}$ point and the head is attached to the corresponding “Most Probable” estimate of $\mathrm{\left(}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}$. This correction should only be applied to $\mathrm{\left(}\overline{\mathit{p}}\mathit{,\varphi }\mathrm{\right)}$ points lying in the blue region given by Eq. (10). See the text for discussion about $\mathrm{\left(}\overline{\mathit{p}}\mathit{,\varphi }\mathrm{\right)}$ points outside the blue region.

Fig. 3 Behavior of the bias correction field (Eq. (9)) on the boundary of the blue valid region (Eq. (10)). The left-most point on the boundary of the blue region is at $\mathrm{(}\sqrt{\mathrm{2}}\mathit{,\pi }\mathit{/}\mathrm{2}\mathrm{)}$.

The “bias correction field” in Fig. 2 shows that the correction applied to $\overline{\mathit{p}}$ gets larger as $\overline{\mathit{p}}$ gets smaller (the actual magnitude is ϕ-dependent). This is consistent with the known behavior in the two-dimensional case. Figure 2 also shows that ϕ itself must be corrected: it should be made larger when ϕ>π/ 2 and smaller when ϕ<π/ 2. The magnitude of this correction also tends to be larger as $\overline{\mathit{p}}$ gets smaller. The exact magnitude of the ϕ-correction at a given $\overline{\mathit{p}}$ depends on ϕ itself with the magnitude being largest in the “mid-latitudes” of the Poincaré sphere and zero on the equator. Generally, $\mathrm{\left(}\overline{\mathit{p}}\mathit{,\varphi }\mathrm{\right)}$-points not on the equator are “attracted” towards the poles, an effect most prominent at low signal-to-noise (between, say, $\sqrt{\mathrm{2}}\mathit{<}\overline{\mathit{p}}\lesssim \mathrm{5}$). The bias correction vector field diminishes as $\overline{\mathit{p}}$ gets larger but for ϕ values near the poles is seen to remain significant even for $\overline{\mathit{p}}\mathit{>}\mathrm{10}$, a regime typically thought of as having signal-to-noise high enough not to have to worry about such effects.

For points infinitesimally close to the boundary defined by Eq. (10), the bias correction for ϕ estimates ϕ₀ = 0 or ϕ₀ = π, whichever pole is nearest. Figure 3 illustrates this behavior. The special case of ϕ = π/ 2 “corrects” to ϕ₀ = π/ 2. More specifically, points on the boundary with π/ 2 <ϕ<π all bias correct to ϕ₀ = π while those with 0 <ϕ<π/ 2 all bias correct to ϕ₀ = 0. The $\overline{\mathit{p}}$ correction on the boundary acts typically and corrects more and more towards the origin as $\overline{\mathit{p}}$ gets smaller.

From Figs. 2 and 3, it can be seen that every possible $\mathrm{\left(}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}$ with ${\overline{\mathit{p}}}_{\mathrm{0}}\mathit{>}\mathrm{0}$ and sinϕ₀ ≠ 0 is associated with a $\mathrm{\left(}\overline{\mathit{p}}\mathit{,\varphi }\mathrm{\right)}$ point within the blue region and vice versa. Unfortunately, a closed form solution for the inverse of Eq. (9) was not found and may not exist.

Intuitively⁵, the $\overline{\mathit{p}}$ correction may be understood by considering a ${\overline{\mathit{p}}}_{\mathrm{0}}$ value of zero. The measured value of $\overline{\mathit{p}}$ cannot be negative and, due to scatter from measurement error, you expect that $\overline{\mathit{p}}$ will almost always be positive (that is, there is a negligibly small chance of zero) even though ${\overline{\mathit{p}}}_{\mathrm{0}}\mathrm{=}\mathrm{0}$. Thus, the measurement process “biases” $\overline{\mathit{p}}$ to larger values by some amount. Therefore, a correction should be made that subtracts a bit from the measured value to obtain a good estimate of ${\overline{\mathit{p}}}_{\mathrm{0}}$. Even for a non-zero ${\overline{\mathit{p}}}_{\mathrm{0}}$, this biasing occurs although with diminishing effect as ${\overline{\mathit{p}}}_{\mathrm{0}}$ gets larger. The same line of reasoning can be used to gain an intuitive understanding of the ϕ bias correction. Consider a source in a state of pure circular polarization. This corresponds to a pole of the Poincaré sphere. Let’s focus on the purely right-circularly polarized case (ϕ₀ = 0). Measurements of this state will be clustered around the pole (including outside the Poincaré sphere) but ϕ is never negative and is almost always positive. Similarly for the left-circularly polarized case (ϕ₀ = π), measurements of ϕ will be clustered near that pole but cannot have ϕ greater than π and almost always will have ϕ less than π. Just as with $\overline{\mathit{p}}$, a correction is needed to remove this effect from ϕ. It is important to remark that both the $\overline{\mathit{p}}$ and ϕ bias corrections are purely coordinate effects. It may seem at first as if the ϕ-correction is inconsistent with the symmetry of the system (from Σ_Q = Σ_U = Σ_V) but it is not. While conventionally the v-axis is always chosen to isolate circular polarization, the effect would arise regardless of the orientation of the v-axis within the Poincaré sphere, which preserves the overall symmetry (i.e., the symmetry is spontaneously broken by introduction of the coordinate system).

It must be stressed that the solution of Eq. (9) only applies to measured points satisfying Eq. (10) (illustrated by the blue region in the figures). Presumably, the area outside the blue region divides into three other regions: one characterized by ${\overline{\mathit{p}}}_{\mathrm{0}\mathit{,}\mathrm{MP}}\mathrm{=}\mathrm{0}$, another by ${\overline{\mathit{p}}}_{\mathrm{0}\mathit{,}\mathrm{MP}}\ne \mathrm{0}$ and ϕ_0,MP = 0, and still another by ${\overline{\mathit{p}}}_{\mathrm{0}\mathit{,}\mathrm{MP}}\ne \mathrm{0}$ and ϕ_0,MP = π. In the latter two cases, there could still be a $\overline{\mathit{p}}$ bias correction at work. It is unclear how to determine mathematically a “MP solution” in the non-blue region of Fig. 2. A formal solution may simply not exist there. The authors’ best guess as to the best practice is that measurements within the “triangle” defined by points (0,0), (0,π), and $\mathrm{(}\sqrt{\mathrm{2}}\mathit{,\pi }\mathit{/}\mathrm{2}\mathrm{)}$ correct to (0,π) if ϕ>π/ 2, (0,0) if ϕ<π/ 2, and (0,π/ 2) if ϕ = π/ 2⁶. The remaining points could correct to the same value as the blue region boundary point whose “correction arrow” passes through the point. Another possibility for these later points is to translate the correction arrow from the boundary of the blue region with the same ϕ to the point and use the intersection of that arrow with the $\overline{\mathit{p}}$-axis as the corrected value. Formal solutions in these regions where not, however, found.

Equation (9) does indeed provide an estimate of the $\mathrm{\left(}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}$ point that would produce a distribution of measured points with maximum at $\mathrm{\left(}\overline{\mathit{p}}\mathit{,\varphi }\mathrm{\right)}$. It is however just a point estimate. The actual $\mathrm{\left(}\overline{\mathit{p}}\mathit{,\varphi }\mathrm{\right)}$ distribution for a small values of ${\overline{\mathit{p}}}_{\mathrm{0}}$ tends to be broad, shallow, and nearly flat, especially in the ϕ dimension. In other words, the errors bars associated with each estimate will be large so it is cautioned that the point estimate should not be taken too literally. Using the ${\overline{\mathit{f}}}^{\mathrm{\prime }}$ distribution, one could construct confidence regions. These regions are however not uniquely determined, even in the one-dimensional case (Simmons & Stewart 1985). Construction of confidence regions is reserved for future work.

4. Marginalized distributions

Experimenters often work with only a partial description of a physical system. This may be necessary when independent variables of a distribution of interest related to the system are not measurable, or desired when there are so-called nuisance variables of little concern. In a polarization study, for instance, it may be that only the degree of linear polarization is the focus and not the angle of linear polarization. In such situations, one can marginalize over “unwanted” independent variables to find a new distribution which no longer depends on them but at the cost of losing some information present in the original distribution.

The marginalization possibilities for f′ are to marginalize over (I) p only; (II) θ only; (III) ϕ only; (IV) p and θ; (V) p and ϕ; and (VI) θ and ϕ. The first and last cases are perhaps the most interesting and shall now be investigated.

4.1. The p-marginalized angular distribution

The p-marginalized angular distribution is $\mathit{M}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right|}{\mathit{p}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathit{,\sigma }\mathrm{)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{f}}^{\mathrm{\prime }}\mathrm{(}\mathit{p,\theta ,\varphi }\mathrm{\left|}{\mathit{p}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathit{,\sigma }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{p}$(11)or, in signal-to-noise variables, $\overline{\mathit{M}}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right|}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathit{,\sigma }\mathrm{)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\overline{\mathit{f}}}^{\mathrm{\prime }}\mathrm{(}\overline{\mathit{p}}\mathit{,\theta ,\varphi }\mathrm{\left|}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathit{,\sigma }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\overline{\mathit{p}}\mathit{.}$(12)This equation can be simplified. Start with $\begin{array}{ccc}\overline{\mathit{M}}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right|}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathit{,\sigma }\mathrm{)}\mathrm{=} {\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\frac{{\overline{\mathit{p}}}^{\mathrm{2}}\sin \mathit{\varphi }}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}\frac{{\overline{\mathit{p}}}^{\mathrm{2}}\mathrm{+}\overline{\mathit{p}}\begin{array}{c}\mathrm{2}\\ \mathrm{0}\end{array}\mathrm{-}\mathrm{2}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\left(\sin \mathit{\varphi }\sin {\mathit{\varphi }}_{\mathrm{0}}\cos \mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{\theta }}_{\mathrm{0}}\mathrm{)}\mathrm{+}\cos \mathit{\varphi }\cos {\mathit{\varphi }}_{\mathrm{0}}\right)}{\mathrm{2}}}\mathrm{d}\overline{\mathit{p}}& & \\ \mathrm{=}\frac{\sin \mathit{\varphi }}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}\frac{\overline{\mathit{p}}\begin{array}{c}\mathrm{2}\\ \mathrm{0}\end{array}}{\mathrm{2}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\overline{\mathit{p}}}^{\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\frac{{\overline{\mathit{p}}}^{\mathrm{2}}\mathrm{-}\mathrm{2}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\left(\sin \mathit{\varphi }\sin {\mathit{\varphi }}_{\mathrm{0}}\cos \mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{\theta }}_{\mathrm{0}}\mathrm{)}\mathrm{+}\cos \mathit{\varphi }\cos {\mathit{\varphi }}_{\mathrm{0}}\right)}{\mathrm{2}}}\mathrm{d}\overline{\mathit{p}}& & \end{array}$and use the definite integral ${\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{x}}^{\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{x}}^{\mathrm{2}}\mathrm{+}\mathit{ax}}{\mathrm{2}}}\mathrm{d}\mathit{x}\mathrm{=}\sqrt{\frac{\mathit{\pi }}{\mathrm{32}}}\mathrm{(}{\mathit{a}}^{\mathrm{2}}\mathrm{+}\mathrm{4}\mathrm{)}{\mathrm{e}}^{{\mathit{a}}^{\mathrm{2}}\mathit{/}\mathrm{8}}\mathrm{erfc}\left(\frac{\mathit{a}}{\mathrm{2}\sqrt{\mathrm{2}}}\right)\mathrm{-}\frac{\mathit{a}}{\mathrm{2}}\mathit{,}$(13)where erfc(x) is the complementary error function defined as⁷$\mathrm{erfc}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{\equiv }\mathrm{1}\mathrm{-}\mathrm{erf}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}}{\sqrt{\mathit{\pi }}}{\mathrm{\int }}_{\mathit{x}}^{\mathrm{\infty }}{\mathrm{e}}^{\mathrm{-}{\mathit{t}}^{\mathrm{2}}}\mathrm{d}\mathit{t}$(14)and erf(x) is the error function to reduce this integral. $\begin{array}{ccc}\overline{\mathit{M}}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right|}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathit{,\sigma }\mathrm{)}\mathrm{=}\frac{\sin \mathit{\varphi }}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}\frac{\overline{\mathit{p}}\begin{array}{c}\mathrm{2}\\ \mathrm{0}\end{array}}{\mathrm{2}}}\left(\sqrt{\frac{\mathit{\pi }}{\mathrm{32}}}\mathrm{(}{\mathit{A}}^{\mathrm{2}}\mathrm{+}\mathrm{4}\mathrm{)}{\mathrm{e}}^{{\mathit{A}}^{\mathrm{2}}\mathit{/}\mathrm{8}}\mathrm{erfc}\left(\frac{\mathit{A}}{\mathrm{2}\sqrt{\mathrm{2}}}\right)\mathrm{-}\frac{\mathit{A}}{\mathrm{2}}\right)& & \end{array}$(15)is obtained where $\mathit{A}\mathrm{\equiv }\mathrm{-}\mathrm{2}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\left(\sin \mathit{\varphi }\sin {\mathit{\varphi }}_{\mathrm{0}}\cos \mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{\theta }}_{\mathrm{0}}\mathrm{)}\mathrm{+}\cos \mathit{\varphi }\cos {\mathit{\varphi }}_{\mathrm{0}}\right)$(16)which may be viewed as a higher-dimensional analog of the angular distribution presented in Vinokur (1965; see also Clarke & Stewart 1986; Naghizadeh-Khouei & Clarke 1993; Quinn 2012).

Examining a few density plots of Eq. (15) for various values of ${\overline{\mathit{p}}}_{\mathrm{0}}$, θ₀, and ϕ₀ reveals that the probability distribution is attenuated near the poles of the coordinate system, an effect that becomes broader for small values of ${\overline{\mathit{p}}}_{\mathrm{0}}$. In particular, one should notice that even for ${\overline{\mathit{p}}}_{\mathrm{0}}\mathrm{=}\mathrm{0}$, the probability density still varies with ϕ. This is due to the pole bias effect discussed earlier.

One would like to go further and find the two one-dimensional marginalized angular distributions (cases IV and V) but the integrations seem unable to be completely performed using elementary or even common special functions.

4.2. The angular-marginalized p-distribution

One is often primarily concerned with the degree of polarization, p, and not the value of the angular variables. Marginalization over the angular variables will therefore produce a new distribution, to be called Q ($\overline{Q}$ in signal-to-noise variables) of high practical importance. In the polar coordinate case, marginalization over the (lone) angular variable, θ, produces the Rice distribution. In the spherical coordinate case, the marginalization must occur over both θ and ϕ. Thus, the function $\overline{Q}$ may be considered a higher-dimensional analog of the Rice distribution. Let $Q\mathrm{\left(}\mathit{p}\mathrm{\right|}{\mathit{p}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathit{,\sigma }\mathrm{)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{-}\mathit{\pi }}^{\mathit{\pi }}{\mathit{f}}^{\mathrm{\prime }}\mathrm{(}\mathit{p,\theta ,\varphi }\mathrm{\left|}{\mathit{p}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathit{,\sigma }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\theta }\hspace{0.17em}\mathrm{d}\mathit{\varphi }$(17)or similarly $\overline{Q}\mathrm{\left(}\overline{\mathit{p}}\mathrm{\right|}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathit{,\sigma }\mathrm{)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{-}\mathit{\pi }}^{\mathit{\pi }}{\overline{\mathit{f}}}^{\mathrm{\prime }}\mathrm{(}\overline{\mathit{p}}\mathit{,\theta ,\varphi }\mathrm{\left|}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathit{,\sigma }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\theta }\hspace{0.17em}\mathrm{d}\mathit{\varphi }\mathit{.}$(18)Let’s calculate. $\begin{array}{ccc}& & \overline{Q}\mathrm{\left(}\overline{\mathit{p}}\mathrm{\right|}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\theta }}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathit{,\sigma }\mathrm{)}\mathrm{=}\frac{{\overline{\mathit{p}}}^{\mathrm{2}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}\frac{{\overline{\mathit{p}}}^{\mathrm{2}}\mathrm{+}\overline{\mathit{p}}\begin{array}{c}\mathrm{2}\\ \mathrm{0}\end{array}}{\mathrm{2}}}\\ & & \mathrm{\times }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{-}\mathit{\pi }}^{\mathit{\pi }}\sin \mathit{\varphi }\hspace{0.17em}{\mathrm{e}}^{\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\mathrm{\left(}\sin \mathit{\varphi }\sin {\mathit{\varphi }}_{\mathrm{0}}\cos \mathrm{\right(}\mathit{\theta }\mathrm{-}{\mathit{\theta }}_{\mathrm{0}}\mathrm{)}\mathrm{+}\cos \mathit{\varphi }\cos {\mathit{\varphi }}_{\mathrm{0}}\mathrm{)}}\hspace{0.17em}\mathrm{d}\mathit{\theta }\hspace{0.17em}\mathrm{d}\mathit{\varphi }\\ & & \mathrm{=}\frac{{\overline{\mathit{p}}}^{\mathrm{2}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}\frac{{\overline{\mathit{p}}}^{\mathrm{2}}\mathrm{+}\overline{\mathit{p}}\begin{array}{c}\mathrm{2}\\ \mathrm{0}\end{array}}{\mathrm{2}}}\\ & & \mathrm{\times }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}\sin \mathit{\varphi }\hspace{0.17em}{\mathrm{e}}^{\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\cos \mathit{\varphi }\cos {\mathit{\varphi }}_{\mathrm{0}}}{\mathrm{\int }}_{\mathrm{-}\mathit{\pi }}^{\mathit{\pi }}{\mathrm{e}}^{\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\sin \mathit{\varphi }\sin {\mathit{\varphi }}_{\mathrm{0}}\cos \mathrm{(}\mathit{\theta }\mathrm{-}{\mathit{\theta }}_{\mathrm{0}}\mathrm{)}}\hspace{0.17em}\mathrm{d}\mathit{\theta }\hspace{0.17em}\mathrm{d}\mathit{\varphi }\mathit{.}\end{array}$(19)The θ-integral can be solved using ${{}^{\mathrm{\int }}}_{\mathrm{-}\mathit{\pi }}^{\mathit{\pi }}{\mathrm{e}}^{\mathit{a}\cos \mathrm{\left(}\mathit{t}\mathrm{\right)}}\mathrm{d}\mathit{t}\mathrm{=}\mathrm{2}\mathit{\pi }{\mathrm{ℐ}}_{\mathrm{0}}\mathrm{\left(}\mathit{a}\mathrm{\right)}$, where ℐ₀(x) is the modified Bessel function of the first kind. ℐ₀ is an even, real function when its argument is real. (Because the range of integration is over a full period of θ, the value of θ₀ is immaterial because it appears as a phase shift and does not affect the integral.) Once the θ-integral is performed, the $\overline{Q}$ distribution no longer explicitly depends on θ₀. We shall therefore drop θ₀ from the conditions. Continuing with the calculation, $\begin{array}{ccc}& & \overline{Q}\mathrm{\left(}\overline{\mathit{p}}\mathrm{\right|}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathit{,\sigma }\mathrm{)}\mathrm{=}\frac{{\overline{\mathit{p}}}^{\mathrm{2}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}\frac{{\overline{\mathit{p}}}^{\mathrm{2}}\mathrm{+}\overline{\mathit{p}}\begin{array}{c}\mathrm{2}\\ \mathrm{0}\end{array}}{\mathrm{2}}}\\ & & \mathrm{\times }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}\sin \mathit{\varphi }\hspace{0.17em}{\mathrm{e}}^{\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\cos \mathit{\varphi }\cos {\mathit{\varphi }}_{\mathrm{0}}}\mathrm{\left(}\mathrm{2}\mathit{\pi }\mathrm{\right)}{\mathrm{ℐ}}_{\mathrm{0}}\mathrm{\left(}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\sin \mathit{\varphi }\sin {\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}\mathrm{d}\mathit{\varphi }\\ & & \mathrm{=}\frac{{\overline{\mathit{p}}}^{\mathrm{2}}}{\sqrt{\mathrm{2}\mathit{\pi }}}{\mathrm{e}}^{\mathrm{-}\frac{{\overline{\mathit{p}}}^{\mathrm{2}}\mathrm{+}\overline{\mathit{p}}\begin{array}{c}\mathrm{2}\\ \mathrm{0}\end{array}}{\mathrm{2}}}\\ & & \end{array}$(20)The ϕ-integral in the last expression is non-trivial but it appears it can be solved. Let $\mathit{\omega }\mathrm{\left(}\mathit{a,\varphi ,}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}\mathrm{\equiv }\sin \mathit{\varphi }\hspace{0.17em}{\mathrm{e}}^{\mathit{a}\cos \mathit{\varphi }\cos {\mathit{\varphi }}_{\mathrm{0}}}{\mathrm{ℐ}}_{\mathrm{0}}\mathrm{\left(}\mathit{a}\sin \mathit{\varphi }\sin {\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}$(21)and $\mathrm{\Omega }\mathrm{\left(}\mathit{a,}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}\mathrm{\equiv }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}\mathit{\omega }\mathrm{\left(}\mathit{a,\varphi ,}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}\mathrm{d}\mathit{\varphi ,}$(22)where $\mathit{a}\mathrm{\equiv }\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\mathrm{\ge }\mathrm{0}$. No computer algebra system tested could solve the integral in Eq. (22) and no applicable integrals were found in tables of integrals. Manual attempts to solve it also failed. From graphs of Ω versus ϕ₀ created numerically⁸ for several values of a, it is noticed that the function may actually be constant in ϕ₀ for a given a. Unfortunately, it is not a simple matter to show $\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial }{\mathit{\varphi }}_{\mathrm{0}}}\mathrm{=}\mathrm{0}$ but if $\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial }{\mathit{\varphi }}_{\mathrm{0}}}$ is graphed for 0 ≤ ϕ₀ ≤ π at several sampled values of a it seems to always show residuals about zero with scatter consistent with numerical noise. The previous numerical observations lead one to suspect very strongly that Ω(a,ϕ₀) is independent of the value of ϕ₀. This may be surprising in light of the bias corrections previously discussed; however, the same marginalized formula seems like it should be found regardless of the rotational orientation of the $\mathrm{\left(}\overline{\mathit{q}}\mathit{,}\overline{\mathit{u}}\mathit{,}\overline{\mathit{v}}\mathrm{\right)}$ coordinate system within the Poincaré sphere. If it is true that Ω(a,ϕ₀) is constant in ϕ₀, any value of ϕ₀ may be chosen and it does not alter Ω(a,ϕ₀)’s value at a given a. In that case, a value of ϕ₀ may be judiciously chosen that simplifies the integral. Good choices are ϕ₀ = 0 or π for which $\mathrm{\Omega }\mathrm{\left(}\mathit{a,}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\sinh \mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{a}}$(23)is found using the common integral ${{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{\pi }}{\mathrm{e}}^{\mathit{x}\cos \mathit{\varphi }}\sin \mathit{\varphi }\hspace{0.17em}\mathrm{d}\mathit{\varphi }\mathrm{=}\mathrm{2}\sinh \mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{/}\mathit{x}$, where sinh(x) is the hyperbolic sine of x. Indeed, simple arguments modifying the factors in the integrand of Eq. (22) can prove this to be a lower limit for Ω(a,ϕ₀) at a given a (see Appendix A). Proof that it is also the upper limit was elusive (except at ϕ₀ = 0, π/ 2, and π). A Taylor series expansion of Ω(a,ϕ₀) in ϕ₀ at ϕ₀ = 0 was found to be consistent with Eq. (23) out to 26-th order using computer algebra systems⁹. Based on physical suspicion corroborated by numerical results, Eq. (23) shall be assumed to be true henceforth.

Returning to the main distribution, $\overline{Q}$, and using Eq. (23), it is proposed that $\overline{Q}\mathrm{\left(}\overline{\mathit{p}}\mathrm{\right|}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,\sigma }\mathrm{)}\mathrm{=}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\frac{\overline{\mathit{p}}}{{\overline{\mathit{p}}}_{\mathrm{0}}}{\mathrm{e}}^{\mathrm{-}\frac{{\overline{\mathit{p}}}^{\mathrm{2}}\mathrm{+}\overline{\mathit{p}}\begin{array}{c}\mathrm{2}\\ \mathrm{0}\end{array}}{\mathrm{2}}}\sinh \mathrm{(}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\mathrm{)}\mathit{.}$(24)Since the value of ϕ₀ is irrelevant, $\overline{Q}$’s dependence on ϕ₀ has been dropped from the notation. This is the higher-dimensional analog of the Rice distribution. Surprisingly, it is a simpler function than the Rice distribution in the sense that it depends on the hyperbolic sine function rather than a modified Bessel function. The $\overline{Q}$ distribution is normalized, ${{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathrm{\infty }}\overline{Q}\hspace{0.17em}\mathrm{d}\overline{\mathit{p}}\mathrm{=}\mathrm{1}$.

Fig. 4 Contour plot of the angular-marginalized $\overline{\mathit{p}}$-distribution given by Eq. (24). The contours are at 0.1 intervals with the lighter shades corresponding to larger values. The dotted and dashed curves (Eqs. (25) and (26)) trace maximums along horizontal and vertical slices and end at the red points at $\mathrm{\left(}\sqrt{\mathrm{2}}\mathit{,}\mathrm{0}\mathrm{\right)}$ and $\mathrm{\left(}\sqrt{\mathrm{3}}\mathit{,}\mathrm{0}\mathrm{\right)}$, respectively. Although the “shape” of the distribution is not affected by the value of σ, the plot assumes that σ ≤ 1/4 so that the plotted values for ${\overline{\mathit{p}}}_{\mathrm{0}}$ up to 4 (=1 /σmax) exist.

Some general properties of $\overline{Q}$ can now be investigated. A contour plot of the function is shown in Fig. 4 with contours at 0.1 intervals. The function achieves a global maximum as $\mathrm{\left(}\overline{\mathit{p}}\mathit{,}{\overline{\mathit{p}}}_{\mathrm{0}}\mathrm{\right)}$ approaches $\mathrm{\left(}\sqrt{\mathrm{2}}\mathit{,}\mathrm{0}\mathrm{\right)}$ of $\frac{\mathrm{2}}{\mathit{e}}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\mathrm{\approx }\mathrm{0.587051}$. Tracing the maximums of horizontal slices (the dotted line in the figure) yields the equation $\mathrm{(}{\overline{\mathit{p}}}^{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}\sinh \mathrm{(}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\mathrm{)}\mathrm{-}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\cosh \mathrm{(}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\mathrm{)}\mathrm{=}\mathrm{0}$ or $\mathrm{(}\overline{\mathit{p}}\begin{array}{c}\mathrm{2}\\ \mathrm{0}\end{array}\mathrm{-}\mathrm{1}\mathrm{)}\tanh \mathrm{(}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\mathrm{)}\mathrm{=}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{.}$(25)This implicit curve intersects the $\overline{\mathit{p}}$-axis at the global maximum at $\mathrm{\left(}\sqrt{\mathrm{2}}\mathit{,}\mathrm{0}\mathrm{\right)}$ and is the “Most Probable” estimator curve for $\overline{Q}$. Tracing the maximums of vertical slices (dashed line) yields $\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\cosh \mathrm{(}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\mathrm{)}\mathrm{-}\mathrm{(}\overline{\mathit{p}}\begin{array}{c}\mathrm{2}\\ \mathrm{0}\end{array}\mathrm{+}\mathrm{1}\mathrm{)}\sinh \mathrm{(}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\mathrm{)}\mathrm{=}\mathrm{0}$ or $\mathrm{(}\overline{\mathit{p}}\begin{array}{c}\mathrm{2}\\ \mathrm{0}\end{array}\mathrm{+}\mathrm{1}\mathrm{)}\tanh \mathrm{(}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\mathrm{)}\mathrm{=}\overline{\mathit{p}}\hspace{0.17em}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{.}$(26)This curve intersects at the $\overline{\mathit{p}}$-axis at $\mathrm{\left(}\sqrt{\mathrm{3}}\mathit{,}\mathrm{0}\mathrm{\right)}$ and is the “Most Likely” estimator curve for $\overline{Q}$.

The MP estimator curve of Eq. (25) (dotted curve in Fig. 4) can be used for all $\overline{\mathit{p}}\mathit{>}\sqrt{\mathrm{2}}$. Yet it was stressed that the full MP solution given Eqs. (9) and (10) is only valid for some values with $\overline{\mathit{p}}\mathit{>}\sqrt{\mathrm{2}}$. Marginalization has thus hidden some of the complexity of the problem. The ML estimator given by Eq. (26) (dashed curve) now has non-trivial behavior unlike the full solution in Eq. (8). It is not easy to anticipate the effect that marginalization will have on various estimators.

It is interesting to compare the long form of these estimator curves and the corresponding curves based on the Rice distribution (Wardle & Kronberg 1974; Simmons & Stewart 1985; Quinn 2012). These curves are based on the hyperbolic trigonometric functions whereas the curves in the two-dimensional case are based on Bessel functions. Their points of intersection with the $\overline{\mathit{p}}$-axis have also been shifted. The differences arise because the Jacobian factor is proportional to p² in spherical coordinates but only p for polar coordinates.

5. Discussion

How do these new results relate to the old results derived using polar coordinates for the “Poincaré disk”? First, the old results are not special cases of the new results when ϕ → π/ 2 and ϕ₀ → π/ 2; for instance, in the special case of ϕ = π/ 2, the $\overline{\mathit{p}}$ equation of Eq. (9) reduces to $\overline{\mathit{p}}\mathrm{-}\mathrm{2}\mathit{/}\overline{\mathit{p}}$ rather than the Wang estimator $\overline{\mathit{p}}\mathrm{-}\mathrm{1}\mathit{/}\overline{\mathit{p}}$ (Wang et al. 1997; Quinn 2012) and $\overline{\mathit{M}}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}$ does not reduce to angular distribution given by Vinokur (1965). Similarly, $\overline{Q}$ is altogether distinct from the Rice distribution. The new results are therefore not direct extensions of the old functions but they do have somewhat similar forms. The three versus two dimensional nature of the parameter space is responsible for this.

Particularly interesting is a comparison of Fig. 4 with Fig. 2a of Quinn (2012). The intercepts of all the estimator curves have undergone shifts of the form $\sqrt{\mathit{n}}\mathrm{\to }\sqrt{\mathit{n}\mathrm{+}\mathrm{1}}$. This again is due to the increased dimensionality of the problem. Thus, when circular polarization has equal measurement error to the linear Stokes axes’ measurement error, it takes a bigger measurement of the magnitude of the polarization vector to be significant than the two-dimensional case suggests. The two-dimensional results can perhaps be recovered by some limiting procedure where σ_v → 0 while σ (=σ_q = σ_u) is held constant. This suggests that decreasing the measurement error on one of σ_q, σ_u, or σ_v should decrease the $\overline{\mathit{p}}$-threshold for a significant measurement. On the other hand, often an instrument measures only linear or only circular polarization. This may correspond to a scenario where the distribution of Stokes parameters is best modeled by assuming a flat (instead of gaussian) distribution for the variables of the undetected type. This would violate our starting assumptions in Eq. (1) and would produce rather different results. As a flat distribution is in an informal sense a “wide” distribution, it is probable that this would raise the detection threshold higher. It should be clear that interpreting polarization data at low signal-to-noise in “human digestible” quantities like degrees of polarization or angle of polarization is difficult.

6. Conclusion

The three-dimensional Stokes sampling distribution in spherical coordinates is given. From it the “Most Likely” and “Most Probable” classical (that is, non-Bayesian) estimators for the “true” Stokes parameters given measured values have been derived. Additionally, a two useful marginalizations of the sampling distribution have been calculated including a higher dimensional analog of the Rice distribution. These results are necessary stepping stones towards a full and proper treatment of polarization measurements.

It is cautioned that a deep understanding of the measurement of Stokes parameters in spherical coordinates like here or the polar coordinates usually used for linear polarization requires a Bayesian approach. The presented results may nonetheless have tolerable accuracy for many experiments.

While the setup of the problem and the presentation of the results has been framed in terms of the study of polarization, the solutions are in fact far more general. They are applicable to the measurement of any “source point” contained in a three-dimensional ball measured with Gaussian error in each Cartesian direction and presented in spherical coordinates. The results could be used to study the location of neutrino production within the Sun or the determination of the epicenter of an earthquake. The results could also influence the interpretation of cosmic microwave background polarization data. Extensions of the work removing the constraint that the error be equal in each direction or that involve correlated variables would be highly desirable.

References

Appendix A: Bounds on Ω(a, ϕ₀)

It helps to have functional bounds on $\mathrm{\Omega }\mathrm{\left(}\mathit{a,}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}\mathrm{\equiv }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}\sin \mathrm{\left(}\mathit{\varphi }\mathrm{\right)}{\mathrm{e}}^{\mathit{a}\cos \mathrm{\left(}\mathit{\varphi }\mathrm{\right)}\cos \mathrm{\left(}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}}{\mathit{I}}_{\mathrm{0}}\mathrm{\left(}\mathit{a}\sin \mathrm{\right(}\mathit{\varphi }\mathrm{\left)}\sin \mathrm{\right(}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\varphi }$(A.1)to show, for instance, that it is at least non-infinite over the ϕ₀-domain of interest. There are three factors (sin(ϕ), e^{acos(ϕ)cos(ϕ₀)}, and I₀(asin(ϕ)sin(ϕ₀))) in the integrand, each of which is simple enough that they allow various bounds to be obtained for the total expression.

Appendix A.1: Lower bound

Under variation of ϕ₀, the e^{acos(ϕ)cos(ϕ₀)} factor is minimized when cos(ϕ₀) is the smallest, which occurs for ϕ₀ = π, so that the factor becomes e^−acos(ϕ). Furthermore, the I₀(asin(ϕ)sin(ϕ₀)) factor is the smallest when its argument is zero, which also occurs at ϕ₀ = π (also ϕ₀ = 0) so that the factor become ℐ₀(0) = 1. After these substitutions, the integral for Ω can be performed explicitly: ${\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}\sin \mathrm{\left(}\mathit{\varphi }\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathit{a}\cos \mathrm{\left(}\mathit{\varphi }\mathrm{\right)}}\mathrm{d}\mathit{\varphi }\mathrm{=}\frac{\mathrm{2}\sinh \mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{a}}\mathrm{·}$(A.2)For any given value of a, this sets a (constant) global lower bound on the integral: the value of integral must be greater than or equal to $\frac{\mathrm{2}\sinh \mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{a}}$ for all values of ϕ₀. It is probably also the greatest lower bound for each ϕ₀.

Appendix A.2: Upper bounds

The e^{acos(ϕ)cos(ϕ₀)} factor is maximized when its exponent is the largest while the I₀(asin(ϕ)sin(ϕ₀)) factor is the largest when the magnitude of its argument is as big as possible.

A constant global upper bound is found by modifying the factors in the obvious way: ${\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}\sin \mathrm{\left(}\mathit{\varphi }\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathit{a}\cos \mathrm{\left(}\mathit{\varphi }\mathrm{\right)}}{\mathrm{ℐ}}_{\mathrm{0}}\mathrm{(}\mathrm{\pm }\mathit{a}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{\varphi }\mathrm{=}\frac{\mathrm{2}\sinh \mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{a}}{\mathrm{ℐ}}_{\mathrm{0}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathit{.}$(A.3)A stricter (but non-constant) global upper bound can be found merely by allowing sin(ϕ₀) → 1 in the Bessel function factor, $\begin{array}{ccc}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}\sin \mathrm{\left(}\mathit{\varphi }\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathit{a}\cos \mathrm{\left(}\mathit{\varphi }\mathrm{\right)}\cos \mathrm{\left(}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}}{\mathrm{ℐ}}_{\mathrm{0}}\mathrm{\left(}\mathit{a}\sin \mathrm{\right(}\mathit{\varphi }\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\varphi }\mathrm{=}\frac{\mathrm{2}\sec \mathrm{\left(}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}\sinh \mathrm{\left(}\mathit{a}\cos \mathrm{\right(}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\left)}\mathrm{\right)}}{\mathit{a}}{\mathrm{ℐ}}_{\mathrm{0}}\mathrm{\left(}\mathit{a}\mathrm{\right)} \mathrm{for}{\mathit{\varphi }}_{\mathrm{0}}\ne \mathit{\pi }\mathit{/}\mathrm{2.}& & \end{array}$(A.4)This last solution has functional dependence on ϕ₀ and is not valid for ϕ₀ = π/ 2 because of the sec(ϕ₀) component. For ϕ₀ = 0 this bound is $\frac{\mathrm{2}\sinh \mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{a}}{\mathrm{ℐ}}_{\mathrm{0}}\mathrm{\left(}\mathit{a}\mathrm{\right)}$ and it decreases smoothly and approaches $\frac{\mathrm{2}\sinh \mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{a}}$ as ϕ₀ → π/ 2⁺. Similarly, for ϕ₀ = π it is $\frac{\mathrm{2}\sinh \mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{a}}{\mathrm{ℐ}}_{\mathrm{0}}\mathrm{\left(}\mathit{a}\mathrm{\right)}$ and it decreases smoothly and approaches $\frac{\mathrm{2}\sinh \mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{a}}$ as ϕ₀ → π/ 2⁻. Of course, if ϕ₀ = π/ 2, Eq. (A.1) can be accomplished directly: ${\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}\sin \mathrm{\left(}\mathit{\varphi }\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{0}}{\mathrm{ℐ}}_{\mathrm{0}}\mathrm{\left(}\mathit{a}\sin \mathrm{\right(}\mathit{\varphi }\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\varphi }\mathrm{=}\frac{\mathrm{2}\sinh \mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{a}} \mathrm{for}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{=}\mathit{\pi }\mathit{/}\mathrm{2}\mathit{,}$(A.5)so the discontinuity in Eq. (A.4) at ϕ₀ = π/ 2 can be “removed”. Equation (A.1) can also be performed for ϕ₀ = 0 and ϕ₀ = π and in both cases $\frac{\mathrm{2}\sinh \mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{a}}$ results, so we have explicit upper bounds at those values too.

Appendix A.3: Best bounds

Combining the previous information, the best bounds are $\frac{\mathrm{2}\sinh \mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{a}}\mathrm{\le }\mathrm{\Omega }\mathrm{\left(}\mathit{a,}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}\mathrm{\le }\{\begin{array}{c}\\ & \mathrm{if}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{=}\mathrm{0}\mathit{,}\frac{\mathit{\pi }}{\mathrm{2}}\mathit{,\pi }\\ & \mathrm{otherwise,}\end{array}$(A.6)which are finite for all ϕ₀ at a given value of a.

All Figures

	Fig. 1 Illustration of the spherical coordinate system used in the paper. The variable θ is used for the azimuthal angle with a range of (−π,π ], and the variable ϕ is measured from the positive v-axis with a range of [ 0,π ].
In the text

Fig. 2 Visualization of the $\mathrm{\left(}\overline{\mathit{p}}\mathit{,\varphi }\mathrm{\right)}$ bias correction field of Eq. (9). The tail of each arrow is located at a measured $\mathrm{\left(}\overline{\mathit{p}}\mathit{,\varphi }\mathrm{\right)}$ point and the head is attached to the corresponding “Most Probable” estimate of $\mathrm{\left(}{\overline{\mathit{p}}}_{\mathrm{0}}\mathit{,}{\mathit{\varphi }}_{\mathrm{0}}\mathrm{\right)}$. This correction should only be applied to $\mathrm{\left(}\overline{\mathit{p}}\mathit{,\varphi }\mathrm{\right)}$ points lying in the blue region given by Eq. (10). See the text for discussion about $\mathrm{\left(}\overline{\mathit{p}}\mathit{,\varphi }\mathrm{\right)}$ points outside the blue region.

In the text

	Fig. 3 Behavior of the bias correction field (Eq. (9)) on the boundary of the blue valid region (Eq. (10)). The left-most point on the boundary of the blue region is at $\mathrm{(}\sqrt{\mathrm{2}}\mathit{,\pi }\mathit{/}\mathrm{2}\mathrm{)}$.
In the text

Fig. 4 Contour plot of the angular-marginalized $\overline{\mathit{p}}$-distribution given by Eq. (24). The contours are at 0.1 intervals with the lighter shades corresponding to larger values. The dotted and dashed curves (Eqs. (25) and (26)) trace maximums along horizontal and vertical slices and end at the red points at $\mathrm{\left(}\sqrt{\mathrm{2}}\mathit{,}\mathrm{0}\mathrm{\right)}$ and $\mathrm{\left(}\sqrt{\mathrm{3}}\mathit{,}\mathrm{0}\mathrm{\right)}$, respectively. Although the “shape” of the distribution is not affected by the value of σ, the plot assumes that σ ≤ 1/4 so that the plotted values for ${\overline{\mathit{p}}}_{\mathrm{0}}$ up to 4 (=1 /σmax) exist.

In the text