Stars Bisected by Relativistic Blades

The governing equations in this setup are the standard special-relativistic hydrodynamic equations:

where $\rho$ is proper fluid density, $[u^{\mu}]=\Gamma(1,\vec{\beta})$ is four-velocity, $\Gamma=(1-\beta^{2})^{-1/2}$ is Lorentz factor, $\beta$ is velocity in units of the speed of light, $c$, which is unity in our setup, $\Theta$ and $\Psi$ are source terms, and $T^{\mu\nu}$ is the stress-energy tensor for a perfect fluid,

where $h=1+\varepsilon+p/\rho$ is total specific enthalpy, $\varepsilon$ is internal energy, $p$ is pressure, and $\eta^{\mu\nu}$ is the Minkowski metric with signature $(-,+,+,+)$. The set of Equations 1 – 3 become closed when choosing an ideal gas equation of state $p(\varepsilon)=(\hat{\gamma}-1)\rho\varepsilon$ where $\hat{\gamma}=4/3$ is the ratio of specific heats at constant pressure and volume.

The source terms are modeled as Dirac delta distributions such that the power density of the engine has the form

where $L_{\rm eng}$ is the engine power integrated over the entire sphere, the Dirac deltas are written as Gaussian approximations:

where $r_{n}$ is the effective radius of the engine nozzle, $\mu=\cos\theta$, $\mu_{n}$ is the direction of the beam, and $f$ is the geometric factor of the jet — i.e., $f=\theta_{0}^{-1}$ for a lamina jet while it is $\theta_{0}^{-2}$ for a classical jet where $\theta_{0}$ is the injection angle. The function $g(t)$ is taken to be a sigmoid decay,

where $\tau$ is the engine duration and $\xi$ is the sharpness of the drop-off. The remaining source terms are constructed from Equation 4, where the radial momentum density source term is

and the baryon loading term is

where we define $\Gamma_{0}$ as the injected Lorentz factor and $\eta_{0}\equiv\dot{E}/\dot{M}_{0}$ as the radiation to baryon ratio where $\dot{E}$ is the energy outflow rate near the engine and $\dot{M}_{0}$ is the initial mass outflow rate​ ¹¹1Some texts call $\eta$ the dimensionless entropy or the initial random Lorentz factor.. The engine duration can be set by requiring $\tau_{\rm bo}<\tau_{*}$ where $\tau_{*}$ is the spin-down time of the magnetar,

$\tau_{\rm bo}$ is the breakout time, $\omega$ is the rotational frequency, $E_{\rm rot}$ is the rotational energy, $L_{*}$ is the spin-down luminosity, $M_{\rm PNS}$ is the proto-neutron star mass, $R_{\rm PNS}$ is the proto-neutron star radius, $B$ is the surface equatorial magnetic field, and $T$ is the rotation period. In reality, it is not as simple as setting $\tau<\tau_{*}$ because the engine must do considerable work to displace enough stellar material to launch a successful jet. The breakout time of the beam is $\tau_{\rm bo}=\Gamma_{\rm ej}R/u_{\rm ej}(r)$, where $R$ is the stellar radius. We can compute $u_{\rm ej}(r)$ by noting the isotropic luminosity of the jet as in (Mészáros & Waxman, 2001),

and balancing the pressure of the jet head with that of sub-relativistic ejecta ahead of it, i.e., $u_{\rm j}^{2}h_{j}\rho_{\rm j}=u_{\rm ej}^{2}\rho_{\rm ej}(r)$, giving

where $q\equiv 4\pi/\Omega$, $\Omega$ is the solid angle, and we’ve made use of the fact that $L_{\rm iso}=qL_{*}$. The breakout time is therefore, $\tau_{\rm bo}\sim 3\Gamma_{\rm ej}u_{\rm ej}RM_{\rm ej}/qrL_{*}$. Now, we revisit $\tau_{\rm bo}<\tau_{*}$:

where the last equality stems from assuming the rotational energy is extracted with perfect efficiency. Equation 13 requires $q^{-1}<\frac{2r}{3R}\beta_{\rm ej}$ for a successful breakout. Quataert & Kasen (2012) compute the half-opening angle constraint for a classical jet (i.e., $q=2/\theta_{j}^{2}$) to be $\theta_{j}<\theta_{j,\rm max}\equiv(\beta_{\rm ej}/2)^{1/2}$. This fixes $r/R=3/8$ in our framework. The constraint on the half-opening angle for the lamina is then

which implies that the lamina must be much more collimated than a classical jet with the same power in order to break out of the star.

To set the engine power, we assume a “split monopole” (Weber & Davis, 1967) field geometry and scale the MSM luminosity calculated by Thompson et al. (2004), $L_{\rm eng}\approx L_{*}|_{t=0}=2.4\times 10^{51}\,\mathrm{erg~{}s^{-1}}(R_{\rm PNS}/12\,\mathrm{km})^{8/3}$. We also assume a relativistic, hot fireball component with $\eta_{0}=1000$, consistent with a radiatively driven engine​ ²²2In reality $\eta$ and other parameters are time variable, but since $L_{\rm*}\propto(1+t/\tau_{*})^{-2}$ we can set the engine duration to $\tau=\epsilon\tau_{*}$ where $\epsilon\ll 1$. This would ensure $|\Delta L_{*}/L_{*}|\sim 2\epsilon$ is small enough to assume rough constancy in our simulations.. In a scenario in which the magnetar is accreting, the spin-down is counteracted by angular momentum transport from the fallback material (e.g., Parfrey et al., 2016; Metzger et al., 2018), so we assume rough constancy within the timescale considered in this work by setting $\xi=100$ in Equation 7. We do not explicitly invoke magnetic fields since we assume all of the magnetic energy outside of the magnetar light cylinder is converted into kinetic energy, e.g., by magnetic dissipation, (see Spruit et al., 2001; Vlahakis & Königl, 2003; Komissarov et al., 2009). Our engine-progenitor model is reminiscent of Duffell & MacFadyen (2015, hereafter DM15) wherein they launch a successful collapsar jet through a star with $\theta_{0,\rm DM15}=0.1$. With this classical jet injection angle as a baseline, we use Equation 14 to set $\theta_{0}=0.005$ to match the power density of the classical jet. For an equatorial outflow, we set $\mu_{n}=0$ in Equation 6. The stellar model is an 18 $M_{\odot}$ Wolf-Rayet star that was originally a 30 $M_{\odot}$ Zero Age Main Sequence (ZAMS) star rotating at 99% breakup. The progenitor was evolved using the Modules for Experiments in Stelar Astrophysics (MESA; Paxton et al., 2011, 2013, 2015, 2018, 2019) code, and we invoke the density profile for this star fitted by DM15,

where $A=A_{*}\dot{M}/4\pi v_{\rm wind}=A_{*}\times 5\times 10^{11}\,\mathrm{g~{}cm^{-1}}$ is the ambient medium mass-loading parameter and the remaining parameters are defined in Table 1. We use an axisymmetric spherical-polar grid with logarithmic radial zones and uniform angular zones. We enforce 1024 radial zones per decade and $\delta\theta=\theta_{0}/N_{\rm beam}$, where $N_{\rm beam}$ is the number of zones within the half-opening angle of the beam. We fix $N_{\rm beam}=10$ in our simulations. The domain range is $r\in[0.001,10]R_{\odot}$ and $\theta\in[0,\pi/2]$, which corresponds to 4096 radial zones by 3142 angular zones. The initial pressure and velocity everywhere are negligible. All variables are made dimensionless through combinations of fundamental constants: $R_{\odot}$, $c$, and $M_{\odot}$. This concludes the initial conditions required to launch the relativistic lamina jet into the stellar progenitor. The problem is simulated using an open source GPU-accelerated second-order Godunov code entitled SIMBI (DuPont, 2023), written by this Letter’s first author. It uses a piecewise linear reconstruction algorithm to achieve second-order accuracy in space and second-order Runge-Kutta is employed for the time integration. $\theta_{\rm PLM}$, a numerical diffusivity parameter for second order schemes, is fixed to 1.5 in our simulations for more aggressive treatment of contact waves.

Our fiducial runs show promising relativistic breakout for highly collimated lamina jets. Figure 1 shows the time-evolved snapshots of the explosion from the initially stationary conditions to the relativistic breakout of the beam. Since the lamina is ultra-thin and radiative, it can easily push aside matter as the beam tunnels through the dense stellar interior, allowing it to get out within four light crossing times of the progenitor. Kelvin-Helmholtz instabilities develop deep in the interior as the beam propagates through the very thin cocoon, and the the jet core is naked once it accelerates down the steep density gradient ahead of it. The maximum Lorentz factor increases monotonically and exceeds the injected value, $\Gamma_{0}$. This hints at the fact that GRBs might not be very sensitive to the intricacies of their complicated central engines outside of a terminal bulk Lorentz factor $\Gamma_{\infty}\sim\eta$ (Zhang et al., 2003; Mészáros, 2006). In just a few percent of a typical MSM spin-down time, the beam breaks out at ultra-relativistic velocities before the cocoon has traversed $\sim 40$% of the star, affecting a clean slice through the progenitor.

Throughout its evolution, the lamina jet appears to experience collimation shocks as evidenced by Figure 2. In Figure 2, we also include another lamina jet with double the opening angle $\theta_{0,\rm wide}=0.01$ to gauge key differences in the cocoon evolution, which depicts a more advanced blast wave for a thinner beam. We believe this is due to the relatively negligible thickness of the blade-like structure of the outflow. The nature of this effect is two-pronged. That is to say, the $\theta_{0}=0.005$ beam’s working surface is is half as small as the $\theta_{0,\rm wide}=0.01$ beam, which means it is impeded by half of the mass and has double the pressure. As the thinner lamina pierces through the star, it interacts with less stellar material that is mixed in the Kelvin-Helmholtz layers of the cocoon and therefore the pressurized cocoon has a lesser impact on the collimation of the flow. Because of this, the beam encounters fewer interactions from rarefaction waves and shocks as the jet-cocoon interface is propagated throughout the star. Another way of putting it is that with the effective working surface of the engine reduced, the wave more easily travels through a star analogous to how a sharper knife more easily cuts through material while keeping the applied force fixed. We suspect that this is also the case for skinny classical jets. However, the knots from the collimation for classical jets are more prominent than for lamina jets as evidenced by previous numerical studies on classical jets (e.g., MacFadyen & Woosley, 1999; Aloy et al., 2000; MacFadyen et al., 2001; Zhang et al., 2003; Tchekhovskoy et al., 2008; Mösta et al., 2014; Mandal et al., 2022). Of course, this would have to be further analyzed in 3D to encapsulate a fuller picture of beam deflection and various instabilities that might arise.

Although the injection angle was $\theta_{0}=0.3^{\circ}$, the lamina at or above $\Gamma_{0}=10$ breaks out with a half-opening angle $\theta_{r}\sim 0.2^{\circ}$ ³³3Measured using $\sin\theta_{r}=E/E_{\rm iso}$ and the carries a maximum Lorentz factor $\Gamma_{\rm core}\sim 30$ at the time of the simulation end. An observer within the stellar equatorial plane and whose line of sight passes through the lamina centroid would see $1/\Gamma=3\%$ of the total structure. If extreme beaming took place — i.e., $\Gamma\gg\theta_{r}^{-1}$ — the lamina would travel through the interstellar medium (ISM) with fluid parcels causally disconnected from their neighbors which would help maintain a non-spreading outflow until it sweeps up a mass $M/\Gamma$ and slows down (Porth & Komissarov, 2015). In Figure 3, we compute the cumulative isotropic-equivalent energy per solid angle at simulation end,

where $\Gamma_{\rm c}\beta_{\rm c}$ is the four-velocity cutoff. Moreover, we estimate the motion of the bulk flow by noting the mean energy-weighted four-velocity,

moving above some fixed value. We are interested in the material which is moving at or above the injected Lorentz factor, $\Gamma\geq\Gamma_{0}=10$. Therefore, the ultra-relativistic component gives $\langle\Gamma\beta\rangle_{E}\sim 15$ as the velocity of the bulk flow. From Figure 3 we find that $E_{k,\rm iso}(>15)=3\times 10^{52}\,\mathrm{erg}$ is focused purely in the equator. This beam has mass $M(>15)|_{\theta=90^{\circ}}=E_{k,\rm iso}/\langle\Gamma\rangle_{E}\approx 10^{-3}M_{\odot}$. Assuming the MSM engine spontaneously shuts off, the lamina will decelerate after sweeping up $7\times 10^{-5}M_{\odot}$ which, for $A_{*}=1$, occurs at a radius $r_{\rm dec}=2\times 10^{16}\,\mathrm{cm}$ or $6\times 10^{5}R$.

Although the lamina breaks out of the star inefficiently (i.e., $\tau_{b}\sim 4R$), the cocoon has little time to traverse the remaining undisturbed star before the jet head outruns the explosion. To get a sense of the remaining energy that might be attributed to a supernova, we estimate this by summing all of the total energy available in the slowest material: $E_{\rm T}(<0.1)\sim 3\times 10^{51}\,\mathrm{erg}$. Only about $30\%$ of this total energy is kinetic at the time of the simulation end. Thus, once the cocoon finishes its journey throughout the remaining stationary star and full conversion of the thermal energy into kinetic energy is complete, the energy liberated in the explosion is of order the canonical supernova explosion energy. About 10% the deposited energy from the engine is shared with the supernova component if it were to shut off spontaneously after four seconds.

The ring will come into causal contact with its edges when the angle of influence, i.e., the proper Mach angle (Konigl, 1980), is comparable to the angular size of the beam, viz., $\theta_{\rm M}/\theta_{r}\geq 1$. The relativistic Mach angle evolves as $\tan\theta_{\rm M}=u_{s}/u\propto 1/\Gamma$ and angular size of the ring is $\theta_{r}=r_{\perp}/r$ where $u_{s}$ is the proper sound speed and $r_{\perp}$ is the size of the ring perpendicular to the flow. Together, these functions imply the causality constraint,

Assuming the thermal energy dominates the rest mast energy of the particles, i.e., $\rho h\propto p$, constancy of energy implies $\Gamma^{2}r^{2}\Delta r\theta_{r}p=\text{constant}$ where $\Delta r$ is the width of the annular blast wave. Energy conservation together with mass conservation, $\Gamma r^{2}\Delta r\theta_{r}\rho=\text{constant}$, leads to the scalings $\Gamma\propto\rho^{1-\hat{\gamma}}$ and $r_{\perp}\propto r^{-2}\rho^{-\hat{\gamma}}$ after utilizing $\Delta r\sim r/\Gamma^{2}$, $r\theta_{r}=r_{\perp}$, and $p\propto\rho^{\hat{\gamma}}$. Equation 18 thus becomes

Equation 19 implies that for an ultra-relativistic ring with $\hat{\gamma}=4/3$, the critical ambient medium density slope is $k=9/5$ where $k<9/5$ implies full causality is reached while $k>9/5$ implies causality is lost. This critical ambient medium slope is different from the $k=2$ classical jet value calculated by Porth & Komissarov (2015). Note that $k=2$ is the value for stellar winds, so Equation 19 implies that small perturbations in the $\phi$-direction might excite interesting instability modes if the ring evolves in a wind-like environment, and we plan to address this in another paper. Moreover, in the scenario in which the ambient medium is not uniform in $\phi$, the ring-like blast wave will become corrugated in the $\phi$ direction like pizza slices of uneven length, which might drive quasi-periodic radiation signatures in the light curves. While interior to the star and in in full 3D, the blade is likely to bend and wobble which might excite further azimuthal instabilities, which we also plan to investigate in a follow up work.

The overall blade carrying a Lorentz factor $\Gamma_{\rm beam}\geq\Gamma_{0}=10$ has half-opening angle $\theta_{\rm beam}=0.2^{\circ}$, and it contains an ultra-relativistic core with $\Gamma_{\rm core}\sim 30$ and half-opening angle $\theta_{\rm core}=0.03^{\circ}$ at simulation end. If the engine was instantaneously shut off from the moment of breakout, the bulk flow, which moves with $\langle\Gamma\beta\rangle_{E}=15$, propagates more than 5 orders of magnitude beyond the stellar radius before slowing down. To gauge when these outflows will become observable on Earth, we compute the photosphere radius assuming a grey atmosphere optical depth of $2/3$,

where $f_{\rm KN}$ are the Klein-Nishina corrections, $\kappa_{\rm T}=0.4/\mu_{e}\,\mathrm{cm^{2}~{}g^{-1}}$, $\mu_{e}=2/(1+X_{H})\simeq 2$ is mean molecular weight per electron, and $A=A_{*}\times 5\times 10^{11}\,\mathrm{g~{}cm^{-1}}$. This implies that the environment becomes optically thin within just a few stellar radii from the source, making the blade visible almost immediately post breakout.

In the scenario where lamina jets are sources of long gamma-ray bursts (lGRBs), their observational implications are immediately realized due to simple geometric arguments (Thompson, 2005; Granot, 2005; DuPont et al., 2023). Once the jet slows down, it will spread in a single dimension, leading to a shallower break in the afterglow light curves that makes it clearly distinguishable from a classical GRB jet configuration. The characteristic afterglows in the low and high frequency bands for such rings were computed by DuPont et al. (2023) and show that, over both frequency ranges, the ring-like blast waves are intermediate between jet-like and spherical outflows. With the afterglow decay distributions being so varied (e.g., Panaitescu, 2005), we believe ring-like explosions are also candidate outflows in this arena. Note that in this work we took the extreme case of a split-monopole axisymmetric engine to achieve such relativistic outflows. A more modest dipolar engine might refocus the blade into bipolar bubbles (e.g., Bucciantini et al., 2007) if the magnetic hoop stress dominates once the magneto-centrifugal wind reaches the termination shock. Another outcome could be that if the engine power still dominates over the magnetic hoop stress, the equatorial wind might produce slower material at breakout due to having a smaller energy dissipation rate, but this likely broadens the types of transients created by the blast wave. In principle, Equations 1 – 7 are scale free, so rescaling opens avenues for providing explanations for transients like super-luminous supernovae, X-ray bursts, fast blue optical transients (FBOTs), and low-luminosity GRBs, to name a few.