Black Hole to Photosphere: 3D GRMHD Simulations of Collapsars Reveal Wobbling and Hybrid Composition Jets

We start our simulation after the collapse of the inner stellar core to a BH and simulate the subsequent collapse of the stellar envelope onto the central BH. A few milliseconds after the onset of the simulation, an accretion disk forms and powerful bipolar relativistic jets are launched. Our choice of the initial stellar progenitor magnetic field and density profiles ensures that the jet power is sufficient to overcome the ram pressure of the infalling gas and sustain a steady outflow from the BH vicinity. As the jet propagates through the stellar envelope, it forms a double-shocked layer cocoon (e.g., Figure 4): The high-density outer part of the cocoon consists of the stellar material is heated up by the forward shock, while the low-density inner cocoon consists of the jet material heated up by the reverse shock (Bromberg et al. 2011b). The pressure of the cocoon confines and collimates the jets. This collimation is effective at large radii; at smaller radii, the jets are confined by the accretion disk winds.

The cocoon regulates the jet power by suppressing the rate at which the stellar envelope feeds the accretion disk. This reduces the accretion rate from the expected scaling of $\dot{M}\propto {t}^{1-2\alpha /3}$ (in the case of a freefalling stellar envelope with a power-law density profile index α) to $\dot{M}\sim {t}^{2(1-\alpha )/3}$ (G22). Figure 1(a) shows that the time evolution of the mass accretion rate at r = 5r_g decays as expected (for our choice of α = 1.5), $\dot{M}\propto {t}^{-1/3}$ during the first few seconds. The jet breaks out from the star after ∼3 s, and when the jet head reaches ∼ 2 R_⋆ at ∼5 s, the jet no longer energizes the cocoon. From this point on, the lateral velocity of the cocoon decreases, and by the time the cocoon shocks the entire star, its velocity has dropped by a factor of the jet-opening angle θ_j (Eisenberg et al. 2022). The decay in the cocoon expansion velocity moderates its effect on the mass accretion rate, which nearly plateaus after ∼5 s.

Zoom In Zoom Out Reset image size

Figures 1(b) and (c) depict the jet luminosity L_j at r = r_g , and the jet power efficiency $\eta ={L}_{j}/\dot{M}{c}^{2}$ on the horizon (as defined in G22), respectively. The jet maintains an order unity energy efficiency of η = 0.7 at all times, from the moment it is launched until the end of the simulation, implying that it is accompanied by a magnetically arrested disk (MAD; Tchekhovskoy et al. 2011). As a result, the jet luminosity on the horizon follows the accretion rate as ${L}_{j}\sim \dot{M}{c}^{2}$, as shown in Figure 1(b). The jet power exhibits intermittency on short timescales (10 ms ≲ Δt ≲ 100 ms), comparable to those observed in GRB lightcurves, with amplitude variations of about half an order of magnitude.

As the jet–cocoon structure expands, it carries with it an increasing amount of stellar material, which will later escape from the star once the cocoon breaks out. Toward the end of the simulation, ∼90% of the stellar mass is affected by the cocoon and becomes hydrodynamically unbound (Figure 1(d)). This suggests that in collapsars, the BH mass at the time of the disk formation M_BH,0 is similar to the final BH mass M_BH, as can be seen in Figure 1(e). We note that in stars with a steeper density profile (α ≳ 2), the accretion rate is highest at early times, before the jet unbinds most of the star. It is therefore possible that in such stars the BH accretes a comparable amount of gas to its initial mass (G22), thereby significantly reducing its dimensionless spin and jet-launching efficiency (Tchekhovskoy et al. 2012). We emphasize that our numerical simulation includes neither self-gravity nor the internal energy of the stellar envelope, and therefore, the amount of bound and unbound mass could be different when the simulation properly accounts for these effects. However, because M_⋆ ∼ 3M_BH,0, the freefall time and escape velocity would only change by up to factors of, respectively, ≲1.5 and 2 in the outermost layers had self-gravity been included. Moreover, in the outermost layers, where most of the stellar mass lies and self-gravity is important, the shocked gas is accelerated to the highest velocities by the relativistic outflow, and thus, the dynamics in the pre-shocked gas is likely to have a weak effect.

G22 showed that when the bound parts of cocoon material fall toward the BH, it hits the polar outflow, gets deflected, and hits the accretion disk, kicking it out of alignment with the BH. Consequently, the accretion disk tilt changes. If these stochastic changes add up constructively, the disk can develop substantial tilt values of up to ∼60° (see the movie at http://www.oregottlieb.com/collapsar.html). Because the jet is launched perpendicularly to the disk (Liska et al. 2018), changes in disk tilt also alter the jet orientation.

The top panel of Figure 2 depicts a 3D rendering of the disk (orange) and jets (red) embedded in the cocoon (gray blue). The jets fly out along the disk rotational axis, which is horizontal in the figure and tilted by ∼45° with respect to the stellar rotation axis, which points from the top-left to the bottom-right corner of the panel. One might expect that a jet that is launched in different directions will be soon choked as it needs to drill a new path for each new direction. However, our simulation shows that as soon as the tilted jets are launched, they are deflected toward the low-density regions drilled by the earlier nontilted jets along the original angular momentum axis of the disk. This behavior is seen in the figure as the tilted jets gradually curve toward an angle of ∼45° as soon as they run into the dense edges of the outer cocoon (gray and blue). While the jets' deviation from the axis is moderated with time, they still feature a ∼0.2 rad tilt with respect to the polar axis upon breakout (Section 4.2). The idea of jittering jets in a star was first proposed by Papish & Soker (2011) as an efficient way for the jets to explode the entire star. We find that although our simulation features jittering jets, they are soon deflected back toward the low-density regions so that the jets successfully break out from the star. Had the jets maintained the altered orientation with which they are launched, they would have failed to pierce through the star and would have been choked. In Section 4.2 and Section 5 we study the consequences of disk tilt on the postbreakout stage and its observational implications, respectively.

Zoom In Zoom Out Reset image size

3D RMHD simulations of highly magnetized jets in a dense medium discovered that they develop current-driven instabilities, primarily the magnetic kink instability, once they run into the dense material and recollimate, forming a magnetic nozzle. This is accompanied by the dissipation of their magnetic energy into heat (Bromberg & Tchekhovskoy 2016). G22 recently found that a similar process also takes place in GRMHD simulations of collapsars: The bipolar jets dissipate most of their magnetic energy and become mildly magnetized above the nozzle. However, in G22, the resolution at the collimation nozzle was marginal, with only ≲6 cells across the jet half-opening angle near the nozzle. Here, our simulation is structured such that at any radius and time, the jet half-opening angle is covered by at least 96 cells, an improvement of more than an order of magnitude compared to G22. Interestingly, as we show next, the magnetic dissipation remains, albeit it is continuously distributed throughout the jet rather than concentrated at the nozzle.

The bottom panel of Figure 2 shows a drop in the jet magnetization from σ₀ > 10 near the BH to σ ∼ 1 within the inner 0.1 R_⋆. Figure 3(a) shows the radial profile of the angle-average magnetization, 〈σ〉. To compute the angle average of a quantity, 〈...〉, we weigh its value by the radial energy flux, excluding the rest-mass energy flux and including only the relativistic outflow, i.e., the matter with asymptotic proper velocity ${u}_{\infty }\equiv {\left({{\rm{\Gamma }}}_{\infty }^{2}-1\right)}^{1/2}\gt 1$, where Γ_∞ is the asymptotic Lorentz factor. Γ_∞ is the maximum Lorentz factor that would be attained if all of the gas thermal and magnetic energy were converted into the kinetic energy: ${{\rm{\Gamma }}}_{\infty }\equiv -{u}_{t}\left(h+\sigma \right)$, where u_t is the covariant time component of the four-velocity, and $h=\tfrac{4p}{\rho {c}^{2}}+1$, where p and ρ are the comoving gas pressure and mass density, respectively. We show the magnetization profile at three times, taken after the jets break out of the star. The magnetization profile in the inner 10⁹ cm is rather smooth, with a similar shape at all times, whereas it shows a more complicated behavior at larger radii. At late times, the jet magnetization outside of the star saturates at values 10⁻² ≲ 〈σ〉 ≲ 10⁻¹.

Zoom In Zoom Out Reset image size

Figure 3(b) depicts the angle-average profiles of σ, u, and u_∞ in the jet with the initial magnetization of σ₀ = 15 (thick lines) at t = 5 s. At r ≲ 10⁹ cm, the proper velocity continuously increases while u_∞ remains constant. This indicates that the magnetic energy is efficiently converted to kinetic energy of the bulk and accelerates the plasma. At r ≳ 10⁹ cm, u_∞ drops and shows a more erratic profile. We attribute this behavior to the entrainment of stellar material from the cocoon into the jet. The mixing of the light jet material with heavy stellar material reduces both 〈σ〉 and 〈u_∞〉. A comparison with Figure 3(a) suggests that the mixing becomes more significant at late times and is likely the main cause of the low 〈σ〉 values outside of the star. Outside the star, the profile of 〈σ〉 flattens, implying that the mixing stops after the ejecta breaks out of the stellar surface. The sharp drop at the farthest radii signifies the jet head.

To examine the dependence of the simulation outcome on the initial magnetization, we show the results of an identical simulation with σ₀ = 200 as thin lines in Figure 3(b). A comparison to the simulation with σ₀ = 15 shows that at the acceleration zone (r < 10⁹ cm), the behavior is similar, but the higher value of σ₀ enables acceleration to higher velocities. In the mixing zone, r ≳ 10⁹ cm, the mixing appears to have a stronger effect on the highly magnetized jet, as both u_∞ and σ drop faster with increasing radius. This result could be influenced by the higher susceptibility to the artificial density that is added by the simulation when the jet plasma hits the minimum density floor value. Nevertheless, we see that the results are qualitatively similar to the lower-magnetization case. Outside of the star, the magnetization and terminal proper velocity saturate at values of 〈σ〉 ∼ 10⁻¹ and 〈u_∞〉 ∼ 10 with peak values of ∼0.3 and ∼30, respectively. This indicates that jets with higher σ₀ can accelerate to higher u_∞ and σ. The decrease of 〈u_∞〉 to ∼0.1σ₀ in both models hints that in order for the jet to reach u_∞ ≳ 100, as expected in GRBs, the initial magnetization should be σ₀ ≳ 10³.

Figure 3(c) depicts the maximum values of σ, u, and u_∞ at each radius at the end of the simulation with σ₀ = 15. This shows that parts in the jet maintain u_∞ ∼ 10 and σ ≳ 0.1 after the breakout. Importantly, ∼0.5% and 20% of the outflow energy is carried by plasma with σ > 0.1 when σ₀ = 15 and σ₀ = 200, respectively (see also Figure 5(c)). Therefore, if σ₀ ≳ 10³ as indicated above, then we anticipate 〈σ〉 ∼ 1 at the emission zone. We conclude that the jet has a hybrid composition at the photosphere, which includes magnetic and thermal components, implying that both play a role in the prompt γ-ray signal (see Section 5.3).

The propagation of the jet in the dense stellar envelope forms a hot cocoon that engulfs the jet. The top panel of Figure 4 shows a 3D rendering of the jet–cocoon structure through the logarithm of the asymptotic proper velocity. While the jet maintains relativistic motion, the shocked jet material in the cocoon moves at mildly relativistic speeds. Some of the shocked stellar material in the outer cocoon is seen in subrelativistic velocities in dark blue. The middle panel depicts the logarithm of the magnetization, which decreases to σ ≲ 0.1 in the jet, about halfway through its journey inside the star, while the cocoon maintains a weak magnetization of σ ∼ 10⁻³. The bottom panel shows the final state of the simulation when the forward shock reaches r ∼ 12 R_⋆ t ∼ 18 s after the beginning of the simulation. Shown in green blue, the cocoon breaks out of the star (yellow) with a tilt angle of ∼15° with respect to the stellar angular momentum axis. The bipolar jets (red) are seen as wobbling blobs. The wobbling originates in the disk tilt, and the blobs emerge due to the intermittent nature of the central engine which increases the mixing (Gottlieb et al. 2020b). The blobs are moving at different velocities, and those that are moving in the same direction may produce internal shocks.

Zoom In Zoom Out Reset image size

The first 3D numerical characterizations of the postbreakout structure of GRB jets were calculated by Gottlieb et al. (2020a, 2021b) for weakly magnetized and hydrodynamic GRB jets, respectively. Remarkably, they found that the angular profile of the isotropic equivalent energy can be modeled by a simple distribution, which consists of a flat core, followed by a power law in the jet–cocoon interface (JCI) and exponential decay in the cocoon. The main difference between hydrodynamic and weakly magnetized jets is the power-law index, where weakly magnetized jets feature a steeper power-law drop in the JCI, owing to suppressed mixing between jet and cocoon material due to the presence of magnetic fields.

These works also considered the energy distribution per logarithmic scale of the asymptotic proper velocity, ${dE}/d\mathrm{log}{u}_{\infty }$, which at the homologous phase reflects the radial distribution of the outflow. They found that the relativistic part of the outflow has more energy if the jet is weakly magnetized, owing to the ability of magnetic fields to stabilize the jet against local hydrodynamic instabilities. The energy distribution in the cocoon (10⁻² ≲ u_∞ ≲ 3) was found to maintain a roughly uniform energy distribution in the proper-velocity space and is independent of the jet magnetization. This result has been recently shown to have important observational implications for collapsars: (i) Barnes et al. (2018) and Shankar et al. (2021) suggested that collapsar jets could be the powering mechanism for Type Ic supernovae (SNe Ic). However, because observations of those SNe indicate that the mildly relativistic component carries orders of magnitudes less energy than the subrelativistic ejecta, the flat distribution found in collapsar simulations cannot explain SN lightcurves, thereby ruling out jets as the sole source of SNe (Eisenberg et al. 2022). (ii) Gottlieb et al. (2022c) considered cocoon cooling emission as a possible source of the optical signal in fast blue optical transients (FBOTs). Because the radial energy distribution translates to the time evolution of the lightcurve, a quasi-uniform distribution suggests that the cooling emission in collapsars decays as L∝t⁻². Indeed, similar values are consistently found in all FBOTs with sufficient sensitivity for measuring the optical lightcurve to support this claim (e.g., Margutti et al. 2019; Ho et al. 2020). Nevertheless, the nature of FBOTs remains an open question, and additional models have been shown to explain the nature of FBOTs using different assumptions (e.g., Margutti et al. 2019; Metzger 2022).

Figure 5 complements the aforementioned works by presenting the first energy distributions at r > R_⋆ of highly magnetized jets launched by the rotation of the central BH. The evolution of the angular distribution of the outflow isotropic equivalent energy in Figure 5(a) indicates that the structure has reached the asymptotic stage in the simulation with σ₀ = 15, such that its shape no longer changes. The overall shape of the structure resembles the structure found for hydrodynamic and weakly magnetized jets (Gottlieb et al. 2020a, 2021b). That is, the jet has a flat core with a characteristic angle of θ_j ≈ 6° at which E_iso drops to half of its value on the polar axis. The cocoon is characterized by an exponential decay, and the JCI maintains a rather moderate power law ≲2, similar to what was found in hydrodynamic jets. This result may come as a surprise because weakly magnetized jets feature a steeper power law than hydrodynamic jets, owing to their stabilization effects. The reason for the similarity between our magnetized jets and hydrodynamic ones may lie in the absence of a self-consistent jet-launching mechanism in previous simulations. The controlled injection of a jet via an artificial boundary condition in other works avoids the stochastic behavior seen in the bottom panel of our Figure 4, which is affected by jet tilt, magnetic dissipation, intermittency of jet launching, etc. These phenomena increase the mixing between the jet and the ambient gas and result in a flatter structure that is more similar to the less stable hydrodynamic jets.

Zoom In Zoom Out Reset image size

Figure 5(b) shows the (radial) distribution of the energy in $\mathrm{log}$ of u_∞ for matter outside of the star. The distribution is in agreement with previous results of a flat energy distribution ⁷ in cocoons of jets with different magnetizations, implying that this is a universal outflow distribution that is insensitive to the underlying physics. Therefore, our results support the conclusions of Eisenberg et al. (2022) by showing that Poynting-flux-driven GRB jets cannot solely account for SNe Ic. It also provides a robust prediction that the expected lightcurve of cocoon cooling emission is similar to that observed in FBOTs (Gottlieb et al. 2022c), irrespective of the jet magnetization. The choice of σ₀ = 15 limits the outflow velocities to u_∞ ≲ 15, inconsistent with those inferred from GRBs. For comparison, we show the energy distribution of jet material with σ₀ = 200 outside the star, when the jet head is at r ≈ 10¹¹ cm (blue line). The distribution remains flat, but the cutoff in the velocity is at u_∞ ∼ 40. We note that the flat energy distribution of the cocoon is only weakly affected by the choice of σ₀ such that it is robust. We conclude that both angular and $\mathrm{log}({u}_{\infty })$ energy distributions are consistent with those of hydrodynamic jets, even to a better extent than with those of weakly magnetized jets.

Finally, Figure 5(c) displays the cumulative energy distribution of the matter that broke out of the star and has magnetization larger than σ, normalized by the total energy. Only a negligible amount of energy is carried by matter with σ ≳ 0.1 at t = 18 s in the jet with σ₀ = 15. However, in the model with σ₀ = 200, about ∼20% of the energy is carried by matter with σ ≳ 0.1. Such a significant amount of energy should have an important contribution to the prompt γ -ray signal (the emission is discussed in Section 5.3).

In Section 3.1 we showed that the jet is launched with a considerable tilt with respect to the stellar rotation axis, but once it encounters the dense shoulders of the enveloping cocoon, it is deflected toward the low-density region along the axis that was paved by the previous jet elements. Therefore, even though the jet's initial orientation can be far from the $\hat{z}$-axis, its final deviation from the axis is moderated by its interaction with the cocoon. To quantify the jet-tilt angle evolution, we compute as a function of radius the jet deviation from the polar axis, θ_t (r), defined as the angle at which the maximum u_∞(r) is measured, which we identify as the location of the jet axis.

The top panel of Figure 6 shows the jet-tilt angle θ_t at r = 0.1 R_⋆ (blue) and at r = R_⋆ (purple). At t ≳ 3 s, the tilt angle at 0.1 R_⋆ deviates significantly from the polar axis with θ_t ≲ 0.5 rad. The jet elements at r = 0.1 R_⋆ reach the stellar surface after ∼2 s, which corresponds to the shift in time between the blue and red peaks in θ_t . The collimation process that the jet orientation undergoes by the cocoon moderates the jet deviation from the stellar rotation axis such that the tilt angle amplitudes drop to θ_t ∼ 0.2 rad upon breaking out from the star.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 6 shows a spacetime diagram of θ_t at each time and radius. The map confirms that θ_t may decrease considerably at r ≲ R_⋆, but outside the star, the interaction between the jet and the dense cocoon is minimal, so the jet maintains a roughly constant tilt angle along light-like streamlines (diagonal paths on the map).

GRB jets are typically considered to have an opening angle θ_j around a fixed axis, with θ_j that is significantly larger than the inverse of the jet Lorentz factor, Γ ≳ 100. After the jet breakout from the star, the blast wave interacts with the interstellar medium, its Lorentz factor decreases and eventually reaches ${\rm{\Gamma }}={\theta }_{j}^{-1}$. At this point emission from the entire jet front reaches the observer and further jet deceleration leads to a jet break in the lightcurve. Adopting this picture, numerous studies (e.g., Frail et al. 2001; Bloom et al. 2003; Guetta et al. 2005; Racusin et al. 2009; Fong et al. 2012) examined the jet-opening angle and obtained similar constraints of 0.07 rad ≲ θ_j ≲ 0.16 rad, based on which the local GRB event rate of ∼ 1 Gpc⁻³yr⁻¹ (e.g., Wanderman & Piran 2010) translates to a total rate of ∼ 100 Gpc⁻³yr⁻¹. The fraction of GRBs in SNe Ib/c also depends on the jet-opening angle. It follows from the aforementioned estimates that GRBs likely exist in ∼0.5%–3% of SNe Ib/c. This result is also consistent with radio surveys, which provide an upper limit of ∼10% (Soderberg et al. 2006).

While the above estimates of θ_j are consistent with our simulations that show θ_j ≈ 0.1 rad, the aforementioned inferred GRB rates assume that the jet orientation does not change during the jet's lifetime. If jets in nature oscillate with θ_t ∼ 0.2, as found in our simulation, then their effective opening angle for detection is ∼ θ_t + θ_j (Figure 7), making them ∼10 times more frequent on the sky than previously estimated. This implies that we are able to observe ∼10% of all GRBs in the universe, and the total GRB rate drops by an order of magnitude to ∼ 10 Gpc⁻³yr⁻¹, which is only ∼0.1% of all SNe Ib/c. One possible explanation for this scarcity is that most jets never manage to break out from the star to generate the GRB signal. This possibility was suggested before to explain the high rate of low-luminosity GRBs (e.g., Eichler & Levinson 1999; Mészáros & Waxman 2001; Norris 2003) and the clustering of many lGRBs duration around ∼10 s (Bromberg et al. 2012). Margutti et al. (2014) and Nakar (2015) proposed that jets fail to successfully drill through massive stars with extended envelopes (e.g., SNe Ib progenitors) and thus explosions of these stars are not coincidentally detected with GRBs. G22 showed that even among GRB progenitors that produce SNe Ic, only certain density, rotation, and magnetic profiles in the star support jet launching and breakout. The mildly relativistic outflow driven by these choked jets could power FBOTs (Gottlieb et al. 2022c), whose estimated rate is ∼10³ Gpc⁻³yr⁻¹ (Coppejans et al. 2020).

Zoom In Zoom Out Reset image size

GRB lightcurves exhibit three characteristic timescales (Nakar & Piran 2002). The first is the signal duration (∼10–100 s), which corresponds to the total work time of the engine. The second timescale is the rapid variability (∼10–100 ms), which corresponds to the stochastic fluctuations in the central engine accretion and launching mechanism and potentially to instabilities that develop during the interaction of the jet with the stellar envelope. Figure 1(b) shows the jet power that varies by about half an order of magnitude over the timescales of ∼10–100 ms, consistent with observations. Thus, we suggest that the rapid variability originates in the central engine, irrespective of the jet–star interaction (Figure 7), in agreement with Gottlieb et al. (2021a).

The third timescale corresponds to quiescent times (∼1–100 s; Ramirez-Ruiz & Merloni 2001; Nakar & Piran 2002), the origin of which is not fully understood. Our simulation shows that quiescent times naturally arise due to the jet tilt. When the jet is not pointing in the direction of the observer, its emission is beamed away, and the lightcurve becomes quiescent, potentially with some periodicity on timescales of ∼1–10 s (e.g., GRB940210). The farther away the observer is from the polar axis, the smaller the ratio between the active signal time to quiescent time, and the lower the observed jet emission efficiency (Figure 7). In addition to the quiescent times from the jet tilt, intrinsic temporal shutoffs of the engine can also produce quiescent times, albeit shorter. Figures 1(d) and 4 (bottom panel) illustrate ∼1 s fluctuations in the jet power, with relativistic jet blobs that are separated by ∼1 lt-s of slow material. If the structure advances homologously to the photosphere, the slow material would produce a quiescent episode in the GRB lightcurve.

The origin of the prompt GRB emission is a matter of active debate, with many fundamental questions that remain unanswered to date, among which (i) is the observed emission subphotospheric or originating in optically thin regions? (ii) What is the underlying mechanism responsible for energizing the electrons that shape the nonthermal spectral tail? The two leading candidates for explaining the latter are particle acceleration in shocks ⁸ (see reviews by, e.g., Piran 1999, 2005; Mészáros 2006; Kumar & Zhang 2015) and acceleration in current sheets formed by magnetic reconnection (e.g., Thompson 1994; Spruit et al. 2001; Lyutikov & Blandford 2003; Giannios 2008; Zhang & Yan 2011; Barniol Duran et al. 2016). These two mechanisms are directly connected to the energy composition in the jet, as the former requires strong shocks that can form in hydrodynamic or mildly magnetized flows (σ ≲ 0.1), whereas the latter requires magnetically dominated flows with σ ≳ 1 (e.g., Sironi et al. 2015). Thus, constraining the jet magnetization at the emission zone is an important step toward solving the prompt emission puzzle.

Early models of steady, Poynting-flux-dominated jets showed that efficient conversion of magnetic energy to bulk motion is possible in narrow jets that are confined by an external medium, such that the jet magnetization can drop to σ ∼ 1 (Tchekhovskoy et al. 2008; Komissarov et al. 2009; Lyubarsky 2009, 2010, 2011). Further acceleration accompanied by a drop in σ is possible if the jets undergo a sudden sideways expansion (Tchekhovskoy et al. 2010; Komissarov et al. 2010), or if the jet is composed of individual pulses that expand in the radial direction (Granot et al. 2011). Such models require the jet plasma to maintain force-free conditions and avoid dissipation of energy via other processes, e.g., magnetic reconnection. Jet collimation often leads to the formation of narrow nozzles in which magnetic fields can become kink unstable and dissipate magnetic energy into heat (e.g., Lyubarsky 2012; Mizuno et al. 2009, 2012; Bromberg et al. 2019). In the case of continuous jets in dense media, Bromberg & Tchekhovskoy (2016) showed that nozzle dissipation is a dominant process that can dissipate half of the jet magnetic energy into heat, rendering the outcome of the ideal acceleration models questionable. Further dissipation may continue above the collimation nozzle mediated by turbulence (Bromberg et al. 2019); however, the state of the plasma that exits the star and when it reaches the emission zone remains an open question.

Our simulations suggest that the jet intermittency and tilt increase the energy dissipation even further by allowing for the entrainment of heavy material into the jet, reducing the jet magnetization to values σ ∼ 10⁻¹ upon breakout from the star (Figure 3), and possibly somewhat higher magnetization if σ₀ ∼ 10³ (see Section 3.2). By assuming the acceleration of a hydrodynamic-dominated jet, the optical depth evolves as τ ∝ t⁻³, we find that the photosphere is located at ⁹ r ∼ 10¹² cm, similar to the photosphere found in hydrodynamic simulations (Gottlieb et al. 2019). If this magnetization persists in the jet through the propagation to the photosphere such that the jet keeps its hybrid composition, as indicated by Figure 3, it renders both shock acceleration and magnetic reconnection as plausible energizing mechanisms. We note that the outflow is yet to reach a self-similar structure and the optical depth is likely to be somewhat different at later times; however, our conclusion does not depend on these details.

In Section 3 we showed that not only do the jets become partly hydrodynamic after breakout, but also their extended structure takes the same shape as that of hydrodynamic jets. It is therefore useful to examine the emission from these jets through previous numerical studies of 3D hydrodynamic jets. Gottlieb et al. (2019) found that owing to the collimation nozzle, which converts the jet energy to thermal, jets reaccelerate at the larger radius of the collimation shock such that their coasting radius is above the photosphere. Consequently, hydrodynamic jets (and initially highly magnetized jets) inevitably produce a substantial photospheric component that can explain the high γ-ray emission efficiency. Our simulations show that the above arguments also apply to highly magnetized jets. Specifically, we find that the postbreakout specific enthalpy allows 30%–80% radiative efficiency (see Figure 7), depending on σ₀.

One important difference that can arise between this work and that of hydrodynamic jets in Gottlieb et al. (2019) is at radii smaller than the dissipation radius, namely before the jet becomes hydrodynamic. In their paper, they found that the mildly relativistic collimation shock at the jet base produces a rest-frame temperature that is maintained at ∼50 keV by pair production, which corresponds to the observed temperatures of a few hundreds keV. If at the collimation shock zone the jet still maintains σ ≫ 1, as might be the case if σ₀ > 100, then the energy dissipation efficiency in the shock is low, and the resulting spectrum would be different.

The hybrid magnetic and thermal composition with a magnetization of σ ∼ 0.1 may broaden the thermal spectrum (see also Thompson & Gill 2014) via reconnection and/or synchrotron emission. In addition to that, our simulation allows the formation of internal shocks that arise from the intermittent jet structure, which can similarly alter the spectrum if their radiative efficiency is high and they occur either above the photosphere or in a moderate optical depth (Parsotan et al. 2018; Ito et al. 2019). We conclude that our simulations suggest that the origin of GRB emission includes multiple components of photospheric, magnetic reconnection, and internal shocks, which can generate the observed Band function.

We summarize the GRB emission in Figure 7, where we show the radiative efficiency ≡ (h − 1)/h as measured at ¹⁰ r = 2 R_⋆ in the simulation where σ₀ = 15. Owing to the jet wobbling, observers at θ = θ_t + θ_j can detect high-efficiency emission. The jet motion also leads to the appearance of long quiescent times in the signal, while the jet intermittency results in a short timescale variability. The efficiency is typically higher close to the polar axis, as more episodes of jet pointing at this direction are expected.