Time-Resolved GRB Polarization with POLAR and GBM

Simultaneous $\gamma$-ray measurements of gamma-ray burst (GRB) spectra and polarization offer a unique way to determine the underlying emission mechanism(s) in these objects as well as probing the particle acceleration mechanism(s) that lead to the observed $\gamma$-ray emission. Herein we examine the jointly-observed data from POLAR and GBM of GRB 170114A to determine its spectral and polarization properties and seek to understand the emission processes that generate these observations. We aim to develop an extensible and statistically sound framework for these types of measurements applicable to other instruments. We leverage the existing 3ML analysis framework to develop a new analysis pipeline for simultaneously modeling the spectral and polarization data. We derive the proper Poisson likelihood for $\gamma$-ray polarization measurements in the presence of background. The developed framework is publicly available for similar measurements with other $\gamma$-ray polarimeters. The data are analyzed within a Bayesian probabilistic context and the spectral data from both instruments are simultaneously modeled with a physical, numerical synchrotron code. The spectral modeling of the data is consistent with a synchrotron photon model as has been found in a majority of similarly analyzed single-pulse GRBs. The polarization results reveal a slight trend of growing polarization in time reaching values of ~30% at the temporal peak of the emission. Additionally, it is observed that the polarization angle evolves with time throughout the emission. These results hint at a synchrotron origin of the emission but further observations of many GRBs are required to verify these evolutionary trends. Furthermore, we encourage the development of time-resolved polarization models for the prompt emission of GRBs as the current models are not predictive enough to enable a full modeling of our current data.


Introduction
Polarization measurements from astrophysical objects are a key piece of information to deciphering the physics and geometry of regions that emit the observed photons. The emission from gamma-ray bursts (GRBs) has been notoriously difficult to understand due to the complexity of modeling their broadband, prompt γ-ray emission. Recent results have provided evidence that the prompt emission is the result of synchrotron radiation from electrons accelerated to ultra high-energies via magnetic reconnection (Burgess et al. 2014;Zhang et al. 2016Zhang et al. , 2018Burgess et al. 2019). Measurements of the optical polarization from a GRB's prompt emission have similarly pointed to a synchrotron origin of the emission (Troja et al. 2017). However, spectral modeling of photospheric based emission has also provided adequate fits to a subset of GRBs (Ryde et al. 2010;Ahlgren et al. 2015;Vianello et al. 2018b). Measurements of polarization can break this degeneracy (Toma et al. 2009;Gill et al. 2018). Photospheric emission will typically produce unpolarized emission although a moderate polarization level is possible in special circumstances (Lundman et al. 2018) and predicts very specific changes of the polarization angle (Lundman et al. 2014). On the other hand, synchrotron emission naturally produces a range of polarized emission depending on the structure of the magnetic field and outflow geometry (Waxman  Lyutikov et al. 2003;Granot 2003). Thus, being able to fit synchrotron emission to the observed spectrum while simultaneously detecting polarization provides a clear view of the true emission process.
Several reports of polarization measurements have been produced by a variety of instruments (Coburn & Boggs 2003;Yonetoku et al. 2011;Hanlon et al. 2008;Chattopadhyay et al. 2017). Measurements by non-dedicated instruments like those reported by BATSE, IBIS, SPI and RHESSI all suffer from problems with instrumental effects or poorly understood systematics (McConnell 2017) making it impossible to draw conclusions based on these results. More recently, AstroSAT has reported polarization measurements of several GRBs (Chattopadhyay et al. 2017). The systematics and procedure related to obtaining these measurements is not immediately clear. The quoted error distributions contain unphysical regions of parameter space (polarization degree > 100%) and are thus questionable. Past measurements of polarization by the first dedicated GRB polarimeter, GAP, provided hints of polarized emission (Yonetoku et al. 2011), while measurements by COSI provided an upper limit on the polarization degree (Lowell et al. 2017). The statistics of these measurements do however not allow constrains on the emission mechanisms. Furthermore, the techniques relied on background subtraction. As both the background and signal counts are Poisson distributed, subtraction is an invalid procedure that destroys statistical information, thus all reported significances are questionable.
The POLAR experiment (Produit et al. 2018) on board the Chinese space laboratory Tiangong-2 observed 55 GRBs and reported polarization measurements for 5 of these GRBs (Zhang et al. 2019). Time-integrated analysis of these GRBs resulted in strict upper limits on the polarization degrees. The most likely polarization degrees found in that analysis are non-zero but remain compatible with an unpolarized emission, leading to the conclusion that GRBs are at most moderately polarized. Using time-resolved analysis it was however found that the polarization of GRB 170114A was most compatible with a constant polarization degree of ∼ 28% with a varying polarization angle. Summing polarized fluxes with varying polarization degrees produces an unpolarized flux. The detection of an evolution in polarization angle within this single pulse GRB could explain the low polarization degrees found for all 5 GRBs. The results presented in (Zhang et al. 2019) do, however, not allow for a detailed time-resolved study of the remaining 4 GRBs, nor do they allow the determination of the nature of the evolution of the polarization angle in GRB 170114A.
Coincidentally, several of the GRBs observed by POLAR were simultaneously observed by the Fermi -GBM. Herein, we present a technically advanced modeling of the polarization and spectral data simultaneously with data from both instruments. This allows the incorporation of information contained in both data sets leading to improved sensitivity and an altogether more robust analysis. This work is divided as follows: The methodology and modeling is described in Sections 2 and 3 and the results are interpreted in Section 4.

Data Analysis and Methodology
For the analysis herein, we have developed a new approach of simultaneously fitting both the spectral data from POLAR and GBM along with the POLAR scattering angle (SA) or polarization data. This alleviates the need for approximate error propagation of the spectral fits into the polarization analysis. Using the abstract data modeling capabilities of 3ML 1 (Vianello et al. 2015), a framework was developed to directly model all data simultaneously with a joint-likelihood in each dataset's appropriate space. Below, we detail each part of the methodology.

Location and Temporal Analysis
Spectral and polarization analysis for both GBM and POLAR rely on knowledge of the sky-position (δ) of the GRB in question. Both being all-sky surveyors, GBM and POLAR lack the ability to image GRBs directly. However, using the BALROG technique (Burgess et al. 2018), we can use the spectral information obtained in the GBM data to locate the GRB. Using a synchrotron photon model (see Section 3), we were able to locate the GRB to RA= 13.10 ± 0.5 deg, DEC=−13.0 ± 0.6 deg. Using this location, spectral and polarization responses were generated for all data types. We note that a standard GBM position 2 exists and, along with their uncertainties, was used for the polarization results presented in (Zhang et al. 2019), however, the standard localization technique has known systematics and now possess arbitrarily inflated error distributions (Connaughton et al. 2015). We found the BALROG derived location much more precise than that of the standard location analysis (see Fig 1), allowing to reduce the systematic errors included in the polarization results presented in (Zhang et al. 2019). Additionally, it has now been shown that the BALROG locations as systematically more accurate (Berlato et al. in prep.).
The chief focus of this analysis is temporal variation in the polarization parameters. We compute the minimum variability time-scale (MVT) (see Vianello et al. 2018b, for details) on the POLAR SA light curve. This yields a MVT of ∼ 0.3 s (Fig. 3). For completeness, the MVTs for both the GBM and POLAR spectral light curves were computed as well. Both analyses yield similar results. Therefore, we are able to analyze data on this time-scale without the concern of summing over evolution of spectral (Burgess & Ryde 2015). However, the raw polarization data do not allow for us to check for variability in the polarization angle prior to fitting. Therefore, it is possible that the angle could change on a time-scale smaller than our selected time-intervals. This could reduce the overall inferred polarization.
With the MVT determined, we utilize the Bayesian blocks algorithm (Scargle et al. 2013) to objectively identify temporal bins for the analysis. The SA light curve is utilized to perform the analysis. The temporal bins created are on the order of the MVT. A total of nine bins were selected and used for spectral and polarization analysis (see Table 4.1). Fig. 1. The BALROG location (red posterior samples) of GRB 170114A derived by fitting the peak of the emission for both the location and spectrum simultaneously. The blue contours display the 1,2, and 3 σ standard GBM catalog location as obtained from the Fermi Science Support Center (FSSC).

Spectral Analysis
The standard γ-ray forward-folding approach to spectral fitting is adopted where we have sky location (δ) dependent responses for both the GBM and POLAR detectors (R γ ) and fold the proposed photon model (n γ ) solution through these responses to produce detector count spectra (n pha ). Thus, for the i th detector in the j th pulse-height amplitude (PHA) channel. Here, δ is the sky location of the GRB. Both POLAR and GBM have their total observed counts Poisson distributed and their backgrounds determined via fitting polynomials in time to off-source regions of the light curves. Thus, Gaussian distributed background counts are estimated by integrating these polynomial models over the source interval of interest. The uncertainty on these estimated counts is derived via standard Gaussian uncertainty propagation. This leads us to use a Poisson-Gaussian likelihood 3 for each detector for the spectral fitting.

Polarization Analysis
To enable performing joint fits of the spectra and the polarization a novel analysis technique was developed. Traditional polarization analysis techniques, such as those employed in (Yonetoku et al. 2011;Chattopadhyay et al. 2017) as well as in (Zhang et al. 2019), rely on fitting data to responses produced for a specific spectrum. This method does not allow joint fits of both the spectrum and polarization parameters nor does it allow naturally including systematic uncertainties from the spectral fits into the systematic uncertainties of the polarization. Here, in order to model the polarization signal seen in the data, we invoke a forward-folding method similar in concept to our approach to spectral analysis. We simulate polarized signals as function of polarization angle, degree and energy to create a matrix of SA distributions (often called modulation curves within the field of polarimetry) which can be compared to the data via the likelihood in data space. For details on the creation of the matrix see Appendix A. Mathematically, where n θ are counts in SA bin k, and R θ is the simulated response of the corresponding scattering bin. In words, we convolve the photon spectrum over the j th photon energy bin with the polarization response to properly weight the number of counts observed in each SA bin. Figure 4 demonstrates how changes in polarization angle and degree appear in the POLAR data space. Hence, our need to simultaneously fit for the photon spectrum which allows for direct accounting of the uncertainties in the weighting. POLAR observed SAs are measured as detector counts and thus Poisson distributed. The pollution of the source signal by background cannot be handled by background subtraction as has been done in previous work. Instead, a temporally off-source measurement of the background polarization is made in order to model the background contribution to the total measurement during the observation intervals. The background measurement is Poisson distributed in each of the k scattering bins. Due to the temporal stability of the background, as presented in (Zhang et al. 2019), we fit a polynomial in time to each of the k scattering bins via an unbinned Poisson likelihood. This allows us to reduce the uncertainty of the background by leveraging the temporal information. We can estimate the on-source background contribution (b k θ ) by integrating the polynomials over time and propagating the temporal fit errors. This implies that the polarization likelihood is also a Poisson-Gaussian likelihood just as with the spectral data. We verify that our approach allows us to identify the latent polarization parameters via simulations in Appendix B. The count rates are corrected for the proper exposure by computing the total dead time fraction associated with each interval. The method employed for dead time calculation is equivalent to that of Zhang et al. (2019).
The full joint likelihood of the data is thus a product over the spectral and polarization likelihoods which is detailed in Appendix C (see also Figure 5). We reemphasize that the spectral model and polarization model communicate with each other through the likelihood. This implies that the posterior density of the model is fully propagated to both datasets without any assumptions such as Gaussian error propagation. As will be seen in the following sections, the resulting parameter distributions can be highly asymmetric.
In a perfect world where all instruments are cross-calibrated over the full energy range, the instruments' various responses would predict similar observed fluxes for each measurement. However, we allow for a normalization constant between GBM and POLAR to account for any unmodeled discrepancies between the instruments. Both POLAR's polarization and spectral data are scaled by these constants which are unity when no correction is required 4 . This constant scale for the effective area by no means accounts for energy-dependent calibration issues.
In order to obtain the posterior parameter distributions, we use MULTINEST (Feroz et al. 2009) to simulate the model's posterior. MULTINEST utilizes nested sampling which is suitable for the multimodal distributions we observe as well as for the non-linear model and high-dimension of our parameter space. For the polarization parameters, we use uninformative priors of appropriate scale. The effective area normalizations are given informative (truncated Gaussians) priors centered at unity with a 10% width. The priors for the spectral modeling are discussed in the Section 3. We run MULTINEST with 1500 live points to achieve a high number of samples for posterior inference. Model comparison is not attempted and thus we do not use the marginal likelihood calculations 5 . As stated, for bothp and φ, we use uninformative priors in each parameters' domain. This is a valid choice for φ, but we note that an informative prior for the expected polarization from synchrotron emission could be used as an assumption. However, as discussed in Section 5, the theoretical predictions for GRB synchrotron models are not mature enough for us to assume such a prior at the current time. Nevertheless, in our work we tested Gaussian priors centered at moderate polarization and found that the data allowed for this assumption. Moreover, we found that our recovered φ was not affected by out choice of prior onp.

Synchrotron Modeling
With the recent finding that synchrotron emission can explain the majority of single-pulse GRBs, we choose to model the time-resolved photon spectrum with a physical synchrotron model. Following Burgess et al. (2019), we set where K is the arbitrary normalization of the flux, B is the magnetic field strength in Gauss, p is the injection index of the electrons, γ cool is the energy to which an electron will cool during a synchrotron cooling time, and Here, K 5/3 is a Bessel function, B crit = 4.13 · 10 13 G, and n e is determined by solving the cooling equation for electrons with the Chang and Cooper method (Chang & Cooper 1970). In mathematical expression, where the injected electrons are defined by a power law of index p where γ inj and γ max are the minimum and maximum injected electron energies respectively and the synchrotron cooling iṡ For our numerical calculations, we created a 300-point grid logarithmically distributed in γ. The linear equations in the implicit scheme form a tridiagonal matrix which is solved numerically with standard methods. The method of Chang & Cooper (1970) is numerically stable and inexpensive as well as shown to conserve particle number in the absence of sources and sinks. Thus, we are able to solve for the synchrotron emission spectrum quickly during each iteration of the fit. The numeric code is implemented in C++ and interfaced with Python into astromodels (Vianello et al. 2018a).
The overall emission is characterized by five parameters: B, γ inj , γ cool , γ max , and p. However, a strong co-linearity exists between B and γ inj as their combination sets the peak of the photon spectrum. Thus, both parameters serve as an energy scaling which forces the setting of one of the parameters. We choose to set γ inj = 10 5 though the choice is arbitrary and does not affect our results. It is therefore important to note that all parameters are determined relatively, i.e., the values of γ cool and γ max are determined as ratios to γ inj . Similarly, the value of B is only meaningful when determining the characteristic energies of γ cool and γ max or hν cool and hν max respectively. In other words, with our parameterization the spectra are scale free. The degeneracies can be eliminated by specifying temporal and radial properties of the GRB outflow which we neglect in this analysis.
Ideally, we would fit for the full set of parameters in the model. However, the already high-dimensionality of the model does not allow us to fit for the cooling regime of the model simultaneously with the polarization due to computational time constraints. Therefore, we first fit the spectral data alone to determine the amount of cooling present in the data. All spectra were found in the slow-cooling regime (Sari et al. 1998;Beniamini & Piran 2013). Thus, we fix the ratio of γ cool to γ inj during the full fits to the slow-cooling regime. Tests revealed that the cooling had no impact on the recovered polarization parameters. Additionally, the lack of high-energy data (via the Fermi -LAT) forces us to fix γ max such that the synchrotron cutoff is above the spectral window. We arrive at a three parameter fit for the spectrum: B, p and the arbitrary spectral normalization (K). B and K are given uninformative scale priors and p a weakly-informative, Gaussian prior centered around p = 3.5. The effective area constants applied to the POLAR response are given truncated Gaussian priors centered at unity with a width of 10% to reflect our prior belief that the instruments are well-calibrated to one another.

Results
In the following two sections, we present the results from the combined polarization and spectral analysis separately. Corner plots of the important (non-nuisance) parameter marginal distributions are displayed in Appendix D.

Polarization
The POLAR polarization data are well described by our modeling of the POLAR instrument. Scattering angle data show good agreement between the data and the model as demonstrated in Figure 6. In order to validate the model's ability to generate the data, we performed posterior predicative checks (PPCs) (Betancourt 2018) of the polarization data for all time intervals. For a subset of posterior samples chosen with appropriate posterior probability, latent polarization and spectral models were generated and subsequent data quantities where sampled from the likelihood. The model was able to sufficiently generate replicated data similar to the observed (see Figure 7) in most cases. It is likely that minor deficiencies still exist in the instrumental responses.
We observe no polarization at the beginning of the pulse and moderate (∼ 30%) polarization as time proceeds, compatible with the results presented in Zhang et al. (2019). Interestingly, we observe a large change in the polarization angle with time (see Figure 8). Although the time-intervals used in this study are different from those used in (Zhang et al. 2019) it can be deduced that the polarization angles found here agree with those in (Zhang et al. 2019). The end of the pulse has a relatively weak signal and thus poorly identified polarization parameters. The 68% credible regions are listed in Table 4.1. Clearly, the level of polarization during the peak of the emission is probabilistically equivalent to both moderate and low polarization whereas during the beginning of the emission the polarization is definitely low even though the ratio of background to total signal is high.
We stress that it is not appropriate to perform model comparison on nested model parameters, e.g., comparing between zero polarization and greater than zero polarization. This includes the use of Bayes factors (Chattopadhyay et al. 2017) which are ill-defined for improper priors and for comparing between discrete values of a continuous parameter (Gelman et al. 2013). Polarization is not a detected quantity, but a parameter. Given that we have detected the GRB, it is important to quote the credible regions of the polarization parameter rather than perform model comparison between discrete values.

Spectra
POLAR and GBM both agree in overall spectral shape and relative normalization of the observed flux. Moreover, the spectral results demonstrate that the synchrotron spectrum is a good, predictive description of the spectral data as displayed in Figure 9. This is both a confirmation that past studies with synchrotron relying on GBM data only are reliable as well as the outstanding calibration between the GBM and POLAR.
As noted above, it is not possible to disentangle the intrinsic parameters of the synchrotron emission without further assumptions. Therefore, we only quote the injection energy in Table 4.1. The evolution of the spectrum is shown in Figure  10. The temporal evolution of the νF ν spectral peak follows a broken power law. We find values between ∼ 3 − 4 for the electron power law injection spectral index. These values are steeper than those of the canonical index expected from shock acceleration (Kirk et al. 2000).
It is possible that other physical spectral models also provide acceptable, predictive fits to the data. However, these models, e.g., subphotospheric dissipation have yet to demonstrate acceptable spectral fits on a large sample of GRBs. Moreover, the numerical schemes (Pe'er & Waxman 2005) required to compute the emission form these models are more complex than that of our synchrotron modeling, require far more computational time, and are not publicly available for replication. Photospheric models also require special geometrical setups to produce polarization. This makes them more predictive and indeed a pertinent set of models to test.

Discussion
For the first time, the polarization and spectrum of GRB prompt γ-ray emission has been fitted simultaneously. Furthermore, the spectral data have been described with a physical synchrotron model consistent with the spectral data of two very distinct spectrometers. We argue that it is unlikely for the spectral and polarization data to conspire to point towards an optically-thin synchrotron origin of the emission. However, the current predictive power of GRB prompt emission polarization theory is not developed enough for our measurements to definitely select synchrotron over other emission mechanisms. Therefore, we speculatively leverage previous spectral results that show that synchrotron emission is dominant mechanism in single-pulse GRBs. Burgess et al. (2019) argued that the observation of synchrotron emission in GRBs invalidates the standard fireball model (Eichler & Levinson 2000). Similar predictions were made before they were supported by data (e.g. Zhang & Pe'er 2009). These results have alluded to a magnetically-dominated jet acceleration mechanism possibly resulting in comoving emission sites or mini-jets (Barniol Duran et al. 2016;Beniamini et al. 2018). These results were arrived at considering spectral analysis alone. The moderate polarization degree observed in this work requires a development in the prediction of the temporal polarization predictions of these models in order to fully interpret their meaning.  While our observations provide broad ranges for the observed polarization degree, the changing polarization angle is easily observed. Figure 11 shows how the both the peak of the synchrotron spectrum and the polarization angle grossly track each other in time. Detailed model predictions for the evolution of the polarization angle during the GRB are not available. We are therefore not able to interpret the change in angle and encourage the community to develop detailed predictions which can be fitted to our data in the future. With more predictive models, appropriate informative priors can be adopted. Even more, spectral parameters can be formulated in terms of polarization parameters making the models stricter and the data more useful. Thus, we hope models develop in the near future.
The combination of POLAR and GBM observations of GRBs enables energy-dependent polarization measurements and is a project currently under development. These measurements will allow us to decipher if polarization increases around the peak of the photon spectrum which would be a signature of synchrotron emission, or if the polarization is higher at low energies as predicted by Lundman et al. (2018). Further multi-messenger studies and missions are encouraged to answer these questions. . The count spectra of POLAR and GBM from the joint spectral/polarization fits. The shaded regions indicate the 2σ credible regions of the fit. Data from a GBM NaI, BGO, and POLAR and displayed in green, black and red respectively.

Software availability
The analysis software utilized in this study are primarily 3ML and astromodels. We have designed a generic, preliminary polarization likelihood for similar X-ray polarization instruments both within 3ML 6 and astromodels 7 . Additionally, the POLAR pipeline 8 we have developed is fully designed to be easily modified for other instruments with polarimetric data. We note that these software distributions are preliminary and encourage the community to participate in their development.  In our situation we do not have a spectral or polarization model for the background process. Thus, we adopt the common procedure of maximizing the probability with respect to b i a priori leading to the profile likelihood referred to as PGSTAT 9 . Thus, with j datasets i.e., spectral or polarization detector's data, the total likelihood for our observations is Previous γ-ray polarization estimates have been achieved via background-subtracted data with the assumption of a Gaussian likelihood. This is improper and can lead to systematically biased results. There have also been attempts to transfer the statistical techniques used in optical polarimetry (Vaillancourt 2006;Quinn 2012) to γ-rays. These techniques are invalid for measurements that infer latent polarization via a secondary measurement such as the Compton scattering angle of a photon. Moreover, these techniques assume that none of the inherent difficulties of γ-ray photon measurement are present, namely, low-counts and dispersion both in energy and scattering angle. We have dealt with this first issue via the proper Poisson-based likelihood. The second issue via our modeling of the responses of our instruments directly in the inference process.