© The Authors 2026

1. Introduction

The evolution of a star is primarily determined by its mass, metallicity, and rotation rate (Meynet & Maeder 1997; Janka et al. 2007). Massive stars, when nuclear fusion stops producing energy, collapse. If their mass exceeds about ∼25 M_⊙, they conclude their lives as stellar black holes (Heger et al. 2003). Gamma-ray bursts (GRBs) are classified based on their duration into two classes: short and long (Kouveliotou et al. 1993). Long GRBs are typically associated with Type Ib and Type Ic supernovae, which implies that their progenitors have lost their outer layers (Woosley & Bloom 2006). Nevertheless, new detections of long-duration GRBs without a supernova (SN) association have questioned the simple duration-defined dichotomy, and motivated multi-parameter class definitions as proposed in Zhang et al. (2009). The collapsar model is currently the most widely acknowledged explanation for the origin of long GRBs (Woosley 1993; MacFadyen & Woosley 1999). The typical energy release in such an event is between 10⁵¹ and 10⁵⁴ erg, with durations ranging from about 2 to 1000 seconds (Piran 2004; Kumar & Zhang 2015).

Variability demonstrated by the prompt phase is one of the key observables used to constrain theoretical models of GRB progenitors. Interaction of shocks emitted with slightly different Lorentz factors within the GRB jet was proposed by Daigne & Mochkovitch (1998). In addition, the collisionless shocks may leave an imprint in the afterglow emission (Spitkovsky 2008). The prompt phase has a rapid onset, and then the resulting light curve can vary over a range of timescales, sometimes showing multiple distinct, spiky episodes spread over the burst duration (Kumar & Zhang 2015). Observations of afterglow light curves suggest a long-lived central engine, a wide distribution of Lorentz factors in the jets, and possibly the deceleration of a Poynting flux-dominated flow (Zhang et al. 2006).

In the collapsar model, the jet drills through the progenitor star, and the jet-star interaction may produce structures in the light curve. This imprint could filter the variability arising primarily from the central engine (Morsony et al. 2010). After breakout, the observed luminosity variations record the central engine variability carried by internal dissipation in the outflow due to shocks or magnetic reconnections (Zhang & Yan 2010; Bromberg & Tchekhovskoy 2016), which are both capable of reproducing the highly variable prompt emission seen in long GRBs (Kobayashi et al. 1997). The collapsar scenario also accommodates longer ‘quiescent’ gaps and late X-ray flares (Margutti et al. 2011) by invoking intermittent engine activity, when the jet can choke. The spectral and timing properties of prompt and X-ray flares point towards dissipation from renewed variable engine activity (Zhang et al. 2006). The alternative to a hyperaccreting black hole in the collapsing star is a magnetar model where different kinds of instabilities operate (Bucciantini et al. 2009). However, the observed properties of some GRBs are in contrast with this scenario and imply a strong evolution of the Lorentz factor and photospheric radius, e.g. Iyyani et al. (2013), that can be reconciled with an evolving black hole central engine.

Numerical simulations of collapsars often simplify the hyperaccretion phase and process of jet ejection by assuming a fixed black hole spin, its constant mass, and uniform flow magnetisation. If those parameters evolve, the recovery of prompt and early afterglow phases phenomenology is more appealing. Even modest changes of the Kerr parameter, or magnetic flux, can map into jet luminosity evolution (Tchekhovskoy et al. 2011). It has been shown that MAD/SANE transitions and flux eruption episodes can intermittently choke accretion, giving rise to multiple pulses and gaps due to engine ‘on-off’ behaviour, which match Swift-era flare phenomenology much better than a fixed central engine schematic (Lloyd-Ronning et al. 2016). Recently, axisymmetric GRMHD models in dynamical spacetime have pointed towards characteristic spin changes on long timescales, relevant for the overall prompt phase duration (Shibata et al. 2024).

In our previous two-dimensional study, we investigated the black hole mass and spin evolution, and we also incorporated effects of self-gravity during the collapse (Janiuk et al. 2023). We demonstrated that for self-gravitating collapsars, the central engine evolves faster and the accretion rate exhibits more rapid variability than in the models without self-gravity. We showed that if the infalling stellar envelope mass rivals or exceeds the seed BH mass, the disc-BH system is strongly perturbed and self-gravity enhances shocks and large-amplitude Ṁ oscillations. In those simulations, we were not able to produce successful jets. In this study, we present three-dimensional simulations with jets emerging from self-gravitating stars.

The article is organised as follows. In Sect. 2.1, we define the initial conditions of our models. Then, we give the equations used to implement self-gravity in the code in Sect. 2.2, and finally, we give the evolution equations and methods used by the GRMHD code in Sect. 2.3. In Sect. 3.1, we present our results about the effects of self-gravity on the jet energetics and properties. In Sect. 3.2, we focus on the effects of self-gravity on the opening angle of relativistic jets. In Sect. 3.3, we discuss the time-dependent evolution of the central black hole. In Sect. 3.4, we discuss accretion flows and the stability of the MAD state. In Sect. 4, we discuss our results in the context of GRB prompt emission, and in Sect. 5, we summarise our conclusions.

2. Simulation setup

2.1. Initial conditions

We used the initial conditions as described in detail in Janiuk et al. (2018), Król & Janiuk (2021), and Janiuk et al. (2023) for the quasi-spherical accretion. Our collapsar models assume that the entire configuration has had enough time to achieve equilibrium. Therefore, we started with the Bondi solution (Bondi 1952), endowed with a small angular momentum. The total mass available for accretion was defined at the start, and no additional matter was added later on (except numerical flooring). The equation of state is polytropic, given by P = Kρ^γ, where P is the pressure, ρ is the density, and K refers to the specific entropy. The relation is parametrised by the polytropic index γ = 4/3. The sonic radius r_s was set at the distance of 80 r_g from the central black hole, where r_g is the gravitational radius. Below this point, matter falls supersonically towards the black hole. The radial velocity $u_s^r$ at critical point was evaluated from

$$\begin{array}{c}\hfill {(u_s^r)}^2=\frac{1}{2r_s}·\end{array}$$(1)

The profile of radial velocity, u^r, as a function of r was obtained by numerically solving the relativistic Bernoulli equation (Michel 1972; Shapiro & Teukolsky 1983), given by

$$\begin{array}{c}\hfill {(1+\frac{\gamma }{\gamma -1}\frac{P}{\rho })}^2(1-\frac{2}{r}+{(u^r)}^2)=\mathrm{constant}.\end{array}$$(2)

The density distribution was evaluated based on the radial velocity, and scaled by mass accretion rate, Ṁ, which was adjusted to satisfy the specified target for the total mass of the evolved stellar core, M_★, in the simulation domain, and is expressed by

$$\begin{array}{c}\hfill {\displaystyle \rho =\frac{\dot{M}}{4\pi r^2{u}^r}·}\end{array}$$(3)

The specific entropy was calculated at the sonic radius from the following formula

$$\begin{array}{c}\hfill \begin{array}{cc}{\displaystyle K={\left(u^{r}4\pi r^2\frac{{\displaystyle c_s^{\frac{2}{\gamma -1}}}}{{\displaystyle {\gamma }^{\frac{1}{\gamma -1}}\dot{M}}}\right)}^{\gamma -1},}\hfill & {\displaystyle c_s^2=\gamma \frac{{\displaystyle P}}{{\displaystyle \rho }},}\hfill \end{array}\end{array}$$(4)

where c_s is the speed of sound.

The angular momentum of the stellar core is scaled with l_isco, representing the angular momentum of a particle orbiting the black hole at the innermost stable circular orbit (Suková & Janiuk 2015; Suková et al. 2017; Palit et al. 2019). In the Kerr metric, the radius of the ISCO, expressed in units of r_g, is given by:

$$\begin{array}{cc}\hfill r_{\mathrm{isco}}& =3+Z_2\mp \sqrt{(3-Z_1)(3+Z_1+2Z_2)},\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill Z_1& \equiv 1+{(1-a^2)}^{1/3}[{(1+a)}^{1/3}+{(1-a)}^{1/3}],\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill Z_2& \equiv \sqrt{3a^2+Z_1^2},\hfill \end{array}$$(7)

where a is the spin of the black hole. For a Kerr black hole, l_isco and the auxiliary energy, ϵ_isco, depend only on the mass and the spin of the black hole (Bardeen et al. 1972), and are given by

$$\begin{array}{cc}\hfill l_{\mathrm{isco}}& =\frac{r_{\mathrm{isco}}^{1/2}-2a/r_{\mathrm{isco}}+a^2/r_{\mathrm{isco}}^{3/2}}{\sqrt{1-3/r_{\mathrm{isco}}+2a/r_{\mathrm{isco}}^{3/2}}},\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill {ϵ}_{\mathrm{isco}}& =\frac{1-2/r_{\mathrm{isco}}+a/r_{\mathrm{isco}}^{3/2}}{\sqrt{1-3/r_{\mathrm{isco}}+2a/r_{\mathrm{isco}}^{3/2}}}·\hfill \end{array}$$(9)

The distribution of the angular velocity u^ϕ is given by

$$\begin{array}{c}\hfill u^{\varphi }(r,\theta )=S(-g^{t\varphi }{ϵ}_{\mathrm{isco}}+g^{\varphi \varphi }{l}_{\mathrm{isco}}){sin}^2\theta ,\end{array}$$(10)

where the g^tϕ and g^ϕϕ are the metric tensor components in the Boyer–Lindquist co-ordinates given by g^tϕ = −2ar/ΣΔ, g^ϕϕ = (Δ − a²sin²θ)/(ΣΔsin²θ), Σ = r² + a²cos²θ, and Δ = r² − 2r + a². The dimensionless parameter S is used to scale the angular momentum with respect to ISCO. For S = 0, the model describes purely spherical accretion. The distribution of the angular momentum was additionally multiplied by sin²θ to ensure the maximum rotational velocity at the equatorial plane (Suková et al. 2017), where an accretion disc should form (MacFadyen & Woosley 1999).

The collapsar models producing relativistic jets should be highly magnetised (Tchekhovskoy & Giannios 2014; Gottlieb et al. 2022). To initialise the magnetic field, we use the plasma parameter β, defined as the ratio of the gas pressure, p = (γ − 1)u, where u is the specific internal energy, to the magnetic pressure, $p_{\mathrm{mag}}=\frac{1}{2}{b}^{\mu }{b}_{\mu }$, where b^μ is the contravariant component of the magnetic four-vector. Our models are scaled to the value of β₀ at 2 r_g. We adopted a uniform magnetic field given by the magnetic potential A_φ, defined as

$$\begin{array}{c}\hfill {\mathbf{A}}_{\phi }\propto (0,0,\frac{1}{2}rsin\theta ).\end{array}$$(11)

The initial mass of the black hole is assumed as M₀ = 3 M_⊙, whereas the mass of the stellar core is M_★ = 25 M_⊙. The initial spin of the black holes is chosen to be a₀ = 0.8. The angular momentum is scaled with S = 2 in all models, similar to Murguia-Berthier et al. (2020). The magnetic field is scaled with β₀ = 1 or β₀ = 10, corresponding approximately to field strengths ranging from B ∼ 10⁸ G in the outer region to B ∼ 10¹⁵ G in the innermost core. This choice was motivated by the magnetic fields of B ∼ 10¹⁵ G observed in magnetars (Turolla et al. 2015; Kaspi & Beloborodov 2017).

In addition, to break the azimuthal symmetry, we introduce the initial random perturbation in the internal energy. The amplitude of the perturbation was set to 5%, as in Mizuta et al. (2018). All models are three-dimensional. The resolution of runs, and their initial parameters are listed in Table 1.

2.2. Self-gravity implementation

Typical GRMHD models of collapsars assume that the central black hole has already formed, and it is surrounded by the accretion disc fuelled via matter fallback from the massive stellar envelope (Mizuno et al. 2004). For the stripped-envelope star, its mass is still several times larger than the initial black hole mass (Heger et al. 2003; Woosley & Heger 2006; Crowther 2007). Therefore, self-gravity cannot be ignored in physically consistent simulations. Although a complete numerical solution of the full Einstein equations in a three-dimensional collapsar scenario is theoretically possible, its computational complexity is very high. In this work, we adopt an alternative technique to model self-gravity, allowing for an efficient run of a set of three-dimensional models. In this way, we aim to investigate the general impact of self-gravity on the magnetised collapsing stars.

Evolution of the mass and spin of the black hole was considered, without self-gravitating force, in Janiuk et al. (2018), Król & Janiuk (2021). The implementation presented in Janiuk et al. (2023), took into account self-gravity as perturbative terms to Kerr metric to simultaneously analyse the black hole mass and spin evolution, as well as self-gravity effects. Our approach is motivated by the Teukolsky equation (Teukolsky 1972). The global vacuum solution is formulated using the CCK formalism (Chrzanowski 1975; Cohen & Kegeles 1974; Wald 1978), enabling us to make the reconstruction of the metric perturbations, following the approach of van de Meent (2017) who proved that the perturbation caused by a particle in any bound orbit around a Kerr black hole affects only two parameters, a and M, measured above the particle’s orbit. Notably, the perturbation disappears inside the orbit.

The modifications of metric mass coefficients in our code are divided into two groups. The first group is responsible for the changes in the central black hole mass, ΔM, due to the mass-energy flux crossing the black hole horizon (primarily mass, but also other forms of energy contained in the radial energy flux $T_t^r$). This flux is integrated over the black hole horizon r_h (note that the size of the event horizon depends on time) at each time step of the simulation using the following formula:

$$\begin{array}{c}\hfill \mathrm{\Delta }M(t)={\int }_0^t{\int }_0^{2\pi }{\int }_0^{\pi }\sqrt{-g}{T}_t^{r}d\theta d\varphi dt.\end{array}$$(12)

The second group comes from the perturbative term δM, which is responsible for the self-gravity and depends on both time and radius. We calculated this perturbation by integrating the energy density, $T_t^t$, over the spherical volume enclosed between the event horizon, r_h, and the given radius, r, as follows:

$$\begin{array}{c}\hfill \delta M(t,r)={\int }_{r_{\mathrm{h}}}^r{\int }_0^{2\pi }{\int }_0^{\pi }\sqrt{-g}{T}_t^{t}d\theta d\varphi dr.\end{array}$$(13)

The evolution of the mass metric coefficient is finally given by the following three terms:

$$\begin{array}{c}\hfill M(t,r)=M_0+\mathrm{\Delta }M(t)+\delta M(t,r).\end{array}$$(14)

Similarly, we introduce the dynamical evolution of the black hole’s angular momentum J(t). The initial angular momentum is J₀, and its evolution ΔJ(t) is governed by the angular momentum flux $T_{\varphi }^r$ crossing the black hole horizon, given by

$$\begin{array}{c}\hfill \mathrm{\Delta }J(t)={\int }_0^t{\int }_0^{2\pi }{\int }_0^{\pi }\sqrt{-g}{T}_{\varphi }^{r}d\theta d\varphi dt.\end{array}$$(15)

The second component is the perturbative term δJ(t.r), which accounts for the total angular momentum density, T_ϕ^t, contained in the volume enclosed between the event horizon, r_h, and the given radius, r. The term was computed as

$$\begin{array}{c}\hfill \delta J(t,r)={\int }_{r_{\mathrm{h}}}^r{\int }_0^{2\pi }{\int }_0^{\pi }\sqrt{-g}{T}_{\varphi }^{t}d\theta d\varphi dr.\end{array}$$(16)

The dynamical evolution of the angular momentum is finally given by

$$\begin{array}{c}\hfill J(t,r)=J_0+\mathrm{\Delta }J(t)+\delta J(t,r).\end{array}$$(17)

To sum up, we have two dynamical metric components in Kerr–Schild co-ordinates, depending on time and distance from the black hole: the mass M(r, t) and the spin parameter a(t, r). Their forms are the same as in Janiuk et al. (2018). The code uses dimensionless units, where c = G = M₀ = 1, so the calculated perturbation to the mass of the black hole needs to be divided by the initial mass of the black hole, M₀, as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}{\displaystyle \xi (t,r)\equiv \frac{M_0+\mathrm{\Delta }M(t)+\delta M(t,r)}{M_0},}& {\displaystyle \alpha (t,r)\equiv \frac{J(t,r)}{\xi (t,r)},}\end{array}\end{array}$$(18)

where ξ and α are auxiliary variables introduced here to simplify metric tensor notation. The six non-zero components of the self-gravitating Kerr–Schild metric are:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill g_{\mathit{tt}}& =-1+\frac{2r\xi }{r^2+{\alpha }^2{cos}^2\theta },\hfill \\ \hfill g_{\mathit{tr}}& =\frac{2r\xi }{r^2+{\alpha }^2{cos}^2\theta },\hfill \\ \hfill g_{t\varphi }& =-\frac{2\alpha r\xi {sin}^2\theta }{r^2+{\alpha }^2{cos}^2\theta },\hfill \\ \hfill g_{\mathit{rr}}& =1+\frac{2r\xi }{r^2+{\alpha }^2{cos}^2\theta },\hfill \\ \hfill g_{r\varphi }& =-\frac{2\alpha r\xi {sin}^2\theta }{r^2+{\alpha }^2{cos}^2\theta },\hfill \\ \hfill g_{\varphi \varphi }& ={sin}^2\theta [(r^2+{\alpha }^2{cos}^2\theta )+\frac{2{\alpha }^{2}r\xi {sin}^2\theta }{r^2+{\alpha }^2{cos}^2\theta }].\hfill \end{array}\end{array}$$(19)

2.3. GRMHD scheme

We used HARM-SELFG, which is our updated version of the general relativistic magnetohydrodynamics code HARM (Gammie et al. 2003; Noble et al. 2006; Sapountzis & Janiuk 2019). The numerical scheme is conservative and shock-capturing, using the HLL method of calculating fluxes (Harten et al. 1983). The code evolves the system of GRMHD equations representing the conservation of rest-mass, energy-momentum, and induction equations, expressed as

$$\begin{array}{c}\hfill \begin{array}{ccc}{\mathrm{\nabla }}_{\mu }\left(\rho u^{\mu }\right)=0,& {\mathrm{\nabla }}_{\mu }\left(T^{\mu \nu }\right)=0,& {\mathrm{\nabla }}_{\mu }(u^{\mu }{b}^{\nu }-u^{\nu }{b}^{\mu })=0,\end{array}\end{array}$$(20)

where u^μ is the contravariant four-velocity, b^μ is the contravariant magnetic four-vector, and ρ is the rest-mass density. The numerical scheme uses the MHD stress-energy tensor, $T_{\mathrm{MHD}}^{\mu \nu }$, which is the sum of the stress-energy tensor for a perfect fluid $T_{\mathrm{fluid}}^{\mu \nu }$ and electromagnetic field $T_{\mathrm{EM}}^{\mu \nu }$, defined as

$$\begin{array}{cc}\hfill T_{\mathrm{fluid}}^{\mu \nu }& =(\rho +u+p)u^{\mu }{u}^{\nu }+p{g}^{\mu \nu },\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill T_{\mathrm{EM}}^{\mu \nu }& =b^2{u}^{\mu }{u}^{\nu }+\frac{1}{2}{b}^2{g}^{\mu \nu }-b^{\mu }{b}^{\nu },\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill T_{\mathrm{MHD}}^{\mu \nu }& =T_{\mathrm{fluid}}^{\mu \nu }+T_{\mathrm{EM}}^{\mu \nu }.\hfill \end{array}$$(23)

The initial conditions are specified using the Boyer–Lindquist co-ordinates (Boyer & Lindquist 1967), and then transformed to the Kerr–Schild co-ordinates. The integration of partial differential equations is done in the modified Kerr–Schild co-ordinates (Noble et al. 2006), with a radial logarithmic grid x^[1], where radius r is replaced by r = e^x[1]. The latitude θ and azimuthal ϕ co-ordinates remain unchanged. The outer radius of the computational grid is set to 1000 r_g, with open boundary conditions, allowing for the free outflow of mass. The boundary conditions in polar and azimuthal directions are reflecting and periodic, respectively.

The final time, t_f, of all simulations was set to 50, 000 t_g. The corresponding conversions from code units to CGS units are as follows: r_g = 4.43 × 10⁵ cm and t_g = 1.48 × 10⁻⁵ s. These units are calculated for M_BH = 3 M_⊙, and are held constant, in contrast to the size of the event horizon, which increases over time. The computational cost of the model with evolving mass and spin of the black hole in high resolution, such as Model-1-NSG is about ≈900 000 CPU hours¹, whereas the model with self-gravity perturbations is about ≈20% more expensive, exceeding one million CPU hours. The scalability is discussed in Janiuk et al. (2025). In general, our implementation of self-gravity is an extra cost of about ≈15 − 25% more CPU hours, which we regard as a reasonable expense for a more physically motivated model of collapsar.

3. Results

3.1. Jet energetics

Astrophysical jets are highly collimated relativistic outflows. The main mechanism responsible for launching Poynting jets is the Blandford–Znajek process (Blandford & Znajek 1977). In these jets, most of the energy is transported via electromagnetic fields rather than the kinetic energy of particles, and the Maxwell stress can even efficiently collimate the outflows provided the jets at launch are narrow (Komissarov 2011). Then, energy is converted into the bulk kinetic energy of particles and emitted as highly energetic radiation (Lyubarsky 2010), observed when the jet axis is pointed towards the observer (Piran 2004). The observed median bulk Lorentz factors of Γ ≈ 300 (Ghirlanda et al. 2018) and values in the range Γ ≈ 130 − 300 (Ravasio et al. 2024), indicate that jets associated with GRBs are highly relativistic. In the previous two-dimensional study of self-gravitating collapsars (Janiuk et al. 2023), no jet was found. In fact, for magnetised collapsars, the parameter space explored in that work was focused on lower black hole spins, slower rotation of the cores and smaller magnetisations. Below, we describe, for the first time, the effects of self-gravity on the properties of jets launched in our three-dimensional simulations.

We measured the jet energetics and variability using a parameter defined as:

$$\begin{array}{c}\hfill \mu ={\mathrm{\Gamma }}_{\infty }=-\frac{T_t^r}{\rho u^r},\end{array}$$(24)

where T_t^r is the component of the energy-momentum tensor, representing the radial energy flux. We used the full stress-energy tensor, which contains both the fluid and electromagnetic parts (Eq. (23)). The quantity μ provides the estimation of the terminal Lorentz factor, Γ_∞, under the assumption that all forms of energy are converted at infinity to the baryon bulk kinetic motion (Vlahakis & Königl 2003; Sapountzis & Janiuk 2019; James et al. 2022).

The method of determining Γ_∞ from our three-dimensional simulations is aimed at finding specific moments in time where jet ejection declines or ceases completely. We measure the terminal Lorentz factor in cones with half-opening angles θ_j of 5°, 10°, and 15°, extending from the horizon radius r_h to a distance of 150 r_g. As a first step, we interpolate the results onto a uniform grid to remove a bias of modified Kerr–Schild co-ordinates, where most of grid cells are located close to the event horizon. Then, we average the 10% most energetic elements in the cone. Our method emphasises contributions from the most energetic elements and is sensitive to physical interruptions in jet ejection to properly probe its variability.

The time evolution of terminal Lorentz factors for all models is shown in Fig. 1. We observe a significant difference between the models with and without self-gravity; however, the initial magnetisation is affecting the whole pattern. First, we note that the choice of the cones half-opening angles θ_j does not affect the emission patterns; however, for smaller angles, maximum spikes are more energetic, which is related to non-uniform jet energy distribution (Tchekhovskoy et al. 2008; Janiuk 2022; Saji et al. 2025). In Model-1-SG and Model-1-NSG, the jets are launched simultaneously, suggesting that self-gravity does not affect the time of jet formation if it is feasible (as supported by a strong magnetic field). Since the jets originate in the vicinity of the black hole, they are initially dominated by the gravitational influence of the central object. The maximum Lorentz factor in both models is of the same order, and in Model-1-SG, it attains the value of Γ_∞ ≈ 190, whereas in Model-1-NSG, the maximum value achieves Γ_∞ ≈ 200.

Fig. 1. Evolution of the terminal Lorentz factor (Γ∞) for Model-1 (left) and Model-2 (right).

In Model-2-SG, the jet ejection does not occur at all, indicating that self-gravity can lead to failed GRBs (in Fig. 1, the terminal Lorentz factor remains Γ_∞ ≈ 1 throughout the evolution). We note that in Model-2-NSG, under identical conditions, but without self-gravity, jet ejection occurs, but it is delayed in comparison to Model-1-NSG, which is possibly a consequence of the weaker magnetic field (Aguilera-Miret et al. 2024). In Model-1-NSG, where β₀ = 1, the jet starts at t ≈ 0.029 s, while in Model-2-NSG at t ≈ 0.158 s, jet formation time is more than 5 times longer in the model with 10 times less magnetic pressure, β₀ = 10. In Model-2-NSG, the terminal Lorentz factor attains at the peak values about Γ_∞ ≈ 500, and then drops to about Γ_∞ ≈ 100. A direct comparison of Γ_∞ values achieved in Model-1 and Model-2 is challenging due to the resolution difference.

For Model-1-SG, and Model-1-NSG, a very sharp initial spike is observed. Then, in the SG case, there is a gradual decline in Γ_∞, until the jet halts completely. We label this break in Fig. 1 as the quiescent interval (the threshold is set to Γ_∞ < 5 for θ_j = 5°) for Model-1-SG. After this break, jet reappears, although it is not as energetic as earlier. At the end of this simulation, the jet in the self-gravitating model ceases, while in the model without self-gravity, it is still present and attains a high Lorentz factor, Γ_∞ = 50 − 150. Our results suggest that in self-gravitating collapsars, a potentially dormant central engine can exist. According to internal shock models (Rees & Meszaros 1994; Kobayashi et al. 1997), it can be responsible for quiescent intervals in the GRB prompt emission. Moreover, we argue that under specific conditions, self-gravity can interrupt the process of jet formation. Different criteria for measuring variability do not qualitatively change the observed variability of the jet, particularly the existence of the quiescent interval, which is caused by the suppression of jet ejection and is independent of the measurement method.

To obtain the electromagnetic power produced by the system, we computed the Blandford–Znajek luminosity (McKinney & Gammie 2004), defined as

$$\begin{array}{c}\hfill L_{\mathrm{BZ}}=-{\int }_0^{2\pi }{\int }_0^{\pi }(b^2{u}^r{u}_t-b^r{b}_t)\sqrt{-g}d\theta d\varphi ,\end{array}$$(25)

and given by the radial energy flux of the electromagnetic stress tensor (Eq. (22)). In practice, we include only the outgoing electromagnetic energy flux, i.e. we integrate only regions where the radial component of the energy flux is directed outwards. The value is integrated over the event horizon r_h, which in our models grows over time due to mass accretion onto the black hole.

The evolution of L_BZ luminosities for all models is shown in Fig. 2. Initially, over the first ∼0.08 s, in two variants of Model-1, we observe the same jet power. Then, about t ≈ 0.1 s in Model-1-SG, there is a small decline. Before the quiescent interval, in Model-1-SG, the power rapidly declines, and during the dormant period, it is about two orders of magnitude weaker than in Model-1-NSG. Afterwards, the luminosity grows again, and finally, about t ≈ 0.4 s, it diminishes. In contrast, the extraction of energy via the Blandford–Znajek process, in Model-1-NSG is more stable over time, keeping the power about L_BZ ≈ (1.57 − 5.27)×10⁵³ erg s⁻¹. The stability of Model-1-NSG is further discussed below.

Fig. 2. Evolution of Blandford–Znajek luminosity for Model-1 (left) and Model-2 (right).

Model-2-NSG shows less stable power than Model-1-NSG, with a peak of L_BZ ≈ 1.36 × 10⁵³ erg s⁻¹. In Model-2-NSG at t ≈ 0.62 s, the Blandford–Znajek power declines, despite still large Lorentz factor, Γ_∞ ∼ 100. This decline is a signature of the gradual jet quenching. In Model-2-SG, we observe the non-zero Blandford–Znajek luminosity, despite the absence of the jet. The pattern resembles that observed in Model-1-SG, with a deep decrease and finally a secondary increase. However, in Model-2-SG, the outflow cannot be called a jet as it is not collimated.

Model-2-SG exhibits the lowest emitted energy E_j, about E_j ≈ 6.54 × 10⁵¹ erg, whereas in Model-2-NSG it is about E_j ≈ 3.07 × 10⁵² erg. The total energy emitted by Model-1-SG is E_j ≈ 3.89 × 10⁵² erg, and for Model-1-NSG, it is about E_j ≈ 2.21 × 10⁵³ erg. We note that the ratio of produced energies, ζ = E_j^SG/E_j^NSG, between the models with and without self-gravity is ζ ≈ 0.176 for Model-1, and for Model-2, it is ζ ≈ 0.213. It is worth noting that the ratio of these emitted energies is close to one fifth, independently of the initial β₀.

We suggest that during the quiescent interval in the central engine, there should be an observable decline in electromagnetic emission from self-gravitating systems. The amount of emitted energy is lower than in the models without self-gravity.

3.2. Jet opening angle

Opening angles of GRB jets can be estimated from their afterglow observations by identifying the achromatic break in the multi-wavelength light curve (Sari et al. 1999; Frail et al. 2001). The analysis of over 100 long GRBs detected by Fermi GBM revealed that the half-opening angles are well described by a log-normal distribution, with approximately 90% of events exhibiting θ_jet < 20° and the median value of about θ_jet ≈ 8° (Goldstein et al. 2016). In general, jet collimation is affected by many physical parameters, such as the black hole spin (Hurtado et al. 2024), magnetic field structure and strength (Fendt 2006), disc wind pressure (Globus & Levinson 2016), and the density of the external medium (Levinson & Begelman 2013).

To measure the opening angles of jets in our simulations, we use two methods. The first one relies on the magnetisation σ, defined as the ratio of energy stored in the magnetic field to the energy stored in the rest-mass density, given by

$$\begin{array}{c}\hfill \sigma =\frac{b^2}{\rho },b^2=b^{\mu }{b}_{\mu }.\end{array}$$(26)

This method is based on identifying the highly magnetised regions where σ > 1. In three-dimensional simulations, to avoid inaccuracies resulting from a potential jet-wobbling, we compute opening angles over N_ϕ cross-sections. In each cross-section, we locate the angular distance from the polar northern axis, where σ > 1. The opening angle θ_jet is then the mean value of all individual half-opening angles measured in all cross-sections. In the second method, similar to the first one, instead of magnetisation, we use criterion Γ_∞ > 5. We compute the opening angle at the radial distance of 100 r_g. Using two measures allows us to analyse the highly magnetised regions independently of the jet launching process to determine the characteristics of the funnel.

In Fig. 3, the top panel shows the evolution of the opening angle for two configurations of Model-1. First, we note that both methods of measuring opening angle give consistent results in Model-1-NSG. Nevertheless, in Model-1-SG, we observe that shortly before the quiescent interval, two measurements start to diverge. This behaviour is related to the decay and subsequent change of sign of the velocity u^r (at the beginning of accretion from the polar regions), which in turn reverses the sign in Eq. (24), as we consider only outflowing elements. The magnetisation arises mostly from the toroidal component B^ϕ of the magnetic field. This indicates that during the break in the jet quenching, while Γ_∞ is very low, the highly magnetised funnel, where previously the jet was ejected, still exists; however, the Poynting flux is unable to transport energy efficiently outwards. Moreover, the opening angle in the model with self-gravity gradually decreases as it approaches the quiescent interval, from θ_jet ≈ 10° to θ_jet ≈ 4°. The ratio of the average gas pressure ⟨P⟩_ϕ at the edge of the jet funnel in both models during the jet launching phase is ⟨P_SG⟩_ϕ/⟨P_NSG⟩_ϕ ≈ 1.27 at t = 0.029 s, increasing to ⟨P_SG⟩_ϕ/⟨P_NSG⟩_ϕ ≈ 7.5 at t = 0.125 s. Then, during the quiescent interval, we observe the opposite relation. Therefore, we find that due to additional pressure from self-gravitating fluid, the opening angle of the jet funnel is smaller.

Fig. 3. Jet opening angle for Model-1 (top), and for Model-1-NSG and Model-2-NSG (bottom), obtained using two methods.

In Fig. 3, the bottom panel shows the evolution of opening angles for non-self-gravitating models, Model-1-NSG and Model-2-NSG. The jet launching is delayed about five times in Model-2, as shown in Sect. 3.1. The chaotic nature of jet is observed in the Model-2-NSG; nevertheless, the opening angle in this case is smaller than in Model-1-NSG, taking the mean value of about ⟨θ_jet⟩_t ≈ 10.2°. This is related to the retarded jet ejection and higher inner pressure, when the jet is formed. In Model-1-NSG, we observe the systematic increase in the opening angle over time, from θ_jet ≈ 11° up to θ_jet ≈ 22°, with the mean value about ⟨θ_jet⟩_t ≈ 17.6°. This gradual increase is related to the declining spin (Hurtado et al. 2024), as discussed below.

To sum up, our results show that self-gravity contributes to jet collimation and keeps the opening angles between θ_jet = 4° −10°, consistently with observations (Goldstein et al. 2016).

3.3. Black hole mass and spin evolution

Our simulations account for the Kerr metric update due to the change in mass and spin of the black hole, allowing us to track the evolution of these parameters and relate them to jet properties. The mass, spin, and accretion rate onto a black hole are shown in Fig. 4, as a function of time.

Fig. 4. Evolution of the black hole accretion rate, black hole spin, and mass for Model-1 (top) and Model-2 (bottom).

The difference between models with and without self-gravity is clearly visible. The mass accretion rate onto a black hole is significantly larger in models with self-gravity, temporarily exceeding Ṁ > 10³ M_⊙ s⁻¹ in Model-1-SG, while in Model-1-NSG the accretion rate only slowly grows over time, reaching at most Ṁ ≈ 8 M_⊙ s⁻¹. The mean value of Ṁ for Model-1-SG and Model-2-SG in the time interval t < 0.4 s is approximately Ṁ ≈ 60 M_⊙ s⁻¹. In the previous study (Janiuk et al. 2023), we observed accretion rates reaching up to Ṁ ≈ 500 M_⊙ s⁻¹ in the self-gravitating models, which is lower than in the current, three-dimensional study. This is explained by the different initial conditions, and the fact that rapid accretion is preceded by mass accumulation, before the rapid collapse. In Model-1-SG, we note that the quiescent interval is prefaced by the rapid mass accretion, and appears when the accretion rate after the peak starts to decrease. On the contrary, Model-1-NSG and Model-2-NSG exhibit less rapid accretion; nevertheless, the rate is increasing over time, more steeply in Model-2-NSG. In both models, Model-2-SG and Model-1-SG, the accretion rate follows a similar pattern, linked to the perturbations imposed by self-gravity, which produce dramatic effects and cannot be reduced by a strong magnetic field. In models without self-gravity, due to various magnetic pressures, the accretion rate in Model-2-NSG is about one order of magnitude larger than in the more magnetised Model-1-NSG.

The spin evolution varies across the models. Initially, in all versions, the evolution proceeds similarly for SG and NSG cases. The spin change is a result of black hole rotational energy extraction, as well as accretion of massive fluid that possesses certain specific angular momentum (Gammie et al. 2004; Sadowski et al. 2011; Shapiro 2005). In Model-1-NSG, the gradual spin decline is related to the Blandford–Znajek process (as outflow of the angular momentum, see Eq. (15)), while in Model-1-SG, the rapid spin decay is not due to extraction, as the Blandford–Znajek luminosity in this model is low (see Fig. 2), but mostly results from accretion and the increasing mass of the black hole. For self-gravitating models, the final spins for Model-1-SG and Model-2-SG are about a ≃ 0.4 (in Model-1-SG slightly lower), whereas in the case of the Model-1-NSG, even though most of the mass is not accreted, the spin is lower, a ≈ 0.35. This result indicates that due to longer and stable jet production, models that do not include self-gravity can extract more rotational energy from the black hole. This trend is visible in both sets of simulations, even though for Model-2-NSG the final spin value was not yet reached. We conclude that self-gravitating collapsars leave black holes with somewhat higher spins, and the equilibrium value is reached earlier than for non-self-gravitating ones, for the same spin at birth.

The same effect is visible in the black hole mass evolution. In the right panels of Fig. 4, we note that black hole masses evolve much more rapidly in self-gravitating models. The quiescent interval appears when the mass of the black hole is almost maximal. In Model-1-SG and Model-2-SG, after the evolution, the final mass of the black hole is M_BH ≃ 28 M_⊙. We note that the accretion rate in Model-1-NSG is relatively low, and after the 50 000 t_g, the mass reaches M_BH ≈ 5 M_⊙, whereas in Model-2-NSG, it is M_BH ≈ 16 M_⊙, showing an immense effect of larger initial magnetisation on the evolution timescale. However, this effect is observed only in the models without self-gravity. The overall pattern of mass and spin evolution is sensitive to the initial magnetisation. The spin in non-self-gravitating highly magnetised collapsars evolves faster than the corresponding weaker magnetised ones. On the other hand, for simulations including self-gravity, the effect of the initial magnetisation on the evolution appears to be significantly weaker than in the corresponding non-self-gravitating models. Self-gravity becomes an important factor determining the system dynamics, increasing the gas pressure and plasma beta, which reduces the role of magnetic field in the system evolution.

3.4. Accretion flow structure and MAD state

Magnetised accretion onto the black hole is mostly driven by the magnetorotational instability (MRI, Balbus & Hawley (1991)). The magnetically arrested disc (MAD) is the model proposed to explain the properties of powerful jets. Early analytic works (e.g. Bisnovatyi-Kogan & Ruzmaikin 1974, 1976), pointed out that a sufficiently strong accumulation of magnetic flux near the event horizon can efficiently suppress accretion. It was noted in GRMHD simulations (Narayan et al. 2003; Proga & Begelman 2003; Igumenshchev et al. 2003) that the effect can become so strong that it almost chokes the accretion flow. In the MAD state, the accumulation of strong magnetic flux near the black hole can suppress the MRI, leading to the dominance of magnetic Rayleigh–Taylor (RTI) instabilities, also known as interchange instabilities (Marshall et al. 2018; White et al. 2019).

The mass inflow strongly interacts with the magnetic barrier at the magnetospheric radius, r_m, where the force of gravity is equal to the magnetic force. The effect can be conveniently quantified under the assumption that all magnetic flux below r_m reaches the black hole horizon (McKinney et al. 2012). However, recent three-dimensional simulations have suggested that in strongly magnetised flows, accretion can be controlled by chaotic magnetic disconnections and local flux eruptions, not by a global MAD-type barrier (Nalewajko et al. 2024).

The dimensionless magnetic flux ϕ_MAD measures the mass-loading of the magnetic field lines and is given by

$$\begin{array}{c}\hfill {\varphi }_{\mathrm{MAD}}=\frac{{\displaystyle \int \int |B^r|\sqrt{-g}d\theta d\varphi }}{\sqrt{\dot{M}{r}_g{(t)}^{2}c}},\end{array}$$(27)

where B^r is the radial component of the magnetic field measured at the horizon, and Ṁ is the mass accretion rate onto the black hole. In this formula, r_g is time-dependent to provide proper length scaling, since the mass of the black hole changes over time. Typically, when ϕ_MAD exceeds a critical value (about ϕ_MAD ≈ 50 in Gaussian units, or about ϕ_MAD ≈ 15 in Lorentz–Heaviside units; the former are used in this study), the magnetic flux near the horizon saturates and suppresses the mass inflow, causing the disc to enter the MAD state (Narayan et al. 2003). In consequence, a force-free magnetosphere forms and the Blandford–Znajek process is activated, so launching of relativistic jets is highly efficient, as demonstrated by Tchekhovskoy et al. (2011). To estimate the jet production efficiency η, we used the following parameter:

$$\begin{array}{c}\hfill \eta =\frac{L_{\mathrm{BZ}}}{\dot{M}{c}^2}·\end{array}$$(28)

The quantities ϕ_MAD and η for Model-1 and Model-2 are shown as a function of time, in Fig. 5. Initially for Model-1, we observe that ϕ_MAD evolves similarly in both models, and in SG and NSG cases, reaching values up to ϕ_MAD ≈ 20. Subsequently, the models start to diverge. In the model with self-gravity, we observe a gradual decline, whereas in the model without self-gravity, this value saturates temporarily at ϕ_MAD ≈ 40, and then slowly decreases over time to ϕ_MAD ≈ 30, still kept within the MAD state. In Model-1-NSG, the jet efficiency is relatively high, reaching up to η ≈ 37% (with the black hole spin at that moment a ≈ 0.65), and then slowly declines to η ≲ 5%. Even though the ϕ_MAD is relatively stable, the efficiency is decreasing due to reduced black hole spin (Tchekhovskoy et al. 2011; Lowell et al. 2023). In Model-1-SG, we note that self-gravity prevents the increase in ϕ_MAD, and reduces its value to about ϕ_MAD ≲ 1, before the peak in accretion rate. This has an impact on the jet production efficiency, which decreases from η ≈ 8% to η ≪ 1%, before the quiescent interval. This suggests that in the self-gravitating collapsars, it is much harder to produce highly efficient jets. Self-gravity acts as the main factor responsible for the increasing mass accretion rate; hence, the magnetic field accumulation in the vicinity of the black hole is not able to choke the inflow effectively. The transition after the quiescent interval, when ϕ_MAD starts to grow rapidly up to ϕ_MAD > 50, is related to the fact that most of the mass from the system is already accreted, and the numerical floor becomes important. The magnetic flux Φ_BH on the event horizon in Model-1-SG and Model-2-SG gradually decreases. In Model-1-SG, when the black hole horizon reaches almost the final size at t = 0.4 s, the flux drops from Φ_BH = 9.3 × 10²⁹ G cm² to Φ_BH = 2.6 × 10²⁹ G cm² at t = 0.7 s. This effect is associated with magnetic reconnection, a process of slowly decaying magnetosphere known as black hole ’balding’ (Bransgrove et al. 2021; Selvi et al. 2024).

Fig. 5. Dimensionless magnetic flux (ϕMAD) and emission efficiency (η) for Model-1 and Model-2.

In Model-2, we observe a similar evolution, as in Model-1. However, the jet efficiency and ϕ_MAD for Model-2-NSG are not as high as in Model-1-NSG. The dimensionless magnetic flux is about ϕ_MAD ≈ 10, and the efficiency oscillates about η ≲ 3%. In contrast to Model-1, in the self-gravity scenario of Model-2, the model is not able to achieve the MAD state. We conclude that, due to self-gravity, a higher accretion rate leads to mass accumulation near the black hole, which in turn suppresses jet formation. The quiescent interval in the models that form a jet can have observational consequences for the jet intermittency. We notice that at late times of the simulation, after t ∼ 0.4 s, most of the mass has been accreted, and Ṁ is very low, so the level of numerical floor can affect the resulting jet efficiency and dimensionless magnetic flux.

In Fig. 6, we present poloidal slices of magnetisation σ, the terminal Lorentz factor Γ_∞, and rest-mass density, ρ, overlaid with magnetic field streamlines taken at the equatorial plane of Model-1-SG. We selected three representative moments in time: (i) when the jet ejection is stable (at t = 0.052 s, with a ≈ 0.79 and M_BH ≈ 3.13 M_⊙); (ii) when the ejection becomes unstable due to episodic mass accretion (at t = 0.111 s, with a ≈ 0.75 and M_BH ≈ 3.5 M_⊙); and (iii) just before the peak accretion rate, when the jet is completely quenched (at t = 0.140 s, with a ≈ 0.74 and M_BH ≈ 6.1 M_⊙). During phase (i), we observe that magnetisation is constant along the jet funnel, and the jet ejection is stable, as shown in the terminal Lorentz factor map. In the density distribution, we do not observe any large mass accumulation. Subsequently, in phase (ii), the launching becomes unstable. We observe that in the southern hemisphere, the jet launching is suppressed, and ejection takes place only in the northern hemisphere. As seen in the L_BZ in Fig. 2 at about t ≈ 0.1 s, as a steep decline in electromagnetic power is a temporary effect lasting about 0.01 s. We note that there is no inherent asymmetry towards the northern hemisphere. This effect is completely random and possibly arises due to the initial perturbation in the internal energy, leading to temporarily greater mass accumulation on one or the other side of the jet. The interesting outcome is that in Model-1-SG, two observers located on the opposite sides of the jet axis would observe different prompt emissions. Finally, during (iii), just before the highest accretion rate peak, we observe that up to 50 r_g the jet is quenched in both hemispheres. On the equatorial plane, the mass density starts to exceed ρ ∼ 10¹² g cm⁻³. It is worth noting that, when rapid accretion begins, the magnetic field streamlines start to be more chaotic. The quiescent interval begins here. This moment is also related to the fastest growth of the size of the event horizon.

Fig. 6. Two-dimensional maps of magnetisation (σ) and terminal Lorentz factor (Γ∞), and equatorial-plane map of rest-mass density (ρ) with overlaid magnetic field lines for Model-1-SG. Columns show t = 0.052, 0.111, and 0.140 s.

In Fig. 7, we show two-dimensional maps of temperature T, magnetisation σ and plasma beta β for Model-1-SG (top row) and Model-1-NSG (bottom row). In the distribution of temperature, we observe that accretion is more chaotic in the model with self-gravity, and temperature achieves slightly lower values, which is related to a larger density of infalling matter. The maps of magnetisation σ show that in Model-1-SG, there are more regions with a low value of magnetisation (σ ≲ 10⁻³), and magnetisation in the vicinity of the black hole is greater in Model-1-NSG. The right panels show the distribution of plasma-β in the poloidal plane co-aligned with the jet axis. The most noticeable effect shown in these figures was discussed in Sect. 3.2, and this is the jet opening angle. In Model-1-NSG, the opening angle is significantly larger, as in Fig. 3. Moreover, the streams of infalling matter with higher values of β are compressed towards the jet axis, and they are more densely packed around the event horizon of the black hole.

Fig. 7. Two-dimensional maps of temperature (T), magnetisation (σ) and plasma beta (β) of Model-1-SG (top row) and Model-1-NSG (bottom row). All snapshots correspond to t = 0.129 s.

To sum up, in our self-gravitating collapsars, the main observable effect is the appearance of the quiescent interval, related to the mass accumulation and choking of the jet funnel, which results in the jet quenching. Moreover, in the models with self-gravity, it is much harder to maintain the stable MAD state.

4. Discussion

In this study, we present the first three-dimensional simulations of jet formation and propagation in collapsing self-gravitating stars, extending our prior studies of self-gravity effects in collapsing stars (Janiuk et al. 2023). We chose initial conditions to simulate the core of a collapsing massive star endowed with low angular momentum. We found that such a self-gravitating star can produce relativistic jets, if supplied with a strong magnetic field with pressure equal to at least 10% of the gas pressure at the core interior. The initial fast rotation of the black hole is slowing down during the collapse, but is not critical to sustain a Blandford–Znajek powered jet. Our implementation of self-gravity is an intermediate step between usual GRMHD simulations of collapsars (Gottlieb et al. 2022, 2024) and future three-dimensional simulations based on the BSSN formulation. Recent axisymmetric dynamical-spacetime GRMHD simulations have demonstrated the important role of dynamical metric evolution in collapsars (Shibata et al. 2024, 2025).

We studied models in the presence and absence of self-gravity under identical initial conditions to obtain a direct probe of jet properties resulting from self-gravitating collapsars. We used here the transonic configuration of a spherical mass distribution around a newly formed black hole, similar to the previously discussed in Janiuk et al. (2018), Król & Janiuk (2021), Janiuk et al. (2023). This distribution simplifies the problem of collapsar initial configuration that would result from a pre-supernova star simulation and stellar evolution models. Notably, smooth initial density profile is used, to avoid imposing density inhomogeneities which appear because of chemical composition changes in the pre-supernova phase (Heger et al. 2000, 2005). These initial inhomogeneities might have an impact on the SG models, because of their interplay with potential hydrodynamical instabilities and disc fragmentation. These effects were earlier discussed in Janiuk et al. (2023) and are further explored in our follow-up paper (Płonka et al., in prep.).

The role of magnetic fields in collapsing stars can be important to halt accretion and limit the black hole spin-up, as discussed previously in Król & Janiuk (2021). In those models, however, only weakly magnetised collapsars and/or slowly rotating black holes were studied. Here we focus on magnetisations high enough to launch bi-polar jets at the cost of the rotation of a spinning black hole. We find the existence of powerful jets in non-self-gravitating collapsars, while the effects of self-gravity are in competition with the gas magnetisation, and have an impact on the MAD state accretion, as already suggested in Janiuk et al. (2023).

We show here that self-gravity can lead to dormant periods in the activity of the central engine. We argue that this occurs due to the mass accumulation in the horizon region and rapid accretion onto the black hole. We observe that in the model without self-gravity, the MAD state remains stable over time, whereas in the model with self-gravity, this state is disrupted. In the self-gravitating case, therefore, we obtain a gradual decay of the jet efficiency, which precedes the quiescent interval. We observe that the opening angle of the jet of self-gravity is significantly smaller, and it gradually decreases until the jet is completely quenched. In particular, the case simulated by Model-1-NSG more efficiently extracts rotational energy from the black hole, resulting in a smaller value of the final black hole spin. We also speculate that when the initial conditions are not able to produce a jet and allow its breakout quickly enough, it potentially can lead to a failed GRB, as in Model-2-SG. This prediction poses a question whether the initial conditions and their parameter space used in standard simulation without including self-gravity can be used to uniquely constrain the favoured scenario for a successful GRB. Our results imply a narrower range of initial conditions describing the pre-collapse star, such as the magnetic field strength, stellar rotation, and black hole spin, that would be capable of producing relativistic jets in collapsars.

The prompt emission in GRBs exhibits variability on three characteristic timescales (Nakar & Piran 2002). The shortest timescales are related to variability in the central engine. The intermediate timescale is a plateau with no activity, and the longest timescale corresponds to the timescale of the burst at jet breakout. The origin of intermediate quiescent times remains unclear (Ramirez-Ruiz & Merloni 2001; Nakar & Piran 2002). Potentially, some of them can be explained by jet-wobbling motion as shown in Gottlieb et al. (2022). The longest timescales are related to the entire time of the engine activity. We argue that self-gravitating collapsars can provide a complementary explanation for intermediate timescales and the appearance of quiescent times. Potentially, the combination of jet-wobbling and quiescent intervals can explain a wide range of long GRB activity patterns.

We highlight that most of our results are not sensitive to the resolution of the simulations, even though Model-2 and Model-1 differ from each other in that respect. The existence of the quiescent interval in Model-1-SG is shaped by the gas dynamics and the resolution does not play a crucial role. The jet opening angles are different in SG and NSG cases because of a different regions of pressure balance. Qualitatively, this result will not change with resolution; however, the quantitative comparison of θ_jet between Model-1-NSG and Model-2-NSG may slightly change for a different number of cells in polar direction (cf. Fig. 3 bottom panel). In this regard, the exact value of time delay for jet launching between these two models may also be affected by resolution; however, the conclusion that self-gravity can stop the process of jet formation remains unchanged.

Finally, our study presents one chosen magnetic field topology, the vertical field. Other possible topologies, such as hybrid fields, have been presented in literature (Mösta et al. 2014; Gottlieb et al. 2022). Because our earlier two-dimensional studies applied dipole fields did not show the appearance of jets, we highlight here that this configuration may not be plausible also in the three-dimensional simulations. However, as shown by Urrutia et al. (2025), the hybrid field resulted in successful jet production. We aim to explore such configurations in our future studies of self-gravitating collapsars.

5. Conclusions

The main findings of this study are as follows:

• Including self-gravity in GRMHD may lead to temporary jet suppression, which could potentially explain some of the quiescent periods observed in the prompt emission of long GRBs.

• The quiescent interval observed in Model-1-SG is caused by the accumulation of mass, and then rapid accretion onto the black hole.

• There is no substantial difference in the jet launching time between models with and without self-gravity, if such ejection is feasible in both models.

• The models can, in the absence of self-gravity, due to more stable jet ejection, extract more rotational energy via the Blandford–Znajek process, leaving black holes with lower spin values.

• The presence of self-gravity is an additional factor that decreases the jet opening angle, which is related to the effect of the extra pressure exerted on the jet funnel.

• In the models including self-gravity, it is much harder to achieve a stable magnetically arrested disc state. Therefore, the jet efficiency is significantly reduced in the models with self-gravity.

• In weaker magnetised case (Model-2-SG), self-gravity disrupts the jet formation, resulting in a failed GRB with no jet ejection. Under the same initial conditions but without self-gravity, jet launching does occur.

• The observation that self-gravity can lead to a failed GRB, suggests that the parameter space of initial conditions capable of producing successful jets is narrower.

• The direct comparison between Model-1-NSG and Model-2-NSG shows that the jet launching time is significantly delayed (about 5 times) in the model with weaker magnetisation.

• Overall, self-gravity reduces the evolutionary timescale of the system, causing the final black hole’s spin and mass to evolve more rapidly. This is consistent with the conclusions of the two-dimensional study (Janiuk et al. 2023).

Acknowledgments

We would like to thank Krzysztof Nalewajko, Gerardo Urrutia Sanchez, and Ore Gottlieb for valuable discussions and insightful comments that helped improve this study. This work has been partially supported by grant No. 2023/50/A/ST9/00527 from the Polish National Science Center. We gratefully acknowledge Polish high-performance computing infrastructure PLGrid (HPC Center: ACK Cyfronet AGH) for providing computer facilities and support within computational grant no. PLG/2025/018232.

References

All Tables

All Figures

	Fig. 1. Evolution of the terminal Lorentz factor (Γ∞) for Model-1 (left) and Model-2 (right).
In the text

	Fig. 2. Evolution of Blandford–Znajek luminosity for Model-1 (left) and Model-2 (right).
In the text

	Fig. 3. Jet opening angle for Model-1 (top), and for Model-1-NSG and Model-2-NSG (bottom), obtained using two methods.
In the text

	Fig. 4. Evolution of the black hole accretion rate, black hole spin, and mass for Model-1 (top) and Model-2 (bottom).
In the text

	Fig. 5. Dimensionless magnetic flux (ϕMAD) and emission efficiency (η) for Model-1 and Model-2.
In the text

	Fig. 6. Two-dimensional maps of magnetisation (σ) and terminal Lorentz factor (Γ∞), and equatorial-plane map of rest-mass density (ρ) with overlaid magnetic field lines for Model-1-SG. Columns show t = 0.052, 0.111, and 0.140 s.
In the text

	Fig. 7. Two-dimensional maps of temperature (T), magnetisation (σ) and plasma beta (β) of Model-1-SG (top row) and Model-1-NSG (bottom row). All snapshots correspond to t = 0.129 s.
In the text