Long-duration gamma-ray bursts (GRBs) accompany the collapse of massive stars and carry important information about the central engine. However, no 3D models have been able to follow these jets from their birth by a rotating black-hole (BH) to the photosphere. We present the first such 3D GRMHD simulations, which span over 6 orders of magnitude in space and time. The collapsing stellar envelope forms an accretion disk, which drags inward the magnetic flux that accumulates around the BH, becomes dynamically-important and launches bipolar jets. The jets reach the photosphere with an opening angle of 6 degrees and a Lorentz factor up to 30, unbind 90% of the star and leave the BH mass essentially unchanged after the initial core-collapse. We find that: (i) The disk-jet system spontaneously develops misalignment relative to the BH rotational axis. As a result, the jet direction wobbles with an angle of 12 degrees, which can naturally explain quiescent times in GRB lightcurves. The effective opening angle for detection, the sum of the tilt and opening angles, suggests that the intrinsic (beaming-corrected) GRB rate is lower by an order of magnitude than standard estimates. This implies that successful GRBs can be much rarer than currently thought and emerge in only ~0.1% of supernovae (SNe) Ib/c. A possible explanation is that jets are either not launched or choked inside most supernova Ib/c progenitors. (ii) The magnetic energy in the jet decreases due dissipation and mixing with the stellar material, resulting in jets with a hybrid composition of magnetic and thermal components at the photosphere, where ~20% of the gas maintains magnetization > 0.1. This indicates that both a photospheric component and magnetic reconnection play a role in the GRB prompt emission.
First collapsar simulation from launching to photosphere

3D visualization of the jet propagation

Download video

3D zoom-in visualization of the tilted disk going MAD

Download video

2D slice of the mass density

Color map is logarithm of mass density
Download video

2D slice of the magnetization and velocity

Color map is logarithm of magnetization (top) and asymptotic proper-velocity (bottom)
Download video

The simulation scale

Color map is logarithm of mass density (top) and magnetization (bottom)
Download video

The simulations have been conducted with the code H-AMR (Liska et al. 2019)