A Generalized Framework for Rapid Polarization Response Generation in Wide Field-of-View Compton Polarimeters

Introduction

The radiation mechanism and the magnetic field structure of gamma-ray burst (GRB) prompt emission remain unsettled, and polarization is the most direct observable that separates the candidate models ([20]). That emission peaks at hard X-ray and gamma-ray energies, from tens of keV to about 1 MeV ([46]), the band this work addresses throughout. The last two decades have seen a rapid growth of polarimetry in high energy astrophysics. This was initiated by a range of polarization measurements of GRBs using instruments like RHESSI ([15, 41]), BATSE ([48]) and those on INTEGRAL ([35, 23, 21, 22]) and followed by several measurements of the Crab and Cygnus X-1 at gamma-ray energies by INTEGRAL SPI and IBIS ([16, 18]). This was followed by a wealth of measurements at X-ray energies thanks to instruments like PoGO+, PoGO+, POLAR-Light, XL-Calibur and most recently IXPE ([47]). In addition, instruments like GAP ([49] (Yonetoku et al. 2011a)), POLAR and AstroSAT CZTI ([14, 13]) have, to date, produced tens of GRB polarization measurements, thereby paving the way for future large scale GRB polarimeters such as LEAP, Daksha ([7]), COSI and finally POLAR-2 which will contain both a gamma-ray polarimeter (the High energy Polarization Detector (HPD)) and the first dedicated X-ray polarimeter for transients (the Low energy Polarization Detector (LPD)) ([30]). In parallel to this growth in measurements, progress has been made regarding the analysis in this field although standard, well tested analysis methods remain elusive. While the initial polarization measurements suffered from significant issues in the analysis, an overview of which is discussed in Gill et al. [20], a first step towards standardized analysis methods has been produced [11] and the first public polarization analysis tools have been published in recent years (see, e.g., Kole et al. [32]). Because the dominant photo-detector interaction is different at different energies, photo-absorption below a few tens of keV and Compton scattering above, the analysis approaches for soft $X$-ray and hard $X$-ray/$\gamma$-ray polarimetry have evolved largely independently. The analysis at $X$-ray energies is now well established and makes use of Stokes parameters. Most $X$-ray polarimeters are pointing instruments observing on-axis sources, for which the Stokes formalism is naturally suited; the off-axis complications central to this work do not arise there. While the use of Stokes parameters for polarimetry has long been used in optical and radio astronomy, at $X$-ray energies the discrete nature of the photons requires the calculation of Stokes parameters for individual photons. A method to deal with this in a mathematically sound manner was presented in Kislat et al. [26] and has since been followed in the majority of X-ray analyses.

At hard X-ray and gamma-ray energies a more wide range of analysis methods has been employed. Initial analysis typically relied on the production of so called modulation curves, which are defined in more detail in the next section. These modulation curves are fitted using single harmonic function with a period of $180^\circ$. Examples of this type of analysis can be found in the AstroSat CZTI (Chattopadhyay et al. [14]; Chand et al. [12]; Sharma et al. [42]) and GAP (Yonetoku et al. [49] (2011b)) GRB analysis. While this method can work well for on-axis observations, as is the case for the analyses presented by the INTEGRAL teams (Bošnjak et al. [8]; Götz et al. [23, 21]), for off-axis polarization observations this method is not correct. The off-axis nature, although appearing trivial at first sight, alters the polarization signal from a simple, single component harmonic, into a complex harmonic function as first discussed in detail in Muleri (2014). The work in Muleri (2014) covers both the low energy X-ray range, where photo-absorption is used to measure the polarization of incoming photons, as well as the hard X-ray range where Compton scattering is used. For the Compton scattering case, on which we will focus here, the dependency of the azimuthal Compton scattering angle on the polarization of the incoming photons is presented in detail. In particular the dependence of the angular distribution on the incoming direction of the photon is derived in general, and the resulting integrals are evaluated in a first-order expansion in the photon energy. That expansion holds up to several tens of keV, the range of X-ray polarimeters, and its approximations fail at the energies of GRB prompt emission. As the field of polarimetry is growing especially in the gamma-ray energy range for wide field of view instruments, we revisit the off-axis case here and evaluate the same description without approximation across the energy range defined above.

In addition, the method is tested on existing data, and the effect of some of the simplifications made in previous analyses are also presented. This regards specifically the standard method where the polarization is retrieved through fitting a harmonic function. In recent years, more advanced methods, such as using forward folding in spectro-polarimetry, have been developed (Burgess et al. [11]; Kole et al. [32]). Such methods rely on comparing the measured modulation curves with simulated ones and retrieve the polarization and spectral parameters through minimization. As a result, they do not rely on fitting the harmonic functions, but rather on a comparison of the harmonic functions in the data with those present in the simulated detector response, so it matters that these complex harmonic functions are well represented in the simulated responses, and we compare them in detail here. This work delivers the off-axis azimuthal scattering angle distribution evaluated numerically and without approximation across that band, for any incidence and polarization angles, verified against GEANT4 simulations of an ideal polarimeter. It adds the effect of the detection threshold on the measured distribution, which every real instrument imposes, with POLAR as the case study. On that description we build the full polarization response of an instrument in two modes that differ in their input: mode 1 takes only the unpolarized simulation at each energy, and mode 2 adds two fully polarized simulations per energy at distinct polarization angles to carry the polarization-dependent effects of the instrument. Either mode needs at least 95 percent fewer Monte Carlo simulations than a full response. The released software builds either response for any planar Compton polarimeter from that small set of its own simulations.

The paper is organized as follows: Section 2 reviews the Compton scattering formalism for polarized photons, sets the notation used throughout, and states the off-axis problem. Section 3 develops the generalized off-axis framework based on numerical integration of the transformed Klein–Nishina cross-section, which holds across the band defined above, and describes the ideal simulations used to verify it; the low-energy expansion it supersedes is recorded in Appendix B. Section 4 extends the framework to realistic instruments through the error-function threshold gate (derived in Appendix C) and the treatment of instrumental effects, and introduces the two response modes of this work. Section 5 describes the current approaches to measuring polarization based on the azimuthal scattering angle distribution and full GEANT4 simulated responses of the instrument [2], and how each response mode is constructed: mode 1 from the unpolarized simulation alone, and mode 2 from the unpolarized simulation and fully polarized simulations at two polarization angles per energy. In § 6, the analytical predictions are compared against the GEANT4-generated responses of an ideal polarimeter and of the POLAR instrument [40] for a range of incidence angles and spectra, and the two modes are compared against the full polarized grids. § 7 then applies both response modes to realistic polarization analyses: the real POLAR observation of GRB 170114A and a grid of synthetic GRBs with known input polarization parameters, quantifying the regime in which the mode 1 assumption holds and showing that mode 2 reproduces the full simulated response at the posterior level. Finally, § 8 summarizes the results and discusses future implementations, enabled by the reduction in computational cost. The software implementing both modes is described in Appendix § 13.

Polarimetry using Compton Scattering

The On-Axis Ideal Case

Photons in the hard X-ray and gamma-ray band of § 1, from tens of keV to about 1 MeV can interact with matter through several processes. In the low-$Z$ materials, typically used in Compton polarimeters, the two dominant processes are photoelectric absorption, in which the photon is fully absorbed and a bound electron is ejected, and Compton scattering, in which the photon scatters inelastically off an electron, transferring part of its energy. Photons can also undergo Rayleigh scattering, an elastic process in which the photon is redirected by the atomic cloud without any energy transfer. Since Rayleigh scattering deposits no energy at the interaction site, it cannot be detected in a Compton polarimeter and is therefore irrelevant to the measurement. In the energy range considered in this work, Compton scattering is the dominant interaction and provides the physical basis for polarization measurements [6] (Knoll 2010).

Compton scattering affects the energy and momentum of both the photon and the emitted electron during interaction. The Compton scattering equation can be derived by applying conservation of relativistic energy and momentum for the collision between a photon and a free electron. This yields the Compton kinematics equation:

Here $E_{\gamma}$ is the energy of the incoming photon, $E'_{\gamma}$ is the scattered photon energy, and $\theta$ is the Compton scattering angle, defined as the angle between the incident momentum $\vec{k}$ and the scattered momentum $\vec{k}'$. It is labeled $\theta_C$ in the following figures; the two symbols denote the same angle. Because the incident photon travels along the $-z_{\gamma}$ axis, $\theta$ is measured from $-z_{\gamma}$ by convention. For readers using standard spherical coordinates, it relates to the polar angle $\theta_s$ (measured from $+z_{\gamma}$) by $\theta= \pi-\theta_s$, equivalently $\cos\theta=-\cos\theta_s$. This convention is used consistently in all text, equations, and figures.

A full three dimensional description of the scattering process is given by the Klein-Nishina cross-section [27]. It provides a fully relativistic treatment of the scattering process, derived from the Dirac equation and Quantum Electrodynamics (QED), which naturally incorporates the polarization state of the incident photon. Consequently, the polarization state of the incident photon directly affects the kinematics of the scattering cross-section. Because the incident photon’s oscillating electric field preferentially accelerates the target electron along the axis of polarization, the resulting scattering probability loses its azimuthal symmetry about the incident photon direction.

The Klein Nishina differential cross section is given by:

Here $r_e$ is the classical electron radius, ($r_e$ is given by $r_e=\frac{\alpha\hbar}{m_ec^2}\approx2.82\times10^{-15}$ m. Here $\alpha$ is the fine structure constant, $\hbar$ is the reduced Planck’s constant, $m_e$ is the electron’s rest mass and $c$ is the speed of light.) and $\theta$ is the same as in equation (1). The angle $\phi$ is the azimuthal scattering angle defined as the angle between incident photon polarization vector $\vec{\xi}$ and the scattered photon direction in a plane perpendicular to incident photon direction as illustrated in Figure 1.

In a physical polarimetry experiment, measurements are recorded relative to a macroscopic instrument frame, typically a large area planar detector defined by an $xy$−plane. Figure 1 shows the simplified case where the intrinsic ”photon frame” is perfectly aligned with this physical detector frame. Here, the photon’s incident momentum direction $\vec{k}$ is aligned with the detector’s negative $z$−axis, and its polarization vector $\vec{\xi}$ is aligned perfectly with the detector’s positive $x$−axis. The scattering angle $\phi$ in this configuration identically matches the azimuthal angle of the vector $\vec{k}^{\prime}$ in the standard spherical coordinate system. However, in practice, a physical detector will measure a more generalized case.

The Klein-Nishina cross-section dictates that polarized photons preferentially scatter in directions perpendicular to their incident polarization vector. For the specific configuration shown in Figure 1, this corresponds to the angles $\phi= \pi/2$, $3\pi/2$, where Equation (2) reaches its maximum. This implies that the azimuthal scattering distribution of the photons contains the necessary information to determine the bulk polarization of the incoming beam. By contrast, an incident unpolarized beam lacks a defined polarization vector and will scatter with equal probability in all azimuthal directions.

In a measurement with a large-area planar detector the polarization vector $\vec{\xi}$ lies at some angle in the detector’s $xy$−plane ($\vec{\xi}$ must be orthogonal to $\vec{k}$, as required by Maxwell’s equations); that angle from the detector’s positive $x$−axis is the polarization angle $\phi_o$ (Figure 2). This is a local polarization angle defined in the detector frame. It should not be confused with the polarization angle defined in the source incidence frame, nor with the IAU North-to-East convention used to report astrophysical polarization angles; the relation between them is fixed by the source geometry introduced in § 3. The Klein-Nishina angle $\phi$ is then measured in the $xy$-plane between the projections of $\vec{k}^{\prime}$ and $\vec{\xi}$, and the azimuthal scattering angle the detector records, $\eta$, measured from the positive $x$-axis, is the rigid rotation $\eta= \phi+ \phi_o$ (to ensure mathematical consistency across all reference frames, all azimuthal angles, $\phi$, $\phi_o$ and $\eta$, are strictly defined as directed angles measured counter-clockwise in the $xy$-plane according to the right-hand rule); strict on-axis incidence is common in the laboratory and rare in orbit, and the off-axis case is treated in the following sections.

A polarimeter records $\eta$ for each scattered photon of a beam incident along the negative $z$-axis. The resulting Azimuthal Scattering Angle Distribution (ASAD) $\mathcal{N}(\eta)$ follows a $\cos^2(\eta-\phi_o)$ modulation, a Fourier second harmonic phase-shifted by $\phi_o$, since the Klein-Nishina cross-section depends on $\cos^2\phi$; its first minimum determines the polarization angle. The polar angle of $\vec{k}^{\prime}$ here is $\pi-\theta$ where $\theta$ is the Compton scattering angle.

The expected shape of the ASAD resulting from a polarized beam incident along the negative $z$-axis of a detector can be seen through simple manipulations of equation (2). By considering only the dependence on the azimuthal scattering angle $\phi=\eta-\phi_o$, we can write

The free parameters of this last equation ($A$, $B$, $\phi_o$) quantify the polarization signature. The parameter B represents the average count rate (offset), and A represents the amplitude of the modulation. A typical ”ideal” shape of the ASAD in terms of these parameters is shown in Figure 3.

The measured ASAD, $\mathcal{N}(\eta)$, is the distribution of recorded scattered photons as a function of the azimuthal scattering angle $\eta$. To quantify the strength of the polarization signal, we define the modulation factor $\mu$ from the maximum and minimum of this distribution, $\mathcal{N}_{\max}$ and $\mathcal{N}_{\min}$:

We can find $N_{max}$ from Equation (3) as $N_{max}=B+A$ and similarly $N_{min}=B-A$. This definition can then be simplified to $\mu=A/B$, which is simply the ratio of the amplitude of the variation to the average value of the ASAD. For a fully polarized beam measured by a perfect detector, this modulation reaches a theoretical maximum known as $\mu_{\max}(E,\theta)$. (It can be referred to as $\mu_{100}$, but the literature reserves $\mu_{100}$ to refer to the measured modulation of a fully polarized beam simulation as a figure-of-merit for the instrument performance, so we label the theoretical limit of modulation as $\mu_{\max}(E,\theta)$. It is a fundamental limit that can be derived directly by comparing Equation (2) (after expanding the $\cos^2 \phi$ term and regrouping) to Equation (3) to yield the amplitude $A$ and the offset $B$, giving:

Real instruments cannot measure this theoretical limit due to instrumental effects that can be captured in simulations as a part of the response in $(E,\theta)$ of that instrument. We define $\operatorname{Response}(E,\theta)$ to denote the instrument’s detection efficiency for a Compton event of incident energy $E$ and scattering angle $\theta$, i.e. the number of recorded coincidences per incident photon. We can then define the maximum measurable modulation factor for an instrument as the weighted average of the theoretical $\mu_{\max}(E,\theta)$ integrated over the response as:

The measured polarization level is then calculated as $\Pi= \mu/\mu_{100}$. The formalism above applies only to on-axis sources, and real wide-FOV instruments observe sources at arbitrary incidence angles.

The Off-Axis Reality (The Problem)

It is evident from equation (2) and 3 that the ASAD (after correcting for instrumental effects) is expected to have only a $180^\circ$ component. This has been the assumption in many, but not all, GRB polarization measurements to date. Those analyses characterized the signal by a single $180^\circ$ modulation and its on-axis modulation factor $\mu_{100}$ regardless of the source position. Examples of this can be found in [14, 13]. This assumption is, however, false when the source location is not along the viewing axis of the polarimeter. An X-ray polarimetry study proved that ASADs from such off-axis sources should include more components in the signal (Muleri 2014). This was also evident from ground calibrations of the POLAR instrument, where additional $360^\circ$ and $90^\circ$ components in the ASAD for off-axis source were measured by Kole et al. [31].

Expanding the ASADs in terms of Fourier harmonics is instructive in showing what components exist in the signal. To illustrate its importance we generated ASADs for different zenith positions $\delta\in[0^\circ, 80^\circ]$ for energies $50$, $100$, and $300\ \mathrm{keV}$ using our model which will be presented in the next section. And we show in Figure 4, the relative amplitudes $a_k$ for the each of the first four harmonics for each ASAD. The relative amplitudes $a_k$ are measured with respect to the constant term in the Fourier expansion. It can be seen that for the polarized ASADs, the $180^\circ$ component is maximum only when we have on-axis incidence $\delta= 0$. And the unpolarized ASADs have harmonics and are not just a constant value once photons are incident from an off-axis angle. In the typical context of classical polarization analysis, as presented in § 5.1, the polarized ASAD is divided by an unpolarized ASAD which is argued to remove the geometric effects in the signal, but as shown in the right column of Figure 4, it does not completely remove the geometric effects due to the projection of Klein-Nishina cross-section which are depicted as extra harmonics. Those additional Fourier components contaminate the $\cos^2\phi$ modulation that traditional analysis fits, and Figure 4 quantifies the part of the off-axis signal that a single modulation factor cannot capture. Assuming a single $180^\circ$ modulation component for off-axis sources therefore biases the recovered polarization, which prompted this study.

Generalized Description of Compton Scattering

To quantify the additional harmonics identified in the previous section, we derive a generalized form of the Klein-Nishina equation valid for an arbitrary source direction and polarization angle. The goal is to obtain a description of the azimuthal scattering angle distribution, in the detector frame and as a function of $\eta$, that accounts for the two rotations defining the source location and the rotation defining the polarization angle. The angles $\theta$ and $\phi$ in the Klein-Nishina equation are defined with respect to the photon frame which we denote by $x_\gamma y_\gamma z_\gamma$. When the coordinates of the incident photon frame align with those of the detector $(x_d y_d z_d)$, only the typical second-order harmonic is observed in the ASAD. However, for a general incident direction the two frames would not align and we expect additional harmonics. This is because the photon polarization vector, by definition, lies in the plane perpendicular to the photon momentum vector. Consequently, if a photon is incident off-axis, projection effects will distort the typical $180^\circ$ component.

Coordinate Transformation

Here, we will utilize the spherical coordinate representation in both the photon frame $x_\gamma y_\gamma z_\gamma$ and the detector frame $xyz$. The description relies on describing the direction of the scattered photon momentum vector $\vec{k}$ in terms of spherical coordinates. In the photon frame, $\theta$ and $\phi$ follow the convention established in § 2.1 (Figure 1). In the detector frame, the corresponding transformed polar and azimuthal angles are $\vartheta$ and $\varphi$, and the measured azimuthal angle is $\eta$. All angles used in the following derivation are summarized in Appendix § 10. An example of a simple transformation from the photon frame to the detector is when the photon is rotated around the $z$-axis with angle $\phi_o$ (the polarization angle). This rotation $R_{\phi_o}^{z}$ makes the angle $\eta$ related to $\phi$ through $\eta= \phi_o + \phi$ as shown in Figure 2. For an arbitrary source location i.e., not on-axis, the transformation from the photon frame to the detector frame requires two rotations as follows: for zenith angle $\delta$ through $R_{\delta}^{y}$ and for azimuth angle $\lambda$ through $R_{\lambda}^{z}$.

We first write the differential Compton scattering distribution in the photon frame, $D_{\gamma}^{\mathrm{pol}}(\theta,\phi,\varepsilon)$ and $D_{\gamma}^{\mathrm{unp}}(\theta,\varepsilon)$ the Klein-Nishina cross-section per unit solid angle (Equation (2)). Here the superscript denotes the polarized or unpolarized beam and the subscript $\gamma$ denotes the photon frame. Integrating over the polar angle and normalizing yields the azimuthal distribution $\Phi_{\gamma}^{\mathrm{pol/unp}}(\phi,\varepsilon)$. Finally, applying the frame rotations re-expresses this in the detector’s azimuthal angle, $\Phi(\varphi) \equiv\Phi(\eta)$, which is the quantity a polarimeter measures. We start with equation (2) and use equation (1) to substitute for $E'_{\gamma}/E_{\gamma}$ as follows:

The ”pol” superscript refers to the polarized portion of the beam. And for the unpolarized portion of the beam, $\cos^2\phi$ is replaced by its azimuthal average, $\langle\cos^2\phi\rangle= \frac{1}{2}$, yielding:

Since $\mathcal{D}$ describes the distribution over a solid angle $\mathrm{d}\Omega$, we integrate over the full range of $\theta$ to obtain the azimuthal angular distributions and to normalize this distribution, it is divided by the integration over the entire solid angle $\Omega$, i.e.

This is now a representation of the normalized distribution $\Phi(\phi)\equiv\Phi(\eta)$ i.e., the modulation curve as viewed in the photon frame when it is anti-aligned with the detector frame. Equation (9) on the right shows that the expected unpolarized distribution for an on axis source and would be a constant for all values of azimuthal scattering angle as outlined earlier. However, for a general photon direction we need to transform the photon frame to anti-align with the detector frame. The transformation would consist of a series of rotations that correspond to the source photon’s degrees-of-freedom, i.e., the azimuthal angle $\lambda$, the zenith angle $\delta$ and the polarization angle $\phi_o$ (Figure 5). We now express the direction of the scattered photon in both frames using spherical coordinates. Following the convention established in § 2.1 ($\theta$ measured from $-z_\gamma$), the photon-frame $z$-component of the scattered photon is $z_\gamma=-\cos\theta$, as used in Equation (10). Throughout, a subscript $\gamma$ denotes a quantity expressed in the incident-photon frame. The incident photon frame unit basis vectors are given by:

Similarly, the detector frame basis unit vectors are labeled by:

The rotations are now applied as follows:

After the first rotation by $\lambda$ around the $z$-axis, the coordinate system is now $x'y'z'$. Then the second rotation by $\delta$ around the $y'$-axis, the coordinate system is now $x''y''z''$. And the final rotation by $\phi_o$ around the $z''$-axis, the coordinate system now matches $x_\gamma y_\gamma z_\gamma$.

Accordingly, the rotation matrices are defined as follows:

The goal remains to obtain the expressions in equations (7) and (8) in terms of the detector frame coordinates, i.e. $\mathcal{D}(\vartheta,\varphi)$. So we need to find the transformation of the relevant terms:

The aforementioned terms are to be transformed according to equation (12) to be expressed in terms of the detector frame coordinate system: $\vartheta$, $\varphi$, $\delta$, $\lambda$, $\phi_o$. The transformed terms are then substituted in the KN cross-section in equations (7) and (8) yielding the expected transformed cross-section. To calculate the transformed normalized angular distribution, one would need to evaluate the transformed version of the integrals in 9 such that:

The above-mentioned equations in 15 would produce the expected modulation curves/distribution for polarized and unpolarized distributions in terms of the measured azimuthal angles in the detector frame $\varphi$ for any source parameters ($\varepsilon$, $\delta$, $\lambda$, $\phi_o$). We evaluate these integrals numerically; the first-order analytical approximation of Muleri (2014), valid at low energies, is recorded in Appendix § 11 and compared with the numerical result in § 6.1.

It is important to note here that, these distributions represent the physical aspect/bias of the ASAD without any consideration for the instrumental effects/biases. Some realistic instrumental consideration affect the physical bias of the ASAD as will be covered later. However, it remains true that this approach thus far is blind to the instrumental effects.

Generalized Framework

The low-energy expansion of Appendix § 11 cannot be carried to the energies of GRB prompt emission, where its expansion parameter is no longer small. No closed form replaces it there, because the transformed Klein-Nishina integrands in Equation (15) contain products of trigonometric functions of $(\vartheta,\varphi)$ coupled through the rotation matrix elements $\mathcal{R}_{ij}(\delta,\lambda,\varphi_o)$. In this work we therefore evaluate the integrals numerically at discretized points $\varphi_i$ with no approximation of the cross-section. Throughout, we call this numerically evaluated description the analytical framework, as opposed to the simulated response of an instrument.

We evaluate the inner integral over $\vartheta$ using fixed-order Gauss-Legendre quadrature, whose roots and weights are computed once Virtanen et al. [45] and reused for every evaluation (32 nodes on each interval of the analytic coincidence window of Equation (17) with no gate or with the sharp gate of Equation (19), and 64 nodes on each of four equal panels of $[0,\pi]$ with the error-function gate of Equation (21)). The outer integral over $\varphi$ is performed on a fixed grid (360 points over $[0,2\pi)$), and the normalization in the denominator is the periodic mean of the same grid. The integrands are evaluated as vectorized array operations in NumPy Harris et al. [24]. With this implementation, a full ASAD for a given set of beam parameters $(\epsilon,\delta,\lambda,\varphi_o)$ is computed in about a millisecond on a modern workstation (the 9,150 mode 1 ASADs of a POLAR response in 11 s), which makes large parameter sweeps practical.

This approach handles the full parameter space at once, unlike the low-energy expansion of Appendix § 11, which is limited to the case where $\lambda=\varphi_o=0$. Even though the source azimuthal angle $\lambda$ produces a pure phase shift, as shown in Equation (24) of Muleri (2014), the polarization angle $\varphi_o$ changes the shape of the ASAD in a non-trivial, energy-dependent way.

Verification Setup

To verify the results of our framework, we utilized a physics-based GEANT4 simulation to produce the ideal detector response Agostinelli et al. [2]. A response that captures purely the physics of the scattering process, without instrumental effects as expected by KN cross-section when projected onto a two-dimensional planar detector as in Figure 6. We set up a dense slab made entirely of hydrogen acting as a perfect scatterer, with a density of 100 g/cm$^3$ to maximize the probability of Compton interaction while keeping essentially all interactions within the top layer of the artificial scatterer. We then fired a plane parallel beam of mono-energetic photons at the detector. A schematic of the simulation geometry is shown in Figure 7.

We were able to fire photons from any source direction denoted by zenith $\delta$ and azimuth $\lambda$, polarized with any polarization angle $\varphi_o$ or unpolarized. We used GEANT4’s GetPostStepPoint() to record the photon’s scattering direction $\vec{k}'$ after its first interaction in the detector. Using $\vec{k}'$ and the predefined $xy$-axes of the detector, we were able to determine the azimuthal scattering angle $\varphi$ whose histogram constitutes the ASAD. We filtered out photons that underwent photoelectric absorption as photons had NaN energy after scattering due to being fully absorbed. We also filtered out photons that had no change in energy after interaction as they constitute Rayleigh scattering. We also minimized Rayleigh scattering cross-section by choosing Hydrogen for our material since it is roughly proportional to $Z^2$, [6]. With these selections, the produced ASADs purely captured the physics of Klein-Nishina cross-section as anticipated without any additional instrumental or physical effects. We generated ASADs for multiple beam parameters to compare with the results of our numerical calculations. The results of these comparisons are presented in § 6.

Realistic Instrument Considerations

Energy Thresholds and Kinematic Limits

The integrals in equation (15) assume a full range of Compton scattering angles can be measured with a detector, but in general this range will be affected by the pre-set energy detection thresholds in the instrument $E_{\mathrm{th}}$.¹ Polarimetry events ideally consist of photons that Compton-scatter off a detector element and are subsequently fully absorbed in a second detector element, typically referred to as two-hit or Compton events. Both detector elements must register energy depositions above the instrument’s trigger threshold for the event to be identified as a Compton coincidence. A photon with a total energy $E_\gamma$ will deposit energy to the recoil electron $E_e=E_\gamma-E'_\gamma$ and deposit its remaining energy $E'_\gamma$ in the second detector element. Thus, the conditions for a detector to record Compton events are $E_e>E_{\mathrm{th}}$ and $E'_\gamma> E_{\mathrm{th}}$. Figure 8 shows the range of valid Compton scattering angles $\theta_{\min} < \theta< \theta_{\max}$ as a function of incident photon energy for three representative threshold energies, illustrating how the window for coincident detection narrows and closes entirely below a threshold-dependent crossover energy $E_{\gamma,\min}$. Using the Compton kinematics equation (1), we can express the scattering angles as:

The condition on the scattered photon energy sets an upper limit on the scattering angle, $\theta_{\max}$, since larger angles correspond to greater energy transfer to the recoil electron, leaving the scattered photon with insufficient energy to trigger the second detector element. Conversely, the condition on the recoil electron energy sets a lower limit $\theta_{\min}$, since smaller angles deposit less energy into the first detector element, i.e., leaving a lower energy recoil electron. The boundary angles are found by setting each condition to equality:

While the angular limits $\theta_{\min}<\theta<\theta_{\max}$ are straightforward to compute in the photon frame, translating them into limits on the transformed angle $\vartheta$ in the detector frame is non-trivial. From Equation (14) where $\cos\theta=-z_{\gamma}(\vartheta,\varphi)$ so that:

This expression cannot be easily inverted to give $\vartheta$ as a function of $\varphi$ alone, since the integral limits would then vary at every point along the $\varphi$ axis, so the adaptive quadrature routine would need to solve a transcendental equation for each $\varphi_i$ to find $\vartheta$ limits, making numerical evaluation impractically slow. A practical solution is to instead introduce a Heaviside gate function:

This function gates the integrand, contributing unity only for physically valid scattering geometries and zero otherwise, without altering the integration limits. The integrals in Equation (15) are then replaced by²:

In the aforementioned equation, it is shown that the solid angle in the photon frame $\mathrm{d}\Omega_{\gamma}=\sin\theta\mathrm{d}\theta\mathrm{d}\varphi$ is equivalent to the transformed solid angle in the instrument frame $\mathrm{d}\Omega=\sin\vartheta\mathrm{d}\vartheta\mathrm{d}\varphi$, since the Jacobian of the rotation matrices is $\lvert\mathcal{J}\rvert=1$.

In a real instrument the detection efficiency does not switch on sharply at a single threshold energy. The trigger efficiency as a function of deposited energy instead follows a gradual turn-on, softened by the combined action of finite energy resolution, channel-to-channel gain and threshold variations, electronic noise, optical crosstalk between detector elements, and the event selection. We describe the turn-on by a Gaussian distribution of effective thresholds with mean $\mu$ and width $\sigma$. Averaging the sharp gate of Equation (19) over that distribution has a closed form (Appendix § 12): because a coincidence survives a threshold exactly when the smaller of the two deposits exceeds it, the Heaviside gate $\mathcal{H}$ is replaced inside the $\vartheta$ integrand by the smooth acceptance weight

where $\Phi_{\mathrm{N}}$ is the standard normal cumulative distribution function, giving an error-function turn-on; $E_e$ and $E'_{\gamma}$ are evaluated at $(\vartheta,\varphi)$ through the scattering kinematics, $\mu$ is the center of the effective threshold distribution, and $\sigma$ is its width. Figure 9 shows the weight for the POLAR values adopted below. This continuum gate is continuous in incident energy by construction and is the threshold model used throughout this work; the sampled mixture of sharp gates from which it is derived remains available in the released code as an option (Appendix § 12). The gate has the form the hardware exhibits: the in-orbit trigger turn-on of each POLAR channel is itself measured to be an error function ((Li et al. 2018)).

For POLAR we adopt $\mu=25$ keV from the measured average detection threshold, about 200 ADC at the typical in-orbit gain of about 8 ADC per keV ((Li et al. 2018)), and $\sigma=75$ keV as an effective width absorbing the optical crosstalk between neighboring channels, the channel-to-channel gain spread, the energy resolution, and the event selection. A response-level scan of $\sigma$ at fixed $\mu$ against the full polarized simulation grids confirms this choice: the optimum lies inside a flat region containing 75 keV (§ 6.3), while the measured per-channel turn-on width is an order of magnitude smaller, so $\sigma$ is an effective parameter rather than a resolution. For other instruments the prescription is to use the instrument’s average detection threshold for $\mu$ and, for $\sigma$, a width representative of the effects that smear the turn-on: for a clean instrument this is simply the energy resolution near the threshold, while for an instrument like POLAR it is an effective parameter to be fitted. We emphasize that this is deliberately an effective single-parameter description, and its limitations, including a mild energy dependence of the fitted width, are discussed in § 6.3.

Geometric and Non-Uniformity Effects

The physical design and imperfections of a real detector further modify the measured ASAD beyond what is predicted by the KN cross-section. We define a detector’s measured ASAD $\mathcal{M}(\varphi)$ with the following product:

Here, $\Phi(\varphi)$ captures the physics-driven distribution and $\mathcal{I}(\varphi)$ encapsulates all instrumental effects. There are two distinct contributions to $\mathcal{I}(\varphi)$: geometric effects and non-uniformity effects.

Typical polarimeters consist of segmented detectors that introduce a bias in the ASAD like POLAR and POLAR-2 ([29]). Photons scattered toward adjacent detectors (i.e., at $\varphi= 0^\circ$, $90^\circ$, $180^\circ$, $270^\circ$) traverse a shorter path length and are therefore detected with artificially higher probability than photons scattered at intermediate angles. This introduces a $90^\circ$ periodicity to the ASAD. Additionally, multiple scattering within a single detector element before absorption in a second element leads to an incorrectly reconstructed azimuthal angle $\varphi$. A further realistic effect is incomplete absorption of the scattered photon in the second detector element, which biases the reconstructed energy and, through the kinematic threshold conditions, the recovered scattering angle. More generally, any departure from azimuthal symmetry of the instrument geometry relative to the incident beam direction introduces additional geometric modulations in the ASAD. On the other hand, the variation between the different detector elements such as gain variations and threshold differences among detector elements would produce complex, non-periodic structures in the ASAD that persist throughout the observation. These non-uniform effects further degrade the measured polarization signal ([20]).

Pointing instruments like PoGOLite [38] and XL-Calibur [1] mitigate geometric effects by rotating during observations, but wide-FOV instruments cannot do so, especially for transient sources such as GRBs. In both cases, Monte Carlo simulations are used to model $I(\varphi)$. These simulations must be calibrated to reproduce the effects of detector geometry, material properties and detection efficiencies.

The framework presented here computes only $\Phi(\varphi)$ and Monte Carlo simulations and calibrations remain essential to extract $I(\varphi)$. Because geometric and efficiency effects are largely independent of the photon’s polarization state, we adopt the approximation:

which implies that unpolarized simulations are sufficient to extract $I(\varphi)$. However, instruments containing passive material surrounding the active detector volume may exhibit polarization-dependent instrumental effects [30]. In that case, this framework would not necessarily hold, depending on the magnitude of the effects, as studied further in § 7.2. Throughout this work we therefore present the response construction in two modes: mode 1, the unpolarized-only mode, which assumes no polarization-dependent instrumental effects, and mode 2, the two-angle mode, which ingests unpolarized simulations like mode 1 and reads the polarization dependence of the instrument from fully polarized simulations at only two polarization angles per energy (§ 5.3.2). The assumption in Equation (23) underlies mode 1 and allows us to generate polarization responses rapidly as presented in § 5.3. In the following section we apply the framework developed here to real polarization analysis.

Application to Polarization Analysis

Classical Analysis

Traditionally, the analysis of polarization data constituted removing geometrical effects and measuring the $\cos^{2}\phi$ modulation in the signal as follows:

• I Histogram the measured azimuthal scattering angle and subtract background to produce the raw uncorrected source ASAD.

• II Using a measured photon energy spectrum and relative source location the ASAD for both 100% polarized and unpolarized emission are created using Monte-Carlo simulations.

• III The instrumental component $I(\varphi)$ is removed from the ASAD to produce the corrected ASAD. This is done by dividing the measured ASAD by its unpolarized (simulated) ASAD counterpart, thereby dividing out $I(\varphi)$ and only leaving $\Phi(\varphi)$.

• IV Equation (3) is then fit to the resulting corrected source ASAD to obtain the fitting parameters $A$, $B$, $\phi_{\mathrm{o}}$. The modulation factor is calculated as $\mu= \frac{A}{B}$.

• V The value of $\mu_{100}$ can be extracted from the simulated $\Phi(\varphi)$ for this specific measurement. To produce this the simulated polarized and unpolarized ASADs are divided and the result is fitted with equation (3) producing the normalizing modulation factor: $\mu_{100}$.

• VI The polarization level can now be calculated using $\Pi= \frac{\mu}{\mu_{100}}$ and the polarization angle is $\phi_{\mathrm{o}}$.

Figure 10 shows a simplified example of the classical analysis approach for a single energy fully polarized ASAD to measure $\mu_{100}$ after correcting with an unpolarized ASAD at the same energy. Here we use the simulated mono-energetic POLAR responses described in § 6.2. The correction removes the geometric effects of the instrument, but the classical analysis remains blind to higher harmonics in the signal due to geometric projection of Klein-Nishina cross-section as discussed in this work. For this particular example these harmonics are relatively small, as discussed later in this work.

Forward-Folding Analysis(threeML)

The above method is relatively straightforward and requires only two simulations for the analysis, one to produce the 100% polarized ASAD and one for the unpolarized one. A prerequisite of this approach is, however, that the spectral shape of the photons is known, as the shape of both the ASADs depends on this. For common sources like the Crab Nebula this is generally not an issue. However, for a source with an unknown spectrum, such as a Gamma-Ray Burst (GRB), one would have to measure the spectrum first before applying this analysis. To ensure the uncertainty on the spectrum carries on into the polarization properly, one can use the forward-folding approach where a full energy dependent polarization response is produced. In this method, developed in Burgess et al. [11], rather than dividing out the $I(\phi)$ component, it is included in the response and the data can be directly fitted against the response³. The analysis, also shown schematically in Figure 11, typically follows the following steps:

• I Generate a full instrument response using Monte-Carlo simulations. This response consists of an array of ASADs simulated for individual energies, polarization angles and polarization levels (0 and 100% suffice, as partial polarization is a linear combination of the two).

• II Similarly, a spectral response is produced using the standard approach used in X-ray and $\gamma$-ray astrophysics.

• III Spectral and polarization parameters are picked from a distribution and forward-folded through both the spectral and polarization response. This produces both a modeled ASAD and a modeled spectrum.

• IV The modeled spectrum and the modeled ASAD (which includes $I(\phi)$ and $\Phi(\phi)$) can be directly fitted against the data and a combined likelihood can be calculated.

• V This method is repeated to produce likelihood distributions for the various spectral and polarization parameters.

This algorithm was originally developed as a plug-in ⁴ for the threeML (Vianello et al. [44]) framework and could be used to analyze the data from POLAR (Kole et al. [32]). Currently, it is being generalized as the polpy package⁵ (Mate et al. [34]) to be usable for other instruments such as AstroSAT CZTI and the upcoming POLAR-2 and Daksha (Bhalerao et al. [7]) missions.

The advantage of this approach is that it allows to analyze sources with unknown spectra such as the GRBs observed only by POLAR (Kole et al. [32]) for which no spectral information was available from other missions. In addition, through its integration in the threeML framework it allows for joint analysis with data from other instruments which observed the GRB, thereby reducing systematic errors on the spectra which further reduce systematic errors on the polarization. A downside of the method is that it is computationally expensive as it requires a large number of simulations to produce the full response. For POLAR, for each GRB (so for a given incoming angle with respect to the instrument), a 3d grid of ASADs was simulated. This response is computationally expensive to produce and each new GRB detected at a particular location requires a new response to be generated.

Rapid Response Generation via Analytical Prediction

While the threeML algorithm presented in § 5.2 is more general than the classical analysis presented in § 5.1 in measuring polarization for off-axis sources simultaneously while also doing spectral fits and joint analysis with other missions, it requires a very large response to be generated. Each response generated to analyze a GRB observed by POLAR took about two weeks on a computer cluster. So for the POLAR catalog of 14 GRBs, it required about 7 months of simulations on a large computer cluster. By utilizing the information presented thus far in this work, a single complete response could be generated much faster as will be shown here.

We can write the ASAD for an unpolarized and a fully polarized beam for any source with parameters $(\varepsilon, \delta, \lambda, \phi_o)$ as:

It was outlined in § 3.2, that $\Phi(\varphi)$ can be calculated for any beam parameters $\varepsilon,\delta,\lambda$ and $\mathcal{I}(\varphi)$ can be taken from an unpolarized ASAD under the assumption that there would be negligible polarization-dependent instrumental effects, as follows. Starting from the unpolarized simulated ASAD $\mathcal{M}^{\mathrm{unp}}(\varphi;\varepsilon,\delta,\lambda)$ of Equation (24), we can peel off $\mathcal{I}(\varphi;\varepsilon,\delta,\lambda)$ by calculating $\Phi^{\mathrm{unp}}(\varphi;\varepsilon,\delta,\lambda)$. Then, we calculate $\Phi^{\mathrm{pol}}(\varphi;\varepsilon,\delta,\lambda)$ which is to be used to calculate a reasonable estimate for the expected simulated fully polarized ASAD $\widetilde{\mathcal{M}}_{\phi_o}^{\mathrm{pol}}(\varphi;\varepsilon,\delta,\lambda)$ at a particular polarization angle. That is:

The last term here is a normalization factor to account for the difference in normalization of $\Phi(\varphi)$ between polarized and unpolarized distributions. We lastly can account for the relative counts in the unpolarized ASAD as follows:

The tilde symbol “$\sim$” is meant to mark the expected-calculated distribution $\widetilde{\mathcal{M}}_{\phi_o}^{\mathrm{pol}}$, such as that calculated through our approach here, not the actual distribution obtained from a GEANT4 simulation $\mathcal{M}_{\phi_o}^{\mathrm{pol}}$. This unpolarized-only construction, requiring no polarized simulations, is the response mode we refer to as mode 1 throughout the remainder of this work; by construction it assumes no polarization-dependent instrumental effects (Equation (23)). As will be shown later, our calculations for the polarized ASAD matched their counterparts generated by GEANT4 simulations for the POLAR instrument response. We checked responses generated for multiple source positions at different energies and polarization angles. Furthermore, the expected-calculated distributions can be used to obtain polarization measurements comparable to those obtained using GEANT4 polarized simulations. This has been verified for real GRB data as well as brighter synthetic GRB data detected by the POLAR instrument, this is presented in § 7.

This approach cuts down the simulation time by a significant factor. A real-life example is the responses typically generated to run polarization measurement analysis using the threeML algorithm outlined in § 5.2. For each GRB (specified by its location), a full POLAR response is generated by simulating 150 mono-energetic beams from 5 to 750 keV in 5 keV steps; at each energy the response stores one unpolarized simulation and 61 fully polarized simulations at 3° steps from $\phi_o=0^\circ$ to $180^\circ$, where the $180^\circ$ entry repeats the configuration of the $0^\circ$ entry to give the fitting software the closed boundary it expects, so that there are 60 independent polarization angles and 62 simulations per energy, 9,300 simulations in all. Mode 1 replaces the 61 polarized simulations per energy by the analytic construction and needs 1 of the 62 simulations per energy, 98 percent fewer, so the simulation time is reduced by a factor of about 60 (Table 1). Mode 2 (§ 5.3.2) needs 3 of the 62, 95 percent fewer, so either mode needs at least 95 percent fewer Monte Carlo simulations than a full response. The polarization-angle step of the generated grid is not tied to the simulated one: both modes refine it in compute only, since the model is analytical. A single full response was typically produced within 2 weeks on a computer cluster, but with this approach, a full response can be produced within a few hours.

Approach	Simulation per energy	Simulation wall time	Construction	$\phi_o$ refinement
Classical full simulation	61 pol + 1 unp	$\sim$2 weeks per GRB	$\cdots$	grows with the number of simulated angles
Mode 1 (unpolarized only)	1 unp	$\sim$5 hr (est.)	13 s	compute only
Mode 2 (two-angle)	1 unp + 2 pol	$\sim$16 hr (est.)	13 s	compute only

Error Propagation on the Analytically Predicted ASAD

The analytical distributions $\Phi^{\mathrm{pol}}(\eta)$ and $\Phi^{\mathrm{unp}}(\eta)$ are obtained by numerical integration of the Klein-Nishina cross section and carry no statistical uncertainty. The only source of error in the predicted ASAD is therefore the unpolarized simulation $\mathcal{M}^{\mathrm{unp}}(\eta_i)$, whose bin counts follow Poisson statistics. The raw bin counts $M_i^{\mathrm{unp}}\cdot f$ are Poisson-distributed with uncertainty $\sqrt{M_i^{\mathrm{unp}}\cdot f}$; rescaling by $1/f$ to return to the normalized histogram gives the per-bin uncertainty:

Combining Equations (25) and (26), the predicted polarized ASAD can be written as:

where $\mathcal{K}_i = \Phi^{\mathrm{pol}}(\eta_i)/\Phi^{\mathrm{unp}}(\eta_i)$ is the analytical polarized-to-unpolarized ratio in the $i^{\mathrm{th}}$ bin, $\mathcal{C}$ is a purely analytical normalization constant, and $\mathcal{N}$ is a data-driven normalization that preserves the total number of counts. To first order, we treat $\mathcal{N}$ as a constant multiplier. This is justified because $\mathcal{N}$ depends on sums over all $N_{\mathrm{bins}}$ bins, so the implicit correlation introduced by varying a single $M_i^{\mathrm{unp}}$ is suppressed by a factor of order $1/N_{\mathrm{bins}}$ relative to the direct dependence. Under this approximation, $\widetilde{M}_i^{\mathrm{pol}}$ is linear in $M_i^{\mathrm{unp}}$, and the uncertainty propagates as:

Equivalently, the fractional uncertainty on the predicted bin equals the fractional uncertainty on the corresponding unpolarized bin:

We note that this treatment neglects the weak bin-to-bin correlation introduced by the global normalization $\mathcal{N}$. A full covariance treatment would introduce off-diagonal terms of order $1/N_{\mathrm{bins}}$ relative to the diagonal, which is negligible for the extents of this work. The uncertainty defined in Equation (29) is the $\sigma_{\mathrm{pred},i}$ used in the goodness-of-fit metric of § 6.

The Two-Angle Response Mode (Mode 2)

The response mode described above requires no polarized simulations: the polarized ASADs are predicted analytically from the unpolarized prior under the assumption of Equation (23). The comparison of that prediction with the full simulated grid of the POLAR responses (§ 6.2) motivates a second mode. When the residual of the mode 1 prediction against the simulated polarized ASADs is mapped over energy and polarization angle (§ 6.2), it is periodic in the polarization angle with period $180^\circ$. Per azimuthal bin, the simulated polarized ASAD as a function of the polarization angle is therefore described by a constant plus one sinusoid in twice the angle: three numbers. The constant is supplied by the unpolarized simulation, and the two sinusoid coefficients by fully polarized simulations at two angles, we chose $\phi_o = 0^\circ$ and $45^\circ$, which isolate the cosine and the sine terms:

That is the whole construction: simulating the unpolarized beam and fully polarized beams at $\phi_o = 0^\circ$ and $45^\circ$ for each energy determines the ASAD at every polarization angle, and at every partial polarization fraction, up to the Poisson noise of the input ASADs. Any two polarized simulations at angles whose difference is not a multiple of $90^\circ$ determine the two coefficients (at a multiple of $90^\circ$ the sinusoid repeats up to sign and the pair carries one number); $0^\circ$ and $45^\circ$ isolate the cosine and the sine terms and condition the solution best, and Equation (31) is the form for that pair, while for another pair the two coefficients are the solution of a two-by-two linear system per bin, and this is followed for the released code. It is validated on every polarized angle not used as input (§ 6.2); as a blind test, every polarized simulation other than the two inputs was removed from the input files and the response rebuilt, and the output was identical. Because the polarized ASADs are simulated rather than predicted, this mode carries the polarization-dependent instrumental effects that mode 1 assumes away, while replacing the 61 polarized simulations per energy of the full response with two. This is the form expected if the instrument response is linear in the polarization state of the beam, as in the on-axis Stokes description of Kislat et al. [26]; whether that linearity holds exactly off axis and through a real instrument is left to future work. The released code accepts any number of input angles from two upward.

**Note—**Abbreviations: pol = fully polarized simulation, unp = unpolarized simulation, per simulated energy. Wall times refer to the POLAR production grid ($3^\circ$ polarization-angle step, 61 simulations per energy): the classical value is the recorded production time for the per-GRB response on a computing cluster, and the mode 1 and mode 2 values are estimates obtained by scaling that time by the required number of simulated ASADs per energy (about 1/60 and 1/20 respectively), assuming simulation time scales linearly with their number. Construction times are measured on a single workstation for the 61-angle by 150-energy grid of the POLAR responses and will vary with hardware; they include reading the input simulation and writing the response files. Refining the polarization-angle step costs the classical approach simulation time in proportion to the number of simulated angles, while both modes refine it in compute only, within seconds, to arbitrarily fine steps; the polarization fraction is likewise a free, compute-only parameter in both modes.

Verification and Systematic Studies

Verification Against Ideal Simulations

We compared the ASADs obtained from the ideal detector simulation, as presented in § 3.3, to the low-energy expansion of Muleri (2014) (Appendix B), and the analytical ASAD derived in this work for a grid of beam energies and zenith locations. We are showing polarized and unpolarized ASADs for each panel. The results are presented in Figure 12. They show that our analytical framework agrees with the simulation over the grid of Figure 12, 50 to 300 keV and zenith angles to $80^\circ$, while the Muleri (2014) approximation diverges at higher energies where the condition $\varepsilon\ll1$ breaks down, most visibly at 300 keV and large zenith angles. The POLAR responses of § 6.2 extend the comparison to 750 keV through a full instrument model, and no approximation enters the integrand at higher energies. The framework therefore applies across the band of § 1, with direct verification to 750 keV.

Verification Against POLAR Instrument Response

For the case study we use POLAR responses generated using GEANT4 that utilize the full mass model of the instrument. Each response was created for a specific GRB observation, defined by its incidence angle with respect to the instrument, to be used in its polarization measurement. The responses consist of polarized and unpolarized ASADs for the grid of $\mathit{E}$ and $\phi_o$ described in § 5.3: 150 energies, one unpolarized simulation and 61 polarized simulations per energy at $3^\circ$ steps from $0^\circ$ to $180^\circ$, of which 60 are independent angles. The POLAR responses and event data used here were obtained from the POLAR collaboration through co-author H.L. For the purpose of this validation, we selected GRB 170101A representing roughly an on-axis case, GRB 170114A representing an intermediate off-axis case, and GRB 170210A representing an extreme off-axis case. The zenith angle of each comparison is taken from the instrument attitude recorded with the data products, and the azimuth from the phase of the simulated polarized distributions, which resolve it to half a degree: GRB 170101A at $\delta= 5.76^\circ$, $\lambda= 72.86^\circ$; GRB 170114A at $\delta= 27.25^\circ$, $\lambda= 6.07^\circ$; and GRB 170210A at $\delta= 82.59^\circ$, $\lambda= 130.9^\circ$. For all comparisons with the POLAR responses, we used the continuum gate of Equation (21) with the adopted POLAR values of § 4.1, $\mu= 25$ keV and $\sigma= 75$ keV, to build each predicted ASAD.

We utilize a modified form of the reduced chi-squared $\chi_\nu^2$ metric as a single summary statistic to assess the agreement between our predicted ASAD and its counterpart from the POLAR response (Gagunashvili [19]). So the per-bin chi-squared $\chi_i^2$ is defined as follows:

Such that: the first term of the denominator is the Poisson variance of the response polarized simulations, evaluated at the model expectation $y_{\mathrm{pred},i}$, with $f$ the normalization factor of Equation (27); and $\sigma_{\mathrm{pred},i}$ represents the propagated errors of predicted counts $y_{\mathrm{pred},i}$ which stem from Poisson error of each bin counts in the response unpolarized simulations. And equivalently, $\Delta_i$ represents the residual per bin and $V_i$ represents the variance per bin. We adopt

the shorthand $\Delta_i$ and $V_i$ for use in the error propagation later-on. We consider the first term from the model rather than from the simulation data, since at finite counts a bin that fluctuates low would otherwise be assigned a smaller variance and therefore a larger weight, which biases the statistic downward.

Two conventions follow from this definition and are applied throughout. First, a bin enters the sum only where the variance is defined and positive, which excludes empty bins. Second, quoted map medians are computed above an energy floor of twice the instrument’s average per-bar threshold, that is 50 keV for POLAR at the adopted gate, and on one set of cells for both response modes: the 58 polarization angles that are not inputs to the mode 2 construction of § 5.3.2. Below that floor the coincidence requirement leaves almost no accepted counts, the model variance follows the predicted counts toward zero, and the ratio in Equation (32) ceases to carry information; the maps are still displayed over the full range, with the sub-floor region marked.

Then, the reduced chi-squared $\chi^2_\nu$ is calculated using:

Where $\nu$ is the number of bins entering the sum. We note here that, this is not a standard reduced chi-squared in the statistical sense because there are no free parameters being fit. We are comparing two fixed distributions, so there is no degrees-of-freedom adjustment for fitted parameters. While it is appropriate for the purpose of gauging agreement, the value $\chi^2_\nu= 1$ does not have the usual interpretation of ”statistically acceptable fit at 1σ”, but it indicates that the residuals between the two distributions are consistent, on average, with the combined uncertainties of both.

The uncertainty on each per-bin $\chi_i^2$ can be obtained by standard error propagation. Since the two distributions (ASAD$_{\mathrm{true}}$, ASAD$_{\mathrm{pred}}$) can be assumed to be statistically independent, the uncertainty on the residual can be derived as $\sigma(\Delta_i) = \sqrt{V_i}$, and the uncertainty on the total variance would be second order so it can be safely dropped. The uncertainty on the $\chi_i^2$ would be dominated by $\sigma(\Delta_i)$ which can be calculated as :

Lastly, for the reduced chi-squared $\chi^2_\nu$, we can similarly propagate the errors through Equation (33) so that:

We show a single ASAD comparison at energy of 150 keV and polarization angle 0 for the mode 1 predicted ASAD versus the GEANT4 POLAR response ASAD of GRB 170114A in Figure 13. The respective errors per each bin are shown in the Figure and the lower panel represents the per-bin chi-squared $\chi_i^2$ with its associated uncertainties. This comparison gives an overall reduced chi-squared $\chi^2_\nu= 1.33 \pm0.12$ over 360 bins.

Now, as Figure 13 shows a decent agreement, we extended the comparisons to the full polarized grid of each response, the 150 energies and 60 independent polarization angles of § 5.3, for the three GRBs, which span the zenith angle range $\delta\in[5^\circ, 83^\circ]$. Before showing the maps of both modes, Figure 14 shows what the mode 1 residual looks like when it is followed across the polarization angle at fixed energy and fixed azimuthal bin for GRB 170114A: the residual is periodic in $\phi_o$ with period 180° and is described by a constant plus one sinusoid in $2\phi_o$. At each energy and azimuthal bin, the residual as a function of $\phi_o$ (not of $\eta$) was expanded as a Fourier series in $\phi_o$ up to the sixth term. In each cell, the squared amplitude of the term with a period of 180° was divided by the summed squared amplitudes of the six sinusoids. The distribution of that ratio over the 50,695 cells above the 50 keV floor has a median of 91.5 percent and quartiles of 80.1 and 96.1 percent, where a residual with no pattern in $\phi_o$ would spread the sum evenly, one sixth per term on average. This is the observation behind the two-angle construction of § 5.3.2.

Figure 15 shows the reduced chi-squared maps of both modes for the three GRBs, on one color scale with white at one. Over the quoted cell set the mode 1 medians are 1.57 for GRB 170101A, 2.55 for GRB 170114A and 1.37 for GRB 170210A (75th percentiles 1.74, 2.86 and 1.64), and the mode 2 medians are 1.01, 1.01 and 1.00 (75th percentiles 1.06, 1.06 and 1.06). The mode 1 maps show the evidence for the presumed polarization dependent instrumental effects. The mode 2 maps carry no such structure; the two input-angle rows are blank because they are ingested rather than predicted. A reduced chi-squared below one means that the residuals between the two distributions are smaller than their combined uncertainty, typically due to low counts.

Systematic Limitations

The mode 1 construction provides a physics-based alternative to dedicated polarized simulations, but it rests on the approximation of Equation (23). We identify several sources of residual discrepancy in the mode 1 construction: (i) the method’s reliance on the unpolarized simulation to capture instrumental effects; (ii) real polarization-dependent instrumental effects, which break the approximation of Equation (23) and are expected when significant passive material surrounds the detector (for POLAR, the entire body of the Tiangong-2 space station); (iii) residual mismatches between the analytical and simulated ASADs; and (iv) the effective threshold distribution is a single fixed approximation, whereas its physically appropriate width is mildly energy-dependent: a response-level scan of $\sigma$ at fixed $\mu$ against the polarized grids of five POLAR responses (the three of this section and those of GRB 170127C at $\delta= 41.82^\circ$ and GRB 170305A at $\delta= 31.54^\circ$) finds the adopted $\sigma= 75\ \mathrm{keV}$ inside the flat optimum of the full energy band. We do not attempt to attribute the observed discrepancies to a single cause. As shown in § 7.2, the net effect of these residuals is negligible for all but the most extreme, exceptionally bright GRBs; the polarization-dependent effects of item (ii), which dominate the mode 1 residual, are modeled directly by the two-angle mode 2, as demonstrated in § 7.3. However, it is noted that this net effect disappears once mode 2 is employed which reproduces the GEANT4 simulation response down to the statistical noise, therefore there is no impact on the measurement even for extremely bright GRBs.

Results and Discussions

In this section we validate the two response modes of this work against the full GEANT4 POLAR response in realistic polarization analyses. We first apply both modes to a real POLAR observation, GRB 170114A, whose incidence angle lies in the region where the mode 1 construction is most strained (§ 6.2), then to a grid of synthetic GRBs with known input polarization parameters that spans the fluence range from a typical bright GRB to the rarest, brightest events, and we close with what the two modes show at the posterior level.

Measuring Polarization for Real GRB Data

To validate the two response modes against the full GEANT4 response in a realistic polarization measurement, we applied all three responses to the POLAR observation of GRB 170114A ($\delta= 27.25^\circ$, $\lambda= 6.07^\circ$) using the threeML forward-folding pipeline as described in § 5.2. This GRB was selected as the baseline test case for three reasons: its incidence angle coincides with the region where higher order harmonics in the ASAD reach maximum amplitude (Figure 4), and the corresponding response exhibits the largest $\chi^2$ residuals in the comparison between GEANT4 and the mode 1 response as shown in Figure 15, and it was bright and long enough while exhibiting a polarization measurement with relatively small credibility intervals as described in Kole et al. [32].

The results are shown in the form of posterior distributions for the polarization parameters $\Pi$ and $\phi_o$ as returned by threeML [44] with the polpy plugin for each response. Although the analysis fits the polarization and spectral parameters simultaneously, for clarity we show only the polarization posteriors and omit the corresponding spectral ones. The three fits share the identical data, spectral response and sampler configuration; only the polarization response differs. Figure 16 shows the result: the GEANT4 response gives $\Pi= 16.6^{+8.2}_{-8.8}$ percent and $\phi_o = 65.4^{+16.2}_{-14.6}$ degrees, the mode 1 response $\Pi= 11.1^{+7.9}_{-6.9}$ percent and $\phi_o = 62.5^{+21.1}_{-18.7}$ degrees, and the mode 2 response $\Pi= 17.1^{+8.3}_{-8.3}$ percent and $\phi_o = 65.7^{+14.9}_{-13.3}$ degrees. At the statistics of this burst the three posteriors overlap, and the mode 2 posterior tracks the GEANT4 one in both parameters.

Effect of Increased Prior Statistics

Because mode 1 relies only on the unpolarized GEANT4 ASADs as a prior, the number of simulations required per response is reduced by a factor of about 60. A natural question is whether redirecting that computational budget into higher-statistics unpolarized simulations would tighten the resulting polarization constraints. To test this, we built the mode 1 response of GRB 170114A a second time from a higher-statistics prior: an independent unpolarized simulation of the same burst geometry with its own spectral response, carrying 30 times the photon count of the nominal unpolarized prior ($1.5 \times10^9$ photons, from the recorded normalization of the simulation). Figure 17 shows the resulting posterior distributions alongside the nominal mode 1 result of Figure 16. Both responses in this comparison are mode 1 constructions, so polarization-dependent instrumental effects are not accounted for in either. The high-statistics prior gives $\Pi= 11.4^{+7.0}_{-6.4}$ percent and $\phi_o = 60.0^{+18.0}_{-16.4}$ degrees against $\Pi= 11.1^{+7.9}_{-6.9}$ percent and $\phi_o = 62.5^{+21.1}_{-18.7}$ degrees for the nominal prior: the credible intervals narrow by 9 percent in degree and 14 percent in angle, and the medians move by +0.3 percentage points and -2.5 degrees, within the intervals. A thirtyfold increase in the prior statistics therefore tightens the measurement only mildly: the width of the posterior is set mainly by the statistics of the data, and the response’s own statistical noise is a minor contribution at this brightness.

Diagnosis of Systematics with Simulated GRBs

The validation in § 7.1 establishes that both response modes recover polarization parameters consistent with the full GEANT4 response for a representative POLAR GRB. However, GRB 170114A is of moderate brightness, and the statistical significance of the measurement is correspondingly limited. To probe the regime in which residual systematics of the two response modes become comparable to the statistical uncertainty, and therefore a dedicated GEANT4 simulation would be required, we constructed a grid of simulated GRBs based on the spectral model and the instrument response for GRB 170114A. Each simulated GRB in the grid was generated with polarization parameters fixed at known input values ($\Pi= 10\%$ or $\Pi= 50\%$, and $\phi_o = 90^\circ$), enabling a direct comparison between measured and true polarization. We adopt the Band-function spectral parameters of GRB 170114A (Kole et al. [32]) as the baseline spectrum, and different fluences. The fluence scan is performed at two input polarization degrees: $\Pi= 10\%$ representative of the value reported in the POLAR GRB catalog Kole et al. [32], and $\Pi= 50\%$ which is near the upper end of values predicted by polarization models of GRB prompt-emission and provides a stress test of the method.

GRBs at the fluence level of GRB 170114A ($1\times$) are detected by Fermi-GBM at a rate of roughly once per week, and GRBs at $5\times$ this fluence are detected at a rate of roughly once every 3 months (Table 2). These fluence levels cover the range populated by the vast majority of GRBs a POLAR-class mission will observe during its operational lifetime. At both $\Pi= 10\%$ and $\Pi= 50\%$, the mode 1 and GEANT4 posteriors at $1\times$ and $5\times$ are statistically consistent, as shown in Figure 18: the median offsets in $\Pi$ stay within the combined $1\sigma$ credible intervals in every cell except the $\Pi= 50\%$, 5× cell, which sits at about $1.0\sigma$. In this fluence range, the two responses are statistically indistinguishable in practice, and the mode 1 method is therefore sufficient for the overwhelming majority of GRBs observed during a mission lifetime.

Fluence multiple	Fluence ($10$–$1000$ keV)	Return time	Comparison
	($\mathrm{erg}\,\mathrm{cm}^{-2}$)	(all-sky)
1×	$2 \times10^{-5}$	$\sim1$ week	GRB 170114A; a typical bright GRB
5×	$1 \times10^{-4}$	$\sim3$ months	brighter than $\sim98\%$ of GBM detected GRBs
50×	$1 \times10^{-3}$	$\sim3$ years	$\sim3$ events per decade of GBM operation
500×	$1 \times10^{-2}$	$\sim100$ years	$\sim20\times$ fainter than BOAT GRB 221009A

For brighter GRBs at 50× and 500× the fluence of GRB 170114A, which are extremely rare, the statistical uncertainty of the recovered $\Pi$, $\phi_o$ collapses significantly, and a small systematic offset between the mode 1 and GEANT4 results becomes resolvable. The mode 1 response underestimates $\Pi$, indicating that the offset reflects a fixed property of the unpolarized-only construction rather than a statistical fluctuation. We attribute it to the residual polarization-dependent instrumental effects that mode 1 assumes away by construction through Equation (23). The offset grows with the degree of linear polarization of the incident beam: at 500× the mode 1 median lies 1.0 percentage points below the GEANT4-response median at $\Pi= 10\%$ and 5.7 percentage points below it at $\Pi= 50\%$ (5.7 at 50×). The recovered polarization angle $\phi_o$ is affected more weakly: an angular offset of $0.3^\circ$ at $\Pi= 50\%$ ($0.2^\circ$ at $\Pi= 10\%$) becomes resolvable only once the statistical precision collapses at 500×, and its absolute scale remains small compared with the credible intervals of GRBs at realistic fluences. This diagnosis defines exactly the regime in which the polarization-dependent effects must be modeled, which is what the mode 2 construction provides.

Mode 2 at the Posterior Level

The mode 2 response of GRB 170114A is built from the unpolarized simulation and the fully polarized simulations at $\phi_o = 0^\circ$ and $45^\circ$ at each of the 150 simulated energies, with every other polarized simulation of the grid held out for validation (§ 6.2); as a blind test, every polarized simulation other than the two inputs was removed from the input files and the response rebuilt, and the output was identical. We repeated every threeML fit of § 7.1 and 7.2 with this response, using the identical data, spectral response and sampler configuration. Figure 16 and Figure 18 show the posteriors and Table 3 collects the recovered medians. Mode 2 reproduces the GEANT4-response fit in every dataset within the credible intervals: the polarization-degree medians agree within 0.06 to 0.5 percentage points and the angle medians within $0.7^\circ$ at the widest and $0.06^\circ$ at the closest. At the brightest cell ($\Pi= 50\%$, 500×), where the posteriors are narrowest, the difference is $0.2^\circ$, at the level set by the statistical noise of the simulated inputs, and the simulated input angle of $90^\circ$ is the angle farthest from both mode 2 inputs. The mode 1 degree offset of § 7.2 is absent in mode 2, which recovers $\Pi= 50.3^{+0.3}_{-0.3}$ percent at $\Pi= 50\%$ and the highest fluence against the input 50 and the GEANT4-response value of 50.5 percent. At $\Pi= 10\%$ the GEANT4-response fit itself departs from the simulated input at the brightest fluences: at 500× its polarization degree sits 1.2 percentage points below the input of 10 percent and its angle $4.1^\circ$ above the input of $90^\circ$, and at 50× the angle sits $7.2^\circ$ above it, three to four posterior widths in each case; at the brightest fluence the counts in the data reach a tenth of the accepted counts of the simulated response at a typical energy and exceed them at and below 100 keV, so the statistical noise of the simulated response itself may no longer be small against that of the data. Mode 2 reproduces that reference fit rather than the input: the appropriate benchmark for a response construction is the fit obtained with the full simulated response, and mode 2 matches it throughout.

Dataset	$\Pi$ GEANT4 [%]	$\Pi$ mode 1 [%]	$\Pi$ mode 2 [%]	$\phi_o$ GEANT4 [$^\circ$]	$\phi_o$ mode 1 [$^\circ$]	$\phi_o$ mode 2 [$^\circ$]
Real 170114A	$16.6^{+8.2}_{-8.8}$	$11.1^{+7.9}_{-6.9}$	$17.1^{+8.3}_{-8.3}$	$65.4^{+16.2}_{-14.6}$	$62.5^{+21.1}_{-18.7}$	$65.7^{+14.9}_{-13.3}$
1x, $\Pi= 10\%$	$22.4^{+13.4}_{-11.9}$	$19.1^{+11.3}_{-11.0}$	$22.2^{+13.5}_{-13.2}$	$127.6^{+16.4}_{-18.8}$	$123.1^{+15.5}_{-18.4}$	$127.4^{+16.3}_{-19.2}$
5x, $\Pi= 10\%$	$8.3^{+3.3}_{-3.6}$	$7.6^{+3.2}_{-3.2}$	$8.5^{+3.7}_{-3.6}$	$69.5^{+12.7}_{-12.0}$	$70.5^{+12.6}_{-11.5}$	$70.2^{+13.2}_{-12.0}$
50x, $\Pi= 10\%$	$10.6^{+0.8}_{-0.8}$	$9.4^{+0.8}_{-0.8}$	$10.5^{+0.8}_{-0.8}$	$97.2^{+1.9}_{-2.0}$	$96.9^{+2.6}_{-2.6}$	$97.3^{+2.4}_{-2.4}$
500x, $\Pi= 10\%$	$8.8^{+0.3}_{-0.3}$	$7.8^{+0.2}_{-0.2}$	$8.7^{+0.3}_{-0.3}$	$94.1^{+0.9}_{-1.4}$	$93.8^{+0.8}_{-0.9}$	$93.5^{+1.0}_{-0.9}$
1x, $\Pi= 50\%$	$48.2^{+11.9}_{-11.8}$	$48.3^{+10.5}_{-10.8}$	$48.6^{+11.8}_{-11.8}$	$95.4^{+7.0}_{-6.9}$	$92.7^{+7.1}_{-6.5}$	$95.3^{+7.0}_{-7.8}$
5x, $\Pi= 50\%$	$48.2^{+3.6}_{-3.4}$	$43.4^{+3.1}_{-3.1}$	$48.0^{+3.4}_{-3.4}$	$86.2^{+2.3}_{-2.0}$	$86.3^{+2.2}_{-2.2}$	$86.7^{+2.0}_{-2.2}$
50x, $\Pi= 50\%$	$51.5^{+0.9}_{-0.9}$	$45.8^{+0.7}_{-0.8}$	$51.5^{+0.8}_{-0.9}$	$91.2^{+0.5}_{-0.5}$	$90.6^{+0.5}_{-0.5}$	$91.1^{+0.5}_{-0.5}$
500x, $\Pi= 50\%$	$50.5^{+0.3}_{-0.3}$	$44.7^{+0.2}_{-0.2}$	$50.3^{+0.3}_{-0.3}$	$90.1^{+0.1}_{-0.1}$	$89.7^{+0.2}_{-0.2}$	$90.2^{+0.2}_{-0.2}$

Future Implementations and Conclusions

The results of this work are the off-axis azimuthal scattering angle distribution described by Muleri (2014) evaluated without approximation across the band of § 1 and extended by the effect of the detection threshold, two response modes built on it, and the software that builds them. The description carries the off-axis projection effects that most GRB polarization analyses to date have neglected. Mode 1 is built from the unpolarized simulations alone and needs about 60 times fewer simulations than the GEANT4-based approach used for the POLAR catalog (Kole et al. [32]), with the small, characterized caveats discussed in § 7.2. Mode 2 adds fully polarized simulations at two polarization angles per energy, models the polarization-dependent instrumental effects from them, and matches the full GEANT4-response measurement from the faintest to the brightest cases tested (§ 7.3). The released package, MENHAG (Modulation ENgine for Harmonic ASAD Generation), builds either response for any planar Compton polarimeter from that instrument’s own simulations (Appendix § 13).

**Note—**All fits share the identical data, spectral response, and sampler configuration; only the polarization response differs. The GEANT4-response fits are the production fits of the responses delivered with the data; the mode 1 and mode 2 fits are seeded, with the seeds recorded in the result files. The responses are built at $\delta= 27.25^\circ$, $\lambda= 6.07^\circ$. Differences quoted in the text are computed from the unrounded medians.

As the polarization response depends on the source incidence angle, the analytical framework makes it practical to generate a dense grid of responses across the full field of view of a detector. For POLAR this was computationally prohibitive as each per-GRB i.e., per incidence location, responses required approximately two weeks on a computing cluster. That is reduced to few hours when employing the analytical method presented here. The method was therefore already leveraged for preliminary analysis of the scientific prospects of the POLAR-2 mission (Kole et al. [30]). This work also opens the door to time-resolved polarization analysis of long GRBs whose incidence angle changes significantly during the prompt emission. Dividing the burst into time bins each with its own response provides a more accurate probe of the polarization angle evolution. This provides more statistically significant constraints on the jet physics and emission mechanism of the prompt emission of GRBs ([11, 32]).

The polarization-dependent systematics identified in § 7.2 are addressed through the two-angle mode which models them from polarized simulations at two angles per energy, removing the polarization-degree bias of the unpolarized-only construction at all fluences probed. For future response-simulation campaigns we therefore recommend including fully polarized simulations at two polarization angles per energy, in addition to the unpolarized simulation, which enables the two-angle construction at any statistics level. Because mode 1 assumes that the instrument modifies polarized and unpolarized ASADs identically (Equation (23)), the difference between the mode 2 and mode 1 responses of a mass model maps the part of the instrumental response that depends on the polarization state, a potential diagnostic for the design and calibration of future instruments.

Beyond the POLAR mission, the framework developed here generalizes naturally to other wide field of view hard $X$-ray and $\gamma$-ray polarimeters, including AstroSAT-CZTI (Chattopadhyay et al. [14]) and the upcoming Daksha (Bhalerao et al. [7]) mission. Compton imagers like COSI (Tomsick et al. [43]) and GRAPE (Oñate Melecio et al. [37]) could also benefit from this approach by considering the polarization model in the frame of the photon without the transformation discussed in § 3.1. Additionally, narrow field-of-view hard $X$-ray polarimeters such as XL-Calibur (Abarr et al. [1]) could benefit from this approach by constructing responses at fine polarization-angle binning at no additional simulation cost, as needed, for polarization measurements.

We have presented a generalized description of the azimuthal Compton scattering angle distribution for off-axis polarized sources, applicable across the hard X-ray and gamma-ray band of Compton polarimeters, from tens of keV to about 1 MeV, together with rapid response-generation methods that enable practical polarization analysis of polarimeter data. The two-angle mode reproduces the results of full Monte-Carlo simulations with two polarized simulations per energy instead of 60. The unpolarized-only mode requires none at all, reproduces them for the overwhelming majority of GRBs observable by Compton polarimeter instruments, and is the path of choice when polarized simulations are not attainable.

The response-generation methods described in this work are implemented in the released software package MENHAG, available at https://github.com/AstrophyJimmy/menhag and described in Appendix § 13.

ACKNOWLEDGMENTS

We thank the organizers and participants of the workshop “Advancing a Universal Analysis Framework for High-Energy Spectro-Polarimetry” (Geneva, February 2025), supported by the Swiss National Science Foundation, for the fruitful discussions that helped refine aspects of this work. M.G.Z. acknowledges the support of the University of New Hampshire Department of Physics & Astronomy through a teaching assistantship during the preparation of this work.

M.G.Z. acknowledges the use of Claude by Anthropic (Opus 5.5 and Fable 5 models) to assist with code packaging, code documentation improvements, and debugging tests for MENHAG. All of which were independently reviewed and scrutinized by the aforementioned author.

Software: MENHAG [52 version 1.0.0], NumPy [24], SciPy [45], Matplotlib [25], Astropy [3, 4, 5], uproot [39], threeML [44], polpy [11], and its repository, [10], MultiNest [17] through PyMultiNest [9].

Angle Definitions

The transformation from the photon frame to the detector frame introduces several angles. Table 4 collects their definitions for reference.

Symbol	Name	Definition
Photon-frame angles ($x_\gamma y_\gamma z_\gamma$)
$\theta$ ($\theta_{\mathrm{C}}$)	Compton scattering angle	Angle between $\vec{k}$ and $\vec{k}'$; polar angle from $-z_\gamma$
$\theta_{\mathrm{s}}$	Spherical polar angle	Standard polar angle from $+z_\gamma$; $\theta= \pi- \theta_{\mathrm{s}}$
$\phi$	Azimuthal scattering angle	Angle between $\vec{\xi}$ and the $xy$-projection of $\vec{k}'$, from $+x_\gamma$
Detector-frame angles ($xyz$)
$\vartheta$	Transformed polar angle	Polar angle of $\vec{k}'$ in the detector frame; integration variable
$\varphi$	Transformed azimuthal angle	Azimuthal angle of $\vec{k}'$ in the detector frame; integration variable
$\eta$	Measured azimuthal angle	Observable ASAD variable; $\eta= \phi+ \phi_o$ for on-axis incidence
Source incidence
$\delta$	Source zenith angle	Incidence angle measured from the detector axis
$\lambda$	Source azimuth angle	Source azimuth in the detector frame
$\phi_o$	Polarization angle	(Local) angle of $\vec{\xi}$ from the detector $+x$-axis

**Note—**A subscript $\gamma$ denotes a quantity in the incident-photon frame. Superscripts “pol” and “unp” denote the polarized and unpolarized beam, respectively.

The Low-Energy Expansion of the Transformed Cross-Section

Following the approach of Muleri (2014), the integrals in equation (15) are too complex to be evaluated analytically, especially in the transformed frame. However, the individual terms in the expanded KN cross-section (equation (7) and 8) show a geometric series with the shape $(1-x)^{-n}$ where $x=-\varepsilon(1-\cos\theta)$. So a first order approximation can be performed where $x \ll1$, which is the condition:

At low energies of a few keV (typical X-ray energies), one can expand each term in equations (7) and (8) geometrically as $1+nx$, i.e.,

Then the transformation in equation (12) is applied with the matrices as defined in equation (13) with $\lambda=\phi_o=0$ so this transformation would be limited to varying zenith angles $\delta$. The effect of varying the zenith angle $\delta$, changes the shape of the modulation curve from a single second harmonic, to multiple harmonics. Now, the transformation of the relevant terms becomes:

Substituting the transformed terms from these equations into the approximated terms in equation (37) we can therefore evaluate the integrals in equation (15) analytically to find:

Which provide expressions for the raw ASAD distribution without normalization for polarized and unpolarized sources as a function of the source energy $\varepsilon$ and zenith angle $\delta$ (the symbol $\mathcal{D}$ is reused here for the unnormalized azimuthal distribution; elsewhere it denotes the cross-section per solid angle). The normalization step in Equation (9) is essential only when constructing an ASAD for a spectrum. Each ASAD energy will correspond to a different normalization factor depending on the value of the integral in the denominator which can also be evaluated analytically under the same approximation. At face value, the distributions of Equation (39) are sufficient for studying mono-energetic sources. These equations clearly show the presence of higher harmonics in both polarized and unpolarized ASADs, since both expressions can be rewritten to resemble the mathematical form of a Fourier series (by expanding all powers of cosines into different harmonics, i.e., $A_n\cos(n\varphi)$). This approach correctly predicted the shape of modulation curves for X-ray energies up to $\sim50$ keV ($\varepsilon\approx0.1$) which was sufficient for X-ray polarimeters Muleri (2014). However, at higher energies, the approximation in Equation (39) breaks down. As a result, energies exceeding $\sim100$ keV require a different approach.

The Threshold Gate as the Limit of Sampled Sharp Gates

A real trigger turns on gradually, so we initially modeled it as a finite set of sharp thresholds. The thresholds are sampled at equal spacing from a Gaussian of mean $\mu$ and width $\sigma$; the gated ASAD of Equation (20) is evaluated at each sampled threshold, each curve is weighted by the Gaussian weight of its sample, and the weighted curves are summed. Samples at negative energies are assigned a sharp threshold of zero and keep their weight. This sampled construction remains available in the released code as a named option.

Letting the number of samples grow turns the weighted sum into an integral of the sharp-gated integrand over the threshold value $t$, with the Gaussian as the weight. A coincidence survives a threshold $t$ exactly when the smaller of the two deposits, $\min(E_e,E'_e)$, exceeds $t$, so the sharp gate is a step function of that smaller deposit, and the integral of a step function against a Gaussian is the Gaussian cumulative distribution evaluated at the smaller deposit. The Heaviside gate inside the $\vartheta$ integrand is therefore replaced by the weight $w(\vartheta,\varphi)$ of Equation (21), which is the error-function gate adopted in § 4.1. The closed form is the limit of the weighted sum of the gated integrands taken before normalization; the normalization of Equation (20) is then applied once to the gated integrals.

The MENHAG Software Package

The framework presented in this work is released as an open-source Python package, MENHAG (Modulation ENgine for Harmonic ASAD Generation), available at https://github.com/AstrophyJimmy/menhag and archived on Zenodo (DOI 10.5281/zenodo.23204599, version 1.0.0). The package computes the fully polarized azimuthal scattering angle distributions of a wide field-of-view planar Compton polarimeter from a small set of simulated inputs, in the two modes of this paper. Mode 1 takes, for each simulated energy, the unpolarized ASAD of the instrument together with the source incidence angles ($\delta$, $\lambda$) and the threshold gate parameters ($\mu$, $\sigma$) of § 4.1, and multiplies the unpolarized ASAD by the analytically computed modulation kernel of Equation (25), renormalized as in Equation (26), for every requested energy, polarization fraction and polarization angle; the gate values are instrument-specific inputs that the user measures or fits, and the sampled mixture of sharp gates of Appendix § 12 remains available as an option. Mode 2 takes the unpolarized ASAD and the fully polarized ASADs at a user-chosen list of polarization angles, two at minimum (default $0^\circ$ and $45^\circ$), and reconstructs the ASAD at any polarization angle and fraction through Equation (31), combining more than two angles by least squares per bin; its output container is identical to that of mode 1, so the two modes are interchangeable in downstream fitting. The statistical uncertainties of the input histograms are propagated to every output bin as described in § 5.3.1.

The package reads and writes the POLAR-style ROOT histogram layout, with FITS and CSV adapters, a documented template for any other instrument format that is adapted by editing two reader and writer functions, and an export to the .prsp polarization-response format read by the polpy plugin for threeML. A configuration file lists the energies, the incidence angles, the gate, the polarization grid, the mode and its input angles, and the file paths; the repository carries the test suite, a minimal configuration, and three executed notebooks that reproduce some of the constructions of this paper. Users of the package should cite this paper as the primary reference for the method and the archived software release through its DOI [52]; the repository carries a citation file that encodes both.