PHOTOSPHERIC EMISSION FROM COLLAPSAR JETS IN 3D RELATIVISTIC HYDRODYNAMICS

There is mounting evidence from recent observations that photospheric emission plays an important role in the prompt phase of gamma-ray bursts (GRBs). Since photons below the photosphere tend to be thermalized due to coupling with the matter, the most direct indication for the existence of photospheric emission is the detection of a thermal-like (blackbody) component in the spectrum. Although rare, such features are reported in the literature (Ryde 2004; Ghirlanda et al. 2013). The most notable example is GRB 090902B (Abdo et al. 2009) in which the overall spectrum was well-fitted by a multi-color blackbody (Ryde et al. 2010). Moreover, it is also worth noting that a non-negligible fraction of GRBs favor a photospheric origin, in the sense that synchrotron emission models have difficulty reproducing their spectral shapes (Preece et al. 1998; Axelsson & Borgonovo 2015).

From a theoretical point of view, in order to properly evaluate the photospheric emission, radiation transfer within the relativistic jet must be considered (Beloborodov 2011; Pe'er and Ryde 2011). Up to now, such calculations have been performed under the one-zone or steady-state approximations (e.g., Pe'er et al. 2005, 2006; Giannios & Spruit 2007; Ioka et al. 2007; Giannios 2008; Beloborodov 2010; Lazzati & Begelman 2010; Vurm et al. 2011; Asano and Mészáros 2013; Ito et al. 2013, 2014; Lundman et al. 2013, 2014; Chhotray & Lazzati 2015). These studies have shown that dissipation and/or multi-dimensional geometry near the photosphere are important to explain the observed spectrum. However, it is obvious that the assumed background hydrodynamics is oversimplified, and that observed temporal variation cannot be explained.

Several studies have investigated the correspondence between the more realistic jet dynamics and photospheric emission by performing two-dimensional (2D) hydrodynamical simulations (Lazzati et al. 2009, 2013; Mizuta et al. 2011; Nagakura et al. 2011; López-Cámara et al. 2014). These studies evaluated emission from a jet breaking out from a massive stellar envelope (i.e., a collapsar jet). They showed that stellar-jet interactions have a significant impact on the resulting light curve. Note, however, that these works assumed that the photosphere itself is a blackbody; radiation transfer, essential for addressing the spectral shape, was not included.

In the present study, we explore the photospheric emission from a collapsar jet based on three-dimensional (3D) hydrodynamical simulations and a post-process radiation transfer calculation. This is the first time that the effects of 3D non-steady hydrodynamics and radiation transfer have been taken into account simultaneously in the context of photospheric emission.¹ We focus on how the 3D dynamics are imprinted on the light curve and spectrum. In particular, we show that precession activity at the jet base can leave a clear imprint in the light curve. Also, we show that the spectrum can have a non-thermal shape, which may account for typical observations of GRBs.

This Letter is organized as follows. In Section 2, we describe our model and numerical procedures used in our calculations. We present the main results in Section 3. Section 4 is devoted to a discussion and conclusions.

We calculate the evolution and propagation of the jet using a relativistic hydrodynamical (RHD) code (for details of the numerical scheme, see Matsumoto et al. 2012; Matsumoto & Masada 2013). We employ a spherical coordinate system ($r,\theta ,\phi$) whose computational domain spans ${10}^{10}\;{\rm{cm}}\leqslant r\leqslant 1.3\times {10}^{13}\;{\rm{cm}}$ and $\pi /4\leqslant \theta ,\phi \leqslant 3\pi /4.$ We use 280 uniformly spaced grid zones in the two angular dimensions, while, in the radial direction, 1280 non-uniform grid zones are used with ${\rm{\Delta }}r=r{\rm{\Delta }}\theta .$

The initial conditions of the simulation include a progenitor star with mass and radius of $\sim 14{M}_{\odot }$ and $\sim 4\times {10}^{10}\;{\rm{cm}},$ respectively (Model 16TI; Woosley & Heger 2006). At the inner boundary ($r={10}^{10}\;{\rm{cm}}$), we inject a relativistic jet with a half-opening angle of ${\theta }_{{\rm{j}}}=5^\circ$ and a kinetic luminosity of ${L}_{{\rm{j}}}={10}^{50}\;\text{erg s}{}^{-1}.$ The initial Lorentz factor and specific enthalpy of the jet are ${{\rm{\Gamma }}}_{{\rm{i}}}=5$ and ${h}_{{\rm{i}}}=100,$ respectively, corresponding to a terminal Lorentz factor of ${{\rm{\Gamma }}}_{{\rm{i}}}{h}_{{\rm{i}}}=500.$ We consider two models for the jet injection. In our fiducial model (Model I), the jet precesses with an inclination angle of ${\theta }_{{\rm{inc}}}=3^\circ$ and period of ${t}_{{\rm{pre}}}=2\;{\rm{s}}.$ Our reference model (Model II) uses steady injection. In both models, a 1% random perturbation is imposed on the pressure at the injection surface and the duration of the injection is $100\;{\rm{s}}.$ We follow the evolution up to the phase when the jet becomes optically thin.

Then, using the output data of the RHD simulation as a background fluid, we treat the photon transfer using a Monte Carlo method. Initially, we inject photons far below the photosphere in a blackbody distribution at local temperature. In the present study, the injection location is the constant-radius surface where the Thomson optical depth $\tau =100$ along a line of sight (LOS) parallel to the central axis of precession (jet axis) for Model I (II).² Since the opacity is dominated by electron scattering, we only take into account the Thomson optical depth.³ Not all of the injection surface is used for injecting photons: to focus on photons related to the jet, injection occurs only where the bulk Lorentz factor exceeds 1.5. After injection, we follow the evolution of the photons until they reach the outer boundary of the calculation, located above the photosphere ($r\lt {10}^{13}\;{\rm{cm}}$) at all times. At each scattering, we take into account the full Klein–Nishina cross section, the recoil effect and the thermal motion of electrons. We assume that the scattering electrons have a Maxwellian distribution at a local temperature obtained from the simulations, typically on the order of $\sim {\rm{keV}}.$

Throughout the paper, the location of the observer is expressed by the observer angle ${\theta }_{{\rm{obs}}},$ which is the angle between the LOS and central axis of precession (jet axis) for Model I (II).⁴

Figure 1 shows snapshots of the RHD simulations. The overall evolution is similar between the two models. In the initial phase, when the jet is penetrating the star, the bulk Lorentz factor of the jets remains relatively low (${\rm{\Gamma }}\lesssim 10$) due to the formation of shocks via strong interaction with the stellar envelope. After breakout, the jet expands in the tenuous medium and accelerates up to ${\rm{\Gamma }}\sim \;\text{a few}\;\times 100.$

Zoom In Zoom Out Reset image size

There are broad similarities in the overall structure of the two models' jets. At the head of both jets, there is a terminal reverse shock. The shocked matter produces a cocoon that envelopes the entire jet in both models. Due to the pressure in the cocoon, the jets remain collimated by forming recollimation shocks. Upstream of the "first" recollimation shock identifiable as the boundary surface of the central red region in the right panels of Figure 1, the jet is causally disconnected with the external region, and freely accelerates by converting its internal energy to kinetic energy. On the other hand, in the downstream region, a complex structure (multiple shocks, turbulence, and mixing) develops due to the strong interaction with the stellar envelope. Note that, while the outflow dynamics in Model II are predominantly determined by the above-mentioned jet–stellar interaction, the precession in Model I induces an additional periodic disturbance in the entire outflow, including the region upstream from the first recollimation shock. To show a clear picture of the 3D internal structure for Model I, we plot an isosurface map of the Lorentz factor in Figure 2. The small wiggling structures produced by the precession can be seen in the figure.

Zoom In Zoom Out Reset image size

The resulting light curves are displayed in the top panels of Figure 3. In both models, the initial rapid increase in the luminosity is produced by the breakout of the jet from the stellar surface, and the photosphere nearly coincides with the position of the forward shock. As the matter dilutes due to the rapid expansion, the photosphere recedes from the forward shock and reveals the inner regions. The subsequent emission strongly reflects the internal structure of the jet at the point of last scattering. Note that the observer time, ${t}_{{\rm{obs}}},$ (arrival time of photons since the first photon reached the observer) is different from the laboratory time t of the simulation.

Zoom In Zoom Out Reset image size

Focusing on Model II, the temporal variation during the initial $\sim 30\;{\rm{s}}$ of the light curve shows few peaks. This behavior is determined by the structures formed in the jet–stellar interaction, since the photosphere is located downstream of the first recollimation shock at this phase. Later, there is a phase of steady (${\theta }_{{\rm{obs}}}\lesssim 4^\circ$), or quasi-steady (${\theta }_{{\rm{obs}}}\gtrsim 4^\circ$), behavior. Note that the two cases are roughly demarcated at an observer angle of ${\theta }_{{\rm{obs}}}\sim 4^\circ$ since the half-opening angle of the recollimation shock is $\sim 4^\circ$ (see the top right panel of Figure 1). The former case does not show variation since the observed photons mainly originate upstream of the shock, where the outflow is largely steady. In the latter case, most of the observed photons are released near the immediate downstream of the shock. Over the period of $30\;{\rm{s}}\lesssim {t}_{{\rm{obs}}}\lesssim 100\;{\rm{s}},$ off-axis observers see the less perturbed material behind the complex jet head, leading to a gradual increase in luminosity. Finally, we see a sudden drop in the luminosity at ${t}_{{\rm{obs}}}\sim 100\;{\rm{s}},$ which corresponds to the time in which the photosphere reaches the tail of the jet and no further injection is occurring. The e-folding time is $\sim 0.2\;{\rm{s}}.$ This means that the light curve can respond to the central engine activity if the variability timescale is longer than the e-folding time.

In Model I, while the long-term behavior ($\sim 10\;{\rm{s}}$ scale) of the light curve can be explained in the same way as above, clear differences are visible at the $\sim {\rm{s}}$ scale, which shows periodic spiky features. This is because the precession induces a small-scale disturbance in the jet that survives up to the emission region. As a result, except for the initial $\sim 10\;{\rm{s}},$ spikes with a period of roughly ${t}_{{\rm{pre}}}=2\;{\rm{s}}$ are produced. Interestingly, this indicates that the central engine activity is not smeared out during the propagation and makes a clear imprint on the emission. This quick response is possible, since the emission is sensitive to any central activity that has a timescale longer than $\sim 0.2\;{\rm{s}}$ as mentioned above. The periodic features are probably washed out in the earliest phase ($\lesssim 10\;{\rm{s}}$) by the strong stellar–jet interaction. However, we note that poor spatial resolution at large radii can also suppress any periodic signal. We note also that the spikes are less prominent for smaller observer angles. This is because ${\theta }_{{\rm{inc}}}\lt {\theta }_{{\rm{j}}},$ leading to mostly steady injection near the precession axis.

In Figure 3 (bottom panels), we show the time integrated spectra, which are broken into a series of instantaneous spectra in Figure 4. Interestingly, non-thermal features are present, even though only thermal photons are injected and only thermal electrons are considered as a scattering medium. There is a significant broadening from a thermal spectrum below and above the peak energies (Ito et al. 2013, 2014).

Zoom In Zoom Out Reset image size

Focusing on energies below the peak (but above $\sim 10\;{\rm{keV}}$), the spectral indices are roughly in the range of $\nu {L}_{\nu }\propto {\nu }^{1.5}\mbox{--}{\nu }^{0.5}$ (Figure 3). This is much softer than the Rayleigh–Jeans part of the blackbody spectrum and is consistent with the observed low energy spectral slopes. This change is produced mainly by the multi-temperature effect. Since there is substantial diversity in both the temperature and the Doppler factor within the jet, we observe various emission components that have different peak energies and luminosities. The temporal evolution of these quantities can be seen in Figure 4. Initially (${t}_{{\rm{obs}}}\lesssim 30\;{\rm{s}}$), the peak energy shows a large variation (up to an order of magnitude change for a fixed observer angle) reflecting the complex structure in the downstream region of the first recollimation shock. From $30\;{\rm{s}}\lesssim {t}_{{\rm{obs}}}\lesssim 100\;{\rm{s}},$ spectral evolution is less extreme, because the photosphere lies in the mostly unperturbed material around the recollimation shock.

The spectrum above the peak is also harder than thermal, but not because of the multi-temperature effect. Instead, the spectrum is formed by a bulk Compton scattering that occurs at regions with sharp velocity shears (Ito et al. 2013, 2014). In our simulations, the most prominent velocity shear is located at the first recollimation shock.⁵ Hence, the non-thermal component is most pronounced when the observer angle is close to the opening angle of the shock (${\theta }_{{\rm{obs}}}\sim 4^\circ$). Also, considering the temporal evolution, the non-thermal components appear at a later phase (${t}_{{\rm{obs}}}\gtrsim 30\;{\rm{s}}$), when the emission region is near, or coincides with, the first recollimation shock.

The trends discussed in the preceding paragraph are confirmed for Model II in Figures 3–4. In the range of ${\theta }_{{\rm{obs}}}\sim 4^\circ \mbox{--}6^\circ ,$ the spectral indices of the non-thermal component ($\nu {L}_{\nu }\propto {\nu }^{-0.3}\mbox{--}{\nu }^{-0.5}$) can reproduce typical observations (Kaneko et al. 2006; Nava et al. 2011; Gruber et al. 2014). In contrast, the non-thermal component is quite weak in Model I because of precession, which makes the structure of the recollimation shock more diffuse. As a result, the efficiency of the bulk Comptonization is lower than that in Model II, due to the smaller velocity gradient at the shock. Note, however, that this result strongly depends on the spatial resolution of the simulation. Due to the prohibitive computational expense, our current setup does not have sufficient resolution to capture structures with length scales of the photon mean free path. Hence, the efficiency of the bulk Comptonization is artificially reduced by numerical diffusion. Our results should therefore be considered a lower bound for the high energy non-thermal components.

We also note that a self-consistent calculation that takes into account the radiation feedback on the dynamics is required for more precise evaluation, since radiation possesses a non-negligible fraction ($\gtrsim 10\%$) of the energy in the jet. While the essential features found in our results are expected to remain unchanged, such calculation is out of the scope of this Letter.

We have evaluated the photospheric emission from a collapsar jet using 3D RHD simulations and post-process radiation transfer calculations. We emphasize that the main difference from previous studies of photospheric emission is that the 3D hydrodynamics and radiation transfer are handled simultaneously. This has enabled us to obtain deeper insight into the emission properties.

Two models with different properties at the jet base are considered: injection with and without a precession motion. We find that both the jet–stellar interaction and precession have a significant impact on the resulting emission. The former effect leads to the presence of a complex structure accompanied by multiple shocks in the downstream region of the first recollimation shock. In addition, the latter effect results in a large number of disturbances and shocks throughout the jet volume. The resulting light curve reflects the structures. One particularly interesting result is that precession imprints a clear signature, manifesting as spikes with the same period as the precession. This implies that, if the central engine activity induces precession motion, it should be directly observable in the light curve. Although not conclusive, the existence of periodicity is indicated in a few bursts (e.g., GRB 970110A and GRB 080319B; Pozanenko et al. 2004; Crider 2006; Beskin et al. 2010). Our results offer a possible explanation for the origin of such features. We suggest a search for additional evidence of periodic behavior in GRB light curves to investigate the possibility.

The above conclusion echoes a previous study by López-Cámara et al. (2014; see also Morsony et al. 2010), which considered the effect on emission of central engine activity with episodic jet injection. They also found that central engine activity can have a direct impact on the light curve. Compared to their study, we find a much clearer response of the light curve to the central engine of the GRB, on timescales as short as $0.2\;{\rm{s}}.$ This is due to the fact that there is no quiescent state in our calculation that would tend to smear out the central engine activity.

Another main result is that the resulting spectra show non-thermal features, which can account for observations of GRBs. As shown in previous studies (Ito et al. 2013, 2014; Lundman et al. 2013, 2014), the broadening from a thermal spectrum is caused by the global structure of the jet. While the low energy slope is produced by the multi-temperature effect, the high energy non-thermal tail is mainly generated by the bulk Comptonization at the recollimation shock. This raises an interesting possibility that the non-thermal spectrum of GRBs can be explained by the multi-D structure of the jet, which inevitably develops during the propagation. Note that this spectral broadening occurs in the absence of relativistic or even non-thermal scattering particles (only thermal particles with non-relativistic temperature are considered). Therefore, it is essentially different from the widely discussed dissipative photosphere models in which the broadening is caused by the population of relativistic electrons or pairs.

Our results show that the low energy spectra of both models are in good agreement with observations, almost independent of observer angle. On the other hand, the typical high energy spectrum is reproduced in only a few cases (Model II with an observer angle in the range of ${\theta }_{{\rm{obs}}}\;4^\circ \mbox{--}6^\circ :$ see Figure 3). However, we again stress that the spectra above the peak energy are artificially reduced due to the lack of spatial resolution in our hydrodynamical calculations. Since shocks and shear flows are inevitably formed within the jet, we expect that calculations with higher spatial resolution, capable of capturing sharper structures, will result in the appearance of a high energy non-thermal tail, or in the enhancement of one already present in our model. Thus, the weak signature or absence of high energy non-thermal tails in our calculations does not imply that the model cannot account for the observed spectra. Such aspects will be addressed in future work.

We thank A. Pozanenko and A. Mizuta for fruitful discussions. This work is supported by the Japan Society for the Promotion of Science (No. 23340069 and No. 25610056). We acknowledge the financial support of a Grant-in-Aid for Young Scientists (B:26800159). Numerical computations and data analysis were carried out on XC30 and PC cluster at Center for Computational Astrophysics, National Astronomical Observatory of Japan. This work was supported in part by the Center for the Promotion of Integrated Sciences (CPIS) of Sokendai.