Photospheric Polarization Signatures From Long Gamma Ray Burst Simulations

A comprehensive understanding of Gamma Ray Bursts (GRBs) has been elusive due to the variety of questions surrounding the radiation mechanism at play in these events. Polarization measurements of GRBs can heavily constrain the relevant radiation mechanisms and the structure of the GRB jet; however, there is a limited number of theoretical predictions that observed GRB polarizations can be compared against. Here, we conduct radiative transfer calculations of a set of two dimensional relativistic hydrodynamic long GRB (LGRB) jet simulations, of a constant and a variable jet, using the Monte Carlo Radiation Transport (MCRaT) code. MCRaT has been enhanced by the inclusion of polarization; it has been first verified by reproducing a variety of results in the literature and then used to obtain the time integrated and resolved polarization degrees and angles of the synthetic LGRBs. While the obtained time-integrated polarization degrees {($\lesssim 1$\%)} are consistent with the constraints from the POLAR experiment, they are lower than other theoretical studies due to the lack of strong gradients in the model jet profiles that we use. The time resolved results suggests that GRBs with wide jets observed on axis will have small polarization degrees ($\lesssim 2\%$) and {constant polarization angles}, during the brightest portion of the light curve. GRBs observed off axis will have larger polarization degrees and polarization angles that change with the temporal structure of radiating shells in the outflow. We then place our results in the context of GRB prompt emission models and future LEAP and POLAR-2 GRB polarimetry detections.


INTRODUCTION
Gamma Ray Bursts have been detected since the late 1960's as a relatively short pulse of gamma ray photons (Klebesadel et al. 1973). These transient events have been categorized based on their duration. Events that last 2 seconds are Short GRBs and are associated with the merger of compact objects (Abbott et al. 2017;Goldstein et al. 2017;Lazzati et al. 2018) while events that last for 2 seconds are denoted Long GRBs (LGRBs) and are associated with core collapse supernovae (Bloom et al. 1999;MacFadyen et al. 2001;Hjorth et al. 2003). Regardless of the type of GRB that is observed, the physical mechanism that produces the high energy X-ray and γ-ray photons that are observed during the first few seconds of these events, known as the prompt emission, is still under investigation.
There are currently two major competing theories: the synchrotron model (SM) (Rees & Mészáros 1994;Zhang & Yan 2010) and the photospheric model (Rees & Mészáros 2005;Pe'er et al. 2006;Beloborodov 2010a;Lazzati et al. 2009). The optically thin SM describes shells of material, which have been launched with varying speeds by a central engine, colliding with one an-other far from the central engine. These collisions produce non-thermal radiation which is able to escape the jet if the opacity τ < 1. This model is able to account for general characteristics of GRBs such as variability and non-thermal spectra, but is in tension with observational relationships such as the Amati, Yonetoku and Golenetskii Correlations (Amati, L. et al. 2002;Yonetoku et al. 2004;Golenetskii et al. 1983;Zhang & Yan 2010). Although Mochkovitch & Nava (2015) showed that the internal shock model can satisfy the Amati relation under certain conditions, there have been other subcategories of synchrotron models developed in an attempt to rectify these discrepancies. These models consider the effects of both globally ordered and random magnetic fields (Toma et al. 2009;Zhang & Yan 2010).
On the other hand, the photospheric model follows photons that have been produced deep in the jet. These photons interact with the matter in the jet until the jet becomes transparent to the radiation. Unlike the SM, this model is able to reproduce most of the observational relationships (Lazzati et al. 2013;López-Cámara et al. 2014). Subphotospheric dissipation events (Chhotray & Lazzati 2015) and the idea of the photospheric region, in which the photosphere is a volume of space in which photons can still be upscattered by sparse interactions with matter in the jet (Parsotan & Lazzati 2018;Parsotan et al. 2018;Ito et al. 2015;Pe'er 2008;Beloborodov 2010b), contribute to the non-thermal nature of the spectra in the photospheric model. Although this model is able to reproduce general characteristics of GRBs, it is not able to fully account for the relatively large amount of low energy photons that are observed in GRB spectra.
The shortcomings of each model have placed them on relatively equal footing; however, polarization measurements of the prompt emission can help break this degeneracy. While there have been a number of polarization measurements made by a variety of instruments, the results have been largely inconclusive due to the difficulty of conducting such an observation with high precision. From the past decade, the largest linear polarization measurement reported was 98±33% from GRB 041219A (Kalemci et al. 2007) and the smallest was 27 ± 11% from GRB 100826A (Yonetoku et al. 2011), where the reported errors are 1σ (also see Gill et al. (2019) for a comprehensive list of detected GRB polarizations measurements). These measurements and many others are time integrated in order to get as much signal as is possible, however the uncertainties are still very large. These relatively high polarizations are typically interpreted under the assumption that only synchrotron radiation can produce such high polarizations (Waxman 2003;Lyutikov et al. 2003;Burgess et al. 2019). A number of studies have shown that GRB jets with ordered magnetic fields can produce high polarizations ranging between 20% and 70% (see e.g. Deng et al. (2016), Lan & Dai (2020), Toma et al. (2009) and Gill et al. (2019), while jets with random magnetic fields produce smaller polarizations. Photospheric emission was originally thought to only produce very small polarizations, however, it has been shown that this model can produce polarization up to ∼ 40% if the jet has significant structure within δθ ∼ Γ −1 , where Γ is the bulk Lorentz factor of the jet (Lundman et al. 2014a(Lundman et al. , 2018Ito et al. 2014), although this configuration may not occur in GRBs. Structure in the jet refers to gradients in the jet profile and/or anisotropy in the outflow of photons (which is related to the expanding outflow) (Lundman et al. 2014b). If the source of soft photons in the photospheric model is due to synchrotron emission then the detected polarization can increase up to ∼ 50% (Lundman et al. 2018).
Building on the polarimetry technology used to acquire past measurements, the POLAR experiment (Produit et al. 2018) recently reported time integrated linear polarizations for 5 GRBs and had enough statistics to conduct a time resolved polarization analysis for one of the GRBs, GRB 170114A (Zhang et al. 2019). Their analysis showed that the GRBs had relatively small upper limits for their time integrated linear polarizations, with the largest upper limit being 68% and the smallest upper limit being 28%; typically they find that the linear polarizations are 10%. The authors claim that although GRB 170114A had a small time integrated linear polarization, the time resolved portions of the GRB show relatively high polarization ( 10%) with a continually changing polarization angle. In acquiring these results, Zhang et al. (2019) make a number of assumptions about the physics of the GRB polarization such as assuming that polarization degree is constant throughout a GRB and that polarization angle can change in time.
The authors interpret the changing polarization angle to be in strain with the photospheric model. Following this analysis, Burgess et al. (2019) analyzed GRB 170114A by combining information from Fermi, which also observed the GRB. They find similar results in that the SSM seems to describe the GRB well despite the other weakness associated with the SSM model. There is ambiguity, however, since the low linear polarization degree of GRB 170114A is easily attainable in the photospheric model. As a result, Burgess et al. (2019) call for more theoretical modeling and predictions of time resolved polarization signatures that will inform future analysis.
There have been no prior analyses of time integrated or time resolved polarization angle or polarization degree under the photospheric model using the realistic profile of a GRB jet. In this work we focus on the photospheric model and show that it is able to account for a changing polarization angle and variable polarization degrees. We use the MCRaT (Monte Carlo Radiation Transfer) code 1 to provide time integrated polarization predictions and the first time resolved analysis of LGRB simulations using the photospheric model. In Section 2 we describe the methods we use to conduct mock observations of our synthetic LGRBs and outline how polarization is handled in MCRaT. Finally, in Section 3 and Section 4, we present the results and discuss the implications to GRB polarimetry missions and understanding the radiation mechanism in GRBs.

METHODS
There are numerous works that have explored polarization in the energy regime applicable to GRBs. These works have focused primarily on the Stokes parameters formalism (see eg McMaster (1961);Ito et al. (2014);Lundman et al. (2014b);Depaola (2003); Krawczynski (2011)). Here, we describe how we conduct mock observations of polarization degree and polarization angle, and then we describe the implementation of polarization in MCRaT and show that it is able to reproduce the results of Depaola (2003), Krawczynski (2011), andLundman et al. (2014b). For more in depth discussions of the Stokes parameters we refer the reader to the aforementioned references. Finally, we outline how we determine equal time of arrival surfaces within the hydrodynamic simulations used in this work.

Mock Observations of Polarization
We produce light curves and spectra in the same manner outlined in Parsotan & Lazzati (2018) and Parsotan et al. (2018), by collecting photons within a given viewing angle and binning them in time and energy. Here, we outline how we calculate the detected polarization degree, the polarization angle and their respective errors.
The Stokes parameters are a vector, S = (I, Q, U, V ) that holds information about the polarization of electromagnetic radiation. I is the intensity of the electromagnetic radiation, Q and U describe the orientation of the polarization ellipse, and V describes the ratio of the principal axis of the polarization ellipse (Rybicki & Lightman 1979). We follow the convention set by McMaster (1961) and Lundman et al. (2014b) where Q = +1 is oriented with the y-axis of the Stokes plane and Q = −1 is oriented with the x-axis of the Stokes plane. The +U axis is rotated 45 • clockwise with respect to the +Q axis and the −U axis is rotated 45 • clockwise with respect to the −Q axis. Furthermore, we normalize the Stokes parameters such that I = 1 at all times giving us s = (1, Q/I, U/I, V /I) = (1, q, u, v). In our simulations, we only consider linear polarization which means that we ignore any contribution by v. This is appropriate since we assume that electron spins, which directly affect v, are isotropically distributed.
The polarization degree, Π, represents the average polarization of the detected photons (Rybicki & Lightman 1979). From the stokes parameters Π is calculated as: where the second portion of the equation takes the normalization by I into account and ignores v since we do not consider circular polarization. Since the photons in MCRaT are weighted to increase computational efficiency (one photon packet in MCRaT represents some number of real photons in the relativistic outflow), we cannot simply add each photon's detected Stokes parameters to calculate q and u; instead, we have to take the photons weight, w, into consideration by averaging the detected photons' Stokes parameters (see Parsotan et al. (2018) for a discussion of the weight). Thus, we calculate q and u as: The error in the polarization degree, σ Π , is given by Kislat et al. (2015) as: where N is the number of photons that were detected and µ is the modulation factor. For a perfect detector µ = 1, which is what we assume in this work. The polarization angle, χ, represents the net direction of the electric field vector once all the detected photons have been summed over (Kislat et al. 2015). χ, the angle between the +q axis and the electric field vector, measured clockwise towards the +u axis of the stokes plane, is given by Kislat et al. (2015) as: The error in the polarization angle, as given by Kislat et al. (2015), is: For the simulations analyzed in this work, where we assume an axisymmetric geometry, χ should be aligned with the positive or negative Stokes Q values due to the sum of the U parameters adding to zero (Lundman et al. 2014b). We verified that the number of photons used in this work are high enough that u ≈ 0. Additionally, χ has π symmetry so we plot it between −90 • and +90 • .

Polarization in the MCRaT code
In MCRaT all photons are initialized to have no polarization; thus, we set s = (1, 0, 0, 0) similar to Lundman et al. (2014b). Each photon becomes 100% polarized from the very first scattering that it undergoes, which does not bias our results. As mentioned before, we only consider linear polarization which means that we ignore any contribution by v.
In order to deal with polarization, we follow the prescription given by Lundman et al. (2014b). We lorentz boost the photons' four momenta from the lab frame to the fluid frame, then from the fluid frame to the electron rest frame. In the electron frame we conduct our scattering using the full Klein Nishina (KN) cross section and then scatter the Stokes parameters using Fanno's Matrix. Afterwards, we deboost the photons from the electron rest frame and the fluid rest frame back to the lab frame.
In the lab frame, the Stokes plane is oriented such that the +Q axis is pointing in the −φ direction and the −Q axis is pointing in the −θ direction. The +U axis is rotated 45 • clockwise with respect to the +Q axis and the −U axis is rotated 45 • clockwise with respect to the −Q axis. χ is measures clockwise from the +Q axis towards the +U axis.φ andθ are the orthonormal azimuthal and polar unit vectors in a spherical coordinate system where the radial unit vector is parallel to the photon's momentum vector. This setup is a natural outcome of using Lundman et al.'s (2014b) definitions.
Each boost to another reference frame entails rotating the Stokes plane using the Muller rotation matrix (McMaster 1961) given by where the angle of rotation, φ, corresponds to the rotation that orients the y-axis of the Stokes plane perpendicular to the photon three momentum and the velocity vector of the frame that the photon will be boosted into.
The equation for φ is given in Appendix B of Lundman et al. (2014b). After each boost, we rotate the Stokes plane again to ensure that the y-axis of the Stokes plane is perpendicular to the plane formed by the reference frame's z-axis and the photon's three momentum. We use the KN cross section to determine if the photon scatters. The photons gets scattered if ξ ≤ σ KN /σ T , where ξ is a random uniformly distributed number between 0 and 1, σ KN , is the KN cross section given by Rybicki & Lightman (1979), and σ T is the Thomson cross section. If ξ ≤ σ KN /σ T then we sample the differential KN cross section to determine the angles, θ sc and φ sc , that the photon will be scattered into. The differential cross section given by Lundman et al.
where r 0 is the classical electron radius, 0 is the incoming photon energy scaled by the electron rest mass, and , is the scattered photons energy. The outgoing photon's θ sc is acquired by rejection sampling the differential cross section integrated over φ sc following the method outlined by Mathews (2013) for maximum efficiency. We acquire φ sc by choosing a random uniformly distributed value between [0, 2π] when q = u = 0, otherwise we apply a rejection sampling method to the differential KN cross section. The case of q = u = 0 removes the KN cross section's dependence on φ sc , which allows us to choose a value between [0, 2π].
Depaola (2003) presents a method of sampling the KN azimuthal distribution, however they use a KN cross section that is in a different form than the one that is employed above. Thus, we outline our method of sampling Equation 7 to acquire φ sc when the incoming photon is polarized. The azimuthal angle that gives the maximum of the KN differential cross section, φ m , can be solved as We then plug φ m and the previously acquired θ sc into the KN differential cross section to get the normalization for the rejection sampling method, dσKN dΩ (θ sc , φ m ). We proceed as follows: (1) Draw a random number ξ uniformly distributed between 0 and 1 (2) Draw a random φ sc uniformly distributed between 0 and 2π 5) If the above condition is met, accept the value of φ sc , otherwise repeat the sampling Once the scattering is completed, we conduct the scattering of the Stokes parameters. First, we use the Muller matrix to rotate the Stokes plane such that the Stokes plane's y-axis is perpendicular to the plane formed by the incoming and outgoing photon three momenta (Mc-Master 1961;Lundman et al. 2014b). Then, we use where, similar to Krawczynski (2011), we have excluded the factors in front since they cancel out when we normalize the scattered Stokes parameters by I as mentioned above. We have also excluded the fourth row and column of the matrix since we only consider linear polarization.
We reproduce Depaola's (2003) result, as shown in Figure A1, which verifies the sampling algorithm of the differential KN cross section. Additionally, reproducing Krawczynski's (2011) results, in Figure A2, tests the Lorentz transform portion of MCRaT.
In order to test the code globally, we reproduce the results of Lundman et al. (2014b). To do so, we imposed the analytic jet structure provided by Lundman et al. (2014b) on the same simulation frames that we will use in Section 3. This is similar to Lazzati (2016) simulating a variety of outflows by imposing an analytic solution onto hydrodynamic simulation files. The domain of this simulation is 2.5 × 10 13 cm along the direction of the jet and 5 × 10 12 cm along the x axis. We used ∼ 6 × 10 5 photons to conduct our code validation for a wide structured jet with θ j = 0.1 radians (∼ 5.7 • ), Γ 0 = 100 and L = 3 × 10 50 erg/s. This is the same case exhibited in Lundman et al.'s (2014b) wide jet with the exception of the value of L that we use, which was chosen to maximize the number of photons that reached the photosphere before they approached the edge of the domain of the hydrodynamic simulation. There are two major differences between the simulation conducted by Lundman et al. For θ v /θ j 1.6 (θ v ≈ 9 • ) the MCRaT polarization begins to decrease again, coming into strain with what is expected. This is due to the fact that the simulation files that we impose the analytic jet equations onto have a finite domain. Even when the photons reach the edge of the domain (2.5 × 10 13 cm) at θ v 9 • , they haven't fully decoupled from the photosphere (which would be located at r ∼ 1×10 14 cm (Lundman et al. 2014b)) thus decreasing the detected polarization.

Equal Time of Arrival Surfaces
In order to relate the structure of the jet in the GRB simulations to the time dependent observables produced by MCRaT (such as the time resolved polarization), we need to calculate the location that photons would be emitting along the observer's line of sight for a given time in the light curve, t detect . We follow the derivation for Parsotan & Lazzati's (2018) Equation 1. The time in which a photon would be emitted from the jet, t j , can be calculated as where t real is the real time acquired from the hydrodynamic simulation. The radius, r j , at which a photon is emitted by the jet is

RESULTS
In this work, we ran MCRaT on two different FLASH 2D special relativistic hydrodynamic (RHD) simulations which both launched a jet into a 16TI progenitor star with a density profile provided by Woosley & Heger (2006). The first simulation which we call the 16TI simulation, has a jet injection radius of 1×10 9 cm, an initial lorentz factor of 5, an opening angle of 10 • , an internal over rest-mass energy ratio, η = 80, and the engine was active for 100 s (Lazzati et al. 2013). The domain of this simulation is 2.5 × 10 13 cm along the direction of the jet. The second simulation, which we denote the 40sp down simulation, has a jet that was injected with similar initial conditions as the 16TI jet with the exception of the jet being pulsed. The 40sp down simulation jet was on for 40 half second pulses, each followed by another half second of quiescence. The luminosity of each pulse was decreased by 5% with respect to the first pulse (López-Cámara et al. 2014). The domain of the 40sp down simulation is 2.56×10 12 cm along the jet axis. These simulations were previously analyzed by Parsotan & Lazzati (2018) and Parsotan et al. (2018) where they focus on various spectral properties and the synthetic light curves. In this paper we focus on the time integrated polarization and the time resolved polarization properties of each simulated GRB.
We ran the MCRaT code during the time in which the central engine of each model was active in order to investigate the effects of the varying central engine on the radiation. The number of photons injected into each simulation was ∼ 1.9 × 10 7 for the 16TI simulation and ∼ 1.4 × 10 6 for the 40sp down simulation. The order of magnitude difference was necessary to maintain a reasonable computation time for the 40sp down simulation. This is a direct result of the jet in the 40sp down simulation moving at Γ ∼ 10 which is an order of magnitude smaller than the 16TI simulation. Figure 2 shows the time integrated polarization degrees, Π, angles, χ, and peak luminosity of the lightcurve, L pk , as a function of observer viewing angles of the synthetic GRBs. These results are acquired by integrating over the time that the jet is active. Here, we plot Π in blue, the light curve peak luminosity, L pk , in black, and χ in purple for the 16TI simulation, in Figure  2(a), and the 40sp down simulation, in Figure 2(b). L pk is the peak of the synthetic light curve when the light curve is binned into 1 second time bins, analogous to the manner in which it is used in the Yonetoku relationship. In general, we find that Π is negatively correlated to L pk , with spearman's rank coefficients r s = −0.6 2 and r s = −0.27 for the 16TI and 40sp down simulations, respectively. Additionally, Π is positively correlated to θ v , where r s = 0.65 and r s = 0.19 for the 40sp down and 16TI simulation respectively. The direct relationship between Π and θ v is easy to see in Figure  2(b) for the 40sp down simulation, however, it becomes more complicated in the 16TI simulation. This is due to the fact that we see a turnover in the polarization in the 16TI simulation at θ v = 8 • , which is consistent with our verification results in Figure 1. The photons at these larger angles are not fully decoupled from the flow, which means that their expected polarization degree is suppressed due to the ongoing scattering that is changing their Stokes parameters. Excluding θ v > 8 • in the analysis of r s in the 16TI simulation changes it to be r s = 0.64, which is consistent with the value acquired from the 40sp down simulation. Accompanying this feature is a switch in the polarization angle being consistent with 0 • to then being consistent with 90 • .

Time Integrated Polarization
The time integrated polarization is much smaller than what is found in other works. Ito et al. (2014) and Lundman et al. (2014b) find Π as low as a few percent and as high as ∼ 40 %, for off axis observers. The difference between our study and theirs can be attributed to the structure of the jet in the RHD simulation and the fact that not all photons are reaching the photosphere (typically located at 10 13 cm); this is shown in Figures 1 and 2. In Figure 3 we show the structure of the 16TI and 40sp down jet as a function of angle at three different times in the jet's evolution, which correspond to various times in the synthetic light curves shown in the next section. The 16TI simulation has a very fast core at all times. Initially, the jet is smaller than the  relatively uniform lorentz factor profile as a function of angle, although it does vary in time, which is expected from a variable jet. Both of these synthetic GRBs have very wide jets that contribute to the extremely low time integrated polarization in Figure 2; these wide jets are relatively uniform and lacking steep gradients, a source of anisotropy in the flow, which means that there is little structure in the flow to produce very high time integrated polarization (Lundman et al. 2014b). Furthermore, the variability in the 40sp down simulation that should produce relatively large polarization (Π 1%) gets washed out as we integrate over the bright and dim portions of the light curve and what is left is the effect of the very wide jet profile that we observe in Figure 3.
Although not all photons decouple from the jet, particularly at large θ v ( 9 • ), there is still enough of an anisotropy in the outflow to produce polarization that is significantly different from Π ≈ 0%. Since the anisotropies in the outflow will only increase as the jet becomes increasingly transparent to the radiation, the polarization is expected to increase. Thus, we consider all of the polarization measurements at large θ v in our results to be lower limits.

Time Resolved Polarization
An advantage of using MCRaT on a time dependent synthetic 2D RHD GRB jet is the ability to produce time resolved polarization predictions. Figure 4 Figure 5 with similar lines. The 16TI simulation does not exhibit much polarization while the polarization angle stays around 0 • during the brightest portion of the light curve. The 40sp down simulation show relatively high polarization degree at larger viewing angles in addition to evolution of χ within the synthetic GRB (Figure 4(d) at t ∼ 11.5 s). We find that a changing χ is indicative of various shells of material coming into the line of sight of the observer due to structure within the jet, while a constant χ indicates that the emission region is staying relatively constant due to lack of structures within the jet.
angles. In the top panel we plot the light curves in black. In the bottom panel we show the time resolved polarization degree in black and the time resolved polarization angle in purple. The dotted purple line is the χ = 0 • reference line. The 16TI simulation is binned uniformly in time with dt = 0.5 s and the 40sp down results are binned non-uniformly with a binning criteria based on the total energy received in each time bin. This energy must be larger than some critical luminosity before the end of the time bin is determined. Figure   4(b) and 4(c) have luminosity cutoffs of 2 × 10 52 erg/s and 2 × 10 51 ergs/s, respectively.
For the case of the 16TI GRB simulation, shown in Figure 4(a), the polarization is very small at all times (Π 1.5 %). This can be attributed to the lack of structure in the jet, as is shown in Figure 3. The brighter portion of the light curve is indicative of low optical depth regions of the outflow (Parsotan et al. 2018); it is during this period that we get the most well constrained polarization measurements (see Section 2.1 where the errors Figure 5. A pseudo-color density plot of a region of the 40sp down GRB simulation. The red dotted line corresponds to the line of sight of an observer located at θv = 7 • . The various lines are surfaces of equal arrival times that are detected by an observer at θv = 7 • at various times. The highlighted regions, colored green, red, and blue, correspond to the start and end of the same colored regions shown in Figure 4(d). We find that the change in χ seen in Figure 4(d) is due to seeing different portions of the jet at various times. Initially, at t = 11 s, the observer sees the core of the jet. Then they observe photons originating from outer regions of the jet at t = 11.5 − 12 s. Finally, the observer sees more of the inner region of the jet again by t = 12.5 s.
are ∝ N − 1 2 ). We see that Π is also at its maximum, at ∼ 2%, and χ ≈ 0 • at all times during the maximum of the light curve. The 16TI simulation does not have multiple shells of material coming into view of the observers line of sight. As a result, we can use it as a control for analyzing the 40sp down simulation.
Things are much different in the case of 40sp down , where the jet is variable and there is lots of time dependent spatial structure on the scale of δθ ∼ Γ −1 , where Γ is ∼ 10 (see Figure 3). At θ v = 2 • , in Figure 4(b), we see that Π is relatively small with a maximum of ∼2%. At this viewing angle we are in the core of the jet where the profile is relatively symmetric and, as a result, do not have much interference with off axis shells of material coming into the observers line of sight. This case is very similar to the 16TI simulation analyzed in Figure  4(a).
Figure 4(c) shows the 40sp down light curve, E pk , Π, and χ for an observer at θ v = 7 • . In this case, we observe larger Π, with the maximum being ∼ 5%. Additionally, the polarization angle changes multiple times during this synthetic GRB observation. On order to emphasize this oscillation in χ, we replot the largest pulse in Figure  4(c) using larger time bins of dt = 0.5 s, in Figure 4(d). We see that within the first few seconds of the GRB, χ changes from 90 • to 0 • and back to 90 • within a δt ∼ 1 s. Then, χ ≈ 0 • during the brightest portion of the light curve, at t ≈ 13 s. The constant value of χ is similar to what we see for the 16TI simulation where there is little structure in the jet and we are seeing the main emission region along the observer's line of sight (LOS). On the other hand, the changing χ suggests that other shells of material within the jet are coming into the observer's line of sight due to variability in the jet.
In order to confirm that a changing χ is representative of various regions of the jet coming to the observers LOS, we relate the times in the light curve shown in Figure 4(d) to equal arrival time surfaces within the jet. In Figure 5 we plot the density in a region of the jet that corresponds to the regions that are emitting during the times of interest in the light curve. These equal time of arrival surfaces are shown as a variety of different colored highlighted regions that correspond to the vertical highlighted regions plotted in Figure 4(d). We find that at t = 11.0 − 11.5 s the dense jet core (ρ ∼ 10 −6 g/cm 3 at y ∼ 2.54 × 10 12 cm and x ∼ 0 cm) is the dominating emission material, scattering photons into the observers LOS (Parsotan et al. 2018) with χ ≈ 90 • . At t = 11.5 − 12 s, the emitting region is between two dense shells in the jet and the observed photons are primarily originating from the denser outer region of the jet, where ρ ∼ 10 −7.5 g/cm 3 located at y ∼ 2.47 × 10 12 cm and x ∼ 5 × 10 11 cm, changing χ to be ≈ 0 • . Finally, by t = 12.5 s, the observer sees emission from the inner region of the jet again, near x ∼ 2 × 10 11 cm and y ∼ 2.48 × 10 12 cm where ρ ∼ 10 −6.5 g/cm 3 , thus changing the polarization angle once more. The width between the two dense shells of material corresponds to the cδt of the changing χ. The outlined effect of a changing χ as related to the lateral and temporal structure of the jet is seen by an observer located at θ v = 3 • − 8 • . The upper limit is related to the fact that photons aren't fully decoupled from the flow at θ v 9 • and the lower limit is due to the core of the jet being the dominating emission region at all times for θ v 3 • .
The changing χ that we observe are consistent with the observed changing polarization angle of GRB 170114A (Zhang et al. 2019). These changing χ from the 40sp down simulation also suggests that the jet that produced GRB 170114A was variable.

SUMMARY AND DISCUSSION
We have used the MCRaT code to conduct Monte Carlo radiation transfer simulations of time resolved and time integrated polarization in LGRBs. MCRaT injects, propagates, and compton scatters photons using the full Klein Nishina cross section. These photons are injected into an outflow described by a FLASH 2D RHD simulation and are subsequently scattered until the end of the simulation. This process of injecting, propagating and scattering photons is repeated until there are no more photons to be injected in the simulation. The implementation of polarization in MCRaT was verified in multiple ways and we were able to recover the time integrated polarization profile of a wide jet as presented by Lundman et al. (2014b). We ran the MCRaT simulations using the RHD simulations of a steady GRB jet, the 16TI simulation, and a variable engine jet, which we denote as the 40sp down simulation.
Primarily, we have found that: • Not all photons at large θ v ( 9 • ) in the MCRaT simulation are decoupled from the synthetic GRB jet which decreases the mock observed Π; as a result, the Π presented in this paper at large θ v are lower limits.
• The time integrated polarizations are generally positively correlated with θ v , with r s ≈ 0.65 for both simulations, and negatively correlated with L pk , with the correlation between Π and L pk being r s = −0.6 and r s = −0.27 for the 16TI and 40sp down simulations respectively.
• The 16TI simulation has very little structure on the scale of Γ −1 which decreases Π and contributes to χ being approximately constant in time.
• The 40sp down simulation shows more structure temporally and spatially which contributes to various shell of materials coming into the line of sight of the observer, as a result, this simulation has larger time resolved Π ∼ 5% and a changing χ.
• The δt for χ to change in a variable jet is indicative of the width of shells of materials within the GRB jet.
• A changing χ also indicates that the observer is seeing various regions of the GRB jet, which may help constrain temporal and lateral variability in the jet structure.
The MCRaT results of the time integrated quantities (Π 1% and χ ≈ 0 • ) are consistent with the results presented by Lundman et al. (2014b), where wider jets produce lower polarization degrees, however our results show the importance of ensuring that all photons are decoupled from the flow. Our simulations have a finite domain which prevents some photons at θ v 9 • from being decoupled from the flow by the time they propagate to the edge of the simulation domain. As Parsotan et al. (2018) mentioned, the peaks in the light curves are composed of photons that are mostly decoupled from the outflow while the quiescent portions are still coupled to the flow. The dimmer regions contribute to lowering the overall time integrated Π, although we do detect some polarization from these dimmer times in the light curves due to structure in the GRB outflow (Lundman et al. 2014b). As the photons become less coupled to the outflow, the asymmetry in the jet will become more pronounced with respect to the radiation and we will expect to detect larger Π. This drives the need to conduct larger domain RHD simulations, to ensure that the photons at larger θ v ( 9 • ) are decoupled from the flow and we have an accurate mock observation of Π at these large angles.
The results presented here are consistent with the results found by Zhang et al. (2019). Our mock detected Π are within the limits that they find for their GRBs (Π 10 %). In particular, Zhang et al. (2019) provide a 99 th percentile upper limit of Π < 28% for GRB 170114A, but quote Π ≈ 4% which is consistent with the Π that we acquire in our study.
Additionally, we are able to show that variable GRB jets can produce changing χ such as is observed in GRB 170114A. Zhang et al. (2019) report the χ of the aforementioned GRB to be 122 • for the first portion of the burst and 17 • for the second portion. While our polarization angles can only change between 0 • and 90 • due to the symmetry of the simulations, we may see a continuous change in χ from 3D MCRaT simulations. Additionally, Zhang et al. (2019) do not specify any errors on their values of χ in GRB 170114A, however, their plots of the fitted Π and χ parameters' confidence contours show that χ is not well constrained. As a result, it is still possible that χ is ∼ 0 • or ∼ ±90 • which is expected for an axis-symmetric jet, such as the ones studies in this work. The fact that the 40sp down simulation produces a changing χ suggests that GRB 170114A had a variable jet, with various parts of the jet coming into the observer's line of sight.
One important assumption that Zhang et al. (2019) make is that Π is constant throughout a GRB and χ can change in time. In our analysis we find that this assumption may not be valid (see for example Figure  4(d) where both Π and χ change). Even if the GRB has a steady jet being injected (which is not known a priori), χ will be approximately constant and Π may change (see Figure 4(a)).
With future GRB polarimetry missions, such as POLAR-2 (Sun 2019) and LEAP (McConnell et al. 2017), expected to observe GRB polarizations with higher precision than that of the POLAR detector, we can use our findings to place constraints on the expected Π and χ of the photospheric model compared to other GRB emission models. Globally, we expect the polarization angle, χ, to be ∼ 0 • or ∼ 90 • . While it is possible that the photospheric model may be able to account for other values of χ in 3D, it is unlikely that the polarization angle will deviate significantly from being aligned perpendicular or parallel to the plane defined by the jet axis and the observer's line of sight. Thus, a measurement of χ that is significantly different from these values combined with a large value of Π may indicate that there is a global magnetic field in the jet (Toma et al. 2009). In the case of a random magnetic field, the observed Π will be low and the observed value of χ would be ∼ 0 • (Toma et al. 2009;Gill et al. 2019), similar to the photospheric case. The ability of the photospheric model to account for very low and high values of time integrated polarization degrees based on the geometry of the jet (Lundman et al. 2014b), makes it very difficult to distinguish from the synchrotron model, with a global or random magnetic field, based solely on that observed parameter. As a result, time resolved Π and χ observations becomes much more important. The results acquired in this work suggests that GRBs observed on axis will have small Π 2% and constant χ, during the brightest portion of the light curve, while GRBs observed off axis will have larger Π and χ that change by ∼ 90 • with the temporal structure of the GRB. If the emission mechanism is synchrotron with a global magnetic field, the Π and χ should not change in time; however with a random magnetic field, Π and χ may vary randomly based on the magnetic field configuration of the GRB at a given time at a given observer line of sight (Gill et al. 2019). There still needs to be more research conducted to pro-duce robust time resolved polarization degree and angle predictions for a variety of jet structures and magnetic field configurations, for each emission model. Unlike the analysis conducted by Lundman et al. (2018), we do not consider the polarization at various energy bands. Lundman et al. (2018) has shown that this can be an interesting area of testing the photospheric model and the capability of testing time resolved polarization of different energy bands can also prove to be fruitful. With this goal in mind, and to ensure that there are enough low energy photons in the outflow (Parsotan & Lazzati 2018;Parsotan et al. 2018), MCRaT will be modified to consider the effects of synchrotron radiation and absorption. This improvement, combined with larger domain RHD simulations will allow us to make stringent time integrated and time resolved predictions of GRB spectra and polarizations.

A. COMPARING MCRAT TO OTHER ANALYSIS
To show that the cross section sampling algorithm and Lorentz boosting of the Stokes parameters are correct, we reproduce some results acquired by Depaola (2003) and Krawczynski (2011). Figure A1 shows the resulting modulation curve when our Klein Nishina (KN) cross section is monte carlo sampled in black. The analytic profile of the cross section as a function of φ sc calculated by Depaola (2003) is shown as the solid blue line. They sample the KN cross section for a photon beam with 100% polarization along the +Q direction, with an energy of 100 keV and 85 • < θ sc < 90 • . Accounting for the differences in the positive stokes parameters in the convention used, we are able to use our method of sampling the stokes parameter to acquire the proper distribution.
In order to verify the algorithm of switching frames of reference, we reproduce the results acquired by Krawczynski (2011). Krawczynski (2011) scattered a beam of photons, 100% polarized along the +Q axis with frequency ω = 10 12 Hz, with a beam of electrons moving at a Lorentz factor of γ = 100. Following their setup, we produce distributions of the resulting lab frame which we show in Figure A2. These distributions are identical, with the exception of normalization, to the distribution acquired by Krawczynski (2011) in their Figure 6. The difference in the signs of Q and U in our distribution with respect to Krawczynski's (2011) Figure 6 is due to the differences in the orientation of the +Q and +U axis of the Stokes plane.