Numerical Simulation of Photospheric Emission in Long Gamma-Ray Bursts: Prompt Correlations, Spectral Shapes, and Polarizations

We explore the properties of photospheric emission in the context of long gamma-ray bursts (GRBs) using three numerical models that combine relativistic hydrodynamical simulations and Monte Carlo radiation transfer calculations in three dimensions. Our simulations confirm that photospheric emission gives rise to correlations between the spectral peak energy and luminosity that agree with the observed Yonetoku, Amati, and Golenetskii correlations. It is also shown that the spectral peak energy and luminosity correlate with the bulk Lorentz factor, as indicated in the literature. On the other hand, synthetic spectral shapes tend to be narrower than those of the observations. This result indicates that an additional physical process that can provide nonthermal broadening is needed to reproduce the spectral features. Furthermore, the polarization analysis finds that, while the degree of polarization is low for the emission from the jet core (Π < 4%), it tends to increase with viewing angle outside of the core and can be as high as Π ∼ 20%–40% in an extreme case. This suggests that the typical GRBs show systematically low polarization compared to softer, dimmer counterparts (X-ray-rich GRBs and X-ray flashes). Interestingly, our simulations indicate that photospheric emission exhibits large temporal variation in the polarization position angle (Δψ ∼ 90°), which may be compatible with those inferred in observations. A notable energy dependence of the polarization property is another characteristic feature found in the current study. Particularly, the difference in the position angle among different energy bands can be as large as ∼90°.


INTRODUCTION
Despite the decades of research since their discovery, the origin of the prompt emission of gamma-ray bursts (GRBs) has been a puzzle.One of the main questions is whether the radiation is optically thin synchrotron or photospheric.Each scenario has advantages and disadvantages; therefore, the issue remains controversial.It may be possible that both emissions play a key role in the prompt phase.
Spectral features have often been invoked as important diagnostics for this issue.While synchrotron emission is favored due to the non-thermal nature of prompt emission, it has been suggested in the literature that it fails to reproduce the low energy part and the width of spectra in many GRBs (e.g., Preece et al. 1998;Axelsson & Borgonovo 2015;Yu et al. 2019;Acuner et al. 2020;Dereli-Bégué et al. 2020;Li et al. 2021).Combined with the fact that it is a natural outcome of the original GRB fireball model (Goodman 1986;Paczynski 1986), these spectral analyses reactivated the study of photospheric emission as a plausible alternative.Various authors have demonstrated that photospheric emission is capable of producing observed non-thermal features by considering dissipation (Pe'er et al. 2005;Ioka et al. 2007;Giannios 2008;Lazzati & Begelman 2010;Beloborodov 2010;Vurm et al. 2011;Chhotray & Lazzati 2015;Ahlgren et al. 2015;Vurm & Beloborodov 2016;Bhattacharya & Kumar 2020) and/or geometrical effects (Ito et al. 2013(Ito et al. , 2014;;Lundman et al. 2013Lundman et al. , 2014) ) in the sub-photospheric region.
However, recently it has been claimed that earlier critiques of synchrotron emission may be unjustified since the analysis is biased by the phenomenological fitting formula employed in the studies (Oganesyan et al. 2017(Oganesyan et al. , 2019;;Burgess 2019;Burgess et al. 2020).Although there are a few exceptions, such as GRB090902B (Ryde et al. 2010) and 220426A (Song et al. 2022), these studies argue that spectra produced by the physics-based calculation of synchrotron emission can fit the great majority of the observed GRB.It should be noted, however, that fitting requires fine-tuning of physical conditions (e.g., magnetic fields, minimum energy of electrons). 10It is Ito et al. unclear whether such conditions can be achieved in reality.
Apart from the issue of the spectral features, polarization properties can provide an independent clue to the prompt emission mechanism (see, e.g., Gill et al. 2021, for a comprehensive overview).Several observational efforts have been devoted to measuring the polarization signature.Early attempts by RHESSI and INTEGRAL have inferred highly polarized emission in few GRBs, but these observations seem to suffer from instrumental systematics (Coburn & Boggs 2003;Rutledge & Fox 2004;Wigger et al. 2004;Kalemci et al. 2007;McGlynn et al. 2007).The first dedicated GRB polarimeter, the GAP instrument on board IKAROS, achieved the first polarization measurement with low systematic uncertainty (Yonetoku et al. 2011a).A high degree of linear polarization (∼ 30−80%) in 3 GRBs was reported in their observations (Yonetoku et al. 2011b(Yonetoku et al. , 2012)).Recent updates were provided by the second dedicated GRB polarimeter POLAR on board Tiangong-2 (Zhang et al. 2019;Kole et al. 2020) and also by CZTI on board AstroSAT (Chattopadhyay et al. 2022), although not dedicated to GRB detection.Contrary to the measurements of GAP, 14 GRBs observed by POLAR are found to be mostly compatible with a low to null polarization.AstroSAT also finds that a majority of their sample (20 GRBs) is consistent with being low to null polarization, whereas a fraction of 5 GRBs is inferred to exhibit high polarization ( 40 %).Time-resolved polarization analysis performed in these observational studies has also suggested that some GRBs (GRB 110826A, 170114A, 160821A: Yonetoku et al. 2011b;Burgess et al. 2019;Sharma et al. 2019) exhibit a significant temporal variation in the polarization position angle.It is discussed in the literature (Burgess et al. 2019) that the low polarization inferred in the majority of GRBs detected by POLAR and As-troSAT may result from the variation of the position angle: the emission may possess intrinsically high polarization, but it is canceled out in the time-integrated measurements as a result of the summation of polarized emissions with various position angles.All polarization measurements that are currently available, however, suffer from low photon statistics and accompanying large error bars.Thus, a clear picture of the polarization properties is yet to be established.Future polarization missions such as POLAR-2 (Kole 2019), LEAP (McConnell et al. 2021), and COSI (Tomsick et al. 2019) are awaited for further constraints.
From a theoretical perspective, many studies have explored polarization properties in the framework of optically thin synchrotron models (e.g., Lyutikov et al. 2003;Granot 2003;Nakar et al. 2003;Waxman 2003;Toma et al. 2009;Lazzati 2006;Zhang & Yan 2011;Gill et al. 2020;Gill & Granot 2021).These studies have shown that, depending on the configuration of magnetic fields, the geometry of outflows, and viewing angles, various degree of polarization ranging from ∼ 0% up to ∼ 70% can be produced.Several studies also explored the polarization in photospheric emissions.The initial study by Beloborodov (2011) finds that, while the photons released at a local emitting region are strongly polarized, the superposition of each emission component suppresses al. (2020).
the net polarization if the outflow is uniform.Following studies have shown up to ∼ 30 − 40% polarization degree can be achieved for an outflow with a sharp velocity gradient in a lateral direction (Ito et al. 2014;Lundman et al. 2014).The study by Lundman et al. (2018) has also shown that if synchrotron emission is induced by sub-photospheric dissipation, the resulting polarization of the photospheric emission can be significant at energies below the spectral peak even when the outflow is uniform.
Combined with spectral features, these predictions of the polarization properties can be used to discriminate between the emission models.It is noted, however, that most studies invoke simplified outflow structures for the theoretical modeling of the GRB emissions.To further test the emission models, it is highly desirable to conduct a study based on realistic jet dynamics.In this respect, an approach to evaluate the photospheric emission based on relativistic hydrodynamical simulation has developed in the last decade.In the early years, a series of studies have explored the properties of long GRBs (LGRBs) by applying a crudely approximated photosphere model to hydrodynamical simulations of a relativistic jet breaking out of massive stellar envelope (Lazzati et al. 2009(Lazzati et al. , 2013;;Mizuta et al. 2011;Nagakura et al. 2011;López-Cámara et al. 2014).Later, a radiation transfer calculation was incorporated to improve the accuracy in evaluating the emission signatures.The initial study by Ito et al. (2015) applied their Monte-Carlo radiation transfer code (Ito et al. 2013(Ito et al. , 2014) ) to compute the light curves and spectra of LGRB based on hydrodynamical simulations.Following the study, several works have carried out numerical simulations of photospheric emission using a similar numerical method (i.e., post-process Monte-Carlo radiation transfer calculation based on a hydrodynamical simulation of LGRB jets; Lazzati 2016; Parsotan & Lazzati 2018;Parsotan et al. 2018;Parsotan & Lazzati 2021;Ito et al. 2019).The latest simulations have updated the method to incorporate the calculation of polarization (Parsotan et al. 2020;Parsotan & Lazzati 2022).Recently, the numerical technique was also applied to explore the photospheric signals in short GRBs (SGRBs) (Ito et al. 2021).
Although possible effects from non-hydrodynamical dissipation processes (e.g., magnetic reconnection and hadronic collisions) are not taken into account, these simulation-based studies have provided an important insight into the imprint of the realistic jet dynamics on the photospheric emission.A notable finding of these studies is that correlations among the spectral peak energy, E p , isotropic energy, E iso , and peak luminosity, L p , that are compatible with the observed GRB prompt correlations arise as a natural outcome of the simulation.While tension from the observation remains in the simulated spectral shapes, the ability to reproduce various empirical prompt correlations is regarded as one of the major advantages of the photospheric emission model (see Parsotan & Ito 2022, for a recent review on the prompt correlations and their theoretical interpretations).
The aim of the current paper is to gain further insight into the properties of photospheric emission in LGRBs.To this end, we extend our previous study (Ito et al. 2019, hereafter, IMN19) and carry out a more comprehensive analysis of the emission properties.In IMN19, we explored light curves and time-integrated spectra for three sets of simulations and mainly focused on the correlation between the spectral peak energy E p and peak luminosity L p (Yonetoku relation).In this paper, we further explore various prompt correlations and conduct a detailed analysis of the spectrum and polarization, which were not investigated in IMN19.While the same set of hydrodynamical simulations is considered as in IMN19, we enlarged the number of photons tracked in the Monte-Carlo radiation transfer calculation by a factor of ∼ 10.The increase in the photon statistics is particularly important for resolving spectral shapes and polarization properties with reasonable accuracy.
This paper is organized as follows.In Section 2, we describe our model and numerical procedures employed in our calculations.In Section 3, we show the transverse structure of the jet obtained from our hydrodynamical simulations, which plays an essential role in the emission properties.Section 4 summarizes the correlations found in the time-integrated emissions.The results of the timeresolved spectral and polarization analysis are presented in Section 5 and Section 6, respectively.Section 7 is devoted to a summary and discussions.

MODEL AND METHODS
As described in IMN19, we analyze the photospheric emissions evaluated in 3 simulations that employ different jet power.The procedure is as follows: First, we compute the evolution of a relativistic jet by performing a 3D hydrodynamical simulation.Then radiation transfer within the jet is solved as a post-process based on the Monte-Carlo method.The details of the numerical method and setup are described in IMN19.The differences from the previous paper are that we have increased the number of photon packets by an order of magnitude in each simulation (from ∼ 5 × 10 8 to ∼ 5 × 10 9 ) and also incorporated the polarization calculation in the radiation transfer calculation.As for the polarization calculation, we follow the method described in Ito et al. (2021) (see also Ito et al. 2014, for further detail).Here we give a brief overview of the numerical setup.
In the hydrodynamical simulations, we inject a relativistic jet inside a massive star which eventually breaks out from the stellar envelope.The subsequent evolution is solved up to the point where it becomes optically thin to radiation. 11We consider three models that employ different jet kinetic powers: L j = 10 49 , 10 50 , and 10 51 erg s −1 .In all simulations, the jet is continuously injected for a duration of t = 100 s with an initial opening angle of θ ini = 5 • and an initial Lorentz factor of Γ ini = 5, at the inner boundary of the simulation, which is located at a radius of r = 10 10 cm.Note that the jet expands to an opening angle of ∼ θ ini + 0.7Γ −1 ini before it encounters the first recollimation shock (Harrison et al. 2018).The initial specific enthalpy is set as h ini = 100 for the models with L j = 10 49 and 10 50 erg/s, while h ini = 180 is adopted in the model with L j = 10 51 erg/s.As for the progenitor star model, we employ model 16TI provided by Woosley & Heger (2006).The breakout of the jet head from the stellar envelope takes place approximately at t ≈ 25, 8, and 4 s in the models with jet powers of L j = 10 49 , 10 50 , and 10 51 erg s −1 , respectively.In line with previous studies on collapsar jets (e.g., Lazzati et al. 2009;Gottlieb et al. 2019), our models also demonstrate that the jet propagates at subrelativistic speeds before the breakout.Nonetheless, in contrast to simulations of similar setups (jet power and stellar radius), our models exhibit a moderately reduced breakout time, which is likely due to the large jet injection radius imposed in the current study.
The Monte-Carlo radiation transfer simulation solves the evolution of the photon packets initially injected at highly optically thick regions.Each packet is treated as an ensemble of photons with the same frequency ν and, as a whole, carry the Stokes parameters: I, the intensity, and Q and U , the two parameters which characterize the linear polarization.Here, the intensity parameter is defined as I = n pack hν, where n pack is the number of photons in the packet.In the current simulation, n pack does not vary throughout the evolution, and an equal number is imposed among all photon packets.Hence, I depends only on ν.From the Stokes parameters, the linear polarization can be quantified: Π = Q 2 + U 2 /I and ψ = 1/2 arctan(U/Q) give the degree and the position angle of linear polarization, respectively.We do not take into account the Stokes parameter V , the circular polarization parameter, because it is irrelevant in the current simulation.12Regarding the coordinate system employed to define Q and U , we use the same prescription as in Ito et al. (2014) (see their Figure 4).Hence, Q > 0 (Q < 0) and U = 0, i.e., ψ = 0 • (90 • ), corresponds to the case when the electric vector of the polarized beam is parallel (perpendicular) to the plane formed by the line of sight (LOS) of the observer and the jet axis.At the injection, photon packets are set to be unpolarized (Q = U = 0), and their frequency is drawn randomly with a probability proportional to the Planck distribution of the local temperature.After the injection, we selfconsistently calculate the evolution of ν (or equivalently I), Q and U due to scatterings by thermal electrons until the packet reaches the outer boundary of the simulation.
Throughout the paper, the observer's location is expressed by the viewing angle θ obs , which is the angle between the LOS and the central axis of the jet.Since our calculation is performed in 3D, the emission depends not only on the viewing angle but also on the azimuthal angle.However, the dependence turns out to be weak in the current models.Therefore, we show the results for a fixed azimuthal angle, i.e., all LOSs are aligned on a single plane.In our current simulations, we adopt a spherical coordinate system (r, Θ, Φ), with the central axis of the jet aligned along the direction of Θ = 90 • and Φ = 90 • .Hence, the plane of the LOS corresponds to the Φ = 90 • plane.Hereafter, we focus on this plane and express the zenith angle of the simulation as θ, with θ = 0 • aligned to the jet axis, i.e., θ = Θ − 90 • .The viewing angle θ obs is defined in exactly the same way as θ, so that θ obs = Θ − 90 • .

TRANSVERSE STRUCTURE OF JET
Before presenting the result of the spectral and polarization analysis, let us give a brief overview of the transverse structure of the jet found in the hydrodynamical simulation since it plays a key role in determining the properties of the emission.In Figure 1, we show the polar distribution of the accumulated isotropic equivalent energy E j,iso (red) and the energy-weighted average value of the Lorentz factor Γ (blue) evaluated at r = 10 13 cm.These quantities are computed as functions of θ: where L j,iso = 4πr 2 Γ 2 [ρc 2 (h − 1/Γ)]β j cosθ β c is the isotropic equivalent jet kinetic power, with c being the speed of light, and ρ and β j denoting the mass density and the velocity normalized by the speed of light, respectively.In the above equations, all physical quantities are evaluated at a given zenith angle θ on a local mesh located at r = 10 13 cm.Here, θ β is the angle between the direction of velocity and radial direction.The time when the forward shock reaches r = 10 13 cm is denoted by t 13 , and the integration is performed from that time for a duration of t max .Hence, the time integration implies the accumulation of the outflow component that passes through the surface at r = 10 13 cm for a duration of t max .This outflow component corresponds to the region emitting radiation from t obs = 0 to ≈ t max .Note that the radius r = 10 13 cm is significantly larger than the stellar radius (≈ 4 × 10 10 cm), implying that the jet is in the phase of free expansion after breaking out from the star.The various lines in the figure correspond to the difference in the duration of the time integration, t max , for the computation of E j,iso and Γ .As found in the previous studies (e.g., Gottlieb et al. 2021a), the outflow structure shows a core region at the center in which E iso,j and Γ are close to uniform that is surrounded by a wing region which shows a rapid decline in these quantities.Such a structure is formed due to the interaction between the jet and the massive stellar envelope.As shown in IMN19 and in the current paper, the transverse structure gives rise to the spectral-luminosity correlation.
As seen in Figure 1, while the opening angle of the core, θ core , remains nearly constant throughout the entire duration for the model with L j = 10 49 erg/s, it rapidly broadens at t ∼ 30 s and ∼ 10 s for models with L j = 10 50 erg/s and 10 51 erg/s, respectively.The sudden spread in these models reflects the time the innermost recollimation shock 13 breaks out from the stellar envelope; therefore, confinement by the cocoon pressure becomes less effective at later time.Before the transition, the opening angle of the core is θ core ∼ 2 − 4 • in all models, which is consistent with the scaling found in previous studies, &   13 As found in the studies involving hydrodynamical simulations of collaspsar jets (e.g., Lazzati et al. 2009;Gottlieb et al. 2019), the jet undergoes the formation of multiple recollimation shocks before it successfully breaks out of the stellar envelope.Here, the term "innermost" recollimation shock refers to the first shock encountered by the jet after its injection at the base.Eventually, this innermost shock also breaks out the stellar envelope, leading to the emergence of the unshocked (uncollimated) jet from the envelope.
Polar distribution of the accumulated isotropic energy E j,iso (red) and the average Lorentz factor Γ (blue) of the outflow evaluated at r = 10 13 cm in a 2D slice taken through the midplane of the hydrodynamical simulation.The dashed, solid, and dotted lines for E j,iso correspond to the values obtained by integrating over the duration tmax = 10, 30 and 100 s, respectively.The dashed, solid, and dotted lines for the Lorentz factor Γ are the average values weighted with the energy of the outflow over the same duration.
Ioka 2013; Gottlieb et al. 2021a).After the breakout of the recollimation shock, the unshocked jet matter that expands following the recollimation shock expells from the stellar envelope and can be visible in the prompt phase as shown in the previous studies (Ito et al. 2015, IMN19) and also in the current paper.

CORRELATIONS IN THE TIME-INTEGRATED EMISSION PROPERTIES
4.1.The Yonetoku (E p − L p ) and Amati (E p − E iso ) correlations In IMN19, we have shown that our simulations reproduce the correlation between the peak energy of the time-integrated spectrum, E p , and the peak luminosity of the light curve, L p , (Yonetoku et al. 2004) quite well.We reexamine the spectral correlation with the current data sets with ∼ 10 times better statistics.The result is displayed in the left panel of Figure 2. As in IMN19, together with the E p and L p that are sampled from the entire emission duration (from t obs = 0 to t obs = t obs,max = 100 s), we also show the cases in which only emission up to a specific time (t obs,max = 10 and 30 s) is considered.This is intended to mimic bursts originating from shorter jet activity since GRBs have diverse durations.As expected, the E p − L p correlation found in IMN19 remains unchanged, confirming the statistical convergence of the result presented in IMN19.It is noted that the slight change in the distribution from that displayed in IMN19 reflects the difference in the choice of the emission durations (t obss,max = 20, 40, and 110 s are employed in Figure 3 of IMN19).In the right panel of Figure 2, we also examine the correlation between E p and the time-integrated energy of the emission, E iso .Our simulation results also agree with the observed E p − E iso correlation (Amati et al. 2002), although there is a slight offset toward higher E iso .
As discussed in IMN19, these correlations arise mainly due to the viewing angle effect.While the viewing angle is within the core of the jet (θ obs θ core ), E p , L p , and E iso do not show significant variation.On the other hand, beyond the core (θ obs θ core ), these quantities simultaneously show rapid decline with the viewing angle.As a result of the viewing dependence at the wing region, a positive correlation is found among E p , L p , and E iso .Regardless of the difference in the jet power spanning two orders of magnitude, all three models agree with the observed correlations.The results are broadly consistent with the independent studies of LGRBs (Parsotan & Lazzati 2018;Parsotan et al. 2018) and also SGRBs (Ito et al. 2021), indicating that the correlations are a robust feature in the photospheric emission arising from a hydrodynamic jet.The theoretical interpretation for this outcome will be explored in future work.
4.2.Prompt emission (E p and L p ) -Lorentz factor (Γ 0 ) correlations As a complementary study, we also compare our results with the empirical correlations between the emission properties (E p , L p , and E iso ) and the bulk Lorentz factor of the outflow, Γ 0 , reported in the literature (e.g., Liang et al. 2010Liang et al. , 2013Liang et al. , 2015;;Ghirlanda et al. 2012;Lü et al. 2012).Figure 3 displays E p −Γ 0 and L p −Γ 0 diagrams of the simulation results together with those taken from the recent observational study by Ghirlanda et al. (2018).
In our simulations, Γ 0 is defined as an energy-weighted average value of the terminal Lorentz factor hΓ along the line of sight.It is computed in the same manner as Equation ( 2) and can be expressed as follows: When depicting the E p (L p ) -Γ 0 relationship in Figure 3, we set θ and t max to coincide with the viewing angle θ obs and the upper limit of t obs , used for determining E p (L p ), respectively.As seen in the figure, our simulation shows a very good agreement with the observed L p − Γ 0 correlation.14While there is a small offset from the center of the observed distribution, the E p − Γ 0 diagram of our simulation also lies within the dispersion of the empirical correlation.
It should be noted, however, that there are caveats for this comparison.One is that Γ 0 is not directly observable, and the values estimated in the literature are deduced by interpreting the hump (onset) in the early afterglow light curve as the deceleration time of the external shock.Therefore, the estimation relies on the correctness of the afterglow model that describes the onset feature.Hence, the indirect measurement has an uncertainty associated with any possible ambiguity in the theoretical modeling (see, e.g., footnote 16 in Parsotan & Ito 2022).The possible error due to ambiguity in the modeling is partially visible in Figure 3.As seen in the figure, different assumptions about the circum-burst medium (CBM), i.e., homogeneous (grey circles) or wind-like (grey stars), result in a systematic shift in the distribution.
Another caveat is on the estimation of Γ 0 based on the simulation.Here we have evaluated the quantity as an energy-weighted average value of the terminal Lorentz factor, hΓ, along the line of sight.Regarding the core region of the jet, the value is likely to result in an overestimate (by a factor of ∼ 2 at most for the current set of simulations) since the photosphere is located in the region at which the flow is still in the acceleration phase. 15n the other hand, the photosphere is located above the acceleration region in the wing.Hence, one can safely determine the terminal Lorentz factor.However, since Γ 0 sharply drops with angle, it is not clear whether the onset time measured by the observer viewing the wing region (θ obs > θ core ) reflects the deceleration time of the jet matter along the line of sight.To remove the theoretical uncertainty, the dynamics and the emission of the structured jet during the onset of the deceleration must be solved.This is beyond the scope of this study.
Nevertheless, keeping the above uncertainties in mind, we can conclude that no severe tension exists between the correlations found in the observations and the simulation.A similar result is obtained in Lazzati et al. (2013), which finds a broad agreement between the observed E iso − Γ 0 correlation and that evaluated based on the simulation.The current study, which employs a  Here, Γ 0 in our simulation is evaluated using Equation (3), and the duration of the integration, tmax, is chosen to coincide with the emission duration.
For comparison, we display observed data on LGRBs from Ghirlanda et al. (2018).This study considered two different assumptions for the density profile of the circum-burst medium (CBM), which follows a form ρ ∝ r −s , to derive Γ 0 .The data is depicted by grey circles for the case of homogeneous density profile (s=0) and stars for the case of a wind-like density profile (s=2).The best-fitting function of the Ep − Γ 0 and Lp − Γ 0 correlations for the two different assumptions are indicated with solid grey lines.
more accurate calculation method, confirms their conclusion that photospheric emission naturally gives rise to the correlation between Γ 0 and the prompt properties.

TIME-RESOLVED SPECTRAL ANALYSIS
In this section, we present the result of the timeresolved analysis.The analysis is similar to those performed by Parsotan & Lazzati (2018); Parsotan et al. (2018).We employ time bins of ∆t = 2s and all spectra are fitted either by Band function (Band et al. 1993) or cutoff power law (CPL) function (Yu et al. 2016) using the χ 2 minimization method.Here, we do not perform a goodness-of-fit evaluation between the two functions.Instead, we adopt the following procedure.In all cases, we first fit the spectrum with Band function, which yields the peak energy, E p , and the photon indices below and above the peak energy, α, and β.When the high energy photon index is smaller than a certain value (β −4.5), the spectral shape above the peak is mostly consistent with an exponential cut-off.For such cases, we employ the fit by the simpler CPL function, which yields E p and α.In Figure 4, we display the results of the spectral fits for the fiducial model (L j = 10 50 erg/s).

Correlation between E p and L
As shown in the previous studies, the peak energy E p and luminosity L tend to decrease with the viewing angle due to the jet structure formed by the interaction with the stellar envelope.The emission at θ obs 3 • does not show a significant difference in E p and L since it originates from the core region of the jet in which the overall jet kinetic power and Lorentz factor are comparable.On the other hand, a sharp decline of these quantities is seen at a larger viewing angle, reflecting the drop in kinetic power and Lorentz factor in the wing region (see middle panel of Figure 1).
As seen in Figure 4, time evolutions of peak energy and luminosity generally show a correlation (E p − L tracking behavior).The clear exception is found at t obs ∼ 20−30 s for θ obs 3 • , where the rapid drop in luminosity is accompanied by the hardening of E p .Thereafter, variability in the emission also vanishes.This sudden transition in the behavior is due to the fact that, while earlier emission originates from a front portion of the jet which is subject to shocks and mixing due to interaction with a dense stellar envelope, the later emission comes from the inner unshocked jet which is causally disconnected to the exterior gas (for detail, see Ito et al. 2015, IMN19).For a wider viewing angle (θ obs 3 • ), the sharp transition is not present since all the emissions are from the shocked region.Nonetheless, the later emission tends to be harder because the jet-stellar interaction is weaker in the later phase.
The other two models also show similar behavior.However, the transition is less prominent even at the core region (θ obs 3 • ) in the lower power jet model (L j = 10 49 erg/s).This is because the entire part of the jet is subject to shocks and mixing.On the other hand, the sharp transition is present for higher jet power (L j = 10 51 erg/s), but the transition occurs at an earlier time (t obs ∼ 10 s).Also, since there is a rapid broadening of the unshocked jet after the transition in the high-power jet model (see bottom panel of Figure 1), a steady hard emission can be observed up to θ obs ∼ 10 • .
It is noted, however, that such a transition may not be present in the observations (see Yoshida et al. 2017, for an attempt to identify the transition feature in the GRB light curves).According to the recent simulation by Gottlieb et al. (2019), which solves the jet dynamics from deeper inside the star with a higher spatial resolution, the recollimation shock does not breakout from the stellar envelope within a reasonable range of parameters.Their result suggests that the time of the transition found in the current simulation underestimates the actual value due to the large injection radius and lack of spacial resolution.Hence, the transition is unlikely to occur within the intrinsic durations in the majority of GRBs (e.g., Zitouni et al. 2018).This suggests that the global hardening should be insignificant; therefore, the E p − L tracking feature is expected to be more robust throughout the emission.
As found in previous studies (Parsotan et al. 2018), the correlation between E p and L shows a good agreement with the observations (Golenetskii et al. 1983;Ghirlanda et al. 2010;Lu et al. 2012). 16In Figure 5, we plot the time-resolved E p and L of the three simulations.It is indicated in Ghirlanda et al. (2010); Lu et al. (2012) that the relation between E p and L of the time-resolved spectra is consistent with that for the time-integrated spectra among the GRBs (Yonetoku relation; Yonetoku et al. 2004Yonetoku et al. , 2010)).The figure shows that our simulation results are in good agreement with the relation.While the viewing angle dependence gives rise to the overall correlation, the E p − L tracking feature at a given viewing angle (single burst) also roughly follows the correlation.Note that the observational data are sparse at L 10 50 erg/s where the simulation results, and possibly also observations, depart slightly from the predictions of the Yonetoku relation.Our simulation data, when fitted to the function L = KE p s , yields a best-fit value of s = 1.87 ± 0.01.The slope is slightly steeper than that of the Yonetoku relation, which is reported as s = 1.6 ± 0.082.When excluding the data with low L or E p values, the slope aligns more closely with the Yonetoku relation.For instance, if we discard the data where L < 10 50 erg/s, the best-fit slope adjusts to s = 1.52 ± 0.03.As discussed in IMN19, this may indicate that the slope of the correlation curve changes at low luminosities.

Spectral shapes
As found in the time-integrated analysis (Ito et al. 2015, IMN19), the time-resolved spectra tend to broaden with viewing angle (Figure 4).Hence, the fitted value of the low energy photon index α tends to decrease with viewing angle.The resulting range of α is softer than that for a blackbody (α = 1).Still, it is harder than the typical observed value (the time-resolved analysis of Yu et al. 2016, found a mean and median value of α ∼ −0.77) (mean and median value of α ∼ −0.77 is found in the time-resolved analysis by Yu et al. 2016) The spectrum above the peak is mostly consistent with exponential cutoff at the core part of the jet (θ obs 3 • ), while a hard energy extension appears at a larger viewing Ito et al.The solid black lines show the light curves with time bins of ∆t = 0.5 s, while the blue, red, and green circles plot the best-fit values for the spectral peak energy Ep, low energy photon index α and high energy photon index β measured in time bins of ∆t = 2 s, respectively.Note that the parameter β is not present in the time intervals where the spectrum is fitted using the CPL function.angle.Therefore, most of the time bins are fitted by the CPL function for the former case, while Band function is used for the latter case.Note that the high-energy photon index β (indicated by the green circle) is absent from the figure for the time intervals where the spectrum is fitted using the CPL function.
For illustrative purposes, we plot the spectrum of the fiducial model (L j = 10 50 erg/s) at a specific time bin for various viewing angles in Figure 6 together with the corresponding best fit analytical function (either CPL or BAND).The shift toward a broader spectrum at a larger viewing angle can be confirmed in the figure.As described in previous work (Ito et al. 2015, IMN19), inhomogeneities in the outflow cause the spectral broadening via the multi-color effect and bulk Comptonization.These effects combine to soften the low-energy spectrum and harden the high-energy spectrum (Lundman et al. 2013;Ito et al. 2013).
We summarize the best-fit spectral parameters for all simulation data in Figure 7.The broad distribution of the peak energy ranging from ∼ keV up to ∼ 10 3 keV mainly reflects the wide range of the viewing angles taken into account in the current study (see Figure 5).Note that the clustering at the highest energy E p ∼ 800 − 900 keV arises from the simulation of the highest jet power L j = 10 51 erg/s.As mentioned earlier (Section 5.1), quasi-steady emission is seen at a later time (t obs 10 s) in a wide range of viewing angles in this simulation.As a result, the data accumulates to produce the peak at the high-energy end of the distribution.
Examining the distribution of spectral indices, we observe that our simulations exhibit a narrow spectrum, characterized by larger values of α and smaller values of β, in contrast to what is typically observed.While the observation shows clustering around ∼ −0.7 for α, our simulations find a sharp peak at around 0.5 above which the number of data drops sharply.The existence of such a cutoff-like feature at α ∼ 0.5 is the expected result for non-dissipative photospheric emission from a spherical outflow (Beloborodov 2010;Bégué et al. 2013;Lundman et al. 2013;Ito et al. 2013).This implies that the multi-color temperature effect is not significant for the population around the cutoff.The high energy photon index is typically in the range −2.5 β −2.0 for observation, while the simulation results show a broader distribution that extends down to smaller values.The broad extension in β arises due to the gradual hardening in the high energy spectra with viewing angle (see Figure 6).
While there is a clear gap between the distributions of α and β for the simulation and those for the observation, some fraction shows an overlap.It is noted, however, that most of the overlapping population in the simulation corresponds to the emissions originating from a wing region that exhibits low E p (< 100keV).Therefore, while this population may be compatible with X-ray-rich (XRR) GRBs or X-ray flashes (XRF) (Katsukura et al. 2020), it does not match with the typical GRB observations.
The tension in the spectral shapes is commonly found in the simulation-based study of photospheric emission.As discussed in our previous studies (Ito et al. 2015, IMN19), calculation with higher spatial resolution may reduce the tension by resolving sharp shear flows, which can enhance the effects of the bulk Comptonization and the multi-color effects.However, while the tension in the high energy photon index may be reduced by introducing intermittent central engine activity (Parsotan et al. 2018), similar results are found in the higher resolution studies in 2D (Lazzati 2016;Parsotan & Lazzati 2018;Ito et al. 2021;Parsotan & Lazzati 2021).Although further 3D study is necessary to confirm this issue, these results suggest that additional physical processes not captured in the global hydrodynamical simulation must be invoked to reproduce the broad spectrum.A potential solution to this problem is a dissipative process that acts around the photosphere.As discussed in Gottlieb et al. (2019Gottlieb et al. ( , 2021)), the spectrum can be broadened by internal shocks that form due to rapid variability that is smeared out in the simulations (Bromberg et al. 2011;Ito et al. 2018;Samuelsson et al. 2022;Samuelsson & Ryde 2022).Moreover, additional internal dissipation that results in a significant non-thermal broadening may arise from other processes, such as magnetic reconnection and hadronic collision (e.g., Vurm & Beloborodov 2016).The quantitative estimation of such effects is currently infeasible and is beyond the scope of this study.

POLARIZATION ANALYSIS
In this section, we present the result of the polarization analysis.It should be emphasized that this is the first study in which the properties of the polarization are examined based on 3D simulations.The imposition of 2D axisymmetry which is often invoked in the literature (Lundman et al. 2014;Ito et al. 2014;Parsotan et al. 2020;Ito et al. 2021;Parsotan & Lazzati 2022), restricts the position angle ψ = 1/2 arct(U/Q) to be one of two orthognal angles, i.e., 0 • (180 • ) or 90 • for the current definition of the Stokes parameter (Section 2), and also forces Π at θ obs = 0 • to vanish.Therefore, a study in 3D is essential to remove these artificial effects on the results of polarization analysis.In Figure 8, we show the time-integrated properties of polarization at three different energy bands (hν = 1−10, 10 − 100, and 100 − 1000 keV) as a function of viewing angle.Following the analysis of the prompt correlations (Section 4), we consider three cases for the emission duration for our time-integrated polarization analysis: t obs = 0 s up to t obs,max = 10, 30, or 100 s.In the figure, the 1σ statistical uncertainty of the degree of polarization (Π = Q 2 + U 2 /I) and position angle (ψ = 1/2 arctan(U/Q)) is determined using error propagation.Specifically, we use the variance and covariance of the Stokes parameters I, Q, and U , denoted by σ 2 I/Q/U and σ IQ/IU/QU , to quantify the variance of Π and ψ as and  (2016).The blue and green lines in the middle and bottom panels distinguish the population of the simulated spectrum with peak energy higher and lower than 100 keV, respectively.
respectively.To visualize the change in the spectral properties with viewing angle, we also display E p in the figure.
As described in the previous section, the viewing angle beyond which E p shows rapid decrease marks the transition from the core to the wing region as the origin of the emission.
As seen in the figure, when the observer LOS is within the core region (θ obs θ core ), the degree of polarization, Π, exhibits a low value: ∼ 4% at most regardless of the energy band.Once the viewing angle enters the wing region (θ obs θ core ), Π tends to increase.A similar trend is found in the 2D simulation of SGRBs conducted in Ito et al. (2021) (see also Parsotan & Lazzati 2022).As discussed in that paper, the sharp lateral gradient of the outflow structure in the wing region results in an enhanced level of anisotropy of emission around the LOS.Hence, the polarization signature at the wing turns out to be larger than that at the core region.The main difference from Ito et al. (2021) is that the current simulations have more complex viewing angle dependence (wiggling features).This presumably reflects the difference between the nature of 3D and 2D simulation, in which the former shows more complex flow structures.It is also worth noting that in the current 3D simulations, Π exhibits a non-zero value at θ obs = 0 • , although small ( 1 %), indicating a departure from axisymmetry.
In Figure 9, we demonstrate how the anisotropy of the emission around the LOS is linked with the degree of polarization Π by considering the model of L j = 10 50 erg/s with t obs,max = 30 s as an illustrative case.The top panel of the figure shows Π as a function of θ obs , which is identical to the right half part of the central panel of Figure 8.On the other hand, the bottom panel shows the difference between the energy-weighted average value of the zenith angle of the last scattering position, Θ sc , and θ obs .If the emission is isotropic around the LOS, Θ sc = θ obs holds.In contrast, the offset between Θ sc and θ obs tends to enlarge as the anisotropy grows.Therefore, the offset of the two angles, Θ sc − θ obs , can be regarded as an indicator of the anisotropy.As shown in the figure, the offset Θ sc − θ obs is modest at the core part of the jet (θ obs 3 • ) since the emission is close to symmetric around the LOS.Note that the slight offset of the average last scattering position toward a larger angle ( Θ sc −θ obs > 0) seen at the low energy band (1-10 keV) reflects the fact that the outer part of the jet tends to show softer emission (i.e., intensity around the LOS at low energies (hν < E p ) tends to increase at larger Θ sc ).Looking at larger viewing angles θ obs 3 • , the offset of the average last scattering position toward the smaller angle ( Θ sc − θ obs < 0) is seen and tends to increase with θ obs .This development of anisotropy arises due to the sharp negative gradient of the intensity in the lateral direction at the wing region.The link between the anisotropy of the emission and the degree of polarization Π can be confirmed in the figure: Π tends to increase as the offset between Θ sc and θ obs becomes larger.
The above description of the anisotropy of the emission and the polarization is not limited to this illustrative case but is a general trend in all cases explored in the current study.Hence, as mentioned while the viewing angle is within the core, a low level of polarization is found, and Π tends to increase with θ at a larger viewing A schematic image for the anisotropy of  4) and ( 5).For reference, the black dotted lines display spectral peak energy Ep.The left, middle, and right column corresponds to the simulation with jet power of L j = 10 49 , 10 50 and 10 51 erg/s, respectively.The top, middle, and the bottom row shows the cases in which emission up to t obs = 10, 30, and 100 s are taken into account, respectively.We discard the viewing angles at which fewer than 100 packets were observed because of the potential for large statistical errors.Hence, the limited range of viewing angles for some curves reflects the paucity of photon packets.As for the position angle, we also discard several plots which show 1σ error bar larger than 40 • .Note that the scale of the vertical axis of Π differs among the panels.  .Schematic picture of the intensity distribution for a given observer's LOS.In the picture, we show two cases: one in which the LOS lies within the jet core (θ obs < θcore: blue) and outside the core (θ obs > θcore: green).The intensity around the given LOS is expressed by a color contour in red with a darker color corresponding to a larger intensity.As drawn in the picture, intensity is close to uniform at θ obs < θcore, while a concentration toward a direction of the jet axis is seen at θ obs > θcore.Note that the picture is a simplified illustration of the overall trend, and the distribution of intensity in the simulation is more complex.For example, anisotropy in the azimuthal direction is also present in the simulation since the flow is not axisymmetric.
We find no general trend regarding the dependence of Π on the energy band.When the viewing angle is small (θ obs θ core ), lower energies (hν = 1 − 10 keV) tend to show the largest Π due to the relatively enhanced anisotropy mentioned earlier.On the other hand, the highest energy (hν = 100 − 1000 keV) tends to show the highest Π as the viewing angle increases.This is because the intensity at higher energy (hν > E p ) shows a sharper negative gradient of intensity in the lateral direction.While Π is in the range of a few % in most viewing angles, it can be larger than 10 % in the wing region, particularly at high energies well above the peak ( E p ) and can be as high as ∼ 20 − 40% in extreme cases.The non-monotonic energy dependence of Π is similar to that found in Ito et al. (2021).
Comparing the results among the different choices of the integration time t obs,max , the change in width of the core can be observed in the viewing angle dependence of Π.As described in Section 3, the core width increases at later times.Hence, the region that shows low polarization tends to extend to a larger θ obs as the integration time increases.This trend is much more pronounced for the high-power jet model (see Section 3 for the physical reason), as shown in Figure 8.We also find that, for a longer t obs,max , Π tends to be larger, particularly at the wing region.This is because the signature of the interaction between the jet and the stellar envelope becomes less significant at later times.As a result, the sharp velocity gradient structure at the wing becomes relatively more stable, which in turn leads to the enhancement of polarization.
It should be again noted, however, that the rapid expansion of the jet core may not take place for a higher spatial resolution simulation with a deeper injection radius (see Section 5.1).A more conservative treatment would be to consider only the emission before the substantial expansion: t obs 30 s for L j = 10 50 erg/s and t obs 10 s for L j = 10 51 erg/s.Even if the affected emission times are excluded, our overall results remain valid: polarization is higher at large viewing angles because the velocity gradient is steep in the wings of the jet.
Looking at the position angle ψ = 1/2 arctan(U/Q) in Figure 8, a clustering around ψ = 0 • (180 • ) and ψ = 90 • is seen: the electric vector of the polarized emission is mostly inclined or perpendicular to the plane formed by the LOS and the jet axis.As mentioned above, ψ is restricted to these two directions when axisymmetry is imposed.Hence, the result indicates that the departure from the axisymmetry in the current 3D simulations is not significant enough to wash out this feature.At the same time, however, a certain deviation from the axisymmetry is also confirmed by the intermediate values of ψ found in the analysis.The departure from axisymmetry is more evident for the case of lower power jet model (L j = 10 49 erg/s) since jet dynamics shows wobbling behavior, which is induced during the jet-stellar interaction due to its lower thrust.It should also be noted that, as in the case of Π, the position angle shows a nonmonotonic energy dependence.While ψ can be aligned in all energy bands in some cases, the difference among the energy band can be as large as 90 • .Such a feature is also indicated in Ito et al. (2021).

Ito et al.
rectly reflect the dynamical changes in the flow structure around the photosphere.In Figure 11, we show the degree and position angle of the polarization as a function of time at different viewing angles for models with jet powers of L j = 10 49 and 10 50 erg/s.The corresponding light curves are also displayed in the figure for comparison.The figure confirms the presence of non-negligible time variability in Π and ψ.It also shows that there is no clear correlation between the temporal behavior of the polarization with the light curves.
Similar to the results of the time-integrated analysis, the time-resolved values of Π tend to increase with θ obs .The non-monotonic dependence of Π on the energy band found in the time-integrated analysis is also confirmed in the time-resolved analysis.It should be noted, however, some level of correlation is seen in the temporal behavior of Π among the different energy bands.For example, at relatively large viewing angles θ obs 4 • , Π at all energy bands tends to increase at later times.As described in Section 6.1, this is because the jet-stellar interactions become less significant at later times.
Regarding the temporal evolution of ψ, it is concentrated around the two orthogonal directions, ψ = 0 • (180 • ) and ψ = 90 • , again indicating that the break of axisymmetry is not so significant.Simultaneously, though, our simulations show temporal change between the two orthogonal directions, i.e. a ∼ 90 • flip in ψ.Two examples of this behavior are visible in Figure 11: in the time period t obs ∼ 20 − 30 s for the fiducial model (L j = 10 50 erg/s), both the θ obs = 6 • (middle right) and θ obs = 8 • (bottom right) show such a flip.Within that period, ψ sharply changes from 180 • (0 • ) to 90 • and, after retaining its value for a few seconds, it again sharply changes back to 180 • (0 • ).Interestingly, while the overall evolution does not suggest a correlation with the luminosity as mentioned above, the rapid rotation is associated with the spike in the light curve.As mentioned earlier, the deviation of axisymmetry is larger in the lower jet power model (L j = 10 49 erg/s).This is evident from the fact that the concentration of ψ in the two orthogonal directions is less prominent in this model (left panels of Figure 11).Accordingly, the rotation of ψ also takes place relatively gradually compared to higher jet power models.
The overall temporal behavior of ψ at different energy bands appears to be somewhat correlated.Its value, however, can vary broadly, and the difference can be as significant as 90 • as already found in the time-integrated analysis.These results suggest that the non-negligible energy dependence of Π and ψ is an inherent feature of the photospheric emission.

Observational implications
To date, observational data of GRB polarization is limited to a sample of relatively bright bursts.Therefore, to facilitate a comparison with these observations, it is necessary to concentrate on the emission emanating from the jet core with E p > 100 keV, since this population is representative of the typical bright bursts.As shown in the previous sections, the core region exhibits low polarization (Π < a few %).This result is consistent with the recent reports by the POLAR collaboration (Kole et al. 2020) and AstroSAT collaboration (Chattopadhyay et al. 2022), which find that the majority of GRBs they have analyzed are consistent with having a low to null polarization across the full burst duration.On the other hand, some of the bursts detected by AstroSAT and previous polarimetry missions (e.g., GAP on board IKAROS; Yonetoku et al. 2012) showed high levels of polarization, which cannot be explained by the emission emanating from the jet core region of our simulations.
Regarding the time-resolved properties, the large temporal change in the position angle (∆ψ ∼ 90 • flip) seen in the current simulation may be in line with the timeresolved analysis of polarizations in a few bursts (GRB 110826A, 170114A, 160821A: Yonetoku et al. 2011b;Burgess et al. 2019;Sharma et al. 2019).However, these analyses suggest an intrinsically high degree of polarization is preferred in these GRBs, contrary to the intrinsically low polarization revealed in the current simulation.Accordingly, our simulations are also incompatible with the possible interpretation for the origin of low timeintegrated polarizations found in the recent observations as a result of summation of intrinsically highly polarized emission that exhibits a large temporal change in the position angle (Burgess et al. 2019;Kole et al. 2020).
To sum up, our simulations predict a low intrinsic polarization for the emission originating from the core region, which is representative of the population of bright bursts.Therefore, while our results are in agreement with recent observations reporting low polarization for the majority of GRBs, they do not reproduce the high polarization inferred in some bursts.It should be noted, however, that all observations so far accompany large error bars, and no robust consensus on the GRB polarimetry has been achieved.Hence, we cannot draw a firm conclusion based on the available observational data.
In the next decade, significantly more sensitive detectors dedicated to GRB polarimetry, POLAR-2 (Kole 2019) and LEAP (McConnell et al. 2021), along with missions such as COSI (Tomsick et al. 2019, which, al-though not exclusively dedicated to GRBs), are likely to provide further insights.Here let us summarize the key characteristics of the polarization found in the current simulations, which might be tested in future missions.One of the robust features in the current results is the inverse correlation of the degree of polarization Π with the spectral peak energy E p or the radiative output (peak luminosity L p and isotropic energy E iso ).In other words, it predicts that dim, soft bursts observed outside the jet core region (E p 100 keV: XRR GRBs or XRFs) tend to show a higher polarization than the above-mentioned typical GRBs viewed within the core (E p 100 keV).To further examine this aspect, Figure 12 illustrates the relationship between Π, E p , and L p .Although the dispersion is large, the negative Π − E p (L p ) correlation can be confirmed from the figure.
Another characteristic found in the time-resolved analysis is the clustering of ψ in the two orthogonal directions associated with the rapid change between the two directions that occurs occasionally.The increased ability to perform a time-resolved polarization measurement in future missions may enable to probe these features.
Lastly, the non-negligible energy dependence of Π and ψ is also a notable feature found in the current study.Observations within the limited energy range of ∼ 10 − 1000 keV alone by POLAR-2 or LEAP, may not be adequate for carrying out an energy-resolved analy-Figure 11.Degree of linear polarization (Π = Q 2 + U 2 /I : upper panels) and position angle (ψ = 1/2 arctan(U/Q) : lower panels) as function of observer time t obs at energy ranges of hν = 1 − 10 keV (red), 10 − 100 keV (blue) and 100 − 1000 keV (green).The time bin employed for the analysis is ∆t = 2 s.The top, middle, and bottom panels in the left (right) column correspond to the results of the model with jet power of L j = 10 49 erg/s (10 50 erg/s) for viewing angles of θ obs = 1 • , 3 • , and 6 • (3 • , 6 • , and 8 • ), respectively.For reference, the corresponding light curve is displayed by the black dotted lines.The error bars represent the 1σ statistical uncertainty, which are evaluated using Equations ( 4) and ( 5).We discard time bins where fewer than 100 packets were observed, due to the large statistical uncertainty.As for the position angle panels, we also omit points that show 1σ error bars larger than 40  sis.Hence, it is important to perform joint observations across different energy bands with other proposed missions, such as utilizing COSI and AMEGO at higher energies (∼ 100 KeV − 5 MeV), along with instruments operating at lower energies like the Low Energy Polarimetry Detector (LPD) (see Section 7.1 of Gill et al. for details), in order to investigate the energy dependence of polarization.

SUMMARY AND DISCUSSIONS
In the present study, we have analyzed the properties of the photospheric emission in the context of LGRBs by conducting 3D relativistic hydrodynamic simulations and post-process Monte-Carlo radiation transfer calculations.This is an update of our previous work (IMN19) in which the number of photon packets tracked in the radiation transfer calculation was smaller by an order of magnitude.While we focused on light curves and timeintegrated spectral properties in IMN19, the improved photon statistics here have enabled us to further explore the time-resolved spectral properties and the polarization properties.
We have considered three models in which the relativistic jet is injected at the interior of the massive stellar envelope with constant power (L j = 10 49 , 10 50 , and 10 51 erg/s) up to 100 s with an opening angle θ ini = 5 • (initially top-hat structure).Under the assumption that the electron distribution is Maxwellian, we have tracked the evolution of the photons subject to numerous electron scatterings before being released at the photosphere.We properly take into account the full Klein-Nishina Compton section in determining the change of the Stokes parameter by the scattering.The summary of the main findings obtained from the current numerical models is listed below.
1. Due to the interaction with the massive stellar envelope, the jet develops a lateral structure in which the core region at the center (θ θ core ) is surrounded by a wing region that shows a rapid decline in power and Lorentz factor with angle (Figure 1).As a result, while the core region exhibits a spectral peak energy E p and luminosity L close to uniform, these quantities sharply decline with the viewing angle for the emission arising from the wing region.This viewing angle dependence at θ obs θ core leads to correlations between E p and the peak luminosity L p as well as between E p and the isotropic energy E iso .The obtained E p -L p correlation agrees well with the Yonetoku relation (Yonetoku et al. 2004), which further validates the key finding of IMN19 (left panel of Figure 2).Additionally, we have demonstrated that the E p -E iso correlation is broadly consistent with the Amati relation (Amati et al. 2002) (right panel of Figure 2).
2. The viewing angle dependence also gives rise to correlations between the prompt emission properties (E p , L p and E iso ) and bulk Lorentz factor of the flow Γ 0 .The correlations are consistent with those inferred from observations (Ghirlanda et al. 2018) (Figure 3).
3. The time-resolved spectral analysis reveals E p -L (luminosity) tracking feature within an individual burst emission (Figure 4), which is also indicated in the observations (Golenetskii et al. 1983).Moreover, the distribution of time-resolved E p -L data extracted from multiple bursts is fairly consistent with the Yonetoku (E p -L p ) relation as found in the study of Lu et al. (2012) (Figure 5).
4. The time-resolved spectral fitting analysis shows that, while a thermal-like spectrum that can be fitted mostly by a cutoff power law (CPL) function is found at the core region, the spectral shape tends to broaden with the viewing angle in the wing region (Figure 6).As a result, while a fraction of low and high energy photon indices found in the analysis can be compatible with those of the observations, it is limited to the emissions arising from the wing region, which represents a population of soft dim bursts (E p < 100 keV: XRR GRBs and XRFs) (Figure 7).Hence, the current numerical models do not reproduce the observed non-thermal spectral features for the emission arising from the core region, which represents a population of the typical GRBs (E p 100 keV).We discuss potential solutions to this problem later in this section.
5. The lateral structure of the jet is also imprinted in the properties of polarization.While the emission from the core region exhibits a low level of polarization (Π 4 %), the degree of polarization tends to increase with the viewing angle for the emission from the wing region and can be as high as Π ∼ 20-40 %, particularly at high energies (hν > E p ) (upper panels of Figure 8).Hence, the current simulations predict that the typical GRBs show lower polarization than their soft dim counterparts (XRR GRBs and XRF) (Figure 12).The low degree of polarization for the typical GRB is consistent with the latest result of POLAR (Kole et al. 2020) and AstroSAT (Chattopadhyay et al. 2022) in which a majority of the observed bright bursts were found to be consistent with having a low to null polarization.On the other hand, it is incompatible with a fraction of GRBs detected by AstroSAT and also with those detected by the earlier polarization mission (e.g., GAP; Yonetoku et al. 2012) in which a high degree of polarization is indicated.However, all measurements have large uncertainty; therefore, we cannot draw a firm conclusion from the comparison.6.While the current simulations are performed in 3D, the resulting jet structure does not show significant deviation from axisymmetry.As a result, the position angle of the polarization ψ tends to be clustered around the two orthogonal directions, which correspond to the restricted angles for axisymmetric outflow (lower panels of Figure 8).The timeresolved analysis finds a presence of ∆ψ ∼ 90 • flip as a result of the transition between the two angles (lower panels of Figure 11).This result may provide a possible origin for the large temporal variation of ψ inferred in a few GRBs (Yonetoku et al. 2011b;Burgess et al. 2019;Sharma et al. 2019), although the intrinsic degree of polarization in the Ito et al.
current models are much lower than those indicated in the observations.7.Both Π and ψ do not exhibit any clear correlation with the photon energy (Figures 8 and 11).Although higher-energy photons tend to show larger Π, this is not always the case.Notably, it is interesting to observe that the difference in ψ between different energy bands can reach up to 90 degrees.Such properties offer an important test for the current numerical models, which can be examined in future polarization missions such as POLAR-2 (Kole 2019) and LEAP (McConnell et al. 2021).
The current and previous studies of photospheric emission based on hydrodynamical simulations (e.g., Lazzati et al. 2013;Parsotan & Lazzati 2018;Parsotan et al. 2018;Ito et al. 2019Ito et al. , 2021) ) have shown that the emission mechanism provides a natural explanation for the origin of empirical correlations found in the prompt emission.The series of studies have also shown, on the other hand, that their spectral shapes turn out to be narrower than those of observed GRBs.The narrow spectral shapes are the main drawback in these global numerical models.As discussed in Section 5.2, the tension may be attributed to the possible sub-photospheric dissipative processes that are not captured in the simulations.Although implementing such effects in global simulations is numerically challenging, further investigation of this issue should be conducted to assess the role of photospheric emission in GRBs properly.
For the first time, we have quantified the GRB polarization signature based on 3D simulations.While the overall properties are similar to those found in the previous 2D simulations (Parsotan et al. 2020;Ito et al. 2021;Parsotan & Lazzati 2022), we also observe some differences due to the lack of axisymmetry (non-zero Π at θ obs = 0 • and ψ not being restricted to two orthogonal angles).The characteristics revealed in the current study may be tested by future polarization missions, which will provide more well-constrained polarization measurements.It should be noted, however, that the inclusion of subphotospheric dissipation mentioned above may influence the resulting polarization properties (Lundman et al. 2018).
Lastly, let us comment on the limitations of the current study (beyond the absence of subphotospheric dissipation).One caveat is the spatial resolution. 17Since the simulations cover a large spatial domain, small-scale 17 Recently, Arita-Escalante et al. (2023) scrutinized the effects of spatial and temporal resolutions (defined by the number of snapshots used per unit time) in hydrodynamical simulations on radiation transfer calculation of MCRaT, a calculation method analogous to the one employed in our current study.As for the temporal resolution, we performed a convergence check using an approach similar to their study.In our simulations, we adopt a temporal resolution of 10 snapshots per second for the initial 300 seconds, which decreases to 1 snapshot per second thereafter.Through comparative analysis with results derived from lower temporal resolutions, we confirm that the radiation properties converge to a consistent solution at the current resolution setting.Similar to the findings reported by Arita-Escalante et al. (2023), our study also indicates that an insufficiently high temporal resolution can induce artificial Comptonization, which consequently hardens the high-energy spectral slope.As for the spatial resolution, we have not undertaken a convergence check.However, while the hydrodynamic properties structures are washed out due to the limited spatial resolution.Consequently, the expected short time-scale variability in the light curve, which arises due to the jetstellar interaction (e.g., Gottlieb et al. 2019), is largely smeared out.Moreover, the jet injection radius employed in the current simulation may not be small enough to achieve numerical convergence, as discussed in Section 5.1 (based on the result of Gottlieb et al. 2019).It is also highly desirable to expand the current models to explore the impact of central engine activity 18 as well as the magnetic field (the jet at its base can be unstable and be chopped for short period of time: Barkov & Baushev 2011).Recent series of 3D hydrodynamical and magnetohydrodynamical simulations by Gottlieb et al. (2020aGottlieb et al. ( ,b, 2021) ) have demonstrated that these ingredients have an important effect on the dynamics and, therefore, on its resulting emission.Improvements to our numerical model to address the issues mentioned above will be the subject of future studies.
have not reached convergence in the current setup, we do not anticipate the emergence of artificial Comptonization reported by Arita-Escalante et al. (2023) in our simulations.This expectation arises from our prescription for radiation transfer calculation: instead of assuming uniform physical quantities within each hydrodynamic mesh (as in their study), we employ linear spatial interpolation from neighboring mesh centers to determine local physical quantities.This approach circumvents abrupt changes in hydrodynamic properties across the mesh which gives rise to the artificial spectral changes.It is worth noting that our test calculations for a relativistic steady spherical outflow confirmed that our interpolation method accurately captures adiabatic cooling effects, even under the constraints of a low spatial resolution. 18The effect of intermittent central engine activity on the photospheric emission has been previously studied by Parsotan et al. (2018) based on 2D hydrodynamical simulation.It should be noted, however, that 3D effects are likely to play a crucial role in the intermittent jet (Gottlieb et al. 2020a).

Figure 3 .
Figure3.Peak energy, Ep (left), and peak luminosity, Lp (right), versus the bulk Lorentz factor, Γ 0 , for different viewing angles (colored symbols).The shape of the symbols represents the difference in the emission duration considered to extract Ep, Lp, and Γ 0 .Here, Γ 0 in our simulation is evaluated using Equation (3), and the duration of the integration, tmax, is chosen to coincide with the emission duration.For comparison, we display observed data on LGRBs fromGhirlanda et al. (2018).This study considered two different assumptions for the density profile of the circum-burst medium (CBM), which follows a form ρ ∝ r −s , to derive Γ 0 .The data is depicted by grey circles for the case of homogeneous density profile (s=0) and stars for the case of a wind-like density profile (s=2).The best-fitting function of the Ep − Γ 0 and Lp − Γ 0 correlations for the two different assumptions are indicated with solid grey lines.

Figure 4 .
Figure 4. Light curves and evolutions of fitted spectral parameters at various viewing angles for the fiducial model (L j = 10 50 erg/s).The solid black lines show the light curves with time bins of ∆t = 0.5 s, while the blue, red, and green circles plot the best-fit values for the spectral peak energy Ep, low energy photon index α and high energy photon index β measured in time bins of ∆t = 2 s, respectively.Note that the parameter β is not present in the time intervals where the spectrum is fitted using the CPL function.

Figure 5 .
Figure5.The relation between spectral peak energy Ep and luminosity L measured in time bins of ∆t = 2 s for various viewing angles.The triangle, circle, and square plot the results for the simulation with jet power of L j = 10 49 , 10 50 and 10 51 erg/s, respectively.

Figure 6 .
Figure 6.Observed spectrum of the fiducial model (L j = 10 50 erg/s) at t obs = 68 − 70 s for various viewing angle.The black circles represent the simulation data, and the error bars display the corresponding 1σ statistical error due to Poisson noise.

Figure 7 .
Figure7.Distributions of the best-fit parameters for the timeresolved spectra.The top, middle, and bottom panels show the spectral peak energy Ep, the low energy photon index α, and the high energy photon index β, respectively.The red histograms summarize the spectra of all simulations.The considered viewing angle spans from θ obs = 0 • to 10 • with an interval of 1 • .The grey histograms are the fitted observed spectral parameters fromYu et al. (2016).The blue and green lines in the middle and bottom panels distinguish the population of the simulated spectrum with peak energy higher and lower than 100 keV, respectively.

Figure 8 .
Figure8.Degree of linear polarization (Π = Q 2 + U 2 /I : upper panels) and position angle (ψ = 1/2 arctan(U/Q) : lower panels) as a function of viewing angle at energy ranges of hν = 1 − 10 keV (red), 10 − 100 keV (blue), and 100 − 1000 keV (green).The error bars for Π and ψ represent the 1σ statistical uncertainty, which is determined using Equations (4) and (5).For reference, the black dotted lines display spectral peak energy Ep.The left, middle, and right column corresponds to the simulation with jet power of L j = 10 49 , 10 50 and 10 51 erg/s, respectively.The top, middle, and the bottom row shows the cases in which emission up to t obs = 10, 30, and 100 s are taken into account, respectively.We discard the viewing angles at which fewer than 100 packets were observed because of the potential for large statistical errors.Hence, the limited range of viewing angles for some curves reflects the paucity of photon packets.As for the position angle, we also discard several plots which show 1σ error bar larger than 40 • .Note that the scale of the vertical axis of Π differs among the panels.

Figure 9 .
Figure9.Degree of linear polarization (Π = Q 2 + U 2 /I: upper panel) and the offset of an energy-weighted average value of the zenith angle of the last scattering position with respect to the viewing angle ( Θsc − θ obs : lower panel) as a function of viewing angle for the model of L j = 10 50 erg/s at t obs ≤ 30 s.As in Figure8, the color represents the different energy ranges, and the error bar in the upper panel indicates the 1σ statistical uncertainty.For reference, the grey horizontal line in the lower panel shows the case in which the offset angle is zero ( Θsc − θ obs = 0).The inset in the bottom panel provides a schematic picture describing how Θsc and θ obs are defined.emissionaround the LOS is displayed in Figure10.The image illustrates how the distribution of the intensity changes with the viewing angle .
Figure10.Schematic picture of the intensity distribution for a given observer's LOS.In the picture, we show two cases: one in which the LOS lies within the jet core (θ obs < θcore: blue) and outside the core (θ obs > θcore: green).The intensity around the given LOS is expressed by a color contour in red with a darker color corresponding to a larger intensity.As drawn in the picture, intensity is close to uniform at θ obs < θcore, while a concentration toward a direction of the jet axis is seen at θ obs > θcore.Note that the picture is a simplified illustration of the overall trend, and the distribution of intensity in the simulation is more complex.For example, anisotropy in the azimuthal direction is also present in the simulation since the flow is not axisymmetric.

Figure 12 .
Figure 12.Degree of linear polarization (Π = Q 2 + U 2 /I) as a function of spectral peak energy Ep (left) and peak luminosity Lp (right) at energy ranges of hν = 1 − 10 keV (top), 10 − 100 keV (middle), and 100 − 1000 keV (bottom).Error bars indicate the 1σ statistical uncertainty evaluated as described in the caption of Figure 8.The color of symbols indicates the difference in the injected jet power: red, blue, and green correspond to L j = 10 49 , 10 50 , and 10 51 erg/s, respectively, The triangle, circle, and square symbols display cases in which emission up to a t obs = 10, 30, and 100 s are taken into account, respectively.