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

Long-duration γ-ray bursts (GRBs) accompany the collapse of massive stars and carry information about the central engine. However, no 3D models have been able to follow these jets from their birth via black hole (BH) to the photosphere. We present the first such 3D general-relativity magnetohydrodynamic simulations, which span over six orders of magnitude in space and time. The collapsing stellar envelope forms an accretion disk, which drags inwardly the magnetic flux that accumulates around the BH, becomes dynamically important, and launches bipolar jets. The jets reach the photosphere at ∼1012 cm with an opening angle θ j ∼ 6° and a Lorentz factor Γ j ≲ 30, unbinding ≳90% of the star. We find that (i) the disk–jet system spontaneously develops misalignment relative to the BH rotational axis. As a result, the jet wobbles with an angle θ t ∼ 12°, which can naturally explain quiescent times in GRB lightcurves. The effective opening angle for detection θ j + θ t suggests that the intrinsic GRB rate is lower by an order of magnitude than standard estimates. This suggests that successful GRBs are rarer than currently thought and emerge in only ∼0.1% of supernovae Ib/c, implying that jets are either not launched or choked inside most supernova Ib/c progenitors. (ii) The magnetic energy in the jet decreases due to mixing with the star, resulting in jets with a hybrid composition of magnetic and thermal components at the photosphere, where ∼10% of the gas maintains magnetization σ ≳ 0.1. This indicates that both a photospheric component and reconnection may play a role in the prompt emission.


INTRODUCTION
Long-duration gamma-ray bursts (lGRBs) are thought to emerge when collimated relativistic jets escape a collapsing massive star (collapsar; Woosley 1993;MacFadyen & Woosley 1999).The prospects of learning from this phenomenon about the properties of GRB progenitor stars and the physics of their central engines inspired analytic studies of the propagation of relativistic jets in stars (e.g., Matzner 2003;Lazzati & Begelman 2005;Bromberg et al. 2011b).The nonlinear behavior of jets implies that numerical tools are essential for modeling the jet evolution before and after it breaks out from the star.This realization motivated the first 2D hydrodynamical simulations of jet propagation in stars (e.g.Aloy et al. 2000;MacFadyen et al. 2001;Zhang et al. 2003Zhang et al. , 2004;;Mizuta et al. 2006;Morsony et al. 2007Morsony et al. , 2010;;Bucciantini et al. 2009;Komissarov & Barkov 2009;Lazzati et al. 2009;Mizuta & Aloy 2009;Nagakura et al. 2011;Mizuta & Ioka 2013;Duffell & MacFadyen 2015).ore@northwestern.eduDespite achieving a significant progress in understanding lGRBs, these models suffer from several major drawbacks: (i) Axisymmetry: 2D models artificially impose axial symmetry, which suppresses nonaxisymmetric modes and leads to numerical artifacts.The growth in computational power and improvements in algorithms has enabled 3D modeling of hydrodynamic jets in collapsars (López-Cámara et al. 2013, 2016;Ito et al. 2015Ito et al. , 2019;;Harrison et al. 2018;Gottlieb et al. 2019Gottlieb et al. , 2020bGottlieb et al. , 2021b)).These studies demonstrated that 3D jets can feature utterly different behavior than 2D, motivating full 3D jet studies.
(ii) Magnetization: Relativistic jet launching from stellar engines can either be driven electromagnetically by the rotation of a compact object -a black hole (BH; Blandford & Znajek 1977) or a neutron star (Goldreich & Julian 1969;Usov 1992;Thompson 1994;Metzger et al. 2011)or thermodynamically by the pair plasma produced by annihilation of neutrinos originating in the accretion disk (e.g.Eichler et al. 1989;Paczynski 1990;MacFadyen & Woosley 1999).It is becoming increasingly clear that the latter energy source falls short of the enormous amounts of energy required to power lGRB jets, lending support to the idea that the outflow is magnetically powered (e.g., Kawanaka et al. 2013;Leng & Giannios 2014;Liu et al. 2015).Numerical simulations of magnetized jets are challenging, and only a few 3D models of magnetized jet evolution have been performed over the years (e.g., Mignone et al. 2010;Porth 2013;Bromberg & Tchekhovskoy 2016;Striani et al. 2016).Simulations of magnetic jets propagating to the stellar surface were performed only recently, although with some limitations: Either the jet was injected with subdominant magnetic energy (e.g.Gottlieb et al. 2020aGottlieb et al. , 2021a) ) or the central engine duration was too short to allow a successful breakout of relativistic material from the star (Gottlieb et al. 2022a, hereafter G22).In both cases, the relativistic outflows were not modeled beyond the breakout phase.Nevertheless, these works highlighted the importance of jet magnetization in jet stability, structure, and propagation.
(iii) All aforementioned studies of jets in stars (apart from G22 and the 2D simulations of Komissarov & Barkov 2009) prescribed the jet launching from a grid boundary rather than including self-consistent launching via the rotation of a magnetized central compact object.Because each chapter in the jet journey influences the following ones, connecting the underlying physics at the central BH with observations can only be achieved through a complete modeling of the entire jet evolution, from a self-consistent jet launching by a spinning central BH to the emission zone.
G22 performed 3D general-relativity magnetohydrodynamic (GRMHD) simulations of collapsars to study the effect of the progenitor structure on the jet launching and breakout.They found that the type of outflow depends on the initial magnetization, rotation, and mass density profiles of the star.If the magnetic field is too weak, a relativistic jet is not launched, and instead a standing accretion shock dominates the outflow.If the stellar rotation is too slow, a quasi-spherical outflow is driven by the engine, rather than a collimated jet.If both the magnetic field is strong and the rotation is fast enough to allow for an accretion disk to form, a relativistic jet is launched.Their simulations showed that the jet is subject to strong magnetic dissipation at the collimation nozzle, which renders the jet essentially hydrodynamic upon breakout from the star.However, the spatial resolution at the narrow collimation nozzle was marginal to verify the robustness of this result.Additionally, it showed that some of the shocked gas remains bound and free-falls toward the BH before it is deflected by the jet to fall onto the accretion disk.This tilts the accretion disk and, subsequently, the jet orientation.
Here, we build on the results of G22 to present a longawaited 3D GRMHD simulation that follows the jets for their entire evolution: from a self-consistent launching near a BH to an unprecedented distance of ∼ 10 6 gravitational radii.To the best of our knowledge, this is the first simulation that features 3D relativistically magnetized jets that reach relativistic Lorentz factors after breakout from the star.With the highest-resolution GRB simulation to date, we investigate the magnetic dissipation at the collimation nozzle and study the observational implications of the disk-jet tilt on the post-breakout outflow at 10 stellar radii.
The structure of the paper is as follows.In §2 we recap the numerical setup of G22 and discuss the modifications required to produce steady relativistic jets that operate for the entire duration of the simulation, ∼ 18 s.In §3 we outline the main physical processes that take place at the BH horizon, the magnetic dissipation, the jet tilt and post-breakout structure.In §4 we discuss the implications of our results for GRB rates, the variability and quiescent times in GRB lightcurves, and the powering mechanism for the GRB prompt emission.In §5 we summarize the results.

SETUP
G22 performed 3D GRMHD simulations of self-consistent jet launching in collapsars, using the 3D GPU-accelerated code H-AMR (Liska et al. 2019).They modeled the conditions for jet launching in collapsars, and connected the progenitor magnetic field, rotation and mass density profiles to the emerging jet properties.However, the central engine in G22 did not launch jets for a sufficiently long time to match the observed durations of lGRBs (∼ 2−100 s) and allow for the jets to break out of the star.G22 suggested that a change in the stellar magnetic field configuration can increase the duration of the active jet phase.Here we use nearly identical initial conditions to G22, but modify the initial magnetic field configuration to attain a longer active jet phase.
The initial conditions of the simulation adopt a Kerr BH of mass M BH,0 = 4M and dimensionless spin a = 0.8, embedded in a cold, compact Wolf-Rayet-like star of radius R = 4 × 10 10 cm and mass M ≈ 14 M .The density profile of the progenitor star at the beginning of the simulation (upon BH formation) is where ρ 0 is set by the requirement that M = R * 0 ρ(r)dV and r g is the BH gravitational radius.G22 found that moderate values of α ∼ 1 can be consistent with many GRB observables: signal duration of 10 s, jet luminosity L j ≈ 10 51 erg s −1 , and the absence of long-term trends in the time evolution of the lightcurve.We adopt α = 1.5 that, as we show here, satisfies all of these observables.
The specific angular momentum profile of the stellar envelope is spherically symmetric, increasing to 70r g , and then plateaus: where ω 0 is constant.We adopt a uniform vertical magnetic field B core = B 0 ẑ inside the stellar core of radius r c = 10 8 cm.While G22 used a magnetic dipole outside the magnetic core, here we employ a shallower magnetic profile with a vector potential.The reason for this choice lies in the need to advect enough magnetic flux to the BH also at late times, so that the central engine remains active throughout the entire duration of the simulation.G22 showed that the jet duration scales with the ratio of the fastest growing magnetorotational instability (MRI) wavelength mode to the disk scale height ∝ B 1/2 , where B is the magnetic field magnitude.Thus, a slow radial decay in the magnetic field allows a longer jet launching duration.Overall, the magnetic vector potential is where µ ≈ B 0 r 2 c is the magnetic moment of the uniformly magnetic core, and the second term determines how fast the magnetic field drops near the stellar edge.The choice of B 0 is such that the maximum magnetization σ ≡ b 2 /(4πρc 2 ) in the star is max σ ∼ 10 −1.5 ; here b is the comoving magnetic field strength.These values correspond to maximal magnetic field in the core of B ∼ 10 12.5 G, which is required for the jet to overcome the ram pressure of the infalling gas, as found by G22.This choice of initially strong magnetic field implicitly assumes that the magnetic field was amplified prior to the BH formation (onset of the simulation), presumably following the rapid free-falling plasma that accumulates in the core after the initial collapse.
The simulation is performed with an ideal equation of state with adiabatic index of 4/3, appropriate for relativistic or radiation dominated gas.For numerical stability purposes we set density floors in the code by setting the maximal magnetization in the simulation to σ 0 = 15.The maximum magnetization is roughly the maximal asymptotic velocity that a fluid element can reach (up to a correction factor of order unity for the thermal energy of launched hot jets).We note that these values are much lower than those of GRBs.Therefore, we also carry out an identical simulation with σ 0 = 200, with which we address the potential effects σ 0 on our results.The plasma magnetization in this case is more sensitive to an artificial injection of floor mass density in regions where σ approaches the limit set by σ 0 .Nevertheless, we can use the outcome both as a consistency check for the results of the simulation with σ 0 = 15, and as a lower limit to the expected magnetization outside of the star.Throughout the paper, we compare the results of the two simulations and address the possible interpretations.
For the numerical integration we use a local adaptive timestep and 4 levels of adaptive mesh refinement (AMR).We use a spherical grid with a logarithmic cell distribution in the radial direction and uniform distributions in θ and φ directions.The radial grid extends from just inside the event horizon, ∼ 6 × 10 5 cm, to ∼ 6 × 10 11 cm with numerical resolution at the base AMR level of N r × N θ × N φ = 384 × 96 × 192 cells, in the r-, θ-, and φ-directions, respectively.We use a novel refinement criterion that at each radius r calculates the jet and cocoon (see its definition below) half-opening angles based on the specific entropy of the fluid, and if either one of them contains less than the desired number of cells, ∆N θ = 96 or ∆N φ = 192, the grid refines to the next AMR level, until it reaches the desired number of cells across each dimension, up to 4 levels of refinement.Therefore, the effective number of cells at the maximum level of refinement is N r × N θ × N φ = 6144 × 1536 × 3072 ≈ 3 × 10 10 .To avoid numerical artifacts associated with the jet propagation along the polar axis, we tilt the metric by 90 • such that the angular momentum orbital plane is the ŷ−ẑ plane, and the stellar rotation is around the x-axis.To avoid confusion, in the text and figures we refer to the ẑ-axis and θ as the stellar rotation axis and the polar angle relative to the rotation axis, respectively.

Early jet evolution
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 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 Ṁ ∝ t 1−2α/3 (in the case of a free falling stellar envelope with a power law density profile index α) to Ṁ ∼ t 2(1−α)/3 (G22).Figure 1(a) shows that the time evolution of the mass accre-  The jet launching efficiency, η, remains of order unity at all times, signaling that the system is accreting in the MAD regime.The jet-cocoon outflow unbinds most of the stellar envelope, ejecting ∼ 90% of the total stellar mass by t ≈ 16 s (panel (d)).As most of the stellar mass becomes unbound early on, the available gas for accretion is limited such that the BH mass M BH increases only by a few percent over the duration of the simulation (panel (e)).The red lines show moving averages of the highly variable quantities.tion rate at r = 5r g decays as expected (for our choice of α = 1.5),Ṁ ∝ t −1/3 during the first few seconds.The jet breaks out from the star after ∼ 3 s, and when the jet head reaches ∼ 2R 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.
Figures 1(b) and (c) depict the jet luminosity L j at r = r g , and the jet power efficiency η = L j / Ṁ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 ∼ Ṁc 2 , as shown in Fig. 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 (Fig. 1d).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 Fig. 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 free-fall 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 a substantial tilt values of up to ∼ 60 • (see 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 Fig. 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 non-tilted 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 ( §3.4).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 §3.4 and §4 we study the consequences of disk tilt on the post-breakout stage and its observational implications, respectively.

Magnetic dissipation
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 Fig. 2 shows a drop in the jet magnetization from σ 0 > 10 near the BH to σ ∼ 1 within the inner 0.1R .Figure 3(a) shows the radial profile of the angleaverage 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 ∞ ≡ (Γ 2 ∞ − 1) 1/2 > 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: Γ ∞ ≡ −u t (h + σ), where u t is the covariant time component of the four-velocity, and h = 4p ρ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 9 cm is rather smooth, with a similar shape at all times, whereas it shows a more complicated behavior at larger radii.At late times, Profiles of the weighted averages of σ (blue), proper velocity u (green), and asymptotic proper velocity u ∞ (black), shown for jets with σ 0 = 15 (thick lines) and σ 0 = 200 (thin lines) when the jet head is at r ≈ 6R (t = 10 s).The jets are launched with high magnetization (σ 1), but the magnetic energy is efficiently converted into kinetic energy until r ≈ 10 9.5 cm, maintaining constant u ∞ .After this point, the mixing between the jets and the star increases, thereby reducing σ further, and u ∞ as well, as seen by their correlated values.The jets become mildly magnetized at r ≈ 10 10 cm.Panel (c): Profiles of the maximum values of the quantities in (b) for the jet with σ 0 = 15 at t = 18 s.It shows that some jet elements reach Γ ∼ 10 after breakout with maximum σ ∼ 0.3.the jet magnetization outside of the star saturates at values 10 −2 σ 10 −1 .Figure 3(b) depicts the angle-average profiles of σ, u and u ∞ in the jet with the initial magnetization of σ 0 = 15 (thick lines) at t = 5 s.At r 10 9 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 9 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 Fig. 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 σ 0 = 200 as thin lines in Fig. 3(b).A comparison to the simulation with σ 0 = 15 shows that at the acceleration zone (r < 10 9 cm), the behavior is similar, but the higher value of σ 0 enables acceleration to higher velocities.In the mixing zone, r 10 9 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 lowermagnetization case.Outside of the star, the magnetization and terminal proper velocity saturate at values of σ ∼ 10 −1 and u ∞ ∼ 10 with peak values of ∼ 0.3 and ∼ 30, respectively.This indicates that jets with higher σ 0 can accelerate to higher u ∞ and σ.The decrease of u ∞ to ∼ 0.1σ 0 in both models hints that in order for the jet to reach u ∞ 100, as expected in GRBs, the initial magnetization should be σ 0 10 3 .
Figure 3(c) depicts the maximum values of σ, u and u ∞ at each radius at the end of the simulation with σ 0 = 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 σ 0 = 15 and σ 0 = 200, respectively (see also Fig. 5(c)).Therefore, if σ 0 10 3 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 §4.3).

Post-breakout structure
The propagation of the jet in the dense stellar envelope forms a hot cocoon that engulfs the jet.The top panel of Fig. 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 −3 .The bottom panel shows the final state of the simulation when the forward shock reaches r ∼ 12R 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.
The first 3D numerical characterizations of the postbreakout structure of GRB jets were calculated by Gottlieb et al. (2020aGottlieb et al. ( , 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 log u ∞ , 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 −2 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); 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 −2 .Indeed, similar values are consistently found in all FBOTs with The energy distribution per logarithmic asymptotic proper velocity is quasi-uniform irrespective to the choice of σ 0 .Both angular and radial distributions are consistent with those of hydrodynamic jets (Gottlieb et al. 2021b).Panel (c): The fraction of energy in matter with magnetization larger than σ out of the total energy.It can be seen that a significant amount (∼ 20%) of the outflow energy is carried by plasma with σ 10 −2 (σ 10 −1 ) when σ 0 = 15 (σ 0 = 200).
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 Fig. 5(a) indicates that the structure has reached the asymptotic stage in the simulation with σ 0 = 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(Gottlieb et al. , 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, since 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 behaviour seen in the bottom panel of our Fig. 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.Fig. 5(b) shows the (radial) distribution of the energy in log of u ∞ for matter outside of the star.The distribution is in agreement with previous results of a flat energy distribution 1 in cocoons of jets with different magnetizations, implying that this is a universal outflow distribution which 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 to the jet magnetization.The choice of σ 0 = 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 σ 0 = 200 outside the star, when the jet head is at r ≈ 10 11 cm (blue line).The distribution remains flat, but the cut-off 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 σ 0 such that it is robust.We conclude that both angular and log(u ∞ ) energy distributions are consistent with those of hydrodynamic jets, even to a better extent than with those of weakly magnetized jets.
Finally, Fig. 5(c) displays the cumulative energy distribution of matter that broke out of the star and has magnetization and upon breakout (purple) shows that θ t is attenuated as the jet propagates in the star and its orientation is focused toward the stellar rotation axis.Bottom: Colormap of the jet tilt angle shows that it varies significantly inside the star, but is conserved after breakout along a light-like streamlines.The asymptotic jet tilt angle is θ t 0.2 rad.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 σ 0 = 15.However, in the model with σ 0 = 200 about ∼ 20% of energy 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 §4.3).

Jet tilt
In §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 initial orientation can be far from the ẑ-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.1R (blue) and at r = R (purple).At t 3 s the tilt angle at 0.1R deviates significantly from the polar axis with θ t 0.5 rad.The jet elements at r = 0.1R 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 breakout from the star.
The panel of Fig. 6 a space-time 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 that the jet maintains a roughly constant tilt angle along light-like streamlines (diagonal paths on the map).

OBSERVATIONAL IMPLICATIONS
We have presented the results of the first 3D GRMHD simulation of collapsars that follow the relativistic jets from their formation at the BH to the photosphere (at ∼ 10 12 cm, see below).We discuss the significance of our findings in the context of observations; however a detailed study that computes the observational signature numerically is required and should be addressed in a future study.

GRB rate
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 Γ = θ −1 j .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 −3 yr −1 (e.g., Wanderman & Piran 2010) translates to a total rate of ∼ 100 Gpc −3 yr −1 .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 (Fig. 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 −3 yr −1 , 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 ∼ 10 3 Gpc −3 yr −1 (Coppejans et al. 2020).

Variability, quiescent times and periodicity
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.Fig. 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 (Fig. 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 (Fig. 7).In addition to the quiescent times from the jet tilt, intrinsic temporal shutoffs of the engine can also produce quiescent times, albeit shorter.Figs.1(d) and 4 (bottom panel) illustrate ∼ 1 s fluctuations in the jet power, with relativistic jet blobs that are separated by ∼ 1 light-second of slow material.If the structure advances homologously to the photosphere, the slow material would produce a quiescent episode in the GRB lightcurve.

Emission mechanism
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 non-thermal spectral tail?The two leading candidates for explaining the latter are particle acceleration in shocks2 (see reviews by e.g.Piran 1999Piran , 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 2009Lyubarsky , 2010Lyubarsky , 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. 2009Mizuno et al. , 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 radiative efficiency at r = 2R for the jet with σ 0 = 15.The jet's wobbling motion is responsible for two features in the signal: i) long quiescent times at all viewing angles and ii) observers far from the polar axis, at θ t + θ j can detect high radiative efficiency.However, we find that in general, the efficiency decreases with angular distance from the axis.the entertainment of heavy material into the jet, reducing the jet magnetization to values σ ∼ 10 −1 upon breakout from the star (Fig. 3), and possibly somewhat higher magnetization if σ 0 ∼ 103 (see §3.2).By assuming the acceleration of a hydrodynamic-dominated jet, the optical depth evolves as τ ∝ t −3 , we find that the photosphere is located at 3 r ∼ 10 12 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 Fig. 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 §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 re-accelerate 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 post-breakout specific enthalpy allows 30% − 80% radiative efficiency (see Fig. 7), depending on σ 0 .
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 σ 0 > 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 at4 r = 2R in the simulation where σ 0 = 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.

CONCLUSIONS
Using a novel GPU-accelerated AMR-enabled GRMHD code H-AMR, we performed a collapsar simulation that is: (i) the first to follow jets from their self-consistent launching near the BH all the way to the emission zone at ∼ 10 12 cm, over 6 orders of magnitude in space and time, (ii) the highest-resolution 3D GRB simulation to date, and (iii) the first collapsar simulation that includes highly magnetized jets that emerge from a star.It is thus a major step forward in un-derstanding the evolution of magnetized jets in general, and in collapsars in particular.
The central engine launches relativistic Poynting-fluxdominated jets intermittently over ∼ 10 − 100 ms timescales, consistent with the observed GRB lightcurve rapid variability.The jet is powered by magnetically saturated, MAD accretion and maintains the energy efficiency of order unity.The extended jet-cocoon structure reaches the photosphere with Lorentz factor Γ ∼ 30 (for σ 0 = 200) and jet opening angle of θ j ∼ 6 • .The powerful outflow (∼ 10 53 erg) unbinds most of the stellar envelope ( 90%), but cannot explain SNe associated with GRBs due to the flat radial energy profile in the logarithm of the proper velocity space.A further investigation of the outflows that includes neutrinos may alter the SN-GRB picture and will be considered in a future work.The large fraction of unbound mass also implies that within the collapsar framework, the BH mass does not significantly change after the initial core-collapse, and therefore, the final BH mass should be similar to the stellar core mass.
The two main findings of our simulation are as follows: 1. Jet tilt: Infalling gas from the bound parts of the cocoon provide non-axisymmetric kicks to the disk and tilt it.Subsequently, the jet axis tilts as well and wobbles over a range of ∼ 0.8 rad.As the tilted jet encounters the dense outer cocoon, it is refocused toward the low-density funnel drilled by the previous jet activity.Closer to breakout, the jet no longer changes its orientation so its post-breakout tilt angle is θ t ≈ 0.2 rad.The jet tilt has important implications for observations of GRBs: (i) The tilt translates into an effective opening angle for detection of ∼ θ t +θ j ≈ 0.3 rad.This suggests that about ∼ 10% of the GRBs in the universe can be observed from Earth, implying that the total (intrinsic) GRB rate is ∼ 10 Gpc −3 yr −1 , about an order of magnitude lower than previous estimates.This raises the possibility that most jets fail to break out from collapsars.This can happen when some progenitor stars keep their outer layers of He and/or H that choke the jets, while others may just lack the conditions (e.g.magnetic and density profiles) essential for jet launching and breakout.(ii) Quiescent times naturally arise due to jet wobbling, when the jet is pointing away from the observer.For an observer who is far from the polar axis, only a few lower-efficiency episodes will be apparent.Shorter (∼ 1 s) quiescent times can also emerge due to the intermittent nature of the central engine.
It is also noteworthy to mention the effect of such wobbling jets on the afterglow lightcurve.For a given observer who is aligned with the axis of one jet, the afterglow lightcurve will take the shape of a regular on-axis jet with a jet break that corresponds to the observed jet opening angle.This is because the emission of other jets is beamed away until their Lorentz factor becomes mildly relativistic, where the precise values depend on the offset of the jet from the observer.It then follows that jets that propagate outside the observer's line of sight may become visible only months after the explosion, and if these jets are comparably strong to the on-axis jet, their contribution to the lightcurve might not be detectable.Jets whose opening angles overlap with each other may lack a clear signature of a jet break in the lightcurve.We leave a full numerical study of the afterglow lightcurve from such structure to a future study.
2. Magnetic dissipation: Whereas the jet is launched as Poynting-flux dominated, its magnetic energy is continuously dissipated during its propagation.At small radii, the jet accelerates by converting magnetic to kinetic energy efficiently, with negligible amount of mixing between the jet and the star.As the jet propagates farther in the star, the mixing increases such that the jet magnetic (and kinetic) energy drops at a faster rate.The jet escapes from the star mildly magnetized, with σ ∼ 10 −1 ; after the jet breakout, the mixing weakens and the magnetization level remains steady.Thus, the jet structure is a hybrid composition of mildly magnetized and thermal parts, whereas its extended (radial and angular) structure is remarkably consistent with those found in hydrodynamic jets.In a companion paper, Gottlieb et al. (2022b), we show that the picture is different for short GRBs jets that propagate in light ejecta.Those jets are subject to weak mixing and thus can retain σ 1 at the photosphere.The consequences of the jet becoming mildly magnetized for the observed emission are the following: (i) A substantial fraction of the magnetic energy is deposited as thermal energy in the jet, thereby enabling the jet to reach the photosphere with enough thermal energy to efficiently generate γ-rays via photospheric emission.(ii) About 20% of the post-breakout plasma maintains σ > 0.1 (when σ 0 = 200), which, together with internal shocks, may transform the spectrum from thermal to the observed Band function.(iii) The need to generate jets with u ∞ implies that σ 0 ∼ 10 3 , in which case we anticipate σ ∼ 1 at the photosphere, such that magnetic reconnection will have an even larger contribution to the prompt signal.
Finally, we outline a few limitations of our numerical setup that might affect our results.First, the simulations do not include any cooling scheme, e.g.via neutrino emission, which may change the disk evolution and alter the jet launching efficiency and duration.Second, our simulation does not model the phase between the core-collapse and formation of the BH.During this stage, a proto-magnetar is anticipated to form and launch magnetized outflows along the axis of rotation.Such outflows may alter the progenitor structure, primarily along the rotation axis, and mitigate the later relativistic jet propagation.Similar low-density regions along the axis may also form by neutrino-antineutrino annihilation or by the fast rotation of the progenitor star.Third, in our simulations, the post-breakout jet's isotropic equivalent lu-minosity is L iso ∼ 10 54 erg s −1 , which translates to prompt γ-ray emission energy L iso,γ ∼ 10 53.5 erg s −1 , about an order of magnitude brighter than the peak luminosity of the observed lGRB distribution function (see e.g., Wanderman & Piran 2010), and corresponds to the brightest known events.In a future study, we will explore the aforementioned caveats by modeling the progenitor structure self-consistently from the time of core-collapse to the BH collapse, including neutrino scheme M1 and examining a variety of jets.We stress that while those effects are important, our main conclusions are anticipated to remain similar, at least qualitatively.

Figure 1 .
Figure 1.Time evolution in the simulation with σ 0 = 15: Panel (a): Accretion onto the BH initially decreases due to the inflation of the cocoon, before plateauing after a few seconds.Panel (b): The jet luminosity, which is calculated as the radial energy flux excluding the rest energy flux, scales proportionally to the mass accretion rate.Panel (c):The jet launching efficiency, η, remains of order unity at all times, signaling that the system is accreting in the MAD regime.The jet-cocoon outflow unbinds most of the stellar envelope, ejecting ∼ 90% of the total stellar mass by t ≈ 16 s (panel (d)).As most of the stellar mass becomes unbound early on, the available gas for accretion is limited such that the BH mass M BH increases only by a few percent over the duration of the simulation (panel (e)).The red lines show moving averages of the highly variable quantities.

Figure 2 .
Figure2.Top: 3D rendering of the σ 0 = 15 jet in the inner 2 × 10 9 cm shows a significant tilt of the disk (orange) and jet (red) axis (horizontal) with respect to the rotation axis of the star at 45 • .Although the jet is launched at different orientations at different times, it is deflected by the heavy outer cocoon material (blue gray) that engulfs the jet region toward the rotation axis at 45 • .Bottom: 3D rendering of the logarithm of jet magnetization shows the deflection of the jet propagation and a drop in the magnetization from σ ∼ 10 to σ ∼ 1.Here, the jet head is located at r = 0.1R .Movies showing the full evolution of the disk-jet tilt are available at http://www.oregottlieb.com/collapsar.html.

Figure 3 .
Figure 3. Panel (a): Radial magnetization profiles of the jets with σ 0 = 15, calculated by taking the weighted-average magnetization over the radial energy flux excluding the rest-mass energy flux and fluid elements with u ∞ < 1.The profiles at different times indicate that the magnetization in the jets does not change over time and remains at σ ∼ 10 −2 after breakout.Panel (b):Profiles of the weighted averages of σ (blue), proper velocity u (green), and asymptotic proper velocity u ∞ (black), shown for jets with σ 0 = 15 (thick lines) and σ 0 = 200 (thin lines) when the jet head is at r ≈ 6R (t = 10 s).The jets are launched with high magnetization (σ 1), but the magnetic energy is efficiently converted into kinetic energy until r ≈ 10 9.5 cm, maintaining constant u ∞ .After this point, the mixing between the jets and the star increases, thereby reducing σ further, and u ∞ as well, as seen by their correlated values.The jets become mildly magnetized at r ≈ 10 10 cm.Panel (c): Profiles of the maximum values of the quantities in (b) for the jet with σ 0 = 15 at t = 18 s.It shows that some jet elements reach Γ ∼ 10 after breakout with maximum σ ∼ 0.3.

Figure 4 .
Figure 4. 3D rendering in the simulation with σ 0 = 15 shows the logarithm of u ∞ (top) and σ (middle) in the jet (yellow)/cocoon (blue) structure when the jet head reaches r = 0.5R at t ∼ 1.5 s.It shows that at this radius the jet magnetization is already σ 0.1.Bottom: The post-breakout outflow after 18 s, when the head of the intermittent jet (red) is at ∼ 12R and the cocoon (green blue) explodes the star (yellow) entirely.The color coding is a combination of the mass density with u ∞ in order to show both the jet and the cocoon.Movies showing the full evolution of the outflow are available at http://www.oregottlieb.com/collapsar.html.

Figure 5 .
Figure5.The energy distribution of matter at r > R , excluding the rest mass energy, is shown according to various parameters.For the simulation with σ 0 = 15 we show the distribution at three different times: t = 6 s (purple), t = 12 s (maroon) and t = 18 s (gray).The blue line represents the distribution of energy in the jet with σ 0 = 200 at t = 10 s.Panel (a): The angular energy distribution integrated over the azimuthal angle of the outflow can be modeled by a flat core and a moderate power law followed by an exponential decay.Panel (b): The energy distribution per logarithmic asymptotic proper velocity is quasi-uniform irrespective to the choice of σ 0 .Both angular and radial distributions are consistent with those of hydrodynamic jets(Gottlieb et al. 2021b).Panel (c): The fraction of energy in matter with magnetization larger than σ out of the total energy.It can be seen that a significant amount (∼ 20%) of the outflow energy is carried by plasma with σ 10 −2 (σ 10 −1 ) when σ 0 = 15 (σ 0 = 200).

Figure 6 .
Figure6.The time evolution of the jet tilt angle θ t in the simulation with σ 0 = 15.Top: The tilt angle deep inside the star (blue) and upon breakout (purple) shows that θ t is attenuated as the jet propagates in the star and its orientation is focused toward the stellar rotation axis.Bottom: Colormap of the jet tilt angle shows that it varies significantly inside the star, but is conserved after breakout along a light-like streamlines.The asymptotic jet tilt angle is θ t 0.2 rad.
Figure7.The radiative efficiency at r = 2R for the jet with σ 0 = 15.The jet's wobbling motion is responsible for two features in the signal: i) long quiescent times at all viewing angles and ii) observers far from the polar axis, at θ t + θ j can detect high radiative efficiency.However, we find that in general, the efficiency decreases with angular distance from the axis.