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

Hirotaka Ito; Jin Matsumoto; Shigehiro Nagataki; Donald C. Warren; Maxim V. Barkov; Daisuke Yonetoku

doi:10.3847/1538-4357/ace775

1. Introduction

Despite the decades of research conducted 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 nonthermal 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 nonthermal features by considering dissipation (Pe'er et al. 2005; Ioka et al. 2007; Giannios 2008; Beloborodov 2010; Lazzati & Begelman 2010; Vurm et al. 2011; Ahlgren et al. 2015; Chhotray & Lazzati 2015; Vurm & Beloborodov 2016; Bhattacharya & Kumar 2020) and/or geometrical effects (Ito et al. 2013, 2014; Lundman et al. 2013, 2014) in the subphotospheric region.

However, it has been recently 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, 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).¹⁰ It is 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 three GRBs was reported in their observations (Yonetoku et al. 2011b, 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 five 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, and 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 AstroSAT 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., Granot 2003; Lyutikov et al. 2003; Nakar et al. 2003; Waxman 2003; Lazzati 2006; Toma et al. 2009; 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) found that, while the photons released at a local emitting region are strongly polarized, the superposition of each emission component suppresses the net polarization if the outflow is uniform. Subsequent 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 subphotospheric 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 explored the properties of long GRBs (LGRBs) by applying a crudely approximated photospheric model to hydrodynamical simulations of a relativistic jet breaking out of massive stellar envelope (Lazzati et al. 2009, 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, 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, 2021; Parsotan et al. 2018; 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 nonhydrodynamical 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 increased 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 time-resolved spectral and polarization analysis are presented in Sections 5 and 6, respectively. Section 7 is devoted to a summary and discussions.

2. Model and Methods

As described in IMN19, we analyze the photospheric emissions evaluated in three 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. 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⁸ to ∼5 × 10⁹) 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 details). 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.¹¹ We consider three models that employ different jet kinetic powers: L_j = 10⁴⁹, 10⁵⁰, and 10⁵¹ erg s⁻¹. 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¹⁰ cm. Note that the jet expands to an opening angle of $\sim {\theta }_{\mathrm{ini}}+0.7{{\rm{\Gamma }}}_{\mathrm{ini}}^{-1}$ 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⁴⁹ and 10⁵⁰ erg s⁻¹, while h_ini = 180 is adopted in the model with L_j = 10⁵¹ 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⁴⁹, 10⁵⁰, and 10⁵¹ erg s⁻¹, 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 that, as a whole, carries the Stokes parameters: I, the intensity, and Q and U, the two parameters that 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: ${\rm{\Pi }}=\sqrt{{Q}^{2}+{U}^{2}}/I$ and $\psi =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.¹² Regarding 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 self-consistently 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 three dimensions, 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°.

3. 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 $\left\langle {\rm{\Gamma }}\right\rangle$ (blue) evaluated at r = 10¹³ cm. These quantities are computed as functions of θ:

$\begin{eqnarray}&&{E}_{j,\mathrm{iso}}(\theta )={\int }_{{t}_{13}}^{{t}_{13}+{t}_{\max }}{L}_{j,\mathrm{iso}}(\theta )\,{dt},\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&\left\langle {\rm{\Gamma }}(\theta )\right\rangle =\displaystyle \frac{1}{{E}_{j,\mathrm{iso}}(\theta )}{\int }_{{t}_{13}}^{{t}_{13}+{t}_{\max }}{\rm{\Gamma }}(\theta )\,{L}_{j,\mathrm{iso}}(\theta ){dt},\end{eqnarray} \tag{ 2 }$

where ${L}_{j,\mathrm{iso}}=4\pi {r}^{2}{{\rm{\Gamma }}}^{2}[\rho {c}^{2}(h-1/{\rm{\Gamma }})]{\beta }_{j}\cos {\theta }_{\beta }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¹³ cm. Here, θ_β is the angle between the direction of velocity and radial direction. The time when the forward shock reaches r = 10¹³ cm is denoted by t₁₃, 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¹³ cm for a duration of ${t}_{\max }$ . This outflow component corresponds to the region emitting radiation from t_obs = 0 to $\approx {t}_{\max }$ . Note that the radius r = 10¹³ cm is significantly larger than the stellar radius (≈4 × 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 $\left\langle {\rm{\Gamma }}\right\rangle$ . As found in the previous studies (e.g., Gottlieb et al. 2021b), the outflow structure shows a core region at the center, in which E_iso,j and $\left\langle {\rm{\Gamma }}\right\rangle$ 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.

$\left\langle {\rm{\Gamma }}\right\rangle $ — **Figure 1.** Polar distribution of the accumulated isotropic energy E_j,iso (red) and the average Lorentz factor $\left\langle {\rm{\Gamma }}\right\rangle$ (blue) of the outflow evaluated at r = 10¹³ 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 ${t}_{\max }=10$ , 30, and 100 s, respectively. The dashed, solid, and dotted lines for the Lorentz factor $\left\langle {\rm{\Gamma }}\right\rangle$ are the average values weighted with the energy of the outflow over the same duration.
Download figure:
Standard image High-resolution image

As seen in Figure 1, while the opening angle of the core, ${\theta }_{\mathrm{core}}$ , remains nearly constant throughout the entire duration for the model with L_j = 10⁴⁹ erg s⁻¹, it rapidly broadens at t ∼ 30 s and ∼10 s for models with L_j = 10⁵⁰ erg s⁻¹ and 10⁵¹ erg s⁻¹, respectively. The sudden spread in these models reflects the time the innermost recollimation shock¹³ breaks out from the stellar envelope; therefore, confinement by the cocoon pressure becomes less effective at later times. Before the transition, the opening angle of the core is ${\theta }_{\mathrm{core}}\sim 2^\circ \mbox{--}4^\circ$ in all models, which is consistent with the scaling found in previous studies, ${\theta }_{\mathrm{core}}\sim (\displaystyle \frac{1}{5}-\tfrac{1}{3})({\theta }_{\mathrm{ini}}+0.7{{\rm{\Gamma }}}_{\mathrm{ini}}^{-1})$ (Mizuta & Ioka 2013; Gottlieb et al. 2021b). After the breakout of the recollimation shock, the unshocked jet matter that expands following the recollimation shock expels 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.

4. 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}_{\mathrm{obs}}={t}_{\mathrm{obs},\max }=100\,{\rm{s}}$ ), we also show the cases in which only emission up to a specific time ( ${t}_{\mathrm{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}_{\mathrm{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.

**Figure 2.** Peak luminosity, L_p (left), and isotropic equivalent energy, E_iso (right), vs the peak energy, E_p, for different viewing angles (colored symbols). The shape of the symbols represents the difference in the emission duration considered to extract L_p, E_iso, and E_p. Observational data of LGRBs taken from Yonetoku et al. (2010; gray points) are added for comparison. The best-fitting function and 3σ systematic error regions of the E_p–L_p and E_p–E_iso correlations are indicated with gray solid and dashed lines, respectively.
Download figure:
Standard image High-resolution image

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 ( ${\theta }_{\mathrm{obs}}\lesssim {\theta }_{\mathrm{core}}$ ), E_p, L_p, and E_iso do not show significant variation. On the other hand, beyond the core ( ${\theta }_{\mathrm{obs}}\gtrsim {\theta }_{\mathrm{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 2 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 (Γ₀) 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, Γ₀, reported in the literature (e.g., Liang et al. 2010, 2013, 2015; Ghirlanda et al. 2012; Lü et al. 2012). Figure 3 displays E_p–Γ₀ and L_p–Γ₀ diagrams of the simulation results together with those taken from the recent observational study by Ghirlanda et al. (2018). In our simulations, Γ₀ is defined as an energy-weighted average value of the terminal Lorentz factor hΓ along the LOS. It is computed in the same manner as Equation (2) and can be expressed as follows:

$\begin{eqnarray}&&{{\rm{\Gamma }}}_{0}(\theta )=\displaystyle \frac{1}{{E}_{j,\mathrm{iso}}(\theta )}{\int }_{{t}_{13}}^{{t}_{13}+{t}_{\max }}h{\rm{\Gamma }}(\theta )\,{L}_{j,\mathrm{iso}}(\theta )\,{dt}.\end{eqnarray} \tag{ 3 }$

When depicting the E_p (L_p)–Γ₀ 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–Γ₀ correlation.¹⁴ While there is a small offset from the center of the observed distribution, the E_p–Γ₀ 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 Γ₀ 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 can be seen in the figure, different assumptions about the circum-burst medium (CBM), i.e., homogeneous (gray circles) or wind-like (gray stars), result in a systematic shift in the distribution.

Another caveat is on the estimation of Γ₀ based on the simulation. Here we have evaluated the quantity as an energy-weighted average value of the terminal Lorentz factor, hΓ, along the LOS. Regarding the core region of the jet, the value is likely to result in an overestimation (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.¹⁵ On 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 Γ₀ sharply drops with angle, it is not clear whether the onset time measured by the observer viewing the wing region ( ${\theta }_{\mathrm{obs}}\gt {\theta }_{\mathrm{core}}$ ) reflects the deceleration time of the jet matter along the LOS. 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–Γ₀ correlation and that evaluated based on the simulation. The current study, which employs a more accurate calculation method, confirms their conclusion that photospheric emission naturally gives rise to the correlation between Γ₀ and the prompt properties.

5. Time-resolved Spectral Analysis

In this section, we present the result of the time-resolved analysis. The analysis is similar to those performed by Parsotan & Lazzati (2018) and Parsotan et al. (2018). We employ time bins of Δt = 2 s, and all spectra are fitted either by the Band function (Band et al. 1993) or a cutoff power law (CPL) function (Yu et al. 2016) using the χ² 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 the 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 cutoff. 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⁵⁰ erg s⁻¹).

**Figure 4.** Light curves and evolutions of fitted spectral parameters at various viewing angles for the fiducial model (L_j = 10⁵⁰ 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 E_p, 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.
Download figure:
Standard image High-resolution image

5.1. 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 that 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 of the emissions are from the shocked region. Nonetheless, 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⁴⁹ 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⁵¹ 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 break out 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 spatial 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 good agreement with the observations (Golenetskii et al. 1983; Ghirlanda et al. 2010; Lu et al. 2012).¹⁶ In Figure 5, we plot the time-resolved E_p and L of the three simulations. It is indicated in Ghirlanda et al. (2010) and 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. 2004, 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⁵⁰ 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=K{{E}_{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⁵⁰ 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.

**Figure 5.** The relation between spectral peak energy E_p 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 powers of L_j = 10⁴⁹, 10⁵⁰, and 10⁵¹ erg s⁻¹, respectively.
Download figure:
Standard image High-resolution image

5.2. 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). 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 angle. Therefore, most of the time bins are fitted by the CPL function for the former case, while the 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⁵⁰ 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 spectral broadening via the multicolor effect and bulk Comptonization. These effects combine to soften the low-energy spectrum and harden the high-energy spectrum (Ito et al. 2013; Lundman 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 ∼1 keV up to ∼10³ 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⁵¹ 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 nondissipative photospheric emission from a spherical outflow (Beloborodov 2010; Bégué et al. 2013; Ito et al. 2013; Lundman et al. 2013). This implies that the multicolor 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 (<100 keV). Therefore, while this population may be compatible with X-ray-rich (XRR) GRBs or X-ray flashes (XRFs; Katsukura et al. 2020), it does not match with the typical GRB observations.

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 multicolor 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 two dimensions (Lazzati 2016; Parsotan & Lazzati 2018, 2021; Ito et al. 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. (2019, 2021a), 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 & Ryde 2022; Samuelsson et al. 2022). Moreover, additional internal dissipation that results in a significant nonthermal 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.

6. 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 (Ito et al. 2014, 2021; Lundman et al. 2014; Parsotan et al. 2020; Parsotan & Lazzati 2022), restricts the position angle ψ = 1/2 arct(U/Q) to be one of two orthogonal 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 three dimensions is essential to remove these artificial effects on the results of polarization analysis.

6.1. Time-integrated Polarization

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}_{\mathrm{obs},\max }=10$ , 30, or 100 s. In the figure, the 1σ statistical uncertainty of the degree of polarization ( ${\rm{\Pi }}=\sqrt{{Q}^{2}+{U}^{2}}/I$ ) and position angle ( $\psi =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 ${\sigma }_{I/Q/U}^{2}$ and σ_IQ/IU/QU, to quantify the variance of Π and ψ as

$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{{\rm{\Pi }}}^{2} & = & {\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial I}\right)}^{2}{\sigma }_{I}^{2}+{\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial Q}\right)}^{2}{\sigma }_{Q}^{2}+{\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial U}\right)}^{2}{\sigma }_{U}^{2}\\ & & +\,2\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial I}\right)\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial Q}\right){\sigma }_{{IQ}}+2\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial I}\right)\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial U}\right){\sigma }_{{IU}}\\ & & +\,2\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial Q}\right)\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial U}\right){\sigma }_{{QU}},\end{array}\end{eqnarray} \tag{ 4 }$

and

$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{\psi }^{2} & = & {\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial I}\right)}^{2}{\sigma }_{I}^{2}+{\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial Q}\right)}^{2}{\sigma }_{Q}^{2}+{\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial U}\right)}^{2}{\sigma }_{U}^{2}\\ & & +\,2\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial I}\right)\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial Q}\right){\sigma }_{{IQ}}+2\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial I}\right)\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial U}\right){\sigma }_{{IU}}\\ & & +\,2\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial Q}\right)\left(\displaystyle \frac{\partial {\rm{\Pi }}}{\partial U}\right){\sigma }_{{QU}},\end{array}\end{eqnarray} \tag{ 5 }$

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 a rapid decrease marks the transition from the core to the wing region as the origin of the emission.

**Figure 8.** Degree of linear polarization ( ${\rm{\Pi }}=\sqrt{{Q}^{2}+{U}^{2}}/I$ : upper panels) and position angle ( $\psi =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 E_p. The left, middle, and right columns correspond to the simulation with jet powers of L_j = 10⁴⁹, 10⁵⁰, and 10⁵¹ erg s⁻¹, respectively. The top, middle, and the bottom rows show 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 that show 1σ error bar larger than 40°. Note that the scale of the vertical axis of Π differs among the panels.
Download figure:
Standard image High-resolution image

**Figure 8.** Degree of linear polarization ( ${\rm{\Pi }}=\sqrt{{Q}^{2}+{U}^{2}}/I$ : upper panels) and position angle ( $\psi =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 E_p. The left, middle, and right columns correspond to the simulation with jet powers of L_j = 10⁴⁹, 10⁵⁰, and 10⁵¹ erg s⁻¹, respectively. The top, middle, and the bottom rows show 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 that show 1σ error bar larger than 40°. Note that the scale of the vertical axis of Π differs among the panels.
Download figure:
Standard image High-resolution image

As can be seen in the figure, when the observer LOS is within the core region ( ${\theta }_{\mathrm{obs}}\lesssim {\theta }_{\mathrm{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 ( ${\theta }_{\mathrm{obs}}\gtrsim {\theta }_{\mathrm{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 simulations, in which the former show more complex flow structures. It is also worth noting that in the current 3D simulations, Π exhibits a nonzero 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⁵⁰ erg s⁻¹ with ${t}_{\mathrm{obs},\max }=30\,{\rm{s}}$ as an illustrative case. The top panel of the figure shows Π as a function of θ_obs, which is identical to the right half 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 found in all cases explored in the current study. Hence, as mentioned earlier, while the viewing angle is within the core, a low level of polarization is found, and Π tends to increase with θ_obs at a larger viewing angle. A schematic image for the anisotropy of emission around the LOS is displayed in Figure 10. The image illustrates how the distribution of the intensity changes with the viewing angle.

**Figure 10.** 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 ( ${\theta }_{\mathrm{obs}}\lt {\theta }_{\mathrm{core}}$ : blue) and outside the core ( ${\theta }_{\mathrm{obs}}\gt {\theta }_{\mathrm{core}}$ : green). The intensity around the given LOS is expressed by a color contour in red with darker color corresponding to a larger intensity. As shown in the picture, intensity is close to uniform at ${\theta }_{\mathrm{obs}}\lt {\theta }_{\mathrm{core}}$ , while a concentration toward the direction of the jet axis is seen at ${\theta }_{\mathrm{obs}}\gt {\theta }_{\mathrm{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.
Download figure:
Standard image High-resolution image

**Figure 10.** 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 ( ${\theta }_{\mathrm{obs}}\lt {\theta }_{\mathrm{core}}$ : blue) and outside the core ( ${\theta }_{\mathrm{obs}}\gt {\theta }_{\mathrm{core}}$ : green). The intensity around the given LOS is expressed by a color contour in red with darker color corresponding to a larger intensity. As shown in the picture, intensity is close to uniform at ${\theta }_{\mathrm{obs}}\lt {\theta }_{\mathrm{core}}$ , while a concentration toward the direction of the jet axis is seen at ${\theta }_{\mathrm{obs}}\gt {\theta }_{\mathrm{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.
Download figure:
Standard image High-resolution image

We find no general trend regarding the dependence of Π on the energy band. When the viewing angle is small ( ${\theta }_{\mathrm{obs}}\lesssim {\theta }_{\mathrm{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 percent 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 nonmonotonic 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}_{\mathrm{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 longer ${t}_{\mathrm{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⁵⁰ erg s⁻¹ and t_obs ≲ 10 s for L_j = 10⁵¹ 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 $\psi =1/2\,\arctan (U/Q)$ in Figure 8, a clustering around ψ = 0° (180°) and ψ = 90° can be 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⁴⁹ 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).

6.2. Time-resolved Polarization

The time-resolved analysis reveals complex temporal behavior in the properties of polarization. As in the case of the light curves and spectra, these properties directly 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⁴⁹ and 10⁵⁰ erg s⁻¹. The corresponding light curves are also displayed in the figure for comparison. The figure confirms the presence of nonnegligible time variability in Π and ψ. It also shows that there is no clear correlation between the temporal behavior of the polarization with the light curves.

**Figure 11.** Degree of linear polarization ( ${\rm{\Pi }}=\sqrt{{Q}^{2}+{U}^{2}}/I$ : upper panels) and position angle ( $\psi =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 a jet power of L_j = 10⁴⁹ erg s⁻¹ (10⁵⁰ 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°.
Download figure:
Standard image High-resolution image

Similar to the results of the time-integrated analysis, the time-resolved values of Π tend to increase with θ_obs. The nonmonotonic 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, that 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.

The temporal evolution of ψ 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⁵⁰ 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⁴⁹ 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 nonnegligible energy dependence of Π and ψ is an inherent feature of the photospheric emission.

6.3. 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 percent). This result is consistent with the recent reports by the POLAR collaboration (Kole et al. 2020) and AstroSAT collaboration (Chattopadhyay et al. 2022), which found that the majority of GRBs they 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 time-resolved analysis of polarizations in a few bursts (GRB 110826A, 170114A, and 160821A; Yonetoku et al. 2011b; Burgess et al. 2019; Sharma et al. 2019). However, these analyses suggest that 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 time-integrated 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.

Over 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, although 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.

**Figure 12.** Degree of linear polarization ( ${\rm{\Pi }}=\sqrt{{Q}^{2}+{U}^{2}}/I$ ) as a function of spectral peak energy E_p (left) and peak luminosity L_p (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 the symbols indicates the difference in the injected jet power: red, blue, and green correspond to L_j = 10⁴⁹, 10⁵⁰, and 10⁵¹ erg s⁻¹, respectively, The triangle, circle, and square symbols display cases in which emission up to t_obs = 10, 30, and 100 s is taken into account, respectively.
Download figure:
Standard image High-resolution image

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 nonnegligible 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 analysis. 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 (see Gill et al. 2021 for details), in order to investigate the energy dependence of polarization.

7. 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 time-integrated spectral properties in IMN19, the improved photon statistics here have enabled us to further explore the time-resolved spectral properties and polarization properties.

We have considered three models in which the relativistic jet is injected in the interior of the massive stellar envelope with constant power (L_j = 10⁴⁹, 10⁵⁰, and 10⁵¹ 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. A 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 ( $\theta \lesssim {\theta }_{\mathrm{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 ${\theta }_{\mathrm{obs}}\gtrsim {\theta }_{\mathrm{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 Γ₀. 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 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 nonthermal 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 three dimensions, 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 time-resolved analysis finds a flip of Δψ ∼ 90° 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 current models is much lower than that 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°. 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. 2019, 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, this tension may be attributed to the possible subphotospheric 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 (nonzero Π 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.¹⁷ Since the simulations cover a large spatial domain, small-scale structures are washed out due to the limited spatial resolution. Consequently, the expected short timescale variability in the light curve, which arises due to the jet–stellar 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¹⁸ 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. (2020b, 2020a, 2021a) 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.

Acknowledgments

We thank B. Zhang, T. Parsotan, and P. Beniamini for fruitful discussions. This work was supported by JSPS KAKENHI grant Nos. JP23H01178, JP20K14473, JP19K03878, and JP19H00693 and 23-12-00220 of the Russian Science Foundation. Numerical computations and data analysis were carried out on Hokusai BigWaterfall system at RIKEN, XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan, and the Yukawa Institute Computer Facility. This work was supported in part by an RIKEN Interdisciplinary Theoretical & Mathematical Science Program (iTHEMS) and an RIKEN pioneering project "Evolution of Matter in the Universe (r-EMU)" and "Extreme precisions to Explore fundamental physics with Exotic particles (E3-Project)."

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

Article metrics

Author e-mails

Author affiliations

ORCID iDs