Pitfalls of statistics-limited X-ray polarization analysis
Victor Mikhalev




KTH Royal Institute of Technology, Department of Physics, 106 91 Stockholm, Sweden.

The Oskar Klein Centre for Cosmoparticle Physics, AlbaNova University Centre, 106 91, Stockholm, Sweden.

mikhalev@kth.se


One of the difficulties with performing polarization analysis is that the mean polarization fraction of sub-divided data sets is larger than the polarization fraction for the integrated measurement. The resulting bias is one of the properties of the generating distribution discussed in this work. The limitations of Gaussian approximations in standard analysis based on Stokes parameters for estimating polarization parameters and their uncertainties are explored by comparing with a Bayesian analysis. The effect of uncertainty on the modulation factor is also shown, since it can have a large impact on the performance of gamma-ray burst polarimeters. Results are related to the minimum detectable polarization (MDP), a common figure of merit, making them easily applicable to any X-ray polarimeter.
The aim of this work is to quantify the systematic errors induced on polarization parameters and their uncertainties when using Gaussian approximations and to show when such effects are non-negligible.
The probability density function is used to deduce the properties of reconstructed polarization parameters. The reconstructed polarization parameters are used as sufficient statistics for finding a simple form of the likelihood. Bayes theorem is used to derive the posterior and to include nuisance parameters.
The systematic errors originating from Gaussian approximations as a function of instrument sensitivity are quantified here. Different signal-to-background scenarios are considered making the analysis relevant for a large variety of observations. Additionally, the change of posterior shape and instrument performance (MDP) due to uncertainties on the polarimeteric response of the instrument is shown.



Pitfalls of statistics-limited X-ray polarization analysis
    V. Mikhalev
    Accepted 23/03/2018
==========================================================





§ INTRODUCTION


The first observations of astrophysical X-ray polarization were made more than forty years ago (Novick et al. <cit.>; Weisskopf et al. <cit.>). The field has been reinvigorated in the past decade by a series of measurements by satellite, including INTEGRAL/SPI (Dean et al. <cit.>; Chauvin et al. <cit.>), INTEGRAL/IBIS (Forot et al. <cit.>; Moran et al. <cit.>), AstroSat/CZTI (Vadawale et al. <cit.>) and IKAROS/GAP (Yonetoku et al.  <cit.>) as well as balloon-borne instruments (PoGOLite (Chauvin et al. <cit.>) and PoGO+ (Chauvin et al. <cit.>)). Results from several on-going missions are expected in the near future (POLAR (Produit et al. <cit.>), X-Calibur (Beilicke et al. <cit.>)) and a dedicated satellite mission is in development (IXPE (Weisskopf et al. <cit.>)). Although some instruments, for example, IXPE, will be able to make polarization measurements of astrophysical sources with high precision, where the statistical analysis becomes relatively trivial, it is always desirable to observe weaker sources or to sub-divide the data. Fine splitting of data may be necessary for understanding physical phenomena, for example, gamma-ray bursts (GRBs) and pulsars require temporal binning, nebulae require spatial binning, and spectral binning is interesting for all objects. It is therefore important to know when the number of photons is sufficient for making an accurate analysis using simple methods and when a more rigorous approach is necessary.

The parameters describing linear polarization are the polarization fraction and the polarization angle. In the frequentist approach, the major challenge in estimating the polarization fraction is that it is a positive definite quantity, meaning that a non-zero fraction is measured even for an unpolarized source, thus introducing a bias. Ways of correcting for this bias have been studied by Simmons & Stewart (<cit.>) and more recently by Maier et al. (<cit.>). A Bayesian approach was first introduced by Vaillancourt (<cit.>) and extensively expanded upon by Quinn (<cit.>), where the shapes of the resulting parameter distributions are described. These preceding works focus on optical measurements of polarization for which some formulae differ from the X-ray counterpart (Kislat et al. <cit.>).

This work quantifies, as a function of measurement sensitivity, the error incurred when using the conventional Stokes parameter analysis. This is shown for both the polarization fraction and angle as well as for their uncertainties. Henceforth the expression “low statistics” is used to indicate poor data quality having a low S/√(N), where S is the number of signal photons and N is the total number of photons (signal and background). Since a point-source polarimeter working in the low-statistics regime is likely to have a low signal-to-background ratio R, whereas a GRB polarimeter is expected to have a higher R (especially if the GRB is bright and short in duration), the study is conducted for different values of R where meaningful. Additionally, the effects of uncertainty on the modulation factor μ_0, defined as the polarization fraction measured for a 100% polarized beam, on the performance of the polarimeter are studied. The uncertainty is typically large for GRB polarimeters which are not optimized for localization of GRBs since μ_0 varies with the photon incidence angle.






§ PARAMETER ESTIMATION


For a Compton scattering (Lei et al. <cit.>) or photo-electric polarimeters (Bellazzini & Muleri <cit.>) with a uniform response, the conditional distribution for a measured polarization angle ψ_i given a polarization fraction p_0, a polarization angle ψ_0 and a modulation factor μ_0 follows


    f(ψ_i|p_0,ψ_0,μ_0)=1/2π(S+B)
    ×(S×(1+μ_0 p_0cos(2(ψ-ψ_0)))+B)
    
    ,

where ψ=ϕ-90^∘ for the Compton process and ψ=ϕ for the photo-electric effect, ϕ is the measured scattering angle, S is the number of signal photons and B is the number of background photons. Treating signal and background as latent variables is beyond the scope of this paper so it is assumed that the background is unpolarized and is known with much higher precision than the polarization parameters. The modulation factor μ_0 is a property of the polarimeter and is derived from calibration. In what follows, the notation with subscript "0" means a physical parameter generating the data, in this case a set of measured scattering angles transformed to the polarization frame {ψ_i}, while subscript "r" denotes reconstructed parameters from this set. The reconstructed parameters are actually sufficient statistics for the data set allowing it to be represented by 2 scalars rather than a set of N angles.

The two most common ways of computing p_r and ψ_r are by performing a χ^2-fit of Eq. <ref> to a histogram of scattering angles or by computing the Stokes parameters. This work only considers the latter as it avoids complicating the analytical form of the likelihood due to binning effects.



 §.§ Stokes parameters


Polarimeters operating in the radio, infra-red, or optical domain measure photon intensities rather than individual photons as opposed to X-ray polarimeters. Clarke et al. (<cit.>) discuss how the Stokes parameter distributions vary depending on the instrumental technique used for measuring these intensities. Since Eq. <ref> is not used in the low-energy domain, not all results presented in that publication can be extrapolated to the X-ray energy band. In particular, the standard deviations and the correlation coefficient of Stokes parameters are affected.

In the X-ray energy band, the Stokes parameters are derived from individual photon events comprising two quantities


    q_i=cos(2ψ_i)
    
    u_i=sin(2ψ_i)
    
    .

For a total number of photons N=S+B the normalized Stokes parameters are written as

    Q_r=1/S∑_i=1^Nq_i
    
    U_r=1/S∑_i=1^Nu_i
    
    .

The normalization is only proportional to the signal because the background is assumed to be unpolarized. Therefore the background contributes only to the variance of (Q_r,U_r).

The Central Limit Theorem (CLT) makes Q_r and U_r Gaussian distributed as long as N is sufficiently large. This is referred to here as the CLT approximation. Thus Q_r and U_r follow the bivariate Gaussian distribution

    B(Q_r,U_r|Q_0,U_0)=1/2πσ_Qσ_U√(1-ρ^2)
    ×exp[-1/2(1-ρ^2)((Q_r-Q_0)^2/σ_Q^2+(U_r-U_0)^2/σ_U^2
    
    -2ρ(Q_r-Q_0)(U_r-U_0)/σ_Qσ_U)]
    
    .

Since the mean, the standard deviation, and the correlation coefficient are sufficient statistics for Gaussian distributed data, so are Q_r and U_r together with the second moments of the data derived in Appendix <ref>. Although not written explicitly in the conditioning, Q_0, U_0, σ_U, σ_Q and ρ are functions of p_0 and ψ_0.



 §.§ Polar coordinates

Stokes parameters (Q,U) are in Cartesian coordinates but can be transformed to polar coordinates to allow the polarization parameters to be expressed in terms of (p,ψ). The sufficient statistics become

    p_r=2√(Q_r^2+U_r^2)/μ_r
    ,


    ψ_r=1/2arctanU_r/Q_r
    ,

as per the derivations in Appendix <ref>.

Now Eq. <ref> can be transformed to polar coordinates by using the general coordinate transformation for a probability density function (p.d.f.) yielding the likelihood

    L(p_r,ψ_r|p_0,ψ_0)=B(Q_r,U_r|Q_0,U_0)×| (J)|
    
    ,

where (J) is the determinant of the Jacobian given by

    ∂ Q_r/∂ p_r   ∂ Q_r/∂ψ_r
    ∂ U_r/∂ p_r   ∂ U_r/∂ψ_r=
    μ_r/2cos(2ψ_r)     -μ_rp_rsin(2ψ_r) 
    μ_r/2sin(2ψ_r)    μ_rp_rcos(2ψ_r)
    =p_rμ_r^2/2
    .


Although the polarization parameters have a very similar form to that of the sufficient statistics (simply replacing the subscript "r" by "0")

    p_0=2√(Q_0^2+U_0^2)/μ_0,


    ψ_0=1/2arctanU_0/Q_0
    ,

p_r does not correspond to the most probable estimate of p_0. This occurs because (σ_Q,σ_U,ρ) depend on (p_0,ψ_0) and the fact that 0≤ p_0 ≤ 1 (since a polarization greater than 100% is unphysical) but there is no upper limit on p_r, that is, 0≤ p_r. Therefore, using the sufficient statistic p_r as the estimator of the polarization fraction p_0 incurs an error.

The error is non-neglible for the one-dimensional likelihood of p_0 (obtained after marginalizing over ψ_0) if the statistical significance of the measurement is low. It results in argmax_p_0 L(p_r,ψ_r|p_0)≠ p_r and p_r≠ p_0, where p_r is the expected value of p_r. Hence p_r, is neither the maximum likelihood nor an unbiased estimator of p_0. This is the case even if μ_r=μ_0 which occurs when there is no uncertainty on μ as is assumed throughout this section.

As the statistical precision of a measurement improves lim_S →∞p_r=p_0 and lim_S →∞argmax_p_0 L(p_r,ψ_r|p_0)= p_r. Conversely, due to symmetry, there is no bias on ψ so ψ_r=ψ_0 and argmax_ψ_0 L(p_r,ψ_r|ψ_0)= ψ_r.




 §.§ Central Limit Theorem approximation of the likelihood


The CLT approximation in B(Q_r,U_r|Q_0,U_0) provides an always easily computable analytical form of the likelihood given by Eq. <ref>. However, if the polarization fraction is high and the number of photons is low, the likelihood is not well described by such an approximation. The full form of the likelihood is given by a product over Eq. <ref>.

    L({ψ_i}|p_0,ψ_0,μ_0)=∏_i=1^N1/2π(S+B)
    ×(S×(1+μ_0 p_0cos(2(ψ_i-ψ_0)))+B)
    
    .

It has N factors and is therefore cumbersome to evaluate when the number of photons is large. This is unlike optical polarimetery where photon intensities of (Q,U) are directly measured instead of individual scattering angles.

Figure <ref> shows the difference between the CLT approximation and L(p_r|p_0,ψ_0,μ_0) as calculated by generating a data set {ψ_i} from Eq. <ref>, computing p_r using the Stokes parameters, repeating for many iterations and making a normalized histogram of p_r. 




In what follows, the concept of minimum detectable polarization (MDP) at 99% confidence level is important. The MDP (Weisskopf et al. <cit.>) is given by 

    MDP=4.29/μ_0 S√(S+B)
    .

Its statistical meaning is that, given an unpolarized source (p_0=0), the probability of measuring p_r>MDP is 1%. This quantity is a standard figure of merit for polarimeter performance and can easily be calculated for any measurement. The number of signal photons in Fig. <ref> is chosen such that MDP=2, since this is the highest ("worst") MDP considered in later sections. As seen from Fig. <ref>, the CLT approximation becomes more accurate as the number of photons increases and the polarization fraction p_0 decreases. In particular, only measurements with high signal-to-background ratio and high polarization fraction require the full likelihood given by Eq. <ref>. Although the error incurred using the CLT approximation is negligible in most cases, it becomes larger when using Eq. <ref> to derive the posterior, as is discussed below. Figure <ref> is for qualitative purposes only; the equivalent figure for ψ_r has been omitted since it leads to the same conclusions.



 §.§ Magnitude and importance of bias


Bias is a frequentist concept which relies on fixing (p_0,ψ_0) and investigating p_r. This approach provides an intuitive understanding for how an unpolarized source can produce a polarized signal (p_r>0). Several previous measurements use the sufficient statistic p_r as an estimate of the polarization fraction p_0 (e.g., Weisskopf et al. <cit.>; Slowikowska et al. <cit.>) and it is therefore necessary to understand when the bias is negligible and when a more sophisticated approach is necessary.

In the frequentist approach, the p.d.f. B(Q_r,U_r|Q_0,U_0) is not a function of ψ_0 (it is a fixed parameter) but of ψ_r. It is assumed, without loss of generality due to the angular symmetry of the problem, that ψ_0=0 so that ρ=0. This can be thought of as rotating the angular coordinate system which does not have any special reference point. The definition of the "zero angle" of such a system has no influence on the polarization fraction. The p.d.f. now simplifies to the product of two normal distributions

    B(Q_r,U_r|Q_0,U_0)_ψ_0=0=1/2πσ_Qσ_U
    ×exp[-1/2((Q-Q)^2/σ_Q^2+(U-U)^2/σ_U^2)]
    
    .

It can now be transformed to polar coordinates similarly to Eq. <ref> yielding

    f(p_r,ψ_r|p_0)_ψ_0=0=B(Q,U)_ψ_0=0| (J)|
    
    .

The relative mean bias β is now given by

    β=p_r-p_0/p_0=∫_0^∞∫_0^π p_r f(p_r,Δψ|p_0)_ψ_0=0dΔψdp_r-p_0/p_0
    ,

where Δψ=ψ_r-ψ_0=ψ_r. 

To understand which parameters have a significant impact on β an approximate analytical expression can be derived by introducing

    σ≡√(2N)/μ_0 S≈2σ_Q/μ_0≈2σ_U/μ_0
    .

It is now possible to write Eq. <ref> as

    f(p_r,Δψ|p_0,0)=p_r/πσ^2
    ×exp(-p_r^2+p_0^2-2p_rp_0cos(2Δψ)/2σ^2)
    .

Integrating over Δψ results in

    f(p_r|p_0)=p_r/σ^2exp(-p_r^2+p_0^2/2σ^2)× I_0(p_rp_0/σ^2)
    
    ,

which is the Rice distribution where I_0 is the modified Bessel function of zeroth order. In the limit of high statistics, the relative mean bias is given by

    lim_p_0/σ→∞β≈σ^2/2p_0^2=N/S^2μ_0^2p_0^2=(MDP/4.29p_0)^2
    
    ,

as shown in Appendix <ref>.

Since p_0 is not known a priori, Eq. <ref> needs to be expressed as a function of p_r by recursion. After some simplification the result is

    β≈1-2x^2-√(1-4x^2)/2x^2
    ,

where x≡MDP/4.29p_r. This shows that MDP/p_r is a good choice of independent variable. Equation <ref> is shown in Fig. <ref> where it is seen that β>0 and increases monotonously, i.e., on average the reconstructed polarization p_r will always be greater than the fixed polarization p_0. Exact numerical integration of Eq. <ref> is also provided for different signal-to-background ratios R yielding similar results to the approximation in Eq. <ref>. Here R=0 is the limit of large N, low S and yet finite MDP. To avoid inaccurately computing β for small S (a problem under the CLT approximation) Eq. <ref> is used for calculating L(p_r|p_0,ψ_0,μ_0) and ultimately its mean when S<200.




To understand when the bias is significant with respect to the statistical uncertainty on p_r, the bias fraction (p_r-p_0)/σ_p_r=β× p_0/σ_p_r is shown in Fig. <ref>. Here σ_p_r is derived (see Appendix <ref>) using standard error propagation yielding

    σ_p_r=2/μ_r√(1/S(N/2S-μ_0^2p_0^2/4))
    .

It becomes clear that this bias cannot be ignored in the low statistics regime. For a measured polarization below MDP, i.e., when MDP/p_r>1, the bias is more than 18% of the statistical uncertainty. Hence, one should not use Eq. <ref> for estimating p_0 in this regime.




A conclusion of this analysis is that significant errors are incurred when dividing the data into several data sets as is done for example by Dean et al. (<cit.>) in order to estimate the statistical uncertainty. The smaller the data set, the bigger the bias, and thus, on average, the result of an integrated measurement will be lower than the mean of its constituent data sets. This is also relevant when fitting models to polarization fraction sub-divided with respect to energy or time. Binning the data will result in a higher reconstructed polarization fraction, thus biasing the fit and therefore physical conclusions should not be drawn from p_r, as is done for example in Vadawale et al. (<cit.>) where the evolution of polarization parameters throughout the Crab pulsar pulse phase is investigated.





§ MAXIMUM A POSTERIORI ESTIMATE


The previous section showed the bias arising when using Eq. <ref> for estimating p_0 independently of the polarization angle. It is not obvious if such a bias exists when using (p_r,ψ_r) for the joint point estimate of (p_0,ψ_0). However, in this case it does not make sense to find a p_r,ψ_r for a fixed (p_0,ψ_0) so here the word "bias" is interpreted as the difference between the Maximum A Posteriori (MAP) estimate, p_MAP, and p_r. Such an analysis requires a Bayesian approach.



 §.§ Central Limit Theorem approximation of the posterior


The posterior P(p_0,Δψ_0|p_r) is derived by applying the Bayes theorem to the likelihood given by Eq. <ref>

    P(p_0,Δψ_0|p_r)=𝒩× P(p_0,Δψ_0) × L(p_r,0|p_0,Δψ_0)
    
    ,

where 𝒩 is the normalization factor, P(p_0,ψ_0) is the prior and Δψ_0=ψ_0-ψ_r. Since in a Bayesian approach the parameters (p_0,ψ_0) vary and the data (p_r,ψ_r) are fixed, L(p_r,0|p_0,Δψ_0) is a function of ψ_0. It follows that σ_Q≠σ_U and ρ≠0. This differs from the typical assumption of σ_Q=σ_U and no correlation between Q and U as done in other works (Simmons & Stewart <cit.>; Vaillancourt <cit.>; Quinn <cit.>; Maier et al. <cit.>).

The Jeffreys prior has the desirable property of being invariant under re-parametrization (Jeffreys <cit.>). In this case it is a uniform prior in the Cartesian coordinates (Q_0,U_0). It will result in a preference of high polarization fractions after transforming to polar coordinates, as discussed in Quinn (<cit.>). This is unphysical because it is more difficult to make a highly polarized photon beam in nature since any disorder in the emission region will lower the polarization fraction. A more realistic approach is to instead take a uniform prior 0≤ p_0 ≤ 1 in polar coordinates. The prior for ψ_0 is also uniform due to symmetry.

Figure <ref> shows the posterior as given by Eq. <ref> for a low-statistics measurement MDP/p_r=1. The asymmetry is apparent as the posterior broadens for low values of p_0. The effect of the prior is illustrated by including the unbound prior scenario where p_0≥0.



Since the posterior is derived using a likelihood which utilizes the CLT approximation, the posterior may, for certain combinations of parameters, be inaccurate. To check the magnitude of the effect, pairs of (p_0,ψ_0) are sampled from the prior and then used to generate a data set. Data sets falling within a narrow window of a chosen p_r and ψ_r (p_r±0.002 and ψ_r±0.36^∘) are selected and the posterior for each data set is found using Eq. <ref> as the likelihood. Finally, the intervals containing 90% of the posteriors marginalized over the polarization angle (for clarity) are shown in Fig. <ref> along with the CLT approximation. The situation is similar for the posterior of the polarization angle marginalized over the polarization fraction.



The most important feature of Fig. <ref> is that the posterior is not the same for fixed (p_r,ψ_r) when the number of signal photons is low. Hence (p_r,ψ_r) are not always sufficient statistics and the full data set {ψ_i} is required for deriving the posterior. However, the posterior converges quickly towards the CLT approximation when R<0.5 since any meaningful measurement with such a low R requires a large number of photons. An effect not shown in the figure is that the difference between the full posterior and the CLT approximation increases as p_r and μ_r increase requiring more photons for good convergence. However, few X-ray polarimeters have μ_0 higher than 0.5 (Krawczynski et al. <cit.>; Kaaret  <cit.>) and synchrotron emission, a proposed mechanism for many astrophysical sources, is not expected to produce a polarization fraction above ∼0.6 (Lyutikov et al. <cit.>).

In conclusion, the full likelihood in Eq. <ref> rather than the CLT approximation in Eq. <ref> is required for sources where the signal-to-background ratio is high but the total number of photons is low (e.g., short-duration GRBs). A guideline is to use the CLT approximation only for data sets where S>10^3. As mentioned before, polarimeters foreseen to be active in the near future (Astrosat/CZTI, POLAR, IXPE, X-Calibur) have μ_0<0.5 so this guideline translates to the CLT approximation being generally valid for .



 §.§ Relative MAP bias


The MAP estimate is at p__MAP> p_r and ψ_0=ψ_r. The inequality occurs because σ_Q and σ_U depend on p_0 but only measurements with high R and p_0^2μ_0^2≫0 are significantly affected. In all other cases σ_Q≈σ_U≈√(2N)/μ_0 S and p__MAP≈ p_r. 

The relative MAP bias is defined as

    β__MAP≡p_r-p__MAP/p__MAP
    .

Table <ref> shows β__MAP for the extreme case MDP/p_r=2 (a monotonously decreasing function) for different polarization parameters. Only scenarios with a sufficient number (S>884 based on the results shown in Fig. <ref>) of signal photons are considered. The results show that β__MAP increases with p_r and μ_0 but its value is negligible for all cases where the CLT approximation is valid. Therefore, when (p_r,ψ_r) are sufficient statistics of the data, they are a good estimate for the MAP.



The broadening of the joint posterior at low values of p_0, as shown in Fig. <ref>, results in p__MAP being a poor estimate for p_0 when marginalizing over ψ_0. More generally, this occurs for any two-dimensional distribution when marginalizing over the nuisance parameter if the estimated parameter has an asymmetric distribution.






§ UNCERTAINTIES ON POLARIZATION PARAMETERS


There are two reasons why it is not always appropriate to use naive Stokes estimates for uncertainties of polarization parameters: non-Gaussianity and the prior. Eq. <ref> gives the correct uncertainty on the polarization fraction when statistics are high and the reconstructed fraction p_r is well below 1. If these conditions are not fulfilled, the marginalized posterior becomes highly asymmetric either due to low statistics or because a part of the likelihood is truncated by the prior at p_0=1 and is thus poorly approximated by a Gaussian.

The asymmetry of the marginalized posterior for the polarization fraction is shown in terms of the skewness in Fig. <ref>. The unbound prior scenario employs the unnormalized uniform prior p_0≥0 and p_r=1 (here all values of p_r produce the same result). It shows that when statistics are low, the marginalized posterior gets a longer tail on the right (positive skewness). This is because the peak moves left, which is equivalent to measuring a higher p_r for a fixed p_0, that is, the same interpretation as for the relative mean bias shown in Fig. <ref>.



Although the marginalized posterior for the polarization angle is always symmetric, it deviates from Gaussianity when statistics are low, for example, approaching the uniform distribution since all angles are equally likely. It is also affected by the prior when p_r is sufficiently high, resulting in a distribution with longer tails because, as shown in Fig. <ref>, the prior only truncates the part of the likelihood contributing to the central part of the marginalized distribution. Therefore, the uncertainty given by

    σ_ψ_r=1/√(2)μ_r p_r√(N)/S

(derived in Appendix <ref>) is not always valid.

In a Bayesian approach, the posterior is used instead of a parameter estimate and its uncertainty since the posterior contains the full information about the parameter. The posterior cannot easily be described in text so a simplified description is necessary such as its peak and the region of Highest Posterior Density (HPD) containing the probability content corresponding to one Gaussian standard deviation.

It is possible to quantify when the HPD region must be derived from the posterior and when Eq. <ref> and Eq. <ref> are good approximations by considering the ratios σ_p__HPD/σ_p_r and σ_ψ__HPD/σ_ψ_r as functions of MDP/p_r, which are shown in Figs. <ref> and <ref>,  respectively. In both figures, the unbound prior scenario uses the unnormalized uniform prior p_0≥0 and p_r=1. Only scenarios with R=0 are considered, meaning that the CLT approximation is valid in the entire range of MDP/p_r for any p_r. Increasing R does not change the results shown in this section but may invalidate the CLT approximation. Such cases cannot be represented because p_r is not a sufficient statistic. However, since the posterior shape does not change drastically (as shown in Fig. <ref>), the following discussion is a good qualitative description of its behavior for any R.

For ease of comparison, σ_p__HPD is defined as half of the region containing 1σ Gaussian probability. Additionally, the Gaussian interval is limited to account for the prior, for example, p_0=0.8±0.3 is the interval [0.5,1.0] which gives an uncertainty of σ_p_r=0.25,  and therefore it does not contain 68.3% (1σ) probability content. 





The unbound prior scenario in Fig. <ref> shows that for MDP/p_r∼1 the marginalized posterior has longer tails than a Gaussian. Additionally, the prior truncates the posterior for high p_r or for low statistics. This can result in either a higher or a lower σ_p__HPD/σ_p_r depending on where exactly the distribution is truncated. Unless p_r is high, σ_p_r is a good approximation for the magnitude of the uncertainty at MDP/p_r∼1 but it is important to remember that the posterior is asymmetric for such low statistics.

For the polarization angle, Fig. <ref> shows that σ_ψ__HPD>σ_ψ_r. The situation is simpler than for p_0 because the angle has a symmetric posterior. The unbound prior shows the deviation from a Gaussian manifesting in longer tails. Figure <ref> shows that when adding the prior 0≤ p_0 ≤ 1, the part of the likelihood which extends past p_0=1 is truncated. This part contributes to the density at the center of the marginalized posterior. Removing it further extends the tails. It follows that Eq. <ref> underestimates the uncertainty on the polarization angle by relative 10%-20% at MDP/p_r∼1.

This analysis shows that in the limit of low statistics, MDP/p_r∼1, the uncertainty on the polarization angle is not well-described by Eq. <ref> and a Bayesian treatment is necessary. For high p_r, such a treatment is required even for high-statistics measurements, MDP/p_r∼0.5, because of the asymmetry in the uncertainty on the polarization fraction.


To illustrate the effects described above, two examples are presented in Table <ref>. The first is a recent measurement of the Crab nebula by PoGO+ (Chauvin et al. <cit.>) and the second is a hypothetical measurement of a high reconstructed polarization fraction highlighting the importance of the prior 0≤ p_0≤ 1. The largest differences are in the polarization fraction and the uncertainty on the polarization angle.







§ UNCERTAINTIES ON THE MODULATION FACTOR


The uncertainty σ_μ on the modulation factor μ_0 can be minimized for most polarimeters designed for measuring point sources by increasing the statistics during calibration tests. However, for polarimeters measuring GRBs it is often impossible to make σ_μ arbitrarily small. The problem is that μ varies depending on the location of the GRB, typically having the highest values for on-axis GRBs but significantly lower values for GRBs located at a large angular separation from the detector axis. Determining the location of the GRB by using a polarimeter involves large uncertainties since it is often not optimized for the task. The uncertainty on the location propagates to a non-negligible σ_μ. Additionally, the primary spectrum of GRBs may not be reconstructed with sufficient precision by the polarimeter, introducing further uncertainties in the simulation required for deducing the μ_0 for a particular GRB. If the GRB is simultaneously observed by a dedicated GRB monitor, the uncertainty on its location and spectrum will be smaller, but no GRB monitor has complete sky coverage or 100% duty-cycle, so it is inevitable that some GRBs will only be seen by the polarimeter.

An example of a GRB polarimeter is POLAR (Produit et al. <cit.>). For a typical bright GRB,  POLAR is expected to have a σ_μ/μ of between 5 and 15% assuming that the burst occurs on-axis (Suarez-Garcia et al. <cit.>). The simplest way to account for this additional uncertainty is to propagate σ_μ through


    p_r=M/μ_r
    ,

where M is the measured modulation. The total symmetric uncertainty is then given by

    σ_tot=p_r√((σ_μ/μ_r)^2+(μ_rσ_p_r/M)^2)
    .


To check for which parameters the symmetry is a good approximation, an additional prior can be added to the posterior Eq. <ref>. For simplicity, this prior is assumed to be Gaussian but it can vary depending on localization sensitivity. The nuisance parameter μ_0 can then be marginalized over, yielding

    P(p_0,Δψ_0|p_r,μ_r)= 𝒩∫^1_0 P(p_0,Δψ_0)
    ×
    exp[-(μ_r-μ_0)^2/2σ_μ^2]× L(p_r,0,μ_r|p_0,Δψ_0)
     dμ_0
    
    ,

where the integral is taken between 0 and 1 because μ_0>1 is unphysical. As an example, the posterior (marginalized over the polarization angle) for a signal-only measurement at MDP=10% (CLT approximation is valid because S is large), p_r=0.4 and μ=0.4 is shown in Fig. <ref> for different σ_μ/μ. Due to the high statistics, MDP/p_r=0.25, the distribution is symmetric for low σ_μ/μ, however, the tail on the right grows rapidly as σ_μ/μ increases. The behavior is governed by the reciprocal Gaussian distribution which the posterior approaches in the limit of high photon statistics


    g(p_0|p_r,μ_r)=1/√(2π)σ_μp_rμ_r/p_0^2exp[-(p_rμ_r/p_0-μ_r)^2/2σ_μ^2]
    ,

shown as an approximation in Fig. <ref>. If not for the prior 0≤ p_0 ≤1, the moments of this distribution would be undefined.

Figure <ref> demonstrates that for σ_μ/μ>10% the shape of the posterior changes significantly and Eq. <ref> is required. In the extreme case of σ_μ/μ=20%, the half-width of the HPD region is 0.092, whereas Eq. <ref> yields 0.086, implying not only that the uncertainty is asymmetric but also that it is significantly larger.



These results can be related to the instrument perfomance by studying the effect σ_μ/μ has on the MDP. A frequentist approach must be followed since MDP is a frequentist concept. MDP is derived by finding the 99% upper limit for the Rayleigh distribution – a special case of the Rice distribution where p_0=0. However, since there is an uncertainty on μ, the Rayleigh distribution must be multiplied by the likelihood for μ (assumed to be a Gaussian for simplicity) and integrated to yield 

    f(p_r|p_0=0,μ_0)=∫_0^∞p_rμ_r/√(2π)πσ^2σ_μ
    ×exp[-μ_r^2 p_r^2/2σ^2-(μ_r-μ_0)^2/2σ_μ^2]dμ_r
    .

Finally, f(p_r|p_0=0,μ_0) can be integrated to find the 99% upper limit so that ∫_0^MDP_σf(p_r|p_0=0,μ_0)dp_r=0.99. Figure <ref> shows the relative increase in the MDP, defined as the ratio MDP_σ/MDP where MDP is given by Eq. <ref>. Although the effect is small for low σ_μ/μ (1% at σ_μ/μ=5%), it becomes significant at larger σ_μ/μ (14% at σ_μ/μ=15%) and deteriorates the instrument performance for extreme values (80% at σ_μ/μ=25%). The result is independent of the intial MDP.





§ CONCLUSIONS

The results presented here provide a means of quantifying the errors incurred when using the simple estimators for polarization parameters as well as for their uncertainties. These errors are related to the well-established figure of merit MDP and the reconstructed polarization fraction, making the results easily applicable to any X-ray polarimeter. Unlike some previous works, this analysis does not ignore the correlation between the Stokes parameters Q and U.

Additionally, the extent to which the reconstructed fraction p_r and polarization angle ψ_r can be used as sufficient statistics is explored. In certain situations, such as for high polarization fraction and high signal-to-background ratio, Q and U are not Gaussian distributed and their likelihood is non-Gaussian, resulting in a significantly different posterior for the polarization parameters when the number of photons is low. In such cases, the Stokes parameters cannot be used at all because (p_r,ψ_r) are not sufficient statistics.

The fact that the reconstructed polarization fraction is biased towards higher values implies that binning a high-statistics data set into smaller subsets to estimate the evolution of the polarization fraction will yield incorrect results when using the simple Gaussian estimator p_r. This work provides a way to decide how coarsely the data must be binned in order to justify the use of Gaussian estimators.

The uncertainties on the estimated parameters are as important as the parameter estimates themselves. When statistics are high (MDP/p_r<0.5) the errors are small and can usually be neglected (justifying Gaussian assumptions made when using simple estimators) unless the reconstructed polarization fraction is high, for example, >0.8. However, in the statistics-limited regime (MDP/p_r>1) the systematic error made from using such simple estimators is non-negligible compared to the statistical uncertainty. Additionally, it is shown how quickly the posterior becomes asymmetric as the statistical power decreases. This effect is strongly dependent on the reconstructed polarization fraction because of the prior 0≤ p_0≤1.

Lastly, the effect of uncertainties on the modulation factor μ is studied. It is shown to be important once the relative uncertainty exceeds 10% and to dominate the performance of an instrument when it is above 20%. This is relevant for the optimization of future GRB polarimeters, since they tend to have large uncertainties on μ due to the difficulty of localizing bursts and measuring the primary GRB spectrum. It also shows that simple Gaussian estimators cannot be used in the high-photon-counts regime for GRB polarimeters when localization uncertainties are high.


This research was supported by the Swedish National Space Board. M. Kiss, M. Pearce, F. Xie and the Referee are thanked for providing constructive feedback on the manuscript.




[2014]Beilicke2014 Beilicke, M., Kislat, F., Zajczyk, A., et al., 2014, JAI, 3, 1440008

[2010]Bellazzini2010 Bellazzini, R., Muleri, F., 2010, NIM A, 623, 766-770

[2013]Chauvin2013SPI Chauvin, M., Roques, J. P., Clark, D. J., et al., 2013, , 769, 137

[2016]Chauvin2016PoGOLite Chauvin, M., Florén, H.-G., Jackson, M., et al., 2016, , 456, L86-L88

[2017]Chauvin2017PoGO+ Chauvin, M., Florén, H.-G., Friis, M., et al., 2017, Sci. Rep., 7, 7816

[1983]Clarke1983 Clarke, D., Stewart, B. G., Schwarz, H. E., et al., 1983, , 126, 260-264

[2008]Dean2008 Dean, A. J., Clark, D. J., Stephen, J. B., et al., 2008, Science, 321, 1183-1185

[2008]Forot2008 Forot, M., Laurent, P., Grenier, I. A., et al., 2008, , 688, L29

[1939]Jeffreys1939 Jeffreys, H., 1939, Theory of probability (Oxford University Press)

[2014]Kaaret2014 Kaaret, P., 2014, arXiv, 1408.5899

[2015]Kislat2015 Kislat, F., Clark, B., Beilicke, M., et al., 2015, Astropart. Phys., 68, 45-51

[2011]Krawczynski2011 Krawczynski, H., Garson, A., Guo, Q., et al., 2011, Astropart. Phys., 34, 550-567

[1997]Lei1997 Lei, F., Dean, A. J., Hills, G. L., 1997, Space Sci. Rev., 82, 309-388

[2003]Lyutikov2003 Lyutikov, M., Pariev, V. I., Blandford, R. D., 2003, , 998-1009

[2014]Maier2014 Maier, D., Tenzer, C., Santangelo, A., 2014, , 126, 459-468

[2016]Moran2016 Moran, P., Kyne, G., Gouiffès C., et al., 2016, , 456, 2974-2981

[1972]Novick1972 Novick, R., Weisskopf, M. C., Berthelsdorf, R., et al., 1972, , 174, L1

[2005]Produit2005 Produit, N., Barao, F., Deluit, S., et al., 2005, NIM A, 550, 616-625

[2012]Quinn2012 Quinn, J. L., 2012, , 538, A65

[1985]Simmons1985 Simmons, J. F. L. and Stewart, B. G., 1985, , 142, 100-106

[2009]Slowikowska2009 Słowikowska, A., Kanbac, G., Kramer, M., et al., 2009, , 397, 103-123

[2010]Suarez-Garcia2010 Suarez-Garcia, E., Haas D., Hajdas W., et al., 2010, NIM A, 624, 624-634

[2018]Vadawale2018 Vadawale S. V., Chattopadhyay, T., Mithun N. P. S., et al., 2018, Nature Astronomy, 2, 50-55

[2006]Vaillancourt2006 Vaillancourt, J. E., 2006, , 118, 1340

[1976]Weisskopf1976 Weisskopf, M. C., Cohen, G. G., Kestenbaum, H. L., Long, K. S., et al., 1976, , 208, L125-L128

[1978]Weisskopf1978 Weisskopf, M. C., Silver, E. H., Kestenbaum, H. L., et al., 1978, , 220, L117-L121

[2010]Weisskopf2010 Weisskopf, M. C., Elsner R. F., O'Dell S. L., 2010, , 7732, 77320E

[2016]Weisskopf2016 Weisskopf, M. C, Ramsey, B., O'Dell, S., et al., 2016, , 9905, 990517

[2011]Yonetoku2011 Yonetoku, D., Murakami, T., Gunji, S., et al., 2011, , L30









§ STOKES PARAMETERS



The Stokes parameters can be constructed from

    q_i=cos(2ψ_i)
    
    u_i=sin(2ψ_i)
    
    .

To find an expression for p_0 , the means q and u must be computed.

    q=∫_0^2πcos(2ψ)f(ψ)dψ
    ,


    q=μ_0 p_0S/2(B+S)cos(2ψ_0)
    ,

where f(ψ) is given by Eq. <ref>. Similarly

    u=μ_0 p_0S/2(B+S)sin(2ψ_0)
    .

The polarization fraction p_0 is then given by

    p_0=2√(Q_0^2+U_0^2)/μ_0
    ψ_0=1/2arctanU_0/Q_0
    ,

where Q_0 and U_0 are the normalized Stokes parameters defined as  

    Q_0=B+S/Sq=μ_0 p_0/2cos(2ψ_0)
    
    U_0=B+S/Su=μ_0 p_0/2sin(2ψ_0)
    
    .

The sufficient statistics for data generated from (Q_0,U_0) are (Q_r,U_r) defined as

    Q_r=1/S∑_i=1^Nq_i=μ_r p_r/2cos(2ψ_r)
    
    U_r=1/S∑_i=1^Nu_i=μ_r p_r/2sin(2ψ_r)
    
    ,

where N=B+S. In polar coordinates this becomes

    p_r=2√(Q_r^2+U_r^2)/μ_r
    ψ_r=1/2arctanU_r/Q_r
    .




§ UNCERTAINTIES ON STOKES PARAMETERS



The uncertainties on Q and U as well as their covariance can be derived by computing their second moments.

    q^2=∫_0^2πcos^2(2ψ)f(ψ)dψ=1/2
    .

Similarly

    u^2=1/2
    ,

and the cross-term is

    qu=1/2π(S+B)∫_0^2πcos(2ψ)sin(2ψ)× f(ψ)dψ=0
    .

Combining Eq. <ref> and Eq. <ref> yields the standard deviation

    σ_q=√(1/2-(μ_0 p_0Scos(2ψ_0)/2(B+S))^2)
    .

Here σ_q describes the dispersion of q but it is σ_Q which is of interest. It is given by

    σ_Q^2=σ_q^2×N/S^2
    ,

where the division by S^2 comes from the definition of normalized Stokes parameters.

    σ_Q=√(1/S(N/2S-μ_0^2p_0^2cos^2(2ψ_0)/4))
    .

Similarly

    σ_U=√(1/S(N/2S-μ_0^2p_0^2sin^2(2ψ_0)/4))
    .

It is finally possible to compute the covariance and the correlation coefficient ρ

    Cov(Q,U)=Cov(q,u)×N^2/S^3=-μ_0^2p_0^2/8Ssin(4ψ_0)
    ,


    ρ=-Sμ_0^2p_0^2sin(4ψ_0)/√(16N^2-8NSμ_0^2p_0^2+S^2μ_0^4p_0^4sin^2(4ψ_0))
    .




§ RELATIVE MEAN BIAS



An approximate expression for the relative mean bias can be derived from the Rice distribution

    f(p_r|p_0)=p_r/σ^2exp(-p_r^2+p_0^2/2σ^2)× I_0(p_rp_0/σ^2)
    ,

where σ≡√(2N)/μ_0 S. The mean p_r is given by

    p_r=√(π/2)σ×exp(-p_0^2/4σ^2)
    ×((1+p_0^2/2σ^2)I_0(p_0^2/4σ^2)+p_0^2/2σ^2I_1(p_0^2/4σ^2))
    
    ,

where I_0 and I_1 are modified Bessel functions of the zeroth and first order, respectively. For high-statistics measurements (where p_0≫2σ) the Bessel functions can be approximated by expanding them to second order.

    lim_x →∞ I_0(x)=e^x/√(2π x)(1+1/8x)
    ,


    lim_x →∞ I_1(x)=e^x/√(2π x)(1-3/8x)
    .

Finally

    lim_p_0/σ→∞p_r=σ^2/2p_0+p_0
    .

So in the limit of high statistics, the relative mean bias β is given by

    lim_p_0/σ→∞β≈σ^2/2p_0^2=N/S^2μ_0^2p_0^2=(MDP/4.29p_0)^2
    .

By using recursion, β can be expressed as a function of p_r instead of p_0 which is not known a priori. This simplifies to

    β≈(MDP/4.29p_0)^2=∑_n=1^∞ C_nx^2n
    
    =∑_n=1^∞1/n+12n  nx^2n=1-2x^2-√(1-4x^2)/2x^2
    ,

where x=MDP/4.29p_r and C_n is the nth Catalan number. 




§ GAUSSIAN UNCERTAINTIES ON POLARIZATION PARAMETERS


The uncertainty on p_r can be derived by using the standard uncertainty propagation formula

    σ_p_r^2=|dp_r/dQ_r|^2σ_Q^2+|dp_r/dU_r|^2σ_U^2+2|dp_r/dQ_r||dp_r/dU_r|Cov(Q,U)
    ,

yielding

    σ_p_r=2/μ_r√(1/S(N/2S-μ_0^2p_0^2/4))
    ,

and similarly for ψ_r

    σ_ψ_r=1/√(2)μ_r p_r√(N)/S
    .






