ABSTRACT

The gravitational wave signature of a binary black hole (BBH) merger is dependent on its component mass and spin. If such black holes originate from rapidly rotating progenitors, the large angular momentum reserve in the star could drive a collapsar-like supernova explosion, hence substantially impacting these characteristics of the black holes in the binary. To examine the effect of stellar rotation on the resulting black hole mass and spin, we conduct a one-dimensional general relativistic study of the end phase of the stellar collapse. We find that the resulting black hole mass at times differs significantly from the previously assumed values. We quantify the dependence of the black hole spin magnitude on the hydrodynamics of the accretion flow, providing analytical relations for calculating the mass and spin based on the progenitor’s pre-collapse properties. Depending on the nature of the accretion flow, our findings have implications for the black hole upper mass gap resulting from pair-instability supernovae, the maximum mass of a maximally rotating stellar black hole, and the maximum effective spin of a BBH formed in a tidally locked helium star–black hole binary.

1 INTRODUCTION

The black hole masses resulting from rapidly rotating pre-collapse stars are not well described by the pre-existing prescriptions that estimate the remnant mass produced at core-collapse (e.g. Fryer & Kalogera 2001; Eldridge & Tout 2004; Fryer et al. 2012). The most widely used remnant mass functions (resulting from the supernova explosion) are based on the work of Fryer et al. (2012). They offer two variations of supernova models to calculate the resulting compact remnant mass, namely the Delayed and Rapid explosion mechanisms. Both models rely on convection-enhanced neutrino heating (Herant et al. 1994; see Fryer, Karpov & Livescu 2021 for a recent review) to drive a successful explosion. The primary difference between them lies in the instability that could drive such an explosion. The Delayed mechanism is driven by a standing accretion shock instability (Blondin, Mezzacappa & DeMarino 2003) and can produce an explosion with a delay of up to one second after the bounce, while the Rapid mechanism utilizes the Rayleigh–Taylor instability (e.g. Sharp 1984), and the explosion occurs within 0.25 s. Under a failed explosion, the Rapid mechanism leads to a direct collapse, while in the Delayed mechanism, some form of explosion may still be possible.

On the other hand, for progenitor stars with a substantial reservoir of angular momentum at core collapse, the supernova explosion might not be determined by convection and therefore may lead to remnant masses different from those predicted by the convection paradigm. This is because the angular momentum reservoir acts as an additional means of centrifugal support and the subsequent viscous dynamics offer an efficient mechanism for energy dissipation during the episode of the explosion.

Fryer et al. (2012) assumed that under the scenario of substantial rotation, the occurrence of relativistic jets during core collapse (e.g. Woosley 1993; MacFadyen & Woosley 1999) would likely disrupt most of the star and only |$\approx$||${\rm M}_\odot$| black hole would be produced. However, a number of subsequent studies, although at times working in the paradigm of rapidly rotating progenitor stars, have either continued to use the Fryer et al. (2012) remnant mass function (e.g. Bavera et al. 2020; Riley et al. 2021) or assumed a more or less direct collapse – i.e. a major mass fraction of the pre-core-collapse star ending up in the black hole, with some (typically the envelope) mass lost during the explosion (e.g. Mandel & de Mink 2016; Marchant et al. 2016; Qin et al. 2018; du Buisson et al. 2020). This is primarily due to the lack of a mass function that captures the effect of stellar rotation on the resulting black hole masses.

Another quality of astrophysical black holes is their spin angular momentum, which, together with their mass, fully determines their property (Israel 1967; Carter 1971; Bekenstein 1972). Although one-dimensional (1D) studies exist that take into account the spin evolution of the black hole during collapse of the star (e.g. Batta & Ramirez-Ruiz 2019; Bavera et al. 2020; Fuller & Lu 2022), they either do not account for the effect of wind outflows (e.g. Bavera et al. 2020), a strong magnetic field, and/or relativistic dynamics on the nature of the accretion flow (Batta & Ramirez-Ruiz 2019; Bavera et al. 2020; Fuller & Lu 2022). In this paper, we aim to address these shortcomings.

The collapse of a sufficiently rapidly rotating star should happen via the formation of an equatorial accretion disc. Such a supernova explosion paradigm, consisting of a newly formed spinning black hole penetrated by a strong magnetic field of the surrounding accretion disc, is also known as a collapsar burst (Woosley 1993; MacFadyen & Woosley 1999) and is believed to be primarily responsible for the observed long-duration gamma-ray bursts (GRBs). The dynamics of the accretion flow under a collapsar scenario has been extensively studied in the past (e.g. Gammie, Shapiro & McKinney 2004; Lee & Ramirez-Ruiz 2006; Kumar, Narayan & Johnson 2008; Tchekhovskoy, Narayan & McKinney 2010, 2011; Sądowski et al. 2014; Sądowski & Narayan 2016; Fujibayashi et al. 2020; Murguia-Berthier et al. 2020; Fujibayashi et al. 2023; Gottlieb et al. 2023; Jacquemin-Ide et al. 2024). The disc evolves viscously via magnetohydrodynamics, resulting in mass outflow and a transient energetic electromagnetic outburst. The resulting dynamics could be highly non-linear, and a detailed study would require a full three-dimensional (3D) magnetohydrodynamic simulation (e.g. Sądowski & Narayan 2016). However, conducting such simulations is computationally expensive and time-consuming. Moreover, presently, most stellar models at core collapse are already 1D in nature, and thus, using these models to conduct a 3D study becomes less effective.

For the current work, we use a simpler approach to study the collapse dynamic, which is based on the earlier work of Kumar et al. (2008) – and later numerically implemented in Fuller & Lu (2022), which is a 1D non-relativistic, height-integrated, time-averaged formalism, but takes into account the effect of vertical wind outflows. We extend their work to a general relativistic framework, hence allowing us to more accurately consider the influence on the accretion flow in the inner region of the accretion disc around the newly formed Kerr black hole. We adopt two variable approaches to treat the accretion flow in the near-horizon region, namely the Novikov & Thorne (1973) accretion disc formalism and the magnetically arrested accretion disc formalism (Tchekhovskoy, Narayan & McKinney 2011).

We use this framework in conjunction with the mesa generated stellar models to calculate the resulting black hole mass and spin relations that would be suitable for black holes born from rotating stars. For the latter, we provide fitting functions that can be readily used by future studies (in particular population synthesis) where the pre-collapse nature of the progenitor star is known. We then compare our results with the earlier work of Fryer et al. (2012), highlighting the crucial difference between the two. We also discuss some implications of our work for the expected birth properties of stellar black holes.

The remainder of this paper is organized as follows. In Section 2, we discuss details of the mesa models used in our analysis. Here, we also develop the formalism used for treating the collapsar dynamics. In Section 3, compare the relativistic formalism with the non-relativitic one. Here, we also present the mass and spin functions for black holes born from rotating stars and discuss their domain of validity. Then, in Section 4, we present some implications of our results, followed by a brief conclusion in Section 5.

2 METHOD

2.1 Models for rotating stars

We use mesa (Paxton et al. 2011, 2013, 2015, 2018, 2019) to generate models of variably rotating helium stars in the mass range 5–80 |${\rm M}_{\odot }$| (separated by |$2\,{\rm M}_\odot$|⁠) with their normalized initial angular velocity |$\Omega _{\rm in}~\in ~[0.01-0.92]$| (separated by 0.04), where the normalization is w.r.t. the critical angular velocity. These models are then evolved from the onset of helium burning to the core carbon depletion.1 No substantial change in the relevant structure (or the angular momentum reservoir) of the star is expected past core carbon depletion, as is evident from the output of such evolved mesa models (e.g. see Fuller & Lu 2022).

Our stellar models exhibit a uniform initial rotation profile at various angular velocities. In mesa, rotation is treated in a shellular approximation (Zahn 1992; Meynet & Maeder 1997), implying that isobaric shells have uniform angular velocity. As the models evolve, they lose mass and hence angular momentum. In addition, the angular momentum is redistributed internally, causing the rotation profile to deviate from that of the initial solid body. In mesa, the transfer of angular momentum is modelled in a diffusive approximation. This is accomplished by contribution from ordinary diffusion, meridional circulations (e.g. Eddington 1925; Sweet 1950), and the magnetic field instability-induced Spruit–Tayler dynamo (Spruit 2002). Although the mapping between the initial and final angular velocities is non-linear, we choose a sufficient number of models (900+) such that it results in stars with different levels of final angular momentum reserves.

For wind mass-loss, we follow the Dutch wind mass-loss prescription as a combination of de Jager, Nieuwenhuijzen & van der Hucht (1988); Vink, de Koter & Lamers (2001), and Nugis & Lamers (2000) with the Dutch_scaling_factor = 0.25 and consider only a single metallicity with |$Z = 10^{-4}$| with the solar scaled relative metal abundances adopted from Asplund et al. (2009). Since the results in Section 3 are calculated so that they remain only a function of the final mass and the angular momentum of the collapsing star, the metallicity of the star is not the variable of interest. However, metallicity would play a direct role in determining the amount of mass and angular momentum retained by the star and also its critical rotation rate (Marchant, Podsiadlowski & Mandel 2023). For certain exceptional cases (e.g. those where angular momentum is sourced from the tidal influence of a binary companion), the angular momentum content can remain large even at larger metallicity; therefore, our approach remains valid for such cases as well. We also include rotation-induced mass-loss as |$\dot{M}(\Omega)=\dot{M}(0) (1- \Omega /\Omega _{\text{crit }})^{- \xi }$| with |$\xi = 0.43$| (Langer 1998). As the angular velocity |$\Omega$| approaches its critical value |$\Omega _{\rm crit}$|⁠, this expression begins to diverge. Hence, we set the variable implicit_mdot_boost = 0.1 in mesa to artificially boost the mass-loss rate until the rotation falls below critical.

2.2 Formalism for collapsar dynamics

Our goal is to determine the final mass and spin of the black hole resulting from the collapse of a rotating star of a given mass and angular momentum distribution. For such a purpose, the detailed nature of the collapse and the dynamics of the evolution of the accretion disc might not be necessary. As such, we focus on the effective (i.e. time-averaged) dynamics of the accretion flow, as detailed below. Unless explicitly stated, we use geometric units with |$G, c = 1$| throughout this section.

2.2.1 Background geometry

As an initial condition, we assume that a black hole of some mass promptly forms during the initial collapse of the stellar core. The spacetime region of this black hole near its horizon (i.e. ignoring the effect of the remaining matter, which will be approximately be distributed in a spherically symmetric fashion) can be modelled using the Kerr metric (⁠|$g_{ab}$|⁠) whose line element in the Boyer–Lindquist coordinates |$(x^a:=~t,r, \theta , \phi)$| reads

(1)

where |$a = J/M; \,\, \Delta = r^2-2Mr+a^2$| and |$\tilde{\rho }^2 = r^2+a^2\cos ^2\theta$|⁠. Here, |$M, J$| are the black hole’s mass and angular momentum, respectively. As the remainder of the stellar matter falls inwards, the metric evolves due to a change in the black hole mass and spin.

For simplicity, we assume that at any given time, there is a negligible amount of matter in the equatorial plane such that the effect of the disc mass on the form of equation (1) can be ignored. Additionally, we assume that the outer shells have a nearly spherical configuration such that the spacetime inside the innermost infalling shell can be assumed to be flat (except for the effect due to the black hole). There are two Killing vector fields |$(\xi , \eta)$| associated with this stationary, axisymmetric spacetime namely |$\xi = (1,0,0,0)\partial _t$| and |$\eta = (0,0,0,1)\partial _{\phi }$|⁠. They correspond to the conservation of energy and the axial component of the angular momentum, respectively.

2.2.2 The rotating fluid

The final state of a massive star evolved through mesa constitutes the initial state of the following analysis. Prior to collapse, we assume (in accordance with the mesa model) that the stellar fluid comprises of matter, with only the axial component of angular momentum being non-zero. Let |$u^a$| be the 4-velocity of a fluid element of this star with components (⁠|$u^t, u^r, u^\theta , u^\phi$|⁠) in Boyer–Lindquist coordinates.2 When the core experiences a collapse, the outer shells will initially be in a free fall. During this phase, we model stellar matter as a fluid with the stress-energy tensor.

(2)

where |$g^{ab}$| (or equivalently |$g^{-1}$|⁠) is the ‘inverse’ metric, P is the pressure and |$\rho , \epsilon$| are the rest mass and internal energy density. We further assume that |$\rho \gg (\epsilon + P)$|⁠. This reduces equation (2) to a dust equation of state. Since the pressure approaches the energy density only in the inner region of the disc, our approximation remains valid while the matter freely falls to form the disc.

2.2.3 Initial formation of the disc

Just prior to infall, the shell is located at a large r value with |$u^a = (\gamma , 0,0,\gamma \Omega)^a \partial _a$| such that |$\Vert u^t \Vert \gg \Vert u^\mu \Vert$| and therefore the shell falls from a near-rest configuration. Here, |$\gamma = {\rm d}t/{\rm d} \lambda$|⁠, |$\lambda$| being the fluid element’s proper time. For simplicity and compatibility with the mesa models, we assume that the stellar matter collapses in a shellular fashion with angular velocity |$\Omega = {\rm d}\phi /{\rm d}t$| (measured at |$r\rightarrow \infty$|⁠). The mass contained within such a shell at time t can be written as

(3)

where |${\rm d}V = \sqrt{g_{\Sigma _t}} {\rm d}r {\rm d}\theta {\rm d}\phi$| and the integral is performed over |$\theta , \phi$| only. Here, |$g_{\Sigma _t} = {\rm det}(g_{\mu \nu })$|⁠, |$g_{\mu \nu }$| being the induced metric on the space-like hypersurface |$\Sigma _t$| at some time t calculated by setting |${\rm d}t = 0$| in the full Kerr metric |$g_{ab}$| (e.g. Paschalidis & Stergioulas 2017). Additionally, |$n_b = (- {\rm d}t/\sqrt{- g^{-1} ({\rm d}t, {\rm d}t)}, 0, 0, 0) = (- {\rm d}t/\sqrt{- g^{tt}}, 0, 0, 0),$| is the unit 1-form normal to |$\Sigma _t$| (where the 1-form |${\rm d}t$| is the dual basis associated with the time-like Killing vector field |$\xi = \partial _t$|⁠). The negative sign in front of |${\rm d}t$| is chosen to keep |$n_b$| future-directed. Thus, under the assumption that |$\Vert u^t \Vert \gg \Vert u^\mu \Vert$|⁠, equation (3) yields

(4)

For |$r \gg M$|⁠, one may approximate the above integral as

(5)

and for |$\gamma = 1$|⁠, we recover the non-relativistic variant of |${\rm d}M$|⁠. The angular momentum contained in the shell with mass |${\rm d}M$| can be written as

(6)

which for |$r \gg M$| can be approximated as

(7)

Again, for |$\gamma = 1$|⁠, we recover the non-relativistic variant of |${\rm d}J$| (cf. Batta & Ramirez-Ruiz 2019). Next, the condition |$g_{ab} u^a u^b = -1$| can be used to estimate the corresponding |$\gamma$| in equations (4) and (6) resulting in two solutions for |$\gamma$| with the positive one being

(8)

In our numerical set-up, we set |$\theta = \pi /4$| in the above equation, owing to the 1D nature of our approach.

Prior to disc formation, we assume that all mass and its associated angular momentum freely falls into the black hole, respecting the condition |$a \lt M$| (see Section 2.3.1). We are interested in finding the stellar radius |$r_0$| at which the infalling shell has sufficient angular momentum to first form an accretion disc. This occurs when the (averaged) specific angular momentum |$\ell$| of the collapsing shell in the equatorial plane equals the specific angular momentum of the innermost stable circular orbit (ISCO) |$\ell _{\rm isco}$| of the Kerr metric. At this point, the matter first circularizes at the ISCO before plunging into the black hole. Using equation (4) and (6), we can calculate the specific angular momentum in the mass shell as

(9)

Additionally, the angular momentum per unit rest mass at radius r of a circular orbit in the equatorial plane of the Kerr metric, assuming prograde orbit is (Bardeen, Press & Teukolsky 1972)

(10)

Choosing |$r = r_{\rm isco}$| in equation (10), the condition

(11)

then gives us the stellar radius |$r_0$| at which the infalling shell has sufficient angular momentum to first form an accretion disc. Since equation (10) is defined per unit rest mass, we accordingly change equation (9) to bring it to a similar form to assist in the comparison in equation (11). The coordinate time interval |$\Delta t~=t(r_0)-t(r_{\rm isco})$| meanwhile gives us the disc formation time post-core-collapse (within an error margin of a sound crossing time as the information of the collapse would take this much time to reach |$r_0$|⁠). As we discuss in Appendix  A, |$\Delta t$| can be well approximated by the corresponding Newtonian free fall time. Therefore, to simplify our calculations, we use the latter approach for calculating |$\Delta t$|⁠.

2.2.4 Evolution of the disc

Once equation (11) is satisfied, the subsequent infalling matter would first be assimilated into the disc before being accreted by the black hole or blown away in the wind. For the accretion phase, we restrict ourselves to the equatorial plane of the black hole, where the accretion disc will reside. The line element in the equatorial plane |$(\theta = \pi /2)$| reads

(12)

where |$A \equiv r^4+r^2 a^2+2 M r a^2$|⁠. At a given time t, let the disc have mass |$M_d$| and angular momentum |$J_d$|⁠. The corresponding effective disc radius |$r_d$| can be calculated by equating

(13)

where the term on the LHS is given in equation (10). The mass and angular momentum evolution of the disc can be written as

(14)

where the dot represents derivative w.r.t. t. Here, |$\dot{M}_{\rm fb}, \dot{J}_{\rm fb}$| represent the mass and the angular momentum assimilation rate into the disc due to fallback and can be calculated using equations (4) and (6). Meanwhile, the mass and angular momentum dissipation rate from the disc (⁠|$\dot{M}_{\rm lost}, \dot{J}_{\rm lost}$|⁠) will be calculated below.

The infalling matter in the disc is accreted on a viscous time-scale. During this process, to account for the effect of wind outflow, we model the local accretion rate as a power law (Blandford & Begelman 1999; Kumar et al. 2008)

(15)

where |$s \in [0.3,0.8]$| is a free parameter (e.g. Kumar et al. 2008; Yuan & Narayan 2014) which we assume to be |$s = 0.5$| (e.g. Kumar et al. 2008). We interpret equation (15) to model the accretion rate as measured from a locally non-rotating frame of reference (⁠|$r, \theta =$| const.). Such observers would have zero angular momentum (apart from the effect of frame dragging) and are also known as zero angular momentum observers (ZAMO). Here, |$\tau _{\rm vis}$| is the viscous time-scale (as measured by ZAMO) at the characteristic radius |$r_d$| and is modelled using the Shakura & Sunyaev (1973) |$\alpha$|-viscosity prescription as (Abramowicz et al. 1996; Gammie & Popham 1998)

(16a)
(16b)

Above |$\Omega _K, \omega$| is the angular velocity of the fluid element (in a circular orbit) and the angular velocity of the frame dragging, respectively, as measured by an observer at radial infinity. In addition, |$\tilde{\Omega }$| is the angular velocity of the fluid as measured by a ZAMO (Bardeen et al. 1972). It is expected that |$\alpha \in [0.01, 0.1]$| with theory supporting the lower value and observation the larger one (King, Pringle & Livio 2007). Here, we take |$\alpha = 0.01$|⁠. A variation in the value of |$\alpha$| would influence the maximal value of |$a^{*} := a/M$| (Sądowski et al. 2011) and therefore may affect the evolution of the black hole (see Sections 2.3.1 and 3.5). To transform equation (15) into the Boyer-Lindquist frame, we convert the ZAMO proper time to the coordinate time as (e.g. see section 10.7 in Gourgoulhon 2018)

(17)

In equation (15), for |$r \le r_t$| – where |$r_t$| is some transitional radius below which the disc undergoes efficient neutrino cooling, resulting in no wind mass-loss and hence efficient accretion – we set |$s = 0$|⁠. We calculate |$r_t$| by setting the RHS of equation (15) to |${ 10^{-2.5} r}/{(2M)} \, M_\odot {\rm s}^{-1}$| (e.g. Kohri, Narayan & Piran 2005; Chen & Beloborodov 2007; Kumar et al. 2008). When |$r \gt r_t$|⁠, the disc enters the advection-dominated accretion flow (ADAF) regime (Narayan & Yi 1994; Narayan & McClintock 2008), where most of the disc mass is blown away in winds. The net mass-loss rate from the disc – owing to outflows and black hole accretion – yields

(18)

where |$\dot{M}_{\rm wind}$| can be evaluated as

(19)

We resort to numerical integration to further evaluate the above equation. Meanwhile, the net angular momentum loss rate yields

(20)

where |$\dot{J}_{\rm BH}$| in defined equation (25) and |$\dot{J}_{\rm wind}$|⁠, represents the angular momentum lost in the wind outflow, assuming that the latter removes the corresponding specific angular momentum from its location of origin. This integral can be evaluated as

(21)

To calculate the last term, we dropped terms |$\mathcal {O}(a/r^2)$| within the integral. This is more accurate than assuming a Schwarzschild geometry with the dropped terms becoming increasingly less unimportant as a decreases or r increases. For the case where |$a = 0$|⁠, equation (21) reduces to the one given in Kumar et al. (2008). Moreover, for the current work, we find equation (21) to also provide a fair match to the numerically integrated value of the |$\dot{J}_{\rm wind}$| term in equation (20).

Finally, to incorporate the influence of the antiparallel relativistic jets on the accretion flow, we assume that at disc formation, the polar region making a half-opening angle of

(22)

with the rotational axis evolves into a relativistic cocoon and is promptly removed from the star. Here |$L_j, r_*, \rho _*$| is the jets’ luminosity, star’s radius, and its mean density, respectively. Additionally, |$\tilde{\eta }$| is the ratio of the speed of light to the jets’ head speed and we assume |$\tilde{\eta } = 10/3$|⁠, e.g. Nakar & Piran 2017 – see Appendix  B for more information on estimating |$\theta _c$|⁠.

2.3 Evolution of black hole with a Novikov–Thorne disc

In the early phase, while the disc is yet to form, we assume that the mass freely falls into the black hole without losing energy or angular momentum. After an initial disc formation, the subsequent disc dynamics is modelled as discussed in Section 2.2.4. This is analogous to modelling the accretion flow as consisting of a relativistic, radiatively efficient, geometrically thin accretion disc that lacks the presence of large-scale magnetic field interactions with the black hole, also known as Novikov & Thorne (1973) disc (NTD).

For such an accretion flow, once the mass reaches |$r_{\rm isco}$|⁠, it is assumed to be directly accreted by the black hole, adding mass at a rate

(23)

where (also note that |$\dot{M}_{\rm acc}$| remains unchanged for |$r \le r_t$|⁠)

(24)

is the energy per unit rest mass of a circular orbit at the radius |$r = r_{\rm isco}$| (Bardeen et al. 1972). This also adds angular momentum to the black hole at a rate3

(25)

where, following Bardeen et al. (1972), |$r_{\rm isco}$| is

(26)

2.3.1 Equilibrium spin of the black hole

In principle, as the black hole accretes mass (and hence angular momentum), its spin parameter will gradually approach the limit |$a = M$|⁠. However, on considering the decelerating impact of disc-emitted photons on the black hole, Thorne (1974) found that the dimensionless spin parameter of the black hole |$a^{*}:= a/M$| should converge to a value of |$a^{*} = 0.9978$|⁠. Since both equations (9) and (24) are functions of a, hence, the equilibrium value of |$a^{*}$| becomes important.4 More recently, Sądowski et al. (2011) recalculated this limit for the case of (advection dominated, optically and geometrical thick) slim accretion discs and found |$a^{*} = 0.9994$| for |$\alpha = 0.01$| and a |$10 \, \times$| super-Eddington mass accretion rate. Since these conditions also manifest in collapsar physics, in this work, we use the latter equilibrium value for |$a^{*}$| when modeling the accretion flow using an NTD. This means that we cap the maximum black hole spin to |$a^{*} = 0.9994$| in our simulation.

2.4 Evolution of black hole with a magnetically arrested disc

Magnetic fields are considered to be crucial in extracting the rotational energy of the Kerr black hole to power the anti-parallel relativistic jets. Of particular interest are the predictions of general relativistic magnetohydrodynamic (GRMHD) simulations of a magnetically arrested disc (MAD) state (e.g. Tchekhovskoy et al. 2011). Here, the accretion disc is threaded by a strong vertical magnetic field that is dragged inwards by the accretion flow, resulting in an accumulation of magnetic flux near the black hole horizon. These magnetic fields can play a central role in extracting the rotational energy of the black hole (Penrose & Floyd 1971; Blandford & Znajek 1977).

To calculate the effect of MAD on the black hole, in conjunction with Section 2.2.4, here we follow the approach of Moderski & Sikora (1996), which accounts for the energy and angular momentum loss of the black hole due to the Blandford–Znajek mechanism as

(27a)
(27b)

where |$\Omega _{\rm H} = r_{\rm H}/(M_{\rm BH}a)$| is the angular frequency of the black hole’s event horizon (measured at |$r \rightarrow \infty$|⁠) and |$r_{\rm H} = M + \sqrt{M^2-a^2}$|⁠. For the MAD state, the variables in equation (27) need to be determined from GRMHD simulations, for which we use the recent analytical fits provided in the 3D GRMHD study of Lowell et al. (2024) as

(28a)
(28b)
(28c)

Here, |$P_{\rm EM}$| is the electromagnetic flux leaving the black hole, which is primarily contained in the jets. Also, |$\Omega _{\rm F}$| is the angular frequency of the magnetic field lines at the horizon. The subscript ‘in’ corresponds to the inner end of the accretion disc in the simulations of Lowell et al. (2024), which is not necessarily at the ISCO.

2.5 Electromagnetic energy of explosion

The energy radiated during a GRB may arise either due to the viscous dynamics of the accretion flow, e.g. via the Blandford & Payne (1982) mechanism (as is the case for an NTD) or a combination of the former plus the Blandford–Znajek mechanism, the latter robbing the black hole of its energy and angular momentum (e.g. for the case of a MAD). This energy would have to climb the gravitational potential well of the black hole in order to be observed by a distant observer.

2.5.1 Gravitational redshift of the emitted energy

We account for the energy lost due to the gravitational redshift of the radiation when measured by an observer at radial infinity as follows. For a stationary observer at radius r, its 4-velocity is|$u^a = \sqrt{-1/g_{tt}} \cdot \xi ^a$|⁠. The specific energy of a radially moving photon as measured by this observer at r is

(29)

Now, |$g_{ab} u^a \xi ^b$| should remain conserved along the photon’s geodesic, given |$\xi$| is Killing. Thus, the observed specific energy at |$r \rightarrow \infty$| is given by |$E_\infty = - g_{ab} u^a \xi ^b$|⁠. Consequently, at a given time, the specific energy of radiation at infinity is related to its value at some radius r as

(30)

where |$g_{tt}$| is calculated at the radius r and |$E(r) = 1-u_t(r)$| is the total specific energy emitted by the fluid element during its decent to radius r from infinity. For computational ease, we evaluate the |$g_{tt}$| term in equation (30) assuming |$r= 2r_{\rm acc}$|⁠, where |$r_{\rm acc}$| is the radius at which the element is accreted (e.g. |$r_{\rm acc} = r_{\rm isco}$| for an NTD). This is a fair approximation as roughly half of the energy would be emitted on each side of this value. Meanwhile, if the energy is sourced from the spin of the black hole (e.g. in the MAD state), one may set |$r = r_{\rm ergo}$|⁠, where |$r_{\rm ergo} = M + \sqrt{M^2-a^2 \cos ^2{\theta }}$| is the outer radius of the ergosphere of the black hole, and assume |$\theta = \pi /2$|⁠. However, under the assumption that the photons move radially outward from the ergosphere (as we did in equation 30), the latter becomes a surface of infinite redshift. Since actual photons would contain angular momentum and hence be able to escape the ergosphere to infinity; thus, to avoid the issue of infinite redshift, for our current purpose, we somewhat ad hocly set |$r=1.01 \times r_{\rm ergo}$| in equation (30) when the energy is sourced from the spin of the black hole.

2.5.2 Isotropic equivalent energy of explosion

We classify the collapsar explosion as a potential source of GRB if the total isotropic equivalent energy |$E_{\rm iso}$| received at radial infinity is greater than |$10^{51}$| erg (Perley et al. 2016). This implies

(31a)
(31b)

where the integral is performed over the total time interval of the explosion, and the subscript of |$E_{\rm iso}$| represents the state of the accretion flow. Additionally, |$\beta = 0.01$| is the fraction of energy that is released in the form of radiation (with the rest lost in neutrinos). In addition, we set |$\beta = 1$| when the radiation energy is sourced from the spin of the black hole. The term |$f:= \hat{\Omega }/4\pi$| represents the fractional probability that the radiation (mainly emitted in the polar jets) aligns with the observer’s line of sight, where |$\hat{\Omega }$| is the total solid angle spanned by the radiation and can be calculated as

(32)

where |$\theta$| is radiation’s half-opening angle w.r.t. the black hole’s rotational axis. Following the calibration of Bavera et al. (2022), we take |$f = 0.05$|⁠.

3 RESULT

By numerically implementing the formalism discussed in Section 2, here, we assess the impact of the progenitor star’s mass and angular momentum reservoir on the properties of the subsequently produced black hole, such as its mass and spin. For the latter, we provide fitting functions that can be readily used by future studies (in particular population synthesis) where the pre-collapse nature of the progenitor star is known. But before that, we first present a brief comparison of the output of the current 1D relativistic implementation to the non-relativistic one.

3.1 Comparison to the non-relativistic approach

A non-relativistic variant of the NTD formalism described in Section 2 is presented in Kumar et al. (2008) and then numerically implemented in Fuller & Lu (2022). To perform a suitable comparison with the above-discussed relativistic treatment, we also modify the non-relativistic formalism to match the requirement set here for the NTD (see Section 2.3.1) and the MAD state (see Section 2.4). To compare the two approaches, we choose a representative stellar model with an initial mass |$35\,{\rm M}_\odot$| that rotates at a range of initial angular velocities |$\Omega _{\rm in}$|⁠. The corresponding angular momentum profile of these models at the moment of collapse is illustrated in Fig. 1.

The angular momentum profile of a helium star evolved from the onset of core-helium burning with an initial mass of $35\,{\rm M}_\odot$ and a range of initial angular velocities (normalized by the critical angular velocity) $\Omega _{\rm in}$. The purple curve with $\Omega = 1$ represents the profile of a star that is rotating at its critical/Keplerian rate (at the point of explosion) and has been plotted for comparison with other models. We note that this does not consider the effect of radiation pressure on the stellar layers and thus overestimates the true value of $\Omega = 1$. $\ell _{\rm isco}$ represents the specific angular momentum at the ISCO if the entire inner mass collapses to form a black hole (respecting $a^{*} \le 0.9994$). Note that models with large $\Omega _{\rm in}$ lose a much larger amount of angular momentum in stellar winds such that their $\ell _{\rm exp}$ converges to a value that is more than an order of magnitude smaller than that of $\Omega = 1$.
Figure 1.

The angular momentum profile of a helium star evolved from the onset of core-helium burning with an initial mass of |$35\,{\rm M}_\odot$| and a range of initial angular velocities (normalized by the critical angular velocity) |$\Omega _{\rm in}$|⁠. The purple curve with |$\Omega = 1$| represents the profile of a star that is rotating at its critical/Keplerian rate (at the point of explosion) and has been plotted for comparison with other models. We note that this does not consider the effect of radiation pressure on the stellar layers and thus overestimates the true value of |$\Omega = 1$|⁠. |$\ell _{\rm isco}$| represents the specific angular momentum at the ISCO if the entire inner mass collapses to form a black hole (respecting |$a^{*} \le 0.9994$|⁠). Note that models with large |$\Omega _{\rm in}$| lose a much larger amount of angular momentum in stellar winds such that their |$\ell _{\rm exp}$| converges to a value that is more than an order of magnitude smaller than that of |$\Omega = 1$|⁠.

As shown in Fig. 2, during the initial stage of collapse, when near-horizon physics is important, the two treatments (i.e. relativistic and non-relativistic) may visibly differ from each other. We find this difference to be evident if the black hole has a lower value of |$a^{*}$| at the onset of disc formation (e.g. Fig. 2a with the spin evolution later shown in Fig. 3). On the other hand, for |$a^{*} \approx 1$| (Fig. 2b), both treatments lead to similar predictions for the location of |$r_t$| and |$r_d$|⁠. This is because a progenitor that can produce a black hole with |$a^{*} \approx 1$| at the onset of disc formation contains a considerable amount of angular momentum in its infalling mass shells. This causes |$r_d$| to assume a large value and thus be located at a relatively far location from the black hole (e.g. see Fig. 2b). For such scenarios, the difference between the relativistic and non-relativistic treatments only becomes apparent if we artificially decrease the value of the transitional radius |$r_t$| (and thus also the value of |$r_d$| as |$r_t$| influences the latter), bringing it closer to the horizon. This is demonstrated in Fig. 2(c) where on decreasing |$r_t$| from its actual location by a factor of 5 (for both the treatments), both |$r_t$| and |$r_d$| begin to differ substantially early on when the disc is near the horizon in the two treatments for the NTD. On the other hand, during later stages, the disc circularizes at a relatively larger distance from the horizon, thus diminishing the effect of the Kerr geometry. Meanwhile, for the MAD state, the inner disc radius is set by equation (28)b, which (for the current model) result in a relatively large inner disc radius. Since |$r_t$| cannot be made smaller than the inner disc radius, in the MAD state, the lowest possible value of |$r_t$| in Fig. 2(c) occurs when |$r_t = r_{\rm in}$| for both the (relativistic and non-relativistic) treatments. As |$r_{\rm in}$| lies relatively far away from the black hole thus, the disc properties for the relativistic and non-relativistic treatments of the MAD state do not show much difference.

The evolution of $r_d$ post-disc formation. The blue (orange) curve represents a relativistic (non-relativistic) model. The solid (dashed) curve represents the NTD (MAD) state. The dotted line represents the transitional radius $r_t$ for both the models. Time axis shows the coordinate (proper) time when a relativistic (non-relativistic) approach is employed. Panel (a) shows the $r_d$ evolution of a star with $M_{\rm exp} = 34.2\,{\rm M}_\odot$ and $\Omega _{\rm in} = 0.25$, and panels (b) and (c) consider $M_{\rm exp} = 33.4\,{\rm M}_\odot$ and $\Omega _{\rm in} = 0.45$. Both stars have same initial mass of $35M_\odot$ (see Fig. 1). Since the blue and orange curves do not visually differ in panel (b), in panel (c), we artificially move $r_t$ inwards by a factor of 5 from its actual location to demonstrate the effect of Kerr geometry on the accretion flow.
Figure 2.

The evolution of |$r_d$| post-disc formation. The blue (orange) curve represents a relativistic (non-relativistic) model. The solid (dashed) curve represents the NTD (MAD) state. The dotted line represents the transitional radius |$r_t$| for both the models. Time axis shows the coordinate (proper) time when a relativistic (non-relativistic) approach is employed. Panel (a) shows the |$r_d$| evolution of a star with |$M_{\rm exp} = 34.2\,{\rm M}_\odot$| and |$\Omega _{\rm in} = 0.25$|⁠, and panels (b) and (c) consider |$M_{\rm exp} = 33.4\,{\rm M}_\odot$| and |$\Omega _{\rm in} = 0.45$|⁠. Both stars have same initial mass of |$35M_\odot$| (see Fig. 1). Since the blue and orange curves do not visually differ in panel (b), in panel (c), we artificially move |$r_t$| inwards by a factor of 5 from its actual location to demonstrate the effect of Kerr geometry on the accretion flow.

Time evolution of the relevant parameters since the onset of disc formation. LHS figure shows the model with $M_{\rm exp} \approx 34.2M_\odot$ and $\Omega _{\rm in} = 0.25$. The RHS figure is the same as LHS but for $M_{\rm exp} = 33.4 {\rm M}_\odot$ and $\Omega _{\rm in} = 0.45$. The legend has the same meaning as in Fig. 2. Here, $M_{\rm BH}, a^{*}$ represent the mass and spin of the black hole, $\dot{M}_{\rm wind}, \dot{J}_{\rm wind}$ the rate of loss of mass and angular momentum in disc winds, $M_d$ the disc’s mass, $r_t$ the transitional radius below which no wind mass-loss occurs, and $r_d$ the characteristic disc radius. Similar to Fig. 2, for the RHS figure, we have artificially shifted $r_t$ inwards by a factor of 5 (for both relativistic and non-relativistic approaches). Note the change in the range of the x-axis for the RHS figures.
Figure 3.

Time evolution of the relevant parameters since the onset of disc formation. LHS figure shows the model with |$M_{\rm exp} \approx 34.2M_\odot$| and |$\Omega _{\rm in} = 0.25$|⁠. The RHS figure is the same as LHS but for |$M_{\rm exp} = 33.4 {\rm M}_\odot$| and |$\Omega _{\rm in} = 0.45$|⁠. The legend has the same meaning as in Fig. 2. Here, |$M_{\rm BH}, a^{*}$| represent the mass and spin of the black hole, |$\dot{M}_{\rm wind}, \dot{J}_{\rm wind}$| the rate of loss of mass and angular momentum in disc winds, |$M_d$| the disc’s mass, |$r_t$| the transitional radius below which no wind mass-loss occurs, and |$r_d$| the characteristic disc radius. Similar to Fig. 2, for the RHS figure, we have artificially shifted |$r_t$| inwards by a factor of 5 (for both relativistic and non-relativistic approaches). Note the change in the range of the x-axis for the RHS figures.

The evolution of the most relevant variables associated with the model discussed in Fig. 2 is presented in Fig. 3. We note that the time axis represents the coordinate (proper) time when a relativistic (non-relativistic) approach is used. Thus, for the variables that involve a time derivative, a meaningful comparison between the two treatments is best achieved at a relatively large distance from the black hole (e.g. see Fig. A1). For the current purpose, this is the case for most part of the collapse. One can check that the tiny but noticeable difference in the black hole properties under the relativistic and non-relativistic approaches in Fig. 3 can mainly be attributed to the relative difference in the location of |$r_t$| for the two treatments. Meanwhile, the relative amplification in the value of disc mass |$M_d$|⁠, and the rate of mass and angular momentum loss in winds |$\dot{M}_{\rm wind}$| and |$\dot{J}_{\rm wind}$| (for the relativistic model) can mainly be attributed to the presence of the |$\gamma$| term in equations (4) and (6).

3.2 Masses of black holes born from rotating progenitors

In this section, we explore the impact of the progenitor star’s mass and angular momentum reservoir on the mass of the subsequently formed black hole.

3.2.1 Masses of black holes as a function of |$M_{\rm exp}$|

To this end, Figs 4 and 5 show the final masses of black holes as a function of the pre-explosion mass |$M_{\rm exp}$| of their progenitor mesa models. For comparison, we use the Fryer Delayed remnant mass function as a source of reference. Additionally, we also overplot the corresponding black hole masses resulting from the non-relativistic approach in Fig. 5.

The masses of black holes as a function of the final mass of their progenitor star (evolved till core carbon depletion). The blue points show the black hole masses resulting from the Fryer Delayed remnant mass function. The label $a^{*}$ represents the dimensionless spin of the black hole. Panel (a) is for the NTD state, while panel (b) is for the MAD state. In both cases (depending on the angular momentum reservoir of the progenitor star), for $M_{\rm exp} \gtrsim 15 \,{\rm M}_\odot$, the black hole mass function experiences a bifurcation into two separate branches. Note that the resulting green branches in panels (a) and (b) show opposite trends in their values of $a^{*}$. The green points (i.e. those with $a^{*} \lesssim 0.12$) on the upper branch in panel (b) are those systems that had low angular momentum already at formation and not because of a GRB-induced spin down.
Figure 4.

The masses of black holes as a function of the final mass of their progenitor star (evolved till core carbon depletion). The blue points show the black hole masses resulting from the Fryer Delayed remnant mass function. The label |$a^{*}$| represents the dimensionless spin of the black hole. Panel (a) is for the NTD state, while panel (b) is for the MAD state. In both cases (depending on the angular momentum reservoir of the progenitor star), for |$M_{\rm exp} \gtrsim 15 \,{\rm M}_\odot$|⁠, the black hole mass function experiences a bifurcation into two separate branches. Note that the resulting green branches in panels (a) and (b) show opposite trends in their values of |$a^{*}$|⁠. The green points (i.e. those with |$a^{*} \lesssim 0.12$|⁠) on the upper branch in panel (b) are those systems that had low angular momentum already at formation and not because of a GRB-induced spin down.

Same as Fig. 4 but now the y-axis shows the ratio of the black hole’s mass to the pre-explosion mass of its progenitor star and the legend shows the isotropic-equivalent energy $E_{\rm iso}$ of the GRB. Additionally, for comparison, we have overplotted the corresponding $M_{\rm BH}/M_{\rm exp}$ values resulting from the non-relativistic approach in lighter colour data points. These are not visible in the LHS figure but can be seen in the RHS one. The grey points are those systems that we do not classify as a potential source of GRB.
Figure 5.

Same as Fig. 4 but now the y-axis shows the ratio of the black hole’s mass to the pre-explosion mass of its progenitor star and the legend shows the isotropic-equivalent energy |$E_{\rm iso}$| of the GRB. Additionally, for comparison, we have overplotted the corresponding |$M_{\rm BH}/M_{\rm exp}$| values resulting from the non-relativistic approach in lighter colour data points. These are not visible in the LHS figure but can be seen in the RHS one. The grey points are those systems that we do not classify as a potential source of GRB.

Figs 4 and 5 show that there is a threshold value of |$a^{*}$| and |$E_{\rm iso}$|⁠, respectively, beyond which the remnant mass trend bifurcates into two separate branches. The magnitude of such an |$a^{*}$| or |$E_{\rm iso}$| and the amount of bifurcation in the resulting black hole masses is a function of the hydrodynamical nature of the accretion flow. In the case of an NTD, the resulting drop in the mass of black holes can at times be larger than 50 per cent (see Fig. 5a). Meanwhile, the MAD state in Fig. 5(b) tends to produce black holes with relatively larger masses compared to the Novikov–Thorne flow. This can be attributed to the difference in the location of |$r_t$| for the two accretion flows. For example, Figs 2(a) and (b) show that once the outer disc transitions to an ADAF (as also discussed in Section 3.1), the location of |$r_t$| for a MAD state is relatively further from the horizon (compared to the NTD state). Thus, in the latter case, the black hole accretes more mass, resulting in larger final masses. Nevertheless, the mass of the black holes is still lower than those resulting from the convection-enhanced neutrino heating explosion, as showcased by the Fryer delayed mechanism.

For stars with relatively smaller reservoirs of angular momentum, the collapse is semi-direct, without much dependence on the nature of the accretion flow. For such stars, we find the masses of the black holes to be typically larger than those predicted by the Fryer delayed mechanism. Although for such slow rotators, one might expect convection-enhanced neutrino heating to additionally contribute towards explosion energetics. However, the scarcity of observations of supernova explosions of stars with mass |$\gtrsim 18 \, {\rm M}_{\odot }$| suggests that direct collapse might be possible (Smartt 2015; Adams et al. 2017). Numerical studies (e.g. Murguia-Berthier et al. 2020) also align with observations, suggesting that massive stars might experience a direct collapse without an ensuing electromagnetic counterpart.

We also observe that, unlike the case of NTD, for the MAD state, the threshold |$E_{\rm iso}$| does not guarantee a downward shift in the black hole masses. For example, in Fig. 5(b), we find very luminous GRBs with |$E_{\rm iso} \gtrsim 5 \times 10^{53}$| erg that also produce relatively massive black holes. Such explosions result from progenitor stars that only develop a black hole accretion disc during the later stage of their collapse (e.g. the scattered green points in Fig. 6), hence minimizing their chances of losing substantial mass in ADAF winds. None the less, the black hole is still able to produce energetic GRB jets as the latter efficiently feeds on its spin angular momentum. Meanwhile, the black hole does not significantly spin down. To achieve the latter, a near maximally rotating black hole needs to accrete |$\approx 20~{{\ \rm per\ cent}}$| of its initial mass (e.g. Jacquemin-Ide et al. 2024), which is not the case for the aforementioned models. A demonstration of this process can also be seen in the bottom LHS of Fig. 3, which shows the spin evolution of a |$34.4M_\odot$| black hole post-disc formation. Here, the black hole’s spin is not significantly depleted as the disc only forms near the end of the accretion phase.

The disc formation time for the various models. The LHS figure is for the NTD state while RHS is for the MAD state. The legend is similar to Fig. 5. Note that the scattered green points (with relatively larger $t_{\rm disk}$ values) on the RHS figure are the same scattered models plotted in green colour in Fig. 5(b).
Figure 6.

The disc formation time for the various models. The LHS figure is for the NTD state while RHS is for the MAD state. The legend is similar to Fig. 5. Note that the scattered green points (with relatively larger |$t_{\rm disk}$| values) on the RHS figure are the same scattered models plotted in green colour in Fig. 5(b).

3.2.2 Masses of black holes as a function of |$M_{\rm exp}$| and |$\ell _{\rm exp}$|

To better quantify the dependence of the black hole properties on the progenitor’s structure, we now consider a broader parameter space constituting the mass of the star prior to explosion |$M_{\rm exp}$| and its (averaged) specific angular momentum |$\ell _{\rm exp}$|⁠.

To this end, as shown in Fig. 7, we find an approximate relation between |$M_{\rm exp}$| and |$\ell _{\rm exp}$| (i.e. the black dashed line)

(33)

where we first expect a progenitor star to experience a GRB. Similarly, the relation

(34)

which is depicted by the black dotted lines in the figure, marks the onset of the region where we first expect to see a strong deviation in the black hole masses |$M_{\rm BH}$|⁠. In other words, equation (34) separates the parameter space belonging to the two branches shown earlier in Figs 4 and 5.

The masses of black holes as a function of their progenitor star’s mass and angular momentum reservoir. The figures are drawn from the fits provided in equation (35). The x-axis represents the mass of the star and the y-axis its specific angular momentum at the point of core carbon depletion. Panel (a) is for the NTD state, while panel (b) is for the MAD state. The region in grey dotted lines in panel (b) corresponds to the very luminous systems in Fig. 5(b) that have $E_{\rm iso} \gtrsim 3 \times 10^{53}$ erg but do not satisfy equation (34) to experience significant ADAF winds.
Figure 7.

The masses of black holes as a function of their progenitor star’s mass and angular momentum reservoir. The figures are drawn from the fits provided in equation (35). The x-axis represents the mass of the star and the y-axis its specific angular momentum at the point of core carbon depletion. Panel (a) is for the NTD state, while panel (b) is for the MAD state. The region in grey dotted lines in panel (b) corresponds to the very luminous systems in Fig. 5(b) that have |$E_{\rm iso} \gtrsim 3 \times 10^{53}$| erg but do not satisfy equation (34) to experience significant ADAF winds.

The parameter space lying above the boundary set by equation (34) is the region where the black hole mass experiences a sudden downward shift in trend. This shift can be attributed to the large reservoir of angular momentum of the progenitor, which causes the accretion disc to form soon after the onset of the collapse. After sufficient time, the accretion flow of such a collapse transitions into the ADAF regime, resulting in a strong mass-loss. Moreover, NTD-hosting black holes in Fig. 7(a) also lose a considerable amount of mass in neutrinos from the inner region (owing to the Urca processes), hence efficiently carrying away energy with them.

Finally, the region plotted in grey dotted lines in Fig. 7(b) corresponds to the previously discussed very luminous systems in Fig. 5(b) that have |$E_{\rm iso} \gtrsim 5 \times 10^{53}$| erg but do not satisfy equation (34). These systems have lower angular momentum reservoirs; however, since, for the MAD state, the GRB jets are primarily powered by the black hole spin, a strong jet is still produced. We also note that above a certain value of |$M_{\rm exp}$| pair-instability supernova might occur (Fowler & Hoyle 1964; Rakavy & Shaviv 1967, see Section 4.1 for more information), thus truncating the black hole masses function; however, we provide the remnant masses for an extended range of |$M_{\rm exp}$| values and leave it on the reader to decide the location of the pair-instability cut-off. Similarly, the effect of pulsation pair-instability (Woosley, Blinnikov & Heger 2007) on the star’s |$M_{\rm exp}$| would need to be determined before using the current results.

In the following, we now provide the fit functions that have been used to plot Fig. 7. A comparison between these fits and the true values is provided in Figs C1 and C2. For readability, let us first define |${\ell }_{\rm exp}:= \bar{\ell }_{\rm exp} \times 10^{14} \text{ [cm}^2{\rm s}^{-1}]$|⁠. Then, the corresponding function for the masses of the black holes (in |${\rm M}_{\odot }$|⁠) used for generating Fig. 7 can be approximated within 8 per cent error (with 90 per cent credibility; see Fig. C1) as

(35a)
(35b)

where the superscripts on |$M_{\rm BH}$| represent the assumption on the nature of the accretion flow during the collapse, and the constant coefficients have appropriate units to make the LHS and RHS dimensionally consistent.

3.3 Spin of black holes born from rotating progenitors

Similar to Fig. 7, Fig. 8 shows the dependence of black hole spin on the progenitor’s |$M_{\rm exp}$| and |$\ell _{\rm exp}$|⁠. Here, Fig. 8(a) shows that in the case of an NTD, for any given value of |$M_{\rm exp}$| the black hole spin increases monotonically as a function of |$\ell _{\rm exp}$|⁠. Moreover, beyond a certain threshold value of |$\ell _{\rm exp}$| (i.e. the dotted red line in Fig. 8a, see equation 34), the black hole typically acquires a spin of |$a^{*} \gtrsim 0.9$|⁠. This coincides with the region where the mass of the black hole experiences a sudden downward bifurcation in Fig. 7(a).

Same as Fig. 7 but now showing the dimensionless spin parameter $a^{*}$ of black holes. Figures are drawn using the fits provided in equation (36).
Figure 8.

Same as Fig. 7 but now showing the dimensionless spin parameter |$a^{*}$| of black holes. Figures are drawn using the fits provided in equation (36).

While in the case of an NTD, the GRB energetics fully result from the accretion dynamics, for the MAD state, the spin angular momentum of the black hole plays a direct role in powering the jets. Consequently, we expect a dramatic variation in the final spin of the black holes depending on the nature of the accretion flow. This can be seen in Fig. 8(b), where the maximum attainable spin value for the MAD state is |$a^{*} \approx 0.59$|⁠. This is because, now, the most rapidly spinning black holes fuel strong GRBs right after their formation, which efficiently spin them down to an equilibrium spin of |$a^{*} \approx 0.12$| (cf. Gottlieb et al. 2023; Jacquemin-Ide et al.2024).

As was done previously, we now provide the fit functions that have been used to plot Fig. 8, with a comparison between these fits and the true values provided in Figs C1 and C3. The spin magnitude of the black hole shown in Fig. 8 can be approximated within 11 per cent error (with 90 per cent credibility) as

(36a)
(36b)

where the superscripts on |$a^{*}$| represent the assumption on the nature of the accretion flow during the collapse. Like equation (35), the constant coefficients have appropriate units to make the LHS and RHS dimensionally consistent. To ensure that the fit values remain positive, one should bound the values of |$a_{\rm NT}^{*}$| and |$a_{\rm MA}^{*}$| from below by zero.

3.4 Domain of validity

Although it is challenging to measure the spin of individual components in a merging BBH, the mass-weighted effective spin parameter |$\chi _{\rm eff}$| of the binary can be relatively well estimated from the observed gravitational wave signal (Abbott et al. 2021, 2023). For a BBH comprising masses |$M_1$| and |$M_2$| with spins |$a^{*}_1$| and |$a^{*}_2$| aligned with their orbital angular momentum,

(37)

Interestingly, recent findings suggest that most stellar mass black holes should be born with a negligible spin (Eggenberger et al. 2019; Fuller, Piro & Jermyn 2019) with values around |$a^{*} \approx 0.01$| (Fuller & Ma 2019). Therefore, binaries with |$\chi _{\rm eff} \gtrsim 0.01$| might have a distinct formation history. The prevalent methods for generating spinning black holes in isolated binaries involve either (i) tidal spin-up of a progenitor star within a tight binary, with a fraction of the acquired angular momentum subsequently inherited by the ensuing black hole (e.g. Mandel & de Mink 2016; Marchant et al. 2016; Qin et al. 2018; Bavera et al. 2020), or (ii) sufficient mass (and hence angular momentum) accretion onto the progenitor star (e.g. Cantiello et al. 2007; Ghodla et al. 2023) or the first-born black hole (e.g. Zevin & Bavera 2022).

3.4.1 Spin-up via tidal interaction

Tidal spin-up can occur while both stars are on their hydrogen main-sequence due to the occurrence of chemically homogeneous evolution (CHE e.g. Mandel & de Mink 2016; Marchant et al. 2016) or when the primary star first evolves into a black hole with now having a tidally locked secondary star that is on its helium main-sequence with negligible hydrogen envelope (e.g. Qin et al. 2018; Bavera et al. 2020; Fuller & Lu 2022). Assuming that the companion star’s spin period promptly synchronizes with its orbital period P, the value of such a star’s (averaged) specific angular momentum |$\ell$| can be estimated as

(38)

where |$I, \omega , M_{\rm star}$| are the moment of inertia, orbital angular velocity, and the mass of the star, respectively. Additionally, |$r_*$| represents the radius and |$r_g r_*$| its radius of gyration.

Using equation (38) one can thus estimate the domain of validity of the fit functions provided in equations (35) and (36). For example, using the helium star models presented in Section 2.1, the blue curve in Fig. 9 shows the threshold maximum orbital period of the binary at the onset of helium main-sequence required to sufficiently spin-up the star for it to satisfy equation (34) by the time the star reaches core carbon depletion. Note that the metallicity of such stars needs to be low to minimize mass and, hence, angular momentum loss in stellar winds while the star evolves to core collapse (our helium star models have |$Z = 10^{-4}$|⁠). The orange curve, on the other hand, considers tidal synchronization till core helium depletion and, therefore, may be applicable to high-metallicity stars as well. Post-core helium depletion, the star with even larger metallicities will undergo relatively less mass-loss and, hence, spin angular momentum loss and thus a strong tidal locking might not be necessary.

The maximum period at which the star would acquire and retain sufficient angular momentum to satisfy equation (34) by core carbon depletion. The blue and orange curves assume tidal locking till the onset of the helium main sequence and core helium depletion, respectively. The x-axis shows the mass of the star at the onset of tidal locking, and the y-axis shows the period of the binary assuming that the other component is a black hole with mass ratio $q = 1$. The grey band of uncertainty results from the variance in the value of $r_g$.
Figure 9.

The maximum period at which the star would acquire and retain sufficient angular momentum to satisfy equation (34) by core carbon depletion. The blue and orange curves assume tidal locking till the onset of the helium main sequence and core helium depletion, respectively. The x-axis shows the mass of the star at the onset of tidal locking, and the y-axis shows the period of the binary assuming that the other component is a black hole with mass ratio |$q = 1$|⁠. The grey band of uncertainty results from the variance in the value of |$r_g$|⁠.

3.4.2 Spin-up via mass accretion

Stellar spin-up via mass accretion might occur while the star is on its hydrogen main-sequence (Cantiello et al. 2007; Ghodla et al. 2023). For such systems, we expect the stars to undergo a GRB if they satisfy equation (33) at core carbon depletion. Additionally, the properties of the subsequently formed black hole could be highly affected if such systems also satisfy equation (34) (see Figs 7 and 8). The latter is likely to be only valid in rapidly rotating low-metallicity stars, i.e. |$Z \lesssim \mathcal {O}(5 \times 10^{-4})$| (Ghodla et al. 2023). However, in contrast to black holes born in tight binaries, black holes born from stars that undergo accretion-induced spin-up might not merge efficiently unless the supernovae kick results in compact and/or significantly eccentric orbits. This is because to efficiently spin up a star; one likely requires the formation of an accretion disc around the mass gainer (e.g. Ghodla et al. 2023). However, for the formation of an accretion disc around the mass gainer, one needs to allow for a relatively larger separation between the binary components, which subsequently results in longer gravitational wave merger time. As such, although produced in relatively large numbers (e.g. Ghodla et al. 2023), such black holes might seldom feature in the BBH mergers detected on Earth by current-generation detectors.

3.5 Variation of the fiducial parameters

There are three main free parameters in this work, namely the power-law index s in equation (15), the parameter |$\tilde{\eta }$| that controls the cocoon half opening angle |$\theta _c$| (see Section 2.2.4), and the disc viscosity parameter |$\alpha$| that can impact the black hole properties. While our choices of the fiducial values for these parameters (i.e. |$s = 0.5, \tilde{\eta } = 10/3$|⁠, and |$\alpha = 0.01$|⁠) are motivated by previous studies, here we discuss the impact on the black hole mass and spin when these values are changed.

In addition to the fiducial values, we consider two variants for s, namely |$s = 0.3, 0.8$|⁠, two variants for |$\tilde{\eta }$|⁠, namely |$\tilde{\eta } = 10/4, 10/2$| and one variant for |$\alpha$|⁠, namely |$\alpha = 0.1$|⁠. When |$\alpha = 0.1$|⁠, we set the maximum spin of the black hole (Section 2.3.1) to |$a^{*} = 0.9846$| (Sądowski et al. 2011) which then influences the location of the corresponding |$r_{\rm isco}$|⁠. In each attempt, we only change one variable while keeping the other two set to their fiducial values. However, we note that |$\theta _c$| depends not only on the magnitude of |$\tilde{\eta }$| but also on that of |$L_j$| at the time of disc formation. The latter can change during the variation of any of the aforementioned parameters, implying that the value of |$\theta _c$| could also change.

The results achieved on variation of these parameters are presented in Appendix  C. Figs C4 and C5 show that a variation in s does not have a significant effect on the masses of the black hole compared to the fiducial case with the black hole masses (especially those corresponding to large |$E_{\rm iso}$|⁠) being slightly smaller for a smaller value of s and vice versa. Meanwhile, a variation in |$\tilde{\eta }$| does have a significant impact on the magnitude |$M_{\rm BH}$|⁠. The corresponding results are presented in Figs C6 and C7. As one might expect, lowering the value of |$\tilde{\eta }$| results in less massive black holes as the jet now remains within the star for a longer period of time, hence efficiently exchanging energy with the stellar matter and displacing a larger volume from the polar region. Finally, Fig. C8 shows that raising the value of the viscosity parameter to |$\alpha = 0 .1$| only impacts the masses of the black holes when the accretion flow is in an NTD state. The inner radius of the accretion disc dictates the amount of mass that would be liberated in energy. Since for |$\alpha = 0.1$|⁠, the maximum value of |$a^{*} = 0.9846$| (in contrast to 0.9994 for |$\alpha = 0.01$|⁠) thus, in the case of NTD, less energy is radiation away. Meanwhile, as discussed in Section 3.1, for the MAD state, the disc may truncate further away from ISCO. Hence, for the latter, changing |$\alpha$| does not strongly influence the location of the inner disc radii.

Similarly, Figs C9 and C10 show the effect of s variation, Figs C11 and C12 show the effect of |$\tilde{\eta }$| variation and Fig. C13 shows the effect of |$\alpha$| variation on the black hole spin, respectively. For all cases, the deviation in the black hole spin from the fiducial case is minimal. In particular, for later use, we note that for the MAD state |$a^{*} \lesssim 0.59$| when |$s = 0.3$|⁠, |$a^{*} \lesssim 0.62$| when |$s = 0.8$|⁠, |$a^{*} \lesssim 0.6$| when |$\tilde{\eta } = 10/4$|⁠, |$a^{*} \lesssim 0.6$| when |$\tilde{\eta } = 10/2$|⁠, and |$a^{*} \lesssim 0.61$| when |$\alpha = 0.1$|⁠.

4 DISCUSSION

In this work, we have adopted two types of accretion flows that could occur near the region surrounding the black hole during collapsar dynamics to predict the black hole’s mass and spin. While a magnetically arrested accretion flow can be invoked to produce the antiparallel GRB jets via the Blandford & Znajek (1977) mechanism, jets can also be produced in astrophysical objects lacking a black hole at the centre (e.g. Metzger et al. 2011; Metzger, Berger & Margalit 2017). Hence, the energy of the jets need not come from the rotational energy of the black hole. For example, one can expect the formation of anti-parallel jets purely from the energy released during the accretion flow via the Blandford & Payne (1982) mechanism. Therefore, the approach of employing a Novikov–Thorne accretion flow for the collapse dynamics is also relevant. However, a hydrodynamical flow (such as in the case of an NTD) would still require the presence of magnetic field lines to channel the outflow (Blandford & Payne 1982). In the following, we discuss some implications of the results detailed in Section 3 depending on whether an NTD or MAD operates to produce the GRB jets followed by some caveats and sources of uncertainty in our work.

4.1 Impact of stellar rotation on the lower end of PISN mass-gap

Stars with helium cores in the mass range of |$\approx 60{\!-\!}130\, {\rm M}_\odot$| may develop the right condition for |$e^-e^+$| pair-production induced pair-instability supernova (PISN). This results in a (local) upper bound on the value of |$M_{\rm exp}$| that could survive during core oxygen fusion (Fowler & Hoyle 1964; Rakavy & Shaviv 1967) without being completely disrupted due to the aforementioned process. Consequently, one also expects the formation of an upper mass-gap in the black hole mass spectrum.

However, rapid rotation can shift the lower boundary of PISN to larger |$M_{\rm exp}$| values (Glatzel, Fricke & El Eid 1985). Marchant & Moriya (2020) found that for such systems, the lower edge of the black hole mass-gap can move upwards by |$4{\!-\!}15~{{\ \rm per\ cent}}$| depending on the efficiency of angular momentum transport within the progenitor stars during their evolution. In contrast, this work suggests that if the collapse of these rapid rotators5 operates under the NTD formalism, such that |$a^{*} \gtrsim 0.9$|⁠, then the maximum black hole mass for such models should converge around |$M_{\rm BH} \approx 35\,{\rm M}_\odot$| for the fiducial choice of parameters. A noticeable deviation from this value occurs when |$\tilde{\eta }$| is changed to a non-fiducial value. However, this also does not result in |$M_{\rm BH}$| that are much larger that |$35 \,{\rm M}_\odot$|⁠. As the latter value is much lower than the expected lower end of the mass-gap near |$M_{\rm BH} \in ~\approx ~[45\, {\rm M}_\odot , 55 \,{\rm M}_\odot ]$|⁠, e.g. Farmer et al. (2019); Marchant & Moriya (2020), consequently, rapid rotators under a Novikov–Thorne accretion flow seem unlikely to extend the lower end of the PISN mass-gap.

4.2 Case for a massive and maximally spinning stellar black hole

We find that in MAD state, the maximum black hole spin accumulates near |$a^{*} \approx 0.62$| (see Section 3.5) while for in NTD state, it takes the value |$a^{*} \approx 1$|⁠. However, for the latter, when |$a^{*} \approx 1$|⁠, the collapsar energetic places an upper limit on the black hole mass as |$M_{\rm BH} \approx 35 {\rm M}_{\odot }$|⁠, which could be even smaller given that the |$M_{\rm exp}$| values corresponding to these massive black holes is large enough for them to experience pair-instability supernovae, e.g. Farmer et al. (2019); Marchant & Moriya (2020). Hence, irrespective of the nature of the accretion flow, this work suggests that a near maximally rotating black hole (i.e. |$a^{*} \approx 1$|⁠) with |$M_{\rm BH} \gtrsim 35 M_\odot$| might not result directly from stellar collapse. Instead, such a black hole could either be a second-generation black hole (i.e. a merger product of an earlier BBH system, e.g. Pretorius 2005; Fishbach, Holz & Farr 2017; Gerosa & Berti 2017) or a product of subsequent mass (and angular momentum) accretion (e.g. Gammie et al. 2004; Sądowski et al. 2011). For the latter, super-Eddington accretion would be the likely scenario, as sub-Eddington accretion would make it difficult to spin up black holes to critical within the lifetime of a massive mass-donor star. Cygnus X-1 with |$a^{*} \gt 0.95$| (Gou et al. 2011, assuming that most of |$a^{*}$| is acquired from the direct collapse of its progenitor star) may appear as a direct counterexample to the above statement, namely black holes with |$a^{*} \approx 1$| and |$M_{\rm BH} \gtrsim 35 \,{\rm M}_\odot$| might not result directly from stellar collapse. However, we note that Cygnus X-1 has a mass of |$\approx 21 \, {\rm M}_\odot$|⁠, which falls below the derived mass limit.

4.3 A possible upper bound on |$\chi _{\rm eff}$| in black hole–helium star binaries

From Fig. 8, we see that for a given |$M_{\rm exp}$|⁠, a higher value of |$\ell _{\rm exp}$| does not always lead to a larger |$a^{*}$|⁠. Instead, once equation (34) is satisfied, the assumption of the nature of the accretion flow becomes crucial. For the case of MAD, a rapid spin-down is observed, resulting in |$a^{*} \approx 0.12$| (cf. Gottlieb et al. 2023; Jacquemin-Ide et al. 2024). This results in an upper bound on the spin with |$a^{*} \lesssim 0.62$|⁠. In the current 1D approach, the latter value is also not very sensitive to a change in the magnitude of |$s, \alpha$| or |$\tilde{\eta }$| (e.g. see Section 3.5). In such a case, under the scenario when only the second-born black hole acquires spin in a tidally locked helium star-black hole binary (e.g. see Section 3.4.1), an upper bound on |$\chi _{\rm eff}$| can be estimated for such systems as

(39)

where |$q:= M_2/M_1 \le 1$| and we set |$q_{\rm max} = 1$|⁠.

On the other hand, the LIGO-Virgo-KAGRA collaboration (Abbott et al. 2021, 2023) has identified more than half a dozen candidate BBHs that do not respect this equality (see e.g. Table C1). If MAD is the mechanism behind GRBs, this hints that an alternative scenario might be behind the formation of such BBHs. Tidal CHE can produce systems with large |$\chi _{\rm eff}$| values since, in this case, both black holes can acquire large spin. However, the latter requires nearly equal mass-ratio binaries (Marchant et al. 2016), which is not the case for the black holes in Table C1. Alternatively, super-Eddington accretion on to the first formed black hole in an ordinary binary system can also generate a rapidly spinning (and massive) black hole. However, the systems in Table C1 likely result from massive stellar progenitors that have a lifespan of only a few million years. This allows a tiny window for mass transfer, which implies the requirement of a rapid mass-transfer rate supported by a very strong super-Eddington accretion. Additionally, the period of the binary should be such that they could subsequently merge within a few billion years. It is uncertain as to how likely these conditions could be met simultaneously.

A more plausible formation pathway for BBHs in Table C1 could be via hierarchical merger (e.g. Fishbach et al. 2017; Gerosa & Berti 2017; Kimball et al. 2021), e.g. in a triple system or dynamical capture in stellar clusters or in the disc of active galactic nuclei (AGNs; e.g. Gerosa & Fishbach 2021; Arca Sedda, Naoz & Kocsis 2023 and reference therein). This is because the merger of even non-rotating black holes can result in second-generation black holes with |$a^{*} \approx 0.7$| (e.g. Pretorius 2005). Additionally, BBHs living in the discs of AGN also have the possibility of accreting mass and angular momentum from the disc.

4.4 Caveats and uncertainties

  • The final structure of the star has been deemed crucial for the following explosion, with some massive stars producing a neutron star instead of a black hole, leading to the so-called islands of explodeability (O’Connor & Ott 2011; Sukhbold et al. 2016; Sukhbold, Woosley & Heger 2018). Although it remains to be seen how this fares in the case of rapidly rotating stars where the dynamics of core collapse are dominated by the collapsar energetics. However, there is a possibility for a magnetar-like explosion (Kasen & Bildsten 2010; Woosley 2010) where the energetics are dictated by the magnetic fields of a proto-neutron star instead of a black hole. We ignore such a scenario in this work. As such, in our analysis, a black hole always forms during the collapse of the core of a massive star.

  • At times, the stochasticity of the core structure can influence the nature of remnant masses. However, any stochasticity would be lost in our approach due to the averaging performed in the analysis. As discussed in Fryer, Olejak & Belczynski (2022), the effect of stochasticity would nevertheless average out for a large enough sample size, which is always true in the context of population studies. This would also be the case with observations in the coming year as the number of gravitational wave detections grows.

  • In the case of non/slowly rotating stars, our formalism results in a direct collapse. Although this might not necessarily be the case, the lack of detections of supernova explosions from high-mass stars (Smartt 2015) as well as the theoretical findings that more massive stars might disappear without leaving any observational feature (Murguia-Berthier et al. 2020), suggest that black holes forming from more massive stars might experience direct collapse (e.g. see Adams et al. 2017 for a failed supernova candidate). However, we note that stochastically varying angular momentum accretion on the proto-neutron star may form jittering jets (Papish & Soker 2011; Soker 2024), which could prevent direct collapse.

  • For the case of an NTD, the accretion flow remains circular with the disc truncating at the ISCO. However, numerical studies find that at a large accretion rate, the inner edge of the disc moves closer to the black hole but cannot be uniquely defined (Abramowicz et al. 2010). The flow also deviates from a Keplerian profile (Abramowicz et al. 2010; Sądowski et al. 2011), becoming super-Keplerain for our choice of |$\alpha = 0.01$|⁠. We do not consider these aspects here. For rapid rotators, the ISCO will approach very close to the black hole, hence reducing our margin of error.

  • In the current work, we ignore the effect of rotation deformation of the star on the mass accretion rate which may affect our results. Additionally, In the case when a non-trivial change in stellar physics results in a substantial change in the star’s angular momentum profile (e.g. shown in Fig. 1) at collapse, our calculations could be affected.

  • Finally, we note that the current results are a 1D approximation of a complex 3D phenomenon. For greater precision, one would need to conduct a full 3D GRMHD study.

5 CONCLUSION

Black holes born from rapidly rotating progenitor stars can experience a very different supernova explosion than those born from non-/slowly rotating stars. This is because a rapidly rotating stars contain a large reservoir of angular momentum that they would need to lose to collapse below a certain characteristic radius. Here, we studied the resulting masses and spins of the black hole when dealing with the collapse of such rotating stars. We showed that for rapid rotators, the black hole mass and spin is a function of the assumption on the accretion flow and can be significantly different from that expected for a non-/slowly rotating progenitor star (Section 3). This discrepancy could influence the merger time, luminosity distance, and gravitational wave properties of the BBH. We find that rapid rotators undergoing a collapsar explosion via a Novikov–Thorne accretion flow might not extend the lower edge of the black hole upper mass-gap (Section 4.1) and that a maximally rotating black hole born directly from stellar collapse might have |$M_{\rm BH} \lesssim 35M_\odot$| as its maximum possible mass (Section 4.2). For the case of a magnetically arrested accretion flow, we find a maximum black hole dimensionless spin of |$a^{*} \approx 0.62$|⁠. For black holes born in helium star black hole binaries, this puts an upper limit on the effective spin of the BBH as |$\chi _{\rm eff} \lesssim 0.31$| (Section 4.3).

ACKNOWLEDGEMENTS

The authors are grateful to Dr Pablo Marchant and the anonymous referee whose feedback significantly improved this paper. SG is supported by the University of Auckland Doctoral Scholarship. JJE acknowledges support of Marsden Fund Council managed through Royal Society Te Apārangi. This work utilized NeSI high performance computing facilities.

DATA AVAILABILITY

The numerical work discussed in Section 2 can be found at https://github.com/SohanGhodla/Collapsar-Formalism. The sage manifolds notebook used to conduct part of the calculation can also be found at the same URL. The simulated data on which the fits were applied have been included in the supplementary material.

Footnotes

1

The inlists for the mesa model, the sage manifolds notebook containing parts of the analytical calculation, and the python script consisting of the numerical implementation of the formalism discussed in Section 2.2 can be found in the data availability statement.

2

The Latin indices run from 0 to 3 and the Greek indices from 1 to 3.

3

Note that we have added the subscript ‘BH’ to the variables to explicitly associate them to the black hole.

4

E.g., for |$a^{*} = 0.9978, 0.9994$| and |$u_t(r_{\rm isco}) \approx 0.68, 0.58$|⁠, respectively.

5

We allow for efficient angular momentum transport within the stellar interior aided by the Spruit (2002) dynamo. Turning this off would preserve more angular momentum within the star and result in favourable conditions for producing rapid rotators.

REFERENCES

Abbott
R.
et al. ,
2021
,
Phys. Rev. X
,
11
,
021053

Abbott
R.
et al. ,
2023
,
Phys. Rev. X
,
13
,
041039

Abramowicz
M.
,
Chen
X.-M.
,
Granath
M.
,
Lasota
J.-P.
,
1996
,
ApJ
,
471
,
762

Abramowicz
M. A.
,
Jaroszyński
M.
,
Kato
S.
,
Lasota
J. P.
,
Różańska
A.
,
Sądowski
A.
,
2010
,
A&A
,
521
,
A15

Adams
S. M.
,
Kochanek
C. S.
,
Gerke
J. R.
,
Stanek
K. Z.
,
Dai
X.
,
2017
,
MNRAS
,
468
,
4968

Arca Sedda
M.
,
Naoz
S.
,
Kocsis
B.
,
2023
,
Universe
,
9
,
138

Asplund
M.
,
Grevesse
N.
,
Sauval
A. J.
,
Scott
P.
,
2009
,
ARA&A
,
47
,
481

Bardeen
J. M.
,
Press
W. H.
,
Teukolsky
S. A.
,
1972
,
ApJ
,
178
,
347

Batta
A.
,
Ramirez-Ruiz
E.
,
2019
,
preprint
()

Bavera
S. S.
et al. ,
2020
,
A&A
,
635
,
A97

Bavera
S. S.
et al. ,
2022
,
A&A
,
657
,
L8

Bekenstein
J. D.
,
1972
,
Phys. Rev. D
,
5
,
1239

Bičák
J.
,
Stuchlik
Z.
,
1976
,
MNRAS
,
175
,
381

Blandford
R. D.
,
Begelman
M. C.
,
1999
,
MNRAS
,
303
,
L1

Blandford
R. D.
,
Payne
D. G.
,
1982
,
MNRAS
,
199
,
883

Blandford
R. D.
,
Znajek
R. L.
,
1977
,
MNRAS
,
179
,
433

Blondin
J. M.
,
Mezzacappa
A.
,
DeMarino
C.
,
2003
,
ApJ
,
584
,
971

Cantiello
M.
,
Yoon
S. C.
,
Langer
N.
,
Livio
M.
,
2007
,
A&A
,
465
,
L29

Carter
B.
,
1971
,
Phys. Rev. Lett.
,
26
,
331

Chen
W.-X.
,
Beloborodov
A. M.
,
2007
,
ApJ
,
657
,
383

Eddington
A. S.
,
1925
,
The Observatory
,
48
,
73

Eggenberger
P.
,
den Hartogh
J. W.
,
Buldgen
G.
,
Meynet
G.
,
Salmon
S. J. A. J.
,
Deheuvels
S.
,
2019
,
A&A
,
631
,
L6

Eldridge
J. J.
,
Tout
C. A.
,
2004
,
MNRAS
,
353
,
87

Farmer
R.
,
Renzo
M.
,
de Mink
S. E.
,
Marchant
P.
,
Justham
S.
,
2019
,
ApJ
,
887
,
53

Fishbach
M.
,
Holz
D. E.
,
Farr
B.
,
2017
,
ApJ
,
840
,
L24

Fowler
W. A.
,
Hoyle
F.
,
1964
,
ApJS
,
9
,
201

Fryer
C. L.
,
Kalogera
V.
,
2001
,
ApJ
,
554
,
548

Fryer
C. L.
,
Belczynski
K.
,
Wiktorowicz
G.
,
Dominik
M.
,
Kalogera
V.
,
Holz
D. E.
,
2012
,
ApJ
,
749
,
91

Fryer
C. L.
,
Karpov
P.
,
Livescu
D.
,
2021
,
Astron. Rep.
,
65
,
937

Fryer
C. L.
,
Olejak
A.
,
Belczynski
K.
,
2022
,
ApJ
,
931
,
94

Fujibayashi
S.
,
Shibata
M.
,
Wanajo
S.
,
Kiuchi
K.
,
Kyutoku
K.
,
Sekiguchi
Y.
,
2020
,
Phys. Rev. D
,
102
,
123014

Fujibayashi
S.
,
Sekiguchi
Y.
,
Shibata
M.
,
Wanajo
S.
,
2023
,
ApJ
,
956
,
100

Fuller
J.
,
Lu
W.
,
2022
,
MNRAS
,
511
,
3951

Fuller
J.
,
Ma
L.
,
2019
,
ApJ
,
881
,
L1

Fuller
J.
,
Piro
A. L.
,
Jermyn
A. S.
,
2019
,
MNRAS
,
485
,
3661

Gammie
C. F.
,
Popham
R.
,
1998
,
ApJ
,
498
,
313

Gammie
C. F.
,
Shapiro
S. L.
,
McKinney
J. C.
,
2004
,
ApJ
,
602
,
312

Gerosa
D.
,
Berti
E.
,
2017
,
Phys. Rev. D
,
95
,
124046

Gerosa
D.
,
Fishbach
M.
,
2021
,
Nat. Astron.
,
5
,
749

Ghodla
S.
,
Eldridge
J. J.
,
Stanway
E. R.
,
Stevance
H. F.
,
2023
,
MNRAS
,
518
,
860

Glatzel
W.
,
Fricke
K. J.
,
El Eid
M. F.
,
1985
,
A&A
,
149
,
413

Gottlieb
O.
,
Jacquemin-Ide
J.
,
Lowell
B.
,
Tchekhovskoy
A.
,
Ramirez-Ruiz
E.
,
2023
,
ApJ
,
952
,
L32

Gou
L.
et al. ,
2011
,
ApJ
,
742
,
85

Gourgoulhon
É.
,
2018
,
Lecture Notes
. p.
350

Herant
M.
,
Benz
W.
,
Hix
W. R.
,
Fryer
C. L.
,
Colgate
S. A.
,
1994
,
ApJ
,
435
,
339

Israel
W.
,
1967
,
Phys. Rev.
,
164
,
1776

Jacquemin-Ide
J.
,
Gottlieb
O.
,
Lowell
B.
,
Tchekhovskoy
A.
,
2024
,
ApJ
,
961
,
212

Kasen
D.
,
Bildsten
L.
,
2010
,
ApJ
,
717
,
245

Kimball
C.
et al. ,
2021
,
ApJ
,
915
,
L35

King
A. R.
,
Pringle
J. E.
,
Livio
M.
,
2007
,
MNRAS
,
376
,
1740

Kohri
K.
,
Narayan
R.
,
Piran
T.
,
2005
,
ApJ
,
629
,
341

Kumar
P.
,
Narayan
R.
,
Johnson
J. L.
,
2008
,
MNRAS
,
388
,
1729

Langer
N.
,
1998
,
A&A
,
329
,
551

Lazzati
D.
,
Begelman
M. C.
,
2005
,
ApJ
,
629
,
903

Lee
W. H.
,
Ramirez-Ruiz
E.
,
2006
,
ApJ
,
641
,
961

Lowell
B.
,
Jacquemin-Ide
J.
,
Tchekhovskoy
A.
,
Duncan
A.
,
2024
,
ApJ
,
960
,
82

MacFadyen
A. I.
,
Woosley
S. E.
,
1999
,
ApJ
,
524
,
262

Mandel
I.
,
de Mink
S. E.
,
2016
,
MNRAS
,
458
,
2634

Marchant
P.
,
Moriya
T. J.
,
2020
,
A&A
,
640
,
L18

Marchant
P.
,
Langer
N.
,
Podsiadlowski
P.
,
Tauris
T. M.
,
Moriya
T. J.
,
2016
,
A&A
,
588
,
A50

Marchant
P.
,
Podsiadlowski
P.
,
Mandel
I.
,
2023
,
preprint
()

Metzger
B. D.
,
Giannios
D.
,
Thompson
T. A.
,
Bucciantini
N.
,
Quataert
E.
,
2011
,
MNRAS
,
413
,
2031

Metzger
B. D.
,
Berger
E.
,
Margalit
B.
,
2017
,
ApJ
,
841
,
14

Meynet
G.
,
Maeder
A.
,
1997
,
A&A
,
321
,
465

Moderski
R.
,
Sikora
M.
,
1996
,
MNRAS
,
283
,
854

Murguia-Berthier
A.
,
Batta
A.
,
Janiuk
A.
,
Ramirez-Ruiz
E.
,
Mandel
I.
,
Noble
S. C.
,
Everson
R. W.
,
2020
,
ApJ
,
901
,
L24

Nakar
E.
,
Piran
T.
,
2017
,
ApJ
,
834
,
28

Narayan
R.
,
McClintock
J. E.
,
2008
,
New A Rev.
,
51
,
733

Narayan
R.
,
Yi
I.
,
1994
,
ApJ
,
428
,
L13

Novikov
I. D.
,
Thorne
K. S.
,
1973
, in
DeWitt
C.
, ed.,
Black Holes (Les Astres Occlus)
.
Gordon and Breach
,
New York
, p.
343

Nugis
T.
,
Lamers
H. J. G. L. M.
,
2000
,
A&A
,
360
,
227

O’Connor
E.
,
Ott
C. D.
,
2011
,
ApJ
,
730
,
70

Papish
O.
,
Soker
N.
,
2011
,
MNRAS
,
416
,
1697

Paschalidis
V.
,
Stergioulas
N.
,
2017
,
Living Rev. Relativ.
,
20
,
1

Paxton
B.
,
Bildsten
L.
,
Dotter
A.
,
Herwig
F.
,
Lesaffre
P.
,
Timmes
F.
,
2011
,
ApJS
,
192
,
3

Paxton
B.
et al. ,
2013
,
ApJS
,
208
,
4

Paxton
B.
et al. ,
2015
,
ApJS
,
220
,
15

Paxton
B.
et al. ,
2018
,
ApJS
,
234
,
34

Paxton
B.
et al. ,
2019
,
ApJS
,
243
,
10

Penrose
R.
,
Floyd
R. M.
,
1971
,
Nature Phys. Sci.
,
229
,
177

Perley
D. A.
et al. ,
2016
,
ApJ
,
817
,
7

Pretorius
F.
,
2005
,
Phys. Rev. Lett.
,
95
,
121101

Qin
Y.
,
Fragos
T.
,
Meynet
G.
,
Andrews
J.
,
Sørensen
M.
,
Song
H. F.
,
2018
,
A&A
,
616
,
A28

Rakavy
G.
,
Shaviv
G.
,
1967
,
ApJ
,
148
,
803

Riley
J.
,
Mandel
I.
,
Marchant
P.
,
Butler
E.
,
Nathaniel
K.
,
Neijssel
C.
,
Shortt
S.
,
Vigna-Gómez
A.
,
2021
,
MNRAS
,
505
,
663

Sądowski
A.
,
Narayan
R.
,
2016
,
MNRAS
,
456
,
3929

Shakura
N. I.
,
Sunyaev
R. A.
,
1973
,
A&A
,
24
,
337

Sharp
D. H.
,
1984
,
Phys. D: Nonlinear Phenom.
,
12
,
3

Sądowski
A.
,
Narayan
R.
,
2016
,
MNRAS
,
456
,
3929

Sądowski
A.
,
Bursa
M.
,
Abramowicz
M.
,
Kluźniak
W.
,
Lasota
J. P.
,
Moderski
R.
,
Safarzadeh
M.
,
2011
,
A&A
,
532
,
A41

Sądowski
A.
,
Narayan
R.
,
McKinney
J. C.
,
Tchekhovskoy
A.
,
2014
,
MNRAS
,
439
,
503

Smartt
S. J.
,
2015
,
PASA
,
32
,
e016

Soker
N.
,
2024
,
Open J. Astrophys.
,
7
,
31

Spruit
H. C.
,
2002
,
A&A
,
381
,
923

Sukhbold
T.
,
Ertl
T.
,
Woosley
S. E.
,
Brown
J. M.
,
Janka
H. T.
,
2016
,
ApJ
,
821
,
38

Sukhbold
T.
,
Woosley
S.
,
Heger
A.
,
2018
,
ApJ
,
860
,
93

Sweet
P. A.
,
1950
,
MNRAS
,
110
,
548

Tchekhovskoy
A.
,
Narayan
R.
,
McKinney
J. C.
,
2010
,
ApJ
,
711
,
50

Tchekhovskoy
A.
,
Narayan
R.
,
McKinney
J. C.
,
2011
,
MNRAS
,
418
,
L79

Thorne
K. S.
,
1974
,
ApJ
,
191
,
507

Vink
J. S.
,
de Koter
A.
,
Lamers
H. J. G. L. M.
,
2001
,
A&A
,
369
,
574

Woosley
S. E.
,
1993
,
ApJ
,
405
,
273

Woosley
S. E.
,
2010
,
ApJ
,
719
,
L204

Woosley
S. E.
,
Blinnikov
S.
,
Heger
A.
,
2007
,
Nature
,
450
,
390

Yuan
F.
,
Narayan
R.
,
2014
,
ARA&A
,
52
,
529

Zahn
J. P.
,
1992
,
A&A
,
265
,
115

Zevin
M.
,
Bavera
S. S.
,
2022
,
ApJ
,
933
,
86

de Jager
C.
,
Nieuwenhuijzen
H.
,
van der Hucht
K. A.
,
1988
,
A&AS
,
72
,
259

du Buisson
L.
et al. ,
2020
,
MNRAS
,
499
,
5941

APPENDIX A: FREEFALL TIME-SCALE

For a shell made of perfect fluid, the collapse is such that the coordinate |$\theta$| remains unchanged. This means that the trajectory of the particles would not cross each other during infall. Solving for the geodesic motion (assuming an infall from rest at radial infinity), one can show that for such a case

(A1)

Equation (A1) depends on the coordinate |$\theta$|⁠, implying that in a Kerr geometry, the equatorial part of a spherical shell falls faster than the polar region (e.g. Bičák & Stuchlik 1976). Integrating equation (A1) when |$a \ne 0$| is difficult, and the resulting level of precision is not required here. For example, in Fig. A1, we compare the infall coordinate time (resulting from integrating equation (A1) with |$a=0$|⁠) with the infall proper time, assuming that both time coordinates are synchronized at some fixed distance |$r_0 =30M$| (which we approximate as the minimum distance from where an infalling shell would first circularize at ISCO in our work). We find that upon reaching ISCO, there is |$\approx 12~{{\ \rm per\ cent}}$| disagreement in the two time intervals. As |$r_0$| moves to larger values, this disagreement becomes smaller. E.g., for |$r_0 = 200M$|⁠, it is reduced to |$\approx 2~{{\ \rm per\ cent}}$|⁠. Finally, when calculating the disc formation time |$\Delta t$| in Section 2.2.3, we note that the matter falls from a finite radius |$r_0$|⁠, from a near-rest configuration, which also needs to be taken into account.

The deviation in the infall coordinate time and the infall proper time (in a Schwarzschild geometry), assuming that both time coordinates are synchronized at some fixed distance $r_0 =30M$, which we approximate as the minimum distance from where an infalling shell would first circularize at ISCO ($r=6M$) in our work.
Figure A1.

The deviation in the infall coordinate time and the infall proper time (in a Schwarzschild geometry), assuming that both time coordinates are synchronized at some fixed distance |$r_0 =30M$|⁠, which we approximate as the minimum distance from where an infalling shell would first circularize at ISCO (⁠|$r=6M$|⁠) in our work.

APPENDIX B: An estimate for |$\theta _c$|

Let us define the breakout time of the jet from the stellar surface as |$t_{\mathrm{br}} = \tilde{\eta } r_*/c$|⁠, where |$\tilde{\eta } \gt 1$| and is a measure of how fast the jet head travels through the stellar material (i.e. |$\tilde{\eta }$| is the ratio of the speed of light to the jets’ head speed). Then, the relativistic cocoon half-opening angle can be defined as

(B1)

Above |$r_*, r_\perp$| is the radius of the star and the transverse radius of the cocoon and |$v_{\rm sh}$| is the velocity of the shock wave in the transverse direction and takes the form |$v_{\rm sh} = \sqrt{ \frac{p_c}{\rho _*}}$| (Lazzati & Begelman 2005), where |$p_c, \rho _*$| are the pressure within the cocoon and the mean density of the star. Lazzati & Begelman (2005) calculates the pressure as

(B2)

where |$L_j$| is the luminosity of the jet. This gives

(B3)

Unfortunately, |$\tilde{\eta }$| remains a free parameter in this work. However, the jet head moves slowly |$(\approx 0.3c)$| within the baryonic medium of the star (Nakar & Piran 2017), implying that |$\tilde{\eta } \approx 10/3$|⁠.

APPENDIX C: ADDITIONAL MATERIAL

See Table C1 and Figs C1C13.

Panels (a) and (b) show the fitting error in the black hole masses given in equation (35), where panel (a) is for NTD and panel (b) for MAD state. Panels (c) and (d) show the fitting error in the black hole spins given in equation (36), where panel (c) is for NTD and panel (d) for MAD state. The intersection of the histogram with the horizonal line shows the percentage error within which 90 per cent of the systems tend to lie – as predicted by equation (35) or (36) – when compared to the true value.
Figure C1.

Panels (a) and (b) show the fitting error in the black hole masses given in equation (35), where panel (a) is for NTD and panel (b) for MAD state. Panels (c) and (d) show the fitting error in the black hole spins given in equation (36), where panel (c) is for NTD and panel (d) for MAD state. The intersection of the histogram with the horizonal line shows the percentage error within which 90 per cent of the systems tend to lie – as predicted by equation (35) or (36) – when compared to the true value.

A comparison between the fit predictions (blue dots) and the simulation data (orange dots) for the black hole mass in the NTD and the MAD state.
Figure C2.

A comparison between the fit predictions (blue dots) and the simulation data (orange dots) for the black hole mass in the NTD and the MAD state.

A comparison between the fit predictions (blue dots) and the simulation data (orange dots) for the black hole spin in NTD and the MAD state.
Figure C3.

A comparison between the fit predictions (blue dots) and the simulation data (orange dots) for the black hole spin in NTD and the MAD state.

The figure is the same as Fig. 5. Panel (a) is for the NTD state, while panel (b) is for the MAD state. The dark colour is for the case when $s = 0.3$, while the light colour for $s = 0.5$. The latter is the same as Fig. 5.
Figure C4.

The figure is the same as Fig. 5. Panel (a) is for the NTD state, while panel (b) is for the MAD state. The dark colour is for the case when |$s = 0.3$|⁠, while the light colour for |$s = 0.5$|⁠. The latter is the same as Fig. 5.

Same as Fig. C4 but with $s = 0.8$.
Figure C5.

Same as Fig. C4 but with |$s = 0.8$|⁠.

Same as Fig. C4 but with instead of changing s here we have changed $\tilde{\eta }$. The dark colour is for the case when $\tilde{\eta } = 10/4$, while the light colour for the fiducial case ($\tilde{\eta } = 10/3$.). The latter is the same as Fig. 5.
Figure C6.

Same as Fig. C4 but with instead of changing s here we have changed |$\tilde{\eta }$|⁠. The dark colour is for the case when |$\tilde{\eta } = 10/4$|⁠, while the light colour for the fiducial case (⁠|$\tilde{\eta } = 10/3$|⁠.). The latter is the same as Fig. 5.

Same as Fig. C6 but with $\tilde{\eta } = 10/2$. Note the change in the lower range of the y-axis.
Figure C7.

Same as Fig. C6 but with |$\tilde{\eta } = 10/2$|⁠. Note the change in the lower range of the y-axis.

Same as Fig. C4 but instead of changing s here we have changed $\alpha$. The dark colour is for the case where $\alpha = 0.1$, while the light colour is for $\alpha = 0.01$. The latter is the same as Fig. 5.
Figure C8.

Same as Fig. C4 but instead of changing s here we have changed |$\alpha$|⁠. The dark colour is for the case where |$\alpha = 0.1$|⁠, while the light colour is for |$\alpha = 0.01$|⁠. The latter is the same as Fig. 5.

Same as Fig. 8 but with $s = 0.3$. The colourbar now shows the residual value of $a^{*}$, i.e. $\Delta a^{*} = a^{*}_{\rm fiducial}-a^{*}_{s = 0.3}$.
Figure C9.

Same as Fig. 8 but with |$s = 0.3$|⁠. The colourbar now shows the residual value of |$a^{*}$|⁠, i.e. |$\Delta a^{*} = a^{*}_{\rm fiducial}-a^{*}_{s = 0.3}$|⁠.

Same as Fig. C9 but with $s = 0.8$.
Figure C10.

Same as Fig. C9 but with |$s = 0.8$|⁠.

Same as Fig. C9 but instead of changing s here we change $\tilde{\eta }$ to $\tilde{\eta } = 10/4$.
Figure C11.

Same as Fig. C9 but instead of changing s here we change |$\tilde{\eta }$| to |$\tilde{\eta } = 10/4$|⁠.

Same as Fig. C11 but with $\tilde{\eta } = 10/2$.
Figure C12.

Same as Fig. C11 but with |$\tilde{\eta } = 10/2$|⁠.

Same as Fig. C9 but instead of changing s here we change $\alpha$ to $\alpha = 0.1$.
Figure C13.

Same as Fig. C9 but instead of changing s here we change |$\alpha$| to |$\alpha = 0.1$|⁠.

Table C1.

LVK inferred parameters for systems with |$\chi _{\rm eff} \ge 0.32$| (Abbott et al. 2021, 2023).

Name|$M_1 \, [{\rm M}_\odot ]$||$M_2 \, [{\rm M}_\odot ]$||$\chi _{\rm eff}$|
GW190517_055101|$39.2^{+13.9}_{-9.2}$||$24.0^{+7.4}_{-7.9}$||$0.49^{+0.21}_{-0.28}$|
GW200208_222617|$51^{+103}_{-30}$||$12.3^{+9.2}_{-5.5}$||$0.45^{+0.42}_{-0.46}$|
GW170729|$50.2^{+16.2}_{-10.2}$||$34.0^{+9.1}_{-10.1}$||$0.37^{+0.21}_{-0.25}$|
GW190805_211137|$46.2^{+15.4}_{-11.2}$||$30.6^{+11.8}_{-11.3}$||$0.37^{+0.29}_{-0.39}$|
GW190620_030421|$58.0^{+19.2}_{-13.3}$||$35.0^{+13.1}_{-14.5}$||$0.34^{+0.22}_{-0.29}$|
GW190519_153544|$65.1^{+10.8}_{-11}$||$40.8^{+11.5}_{-12.7}$||$0.33^{+0.20}_{-0.24}$|
GW200306_093714|$28.3^{+17.1}_{-7.7}$||$14.8^{+6.5}_{-6.4}$||$0.32^{+0.28}_{-0.46}$|
Name|$M_1 \, [{\rm M}_\odot ]$||$M_2 \, [{\rm M}_\odot ]$||$\chi _{\rm eff}$|
GW190517_055101|$39.2^{+13.9}_{-9.2}$||$24.0^{+7.4}_{-7.9}$||$0.49^{+0.21}_{-0.28}$|
GW200208_222617|$51^{+103}_{-30}$||$12.3^{+9.2}_{-5.5}$||$0.45^{+0.42}_{-0.46}$|
GW170729|$50.2^{+16.2}_{-10.2}$||$34.0^{+9.1}_{-10.1}$||$0.37^{+0.21}_{-0.25}$|
GW190805_211137|$46.2^{+15.4}_{-11.2}$||$30.6^{+11.8}_{-11.3}$||$0.37^{+0.29}_{-0.39}$|
GW190620_030421|$58.0^{+19.2}_{-13.3}$||$35.0^{+13.1}_{-14.5}$||$0.34^{+0.22}_{-0.29}$|
GW190519_153544|$65.1^{+10.8}_{-11}$||$40.8^{+11.5}_{-12.7}$||$0.33^{+0.20}_{-0.24}$|
GW200306_093714|$28.3^{+17.1}_{-7.7}$||$14.8^{+6.5}_{-6.4}$||$0.32^{+0.28}_{-0.46}$|
Table C1.

LVK inferred parameters for systems with |$\chi _{\rm eff} \ge 0.32$| (Abbott et al. 2021, 2023).

Name|$M_1 \, [{\rm M}_\odot ]$||$M_2 \, [{\rm M}_\odot ]$||$\chi _{\rm eff}$|
GW190517_055101|$39.2^{+13.9}_{-9.2}$||$24.0^{+7.4}_{-7.9}$||$0.49^{+0.21}_{-0.28}$|
GW200208_222617|$51^{+103}_{-30}$||$12.3^{+9.2}_{-5.5}$||$0.45^{+0.42}_{-0.46}$|
GW170729|$50.2^{+16.2}_{-10.2}$||$34.0^{+9.1}_{-10.1}$||$0.37^{+0.21}_{-0.25}$|
GW190805_211137|$46.2^{+15.4}_{-11.2}$||$30.6^{+11.8}_{-11.3}$||$0.37^{+0.29}_{-0.39}$|
GW190620_030421|$58.0^{+19.2}_{-13.3}$||$35.0^{+13.1}_{-14.5}$||$0.34^{+0.22}_{-0.29}$|
GW190519_153544|$65.1^{+10.8}_{-11}$||$40.8^{+11.5}_{-12.7}$||$0.33^{+0.20}_{-0.24}$|
GW200306_093714|$28.3^{+17.1}_{-7.7}$||$14.8^{+6.5}_{-6.4}$||$0.32^{+0.28}_{-0.46}$|
Name|$M_1 \, [{\rm M}_\odot ]$||$M_2 \, [{\rm M}_\odot ]$||$\chi _{\rm eff}$|
GW190517_055101|$39.2^{+13.9}_{-9.2}$||$24.0^{+7.4}_{-7.9}$||$0.49^{+0.21}_{-0.28}$|
GW200208_222617|$51^{+103}_{-30}$||$12.3^{+9.2}_{-5.5}$||$0.45^{+0.42}_{-0.46}$|
GW170729|$50.2^{+16.2}_{-10.2}$||$34.0^{+9.1}_{-10.1}$||$0.37^{+0.21}_{-0.25}$|
GW190805_211137|$46.2^{+15.4}_{-11.2}$||$30.6^{+11.8}_{-11.3}$||$0.37^{+0.29}_{-0.39}$|
GW190620_030421|$58.0^{+19.2}_{-13.3}$||$35.0^{+13.1}_{-14.5}$||$0.34^{+0.22}_{-0.29}$|
GW190519_153544|$65.1^{+10.8}_{-11}$||$40.8^{+11.5}_{-12.7}$||$0.33^{+0.20}_{-0.24}$|
GW200306_093714|$28.3^{+17.1}_{-7.7}$||$14.8^{+6.5}_{-6.4}$||$0.32^{+0.28}_{-0.46}$|
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Supplementary data