Time-varying Polarized Gamma-Rays from GRB 160821A: Evidence for Ordered Magnetic Fields

Gamma-ray bursts (GRBs) are the most intense astrophysical outbursts in the universe. In the last several decades, the spectra of prompt gamma-ray emission of GRBs have been extensively studied using various space observatories such as the Burst and Transient Source Experiment (BATSE) on board the Compton Gamma-Ray Observatory (CGRO; Fishman 2013), Niel Gehrels Swift (Gehrels & Swift 2004), Fermi (Atwood et al. 2009; Meegan et al. 2009), etc. The radiation process producing the prompt gamma-ray emission, however, still remains a mystery. Polarization along with spectrum measurements can provide an insight into this long-standing enigma. Performing polarization measurements of prompt emission is highly challenging, mainly because of the scarcity of incident photons and the brevity of the event. Previously, polarization measurements of prompt gamma-ray emission were attempted only for a few cases by the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI; Coburn & Boggs 2003; Wigger et al. 2004), the INTErnational Gamma-Ray Astrophysics Laboratory (INTEGRAL; Götz et al. 2009, 2013, 2014; McGlynn et al. 2009), the GAmma-ray burst Polarimeter (GAP; Yonetoku et al. 2011, 2012), the Cadmium Zinc Telluride Imager (CZTI; Rao et al. 2016; Chattopadhyay et al. 2017; Chand et al. 2018, 2019), POLAR (Burgess et al. 2019; Zhang et al. 2019), etc. but the results were statistically less significant and sometimes unconvincing (for a recent review please refer to McConnell 2017). In this Letter, for the first time, we present conclusive evidence of polarization across the GRB 160821A in the energy range 100–300 keV using the CZTI instrument on board AstroSat (Singh et al. 2014).

On 2016 August 21, the Burst Alert Telescope (BAT) on board the Neil Gehrels Swift observatory (Gehrels & Swift 2004) triggered and located the GRB 160821A at (R.A., decl.) =(171.248, 42.343) with 3' uncertainty (Markwardt et al. 2016) at 20:34:30 UT, along with other space observatories such as Konus-Wind and the Gamma-Ray Burst monitor on board the CALorimetric Electron Telescope (CALET). However, due to solar observing constraints Swift could not slew to the BAT position until a week later. Hence, there were no X-ray or ultraviolet (UV)/optical telescope observations for this burst. Half an hour after the trigger time, an optical transient was detected by ground-based telescopes, but no redshift measurements could be made.

The Fermi Gamma-Ray Burst Monitor (GBM; Meegan et al.2009) also triggered on the burst at 20:34:30.04 UT (Stanbro & Meegan 2016). The GBM light curve included a precursor emission starting from trigger time, T₀ until T₀+112 s, followed by a bright emission episode with a duration T90¹¹ of 43 s. For the time interval T₀-4.1 s–T₀+194.6 s, an energy flux of 2.86 ±0.007 × 10⁻⁶ erg cm⁻² s⁻¹ is obtained in 10–1000 keV band. This makes the burst the third brightest GRB observed by Fermi to date in terms of energy flux. The Fermi-Large Area Telescope (LAT; Atwood et al. 2009) detected emission in the LAT Low Energy (LLE) data (30–100 MeV), starting at T₀+116 s and emission above 100 MeV starting at T₀+130 s (Arimoto et al. 2016), with LAT emission extending up to ∼2000 s beyond the duration of GBM emission. AstroSat-CZTI also detected the burst for a duration T90 of 42 s (Bhalerao et al. 2016) and captured only the main episode of the burst. GRB 160821A was incident on CZTI from the direction, θ = 1562 and ϕ = 592, thus coming through the side veto detector. Polarization measurement was attempted in the energy range 100–300 keV using ∼2549 detected Compton events (Chattopadhyay et al. 2014, 2017). In Figure 1(a), a composite 1 s binned light curve obtained from various detectors on board Fermi, AstroSat, and Swift satellites are shown. The present work focuses only on the results of the study of the main episode of the burst observed by both Fermi and AstroSat.

Zoom In Zoom Out Reset image size

The measurement of polarization is obtained from the azimuthal distribution of Compton-scattered photons, which lie preferentially in a direction orthogonal to the incident electric field vector (Covino & Gotz 2016; McConnell 2017). The azimuthal distribution is fitted with the cosine function of the form

$$\begin{eqnarray}&&C(\eta )=A\,\cos (2(\eta -{\phi }_{0}+\pi /2))+B\end{eqnarray} \tag{ 1 }$$

where ϕ₀ is the polarization angle (PA) of the incident photons as measured in the CZTI detector plane, A/B is the modulation factor (μ), and η is the azimuthal angle (see also Zhang et al. 2019). The measured detector plane PA is converted into the celestial/sky reference frame by taking into account the satellite orientation at the time of observation, and these values are reported throughout the Letter unless otherwise mentioned. The polarization fraction (PF) is calculated by normalizing the modulation with μ₁₀₀ (the modulation factor obtained for 100% polarized emission coming from the same direction of the GRB with the same spectrum) for the respective detector PA.

A low polarization ($\mathrm{PF}={21}_{-19}^{+24} \%$) was found at 90% (2.2σ) confidence level for two parameters of interest, when the entire pulse constituting a time interval of 115–155 s was analyzed (Figure 1(b)). All the errors quoted for polarization measurements in this Letter are at a 68% confidence interval for two parameters of interest including the systematic errors (Appendix A), unless otherwise mentioned. We find this result to be in agreement with the recent polarization observations of GRBs made by POLAR (Zhang et al. 2019), where they find the time-integrated emission of GRBs to possess a low PF. They suspect that such low polarization could be due to a change in PA within the pulses/across different pulses.

The evolution of the light curves observed by Fermi in different energy bands were characterized by studying the ratio of photon counts observed in high and low energies, parameterized as the hardness ratio (HR); see the uppermost and second from the top sections in Figure 1(c). We found that the emission above 30 MeV changed distinctly with respect to the rest of the burst after T₀+129 s and T₀+140 s. In addition to this, a fine time-resolved polarization analysis of the main episode was conducted. PF and PA for successive 10 s intervals shifted by 2.5 s were measured, thereby studying the evolutionary trend in PF and PA (Figure 1(c)). Such a methodology including overlapping time intervals was adopted because of the limited number of photons available in the smaller time intervals. Therefore, during the analysis we tried to constrain the PF such that at least the lower limit of 68% confidence interval of one parameter of interest was greater than zero. This is because if the time interval is unpolarized i.e., PF is consistent with zero, then its PA has no physical meaning. A change in PA was observed to occur twice as the burst transited from its rise to peak and later to its decay phases while PF was greater than zero.

We note that during the times when the PA angle makes a change, a decrease in PF is expected. Thus, based on the clear change in PA where we could constrain the polarization at a higher significance (i.e., the lower limit of a 99% confidence interval of two parameters of interest of PF are greater than zero) and the observed change in HR at high energies, we performed a relatively coarser time-resolved polarization analysis of the GRB by dividing the main episode into three non-overlapping time intervals: 115–129 s, 131–139 s, and 142–155 s, which correspond to the rise, peak, and decay phase of the burst, respectively.

2.1. Results

During the first time interval, the burst emission has a PF of ${71}_{-41}^{+29} \%$ and a PA of ${110}_{-15}^{+14\,\circ }$. In the second interval, a $\mathrm{PF}\,={58}_{-30}^{+29} \%$ with a corresponding $\mathrm{PA}={31}_{-10}^{+12\,\circ }$ and in the third interval, a $\mathrm{PF}={61}_{-46}^{+39} \%$ with a corresponding $\mathrm{PA}={110}_{-26}^{+25\,\circ }$, are determined (Figure 1(c)). By performing Monte Carlo simulations (Appendix B) of the data set of each interval, the posterior distributions of PF and PA of these intervals were also obtained. Intervals 1, 2, and 3 were found to be polarized at confidence levels of 99.8% (3.5σ), 99.97% (4σ), and 99.1% (3.1σ), respectively, for two parameters of interest (Figure 2). As the burst transits from its rise to the peak phase and then into its decay phase, the PA shifts by $81^\circ \pm 13^\circ$ and $80^\circ \pm 19^\circ$, respectively, (Figure 1(c)). The statistical significance of the change in PAs, Δ PA₁, and Δ PA₂ are determined at 3.5σ and 3.1σ, respectively, which are the minimum of the obtained statistical significance of the two intervals to be polarized. Other cases of varying polarization that were reported earlier are GRB 041219A (Götz et al. 2009), GRB 100826A (Yonetoku et al. 2011), and GRB 170114A (Zhang et al. 2019) observed by INTEGRAL, GAP, and POLAR, respectively. Consequently, the low PF found across the burst can now be understood as an artifact of the temporal change of PA occurring within the burst. The results of the polarization analysis of the three time intervals are listed in Table 1.

Zoom In Zoom Out Reset image size

Table 1.  Results of the Temporal Polarization Analysis

Time Intervals (s)	115–129	131–139	142–155
Compton events	896	1124	523
Background events (s−1)	20.7 ± 4.5	20.7 ± 4.5	20.7 ± 4.5
MDP (99% confidence)	40%	32%	57%
μa	0.289	0.250	0.243
μ100	0.409	0.435	0.403
PFa	${71}_{-41}^{+29} \%$	${58}_{-30}^{+29} \%$	${61}_{-46}^{+39} \%$
PAa	${110}_{-15}^{+14\,\circ }$	${31}_{-10}^{+12\,\circ }$	${110}_{-26}^{+25\,\circ }$
Confidence level	99.8%	99.97%	99.1%

Note.

^aThe average values of μ, PF and PA from the posterior distribution are reported.

Download table as:  ASCIITypeset image

In order to estimate the average PF across the burst, (i) the first and the third intervals were combined because they had nearly same PAs (fourth section from the top of Figure 1(c)); and (ii) Monte Carlo simulations involving combined fits of this new interval and the second interval with the cosine function were performed. The PF related parameters A and B of the cosine functions were linked across the intervals, while the PAs were kept free. Thus, we found the average PF and the PAs for the new interval and the second interval. In Figure 3, the posterior distributions of the average PF across the burst and the corresponding change in PA estimated by taking the difference of the PAs of the new interval and interval 2 are shown. The average PF across the burst is estimated to be ${66}_{-27}^{+26} \%$ at $99.99992 \% \,(5.3\sigma )$ confidence for two parameters of interest as shown in Figure 3. Also, we note that the change in PA estimate (${80}_{-18}^{+17\,\circ }$) is consistent with the average of the change in PAs that were found occurring within the burst. Previously, such a statistically significant polarization was reported for the burst GRB 021206 by Coburn & Boggs (2003); however, the claim was contested by the analyses performed by Rutledge & Fox (2004) and, subsequently, Wigger et al. (2004). Recently, POLAR found that the time-integrated emission of five bright GRBs possess the most probable PF, between 4% and 11% at a confidence level of 99.9% (Zhang et al. 2019). To date no other polarization measurement of GRBs reported by BATSE, INTEGRAL, GAP, AstroSat, or POLAR have obtained a statistical significance greater than ∼99.9% (Covino & Gotz 2016; McConnell 2017; Zhang et al. 2019).

Zoom In Zoom Out Reset image size

Traditionally, the GRB prompt emission spectrum is modeled using the phenomenological Band function¹² (Band et al. 1993). The time-resolved spectral analysis of the main episode of the burst, however, shows significant deviations from the pure Band function in the brightest bins (Preece et al. 2014; Vianello et al. 2018). The deviation in lower energies is modeled by adding a blackbody (BB) function at kT ∼ 30 keV and, at higher energies by adding a cutoff at E_c ∼ 2–50 MeV (Appendix C). Thus, the spectrum is best modeled using a BB + Band × Highecut (Figure 4(a)), where the BB can be related to the thermal emission (a relic of the dense fireball that formed at the central engine after the explosion), and the rest to the non-thermal emission coming from the optically thin region of the outflow (Guiriec et al. 2011; Axelsson et al. 2012; Iyyani et al. 2013, 2016; Burgess et al. 2014; Vianello et al. 2018). The evolution of the spectral-fit parameters are shown in Figure 4(b). We note that the low-energy part of the spectrum characterized by the low-energy power-law index, α and the spectral peak, E_p are nearly steady throughout at ∼−0.97 and 800 keV, respectively. However, the high-energy part of the spectrum characterized by the high-energy power-law index β and cutoff E_c vary significantly such that β decreases, whereas E_c increases with time. After T₀+140 s, these trends are reversed. This spectral behavior is consistent with the trend observed in HR reported above.

Zoom In Zoom Out Reset image size

GRBs with long durations, >2 s, are associated with the death of massive stars. A highly collimated outflow (jet) of opening angle ${\theta }_{j}$, moving at a relativistic speed (parameterized by Lorentz factor,¹³ Γ) is produced after the core of the massive star collapses to a black hole (or a magnetar) and begins to accrete the surrounding stellar matter. The radiation emitted from this relativistic outflow is beamed toward the observer within a cone of 1/Γ, which is the visible region around the line of sight. This is referred to as the relativistic aberration/beaming.

In a classical fireball model (Mészáros 2006; Pe'er 2015; Iyyani 2018), where the outflow is non-magnetized, the non-thermal emission is generally expected to be produced in shocks created in the optically thin region above the photosphere from where the thermal emission is expected. Energetic electrons produced in the shocks then cool by processes such as synchrotron emission in random magnetic fields generated in the shocks (Ghisellini & Lazzati 1999), or inverse Compton scattering (Lazzati et al. 2004), both being inherently locally axisymmetric around the outflow's velocity vector. For a fixed viewing angle,¹⁴ ${\theta }_{v}\gt 0$, in an axially symmetric jet, the polarization vector integrated over the spatially unresolved emitting region should point either perpendicular to or in the plane containing the axis of the jet and the line of sight of the observer. Thus, the PA is expected to change by 90° when the width of the jet parameterized by ${\rm{\Gamma }}\,{\theta }_{j}$ changes. In such a case, when the orientation of the polarization vector becomes perpendicular to the observer plane, the PF is expected to be low (<10%; Ghisellini & Lazzati 1999; Granot 2003; Toma et al. 2009). However, here we observe that during the entire burst, the PF values are >15%. Thus, a change in width of the visible portion of the jet cannot account for the observed change in PA. Inherently axisymmetric emission models referred to above are ruled out.

The total radiative energy¹⁵ released in the prompt emission of GRB 160821A is estimated to be E_γ,iso = 6.9 × 10⁵³ (3.6 × 10⁵⁵) erg, for a redshift, z = 0.4 (2), if the radiation were isotropic. These are among the highest known values for long GRBs (Racusin et al. 2011). Such a high E_γ,iso suggests that the emission is strongly collimated, with the jet pointing toward the observer such that the line of sight lies within the jet cone or just outside the edge of the jet (${\theta }_{j}+1/{\rm{\Gamma }}$). In the above scenario, the strong observed polarization can be explained only by synchrotron emission produced in magnetic fields that are highly ordered within the viewing cone of 1/Γ.

If we assume that the observed burst emission is a single emission episode, then it is difficult to explain, using any known physical model, the observed high polarization along with a change in PA. On the other hand, the burst emission consists of multiple emission episodes and, depending on the dominance of the synchrotron emissions coming from the different regions, a change in PA can happen with time (also see Lazzati & Begelman 2009).

Thus, for the first time conclusive evidence of high and varying linear polarization is detected in a GRB. The observations are extremely constraining and challenging within the framework of currently proposed physical models for GRBs. This motivates further research into the development of a physical scenario that can explain the observations self consistently.

We would like to thank Dr. Christoffer Lundman, Prof. Pawan Kumar, and Dr. Santosh Roy for enlightening discussions. This publication uses data from the AstroSat mission of the Indian Space Research Organisation (ISRO), archived at the Indian Space Science Data Centre (ISSDC). CZT-Imager is built by a consortium of institutes across India, including the Tata Institute of Fundamental Research (TIFR), Mumbai, the Vikram Sarabhai Space Centre, Thiruvananthapuram, ISRO Satellite Centre (ISAC), Bengaluru, Inter University Centre for Astronomy and Astrophysics, Pune, Physical Research Laboratory, Ahmedabad, Space Application Centre, Ahmedabad. This research has made use of Fermi data obtained through High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center. The Geant4 simulations for this paper were performed using the HPC resources at The Inter-University Centre for Astronomy and Astrophysics (IUCAA) and Physical Research Laboratory (PRL).

CZTI on board AstroSat is an X-ray spectroscopic instrument and is experimentally verified for polarization measurement capability in 100–400 keV energy range for on-axis sources (Vadawale et al. 2015). CZTI functions as a wide field monitor at energies >100 keV because the CZTI collimators and other supporting structures become largely transparent above this energy. This enables it to detect GRBs and measure their polarization. In order to calculate the modulation factor for off-axis sources, a mass model of AstroSat was developed in Geant4 (version 4.10.02.p02). There are several possible sources of systematic errors in the measurement of polarization (Chattopadhyay et al. 2017). Below we list these sources and the error estimates are quoted for the case of GRB 160821A:

(i) An uncertainty can be induced in the observed μ due to different photon propagation paths through the spacecraft structures. This uncertainty is estimated by conducting $\sim {10}^{4}$ Geant4 simulations of this burst with the same spectra and incident direction to produce the same number of observed Compton events. The uncertainty on μ is thus estimated to be ∼8%–16% according to the different number of Compton events. (ii) The selection of background is also expected to induce some systematic errors. This was investigated by taking both pre- and post-GRB background independently, as well as in combination, in order to determine the modulation amplitude. The uncertainty on μ due to this is found to be <1%. (iii) The effect of localization uncertainty is studied. GRB 160821A is localized at R.A. = 171.25 and decl. = +42.33 at ±3' accuracy (Markwardt et al. 2016). The contribution of localization uncertainty on the obtained modulation amplitude and PA is found to <1%. (iv) The uncertainty associated with the spectral model of the GRB is investigated by varying the power-law index within its estimated 1 sigma error and we find the associated uncertainty on the observed modulation amplitude to be <1%. (v) Another source of systematic error could be the unequal quantum efficiency of the CZTI pixels. The relative pixel efficiency across the CZTI plane is found to vary within ≤5%, which produces negligible false-modulation amplitude.

In the μ₁₀₀ estimations, the statistical error is quite small as the simulations are done for a large number (10⁸) of incoming photons. The systematic errors are those associated with the sources (iii) and (iv), while the uncertainty associated with source (i) is estimated to be ≤1%. The value of μ₁₀₀ strongly depends on the fitted PA. The uncertainty associated with this in the estimation of PF is taken care of in the Monte Carlo process to obtain the posterior distribution of PF (described in Appendix B), where for each fitted PA, the corresponding μ₁₀₀ is used.

All of these uncertainties are propagated into the reported values of the limits of the 68% confidence interval of two parameters of interest, namely the observed PF and angle.

The polarization signature in the GRB is estimated through the non-uniform azimuthal distribution of GRB counts. We calculate the normalized counts for eight bins whose mid-values correspond to azimuthal angles: 0°, 45°, 90°, 135°, 180°, 225°, 270°, and 315°. These angles are estimated in an anti-clockwise direction when CZTI is viewed from top. The corrected modulation curves are fitted for large number of iterations (10 million) using the Monte Carlo method for estimating the modulation factor and PA as follows.

• (a)  
  For each azimuthal angle (η_i), and its corresponding normalized counts (${y}_{{\rm{i}}}$) and error (y_err,i), we pick a random value (y_ran,i), from a normal distribution, which is defined by the mean value, ${\mu }_{\mathrm{mean}}={y}_{{\rm{i}}}$ and the standard deviation, ${\sigma }_{\mathrm{SD}}={y}_{\mathrm{err},{\rm{i}}}$. A normal distribution is assumed because, in the case of GRB 160821A, each azimuthal angular bin has over 20 Compton events (Cash 1979).
• (b)  
  The A, B, and ϕ₀ parameters of the fitting cosine function given in Equation (1) are estimated for each set of randomly drawn y_ran,i (where i = 1–8) values via least square curve fitting method. For the fitted PA, the corresponding interpolated value of μ₁₀₀, from a table of μ₁₀₀ values that were generated for a discrete grid of PAs via the Geant4 simulations of 100% polarized emission for the same GRB spectra and incoming direction, is chosen. Thereby, estimating the PF values.
• (c)  
  The above two steps are repeated for a large number of times (10⁷), thereby obtaining the various required likelihood distributions of PF and PA.
• (d)  
  The obtained likelihood is then filtered through the prior condition so that the PF has to lie between 0% and 100%. The simulation runs that satisfied the prior condition were then used to generate the posterior distribution.
• (e)  
  Finally, the respective 2D histograms of posterior distributions of PA and PF are made.

For analysis of the Fermi GBM data, three NaI detectors with the highest count rates were chosen, namely n6, n7 and n9. These detectors observed the GRB at an off-axis angle less than 40°. In addition, the BGO detector BGO 1, which had the strongest detection, was chosen for analysis. LAT data (>100 MeV) belonging to the P8_TRANSIENT020E class and its corresponding instrument response were used, whereas the LAT-LLE spectra (30–130 MeV) were obtained using the same method as that for GBM data. The spectrum for each Fermi detector was extracted using the software Fermi Burst Analysis GUI v 02-01-00p1 (gtburst3¹⁶ ). The background was obtained by fitting a polynomial function to the data regions before and after the GRB for time intervals, T₀ − 410 s to T₀ − 10 s and T₀ + 210 s to T₀ + 250 s, respectively.

A joint spectral analysis of Fermi GBM along with LAT-LLE and LAT data was carried out in Xspec (Arnaud 1996) 12.9.0n software and Pg-stat (Arnaud 2013) statistic was used. For NaI data, energies between 30 and 40 keV corresponding to iodine K-edge and extreme edges such as energies below 10 keV and above 850 keV were removed. In case of BGO, LAT-LLE, and LAT, data between energies 300 keV–10 MeV, 30 MeV–130 MeV, and 100 MeV–5 GeV were used, respectively.

For time-resolved spectroscopy of the burst, time intervals were defined using Bayesian Blocks algorithm on the detector (n6) with largest number of counts. In the tails of the emission episode, due to low signal-to-noise ratio, certain blocks were combined to get a larger time interval so that good spectral constraints could be obtained.

For estimating the effective area correction¹⁷ for NaI and LAT detectors, the bright time bins were simultaneously fitted with the best-fit model; i.e., BB + Band × Highecut multiplied with a constant of normalization, while the constant of normalization of the BGO 1 detector was fixed to unity. The relative normalization constants ${0.97}_{-0.01}^{+0.01}$ for n6, ${0.92}_{-0.01}^{+0.01}$ for n7, ${0.94}_{-0.01}^{+0.01}$, and for n9 and ${0.84}_{-0.06}^{+0.06}$ for LAT were obtained. Throughout the time-resolved spectral analysis, the effective area corrections for the detectors were frozen to these obtained values.

C.1. Statistical Significance Test

We conduct Monte Carlo simulations in order to ascertain the statistical significance of the deviations in the spectrum from a pure Band function, which have been modeled using a BB and cutoff at lower and higher energies, respectively (see the fit residuals in Figure 5). Because this process is computationally intensive, we present here the simulation study done for the brightest time bin (134.59–135.71 s) only, and adopt the obtained Δ Pg-stat distribution as a reference to assess the statistical significance of these components in other time bins.

• (a)  
  BB: the difference in statistics, i.e., Δ Pg-stat, obtained for the model BB + Band ×  Highecut (f_BHec) from f_BHec is 24.5. We assume the model f_BHec with the best-fit model parameters as the null hypothesis (H0) and generate nearly 40,000 realizations and its corresponding background spectra using the fakeit command in Xspec. Each of these realizations are then fit with the null hypothesis model f_BHec and the alternate hypothesis (H1) model BB + f_BHec, and the respective Δ Pg-stat values are recorded. The probability of observing any Δ Pg-stat of >24.5 is found to be 10^−4.3 (Figure 6(a)), which corresponds to a significance level of ∼4σ.
• (b)  
  Highecut: the Δ Pg-stat obtained for the model BB +f_BHec from BB + Band is 201. In this case, we assume the model BB + Band with the best-fit model parameters as the H0 and generate nearly 54,000 realizations and its corresponding background spectra. The model BB + f_BHec is the H1, and the respective complementary cumulative distribution of Δ Pg-stat is obtained. The probability to obtain Δ Pg-stat = 14 is ${10}^{-4.5}$ (Figure 6(b)) which indicates that the probability to obtain any Δ Pg-stat > 201 is ≪10^−4.5, which corresponds to a significance level of >4.2σ. Thus, the addition of Highecut to the Band function is highly significant.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size