Variety of disk wind-driven explosions in massive rotating stars

We perform a set of two-dimensional, non-relativistic, hydrodynamics simulations for supernova-like explosion associated with stellar core collapse of rotating massive stars to a system of a black hole and a disk connected by the transfer of matter and angular momentum. Our model of the central engine also includes the contribution of the disk wind. In this work, we specifically investigate the wind-driven explosion of rotating, large-mass progenitor stars with the zero-age main-sequence mass of $M_\mathrm{ZAMS}=20\,M_\odot$ from arXiv:2008.09132 . This study is carried out using the open-source hydrodynamic code Athena++, for which we implement a method to calculate self-gravity for axially symmetric density distributions. We, then, investigate the explosion properties and the $^{56}$Ni production as a function of (varying) some features of the wind injection. We find a large variety of explosion energy with $E_\mathrm{expl}$ ranging from $\sim 0.049\times10^{51}$~erg to $\sim 34\times10^{51}$~erg and ejecta mass $M_\mathrm{ej}$ from 0.58 to 6 $M_\odot$, which shows a bimodal distribution in high- and low-energy branches. We demonstrate that the resulting outcome of a highly- or sub-energetic explosion for a certain stellar structure is mainly determined by the competition between the ram pressure of the injected matter and that of the infalling envelope. In the nucleosynthesis analysis the $^{56}$Ni mass produced in our models goes from $<0.2~M_\odot$ in the sub-energetic explosions to $2.1~M_\odot$ in the highly-energetic ones. These results are consistent with the observational data of stripped-envelope and high-energy SNe such as broad-lined type Ic SNe. However, we find a tighter correlation between the explosion energy and the ejecta mass than that observationally measured.


INTRODUCTION
Gamma-ray bursts (GRBs) are extragalactic cosmological sources of gamma-rays and among the most energetic events that we can observe in the Universe.They are short (GRBs typically last from few milliseconds up to a few minutes) and very intense flashes of gamma-rays of variable intensity, with fluxes up to ∼ 100 photons cm −2 s −1 usually ranging from hundreds of keV up to ∼1 MeV.GRBs release gamma-rays reaching a total isotropic equivalent radiation energy (i.e. the radiated energy if the GRB was equally bright in all directions) of ∼ 10 53 − 10 54 erg.GRBs properties, such as the total energy, spectra, and duration, can be useful source of information about their progenitor (Mészáros 2006, Woosley & Bloom 2006).
Over the years evidence has showed that GRBs of the "long-soft" variety (lGRBs) are likely to originate from the deaths of massive stars (Woosley 1993, Woosley & Heger 2006, Woosley & Bloom 2006, Janiuk, Agnieszka et al. 2013) and many gamma ray bursts have been now associated with bright supernovae (SNe) (Woosley & Heger 2006).
The Photometric and spectroscopic observations suggest that GRBs and their SNe progenitor have aspherical features.The signature of a conical geometry of the bursts manifests itself as a broadband break in the power-law decay of the GRB afterglow, known as "jet break".This break can be explained by relativistic beaming of light emitted by a decelerated relativistic jet (Frail et al. 2001, Piran 2004) and it is predicted to be achromatic.
Several scenarios have been proposed to explain the GRBs and associated SNe (Woosley (1993), Piran (2004), Woosley & Heger (2006), Trigo-Rodríguez et al. (2017), Obergaulinger & Aloy (2022)).One of the most promising scenarios is the collapsar scenario.The collapse of the core of a massive star (≳ 8  ⊙ ) at the end of its hydrostatic evolution is the starting point for a complex sequence of events with many possible outcomes.Specifically, progenitors with an even higher mass (> 16  ⊙ ), as shown by Woosley & Heger (2006), are likely to undergo a failed supernova and form a black hole (BH) with an accreting disk.It has been shown that in failed supernova the disk wind generated by viscous dissipation inside the accretion disk may naturally be a source of the SN energy (Woosley (1993), MacFadyen & Woosley (1999), Popham et al. (1999)) with an explosion energy  expl > 10 52 erg and it has been found to be rich in 56 Ni (as shown by Hayakawa & Maeda (2018)).Also recent numerical studies based on this scenario have confirmed that a large amount of 56 Ni (≥ 0.1  ⊙ ) can be synthesized in the outflow from the disk (e.g., Just et al. 2022 andFujibayashi et al. 2023a).In this scenario, the interaction between the new-born BH and the still accreting stellar material is the engine for the relativistic jets.
Another promising scenario for the GRBs and associated Type Ic-BL SNe is the so-called proto-magnetar scenario, in which highly magnetized and fast-rotating proto-neutron star generates the relativistic outflow.In this scenario, rotation leads to global asymmetries of the shock wave, which translates into the formation of highly collimated, mildly relativistic bipolar outflow (known as MHD-driven supernova) as shown by Obergaulinger & Aloy (2017), Obergaulinger & Aloy (2020), Aloy & Obergaulinger (2021), Obergaulinger & Aloy (2021).In their study Grimmett et al. (2021) investigate the produc-tion of 56 Ni by performing hydrodynamics simulations based on this scenario.In their most energetic model (in which they measure an energy deposition rate > 10 52 erg s −1 ), they find large masses of ejected 56 Ni (> 0.05 − 0.45  ⊙ ) which is in good agreement with the ranges inferred from the light curves of SNe Ic-BL (0.12 − 0.8  ⊙ with median at 0.28  ⊙ as measured by Taddia et al. 2019a).Therefore, both GRB formation models, the proto-magnetar and the collapsar scenario, seem to be equally plausible at the current moment.
The different properties of the ejecta such as mass, composition, velocity, and geometry strongly depend on the explosion mechanism of the SN.A key to investigate the ejecta properties, is to study its asphericity.Since the progenitor stars have to be rapidly rotating, when they collapse in both previously mentioned scenarios, the resulting ejecta may naturally have aspherical features, as seen in observations of SN 1998bw.When the injection of the energy in the stellar envelope is aspherical, the matter can keep infalling also after a successful explosion.As the energy source of the ejecta may be the infalling mass to the central engine, the feedback of the injected outflow on the infall stellar envelope is an important effect in the scenario.In the case of relativistic bipolar jets, the feedback effect have been extensively studied (Papish & Soker 2013, Liu et al. 2019).However, such an effect by sub-relativistic outflow has not yet been studied in a systematic manner.
The motivation for this work is, therefore, to explore the properties of the ejecta based on the collapsar scenario, with a focus on the late-phase mass ejection after BH formation.We perform a set of two-dimensional hydrodynamics simulations of axisymmetric models of the ejecta generated by the collapse of rotating massive star.Based on the collapsar scenario, we assume that the explosion is powered through a BH-accretion disk system.We vary several parameters controlling the properties of the mass and energy injection to investigate their impact on the final ejecta.
The paper is structured as follows.We begin by explaining the hydrodynamic code we utilize in this work (the hydrodynamic equations it solves, the model for the central engine, the equation of state), the characteristics of the progenitor star we employ (taken from Aguilera-Dena et al. 2020) and the setup for our simulations (inner boundary conditions and the free parameters of our models), in Section §2.In Section §3 we present the results of the simulations focusing on the hydrodynamics of the explosion, the ejecta property and the 56 Ni production with a systematic variation of the initial parameters.Here we also compare our results with observational data and a general relativistic neutrino-radiation viscous-hydrodynamics simulation performed using the same progenitor from the literature.We discuss the implications of our results also considering the observational counterpart.We summarize this work in Section §5.The Appendixes provide the description of the multipole expansion of the gravitational potential we implement in our code and an insight to the model of the disk wind we used.

METHOD
We study the explosion of rotating massive stars ( ∼ 20  ⊙ ) by performing a set of 2D non-relativistic simulations using the opensource multi-dimensional hydrodynamics code Athena++ (Stone et al. 2020).The nucleosynthesis calculation is performed at posteriori using the reaction network torch (Timmes et al. 2000) on tracer particles.

The scenario
For this study, we consider the case of a failed core collapse supernova CCSN, in which the neutrino-driven explosion in the proto-neutron star (PNS) phase does not occur, leading the PNS to collapse into a BH.In this collapsar scenario we model the explosion of a compact progenitor star after the formation of a BH (but see Burrows et al. (2019) for a different scenario).As a progenitor we employ the model provided by Aguilera-Dena et al. (2020) of a rapidly rotating, rotationally mixed star with a 20 ⊙ ZAMS mass.We then build a semi-analytical model for the central engine by taking into account the BH and disk evolution, which in this scenario are governed by the transfer of matter and angular momentum.This method is based on the prescriptions provided by Kumar et al. (2008) on which we add the contribution of the disk following Hayakawa & Maeda (2018).
We chose this progenitor because this kind of stars, with masses ranging from 4 to 45  ⊙ were proposed to be progenitors of both superluminous SNe and long gamma-ray bursts (Japelj et al. 2016, Margalit et al. 2018, Aguilera-Dena et al. 2020).We specifically employ the model with  ZAMS = 20  ⊙ because it is supposed to fail the explosion (Ertl et al. 2016 andMüller et al. 2016) as it has a very compact core with a core compactness of  2.5 > 0.6.Here, the core compactness is calculated by following O' Connor & Ott (2011): where () is the radius at which its enclosed mass is , and it is measured at a mass coordinate of 2.5 ⊙ at the core collapse.This quantity measures the gravitational binding energy near the core of pre-SN stars and is considered as an indicator of whether the collapse of a non-rotating stellar core leads to a successful explosion, or ends up with the formation of a BH instead.Sukhbold & Woosley (2014) found that, if  2.5 > 0.45 at the core collapse, the core collapse is likely to fail the explosion and form a BH.Therefore this  ZAMS = 20  ⊙ progenitor well suits our central engine model in this sense.

Hydrodynamic equations
We perform two-dimensional non-relativistic hydrodynamic simulations using an open-source code Athena++.In addition to the original functions, we newly implement the gravitational potential Φ by solving the Poisson's equation under the cylindrical symmetry.The set of equations solved in this work is as follows: where ,   (with  = ,  and ),  and  t are the density, the velocity components, the pressure and the total energy density of the fluid respectively. t =  kin +  int is the sum of the kinetic energy density  kin = (1/2)  2 and the internal energy density  int .Φ is the gravitational potential which satisfies the Poisson's equation: where  is the gravitational constant,  is the density, and Δ is the Laplacian.This system of equations shows the continuity equation (2), the Euler's equation for the radial, latitudinal, and longitudinal components of the momentum (respectively equation (3), equation (4), and equation ( 5)), and the energy equation ( 6).We compute these equations using a finite volume method on a spherical grid.

The gravity solver
To evaluate the self-gravity in the spherical-polar coordinates, we implemented a gravitational potential solver in our code.In the solver we first use the method of Green's function to obtain the integrated form of the gravitational potential Φ(r) as: where r is a position vector.The potential Φ(r) also satisfies the Poisson's equation (equation ( 7)).To perform the integration, we use a multipole expansion described in Hachisu (1986).We provide more details on the implementation of this solver in Appendix A.

The computational setup
In this work we use an axisymmetric grid with spherical-polar coordinate.Our domain extends from 0 to  for the polar dimension, and from 10 8 cm ( in ) to 3.3 × 10 10 cm ( out ) for the radial dimension.The inner radius determines the inner boundary inside which the enclosed mass is 1.28  ⊙ and it roughly corresponds to the dimension of the iron core at core-collapse which we cut out from the computational domain.The outer radius extends over the stellar surface ( star = 2.7 × 10 10 cm).The initial mass in the computational domain is  domain = 14.2  ⊙ .
The computational domain is discretized by 128 grid points uniformly in the -direction and 220 grid points with geometric spacing in the -direction, in which the mesh size increases with a constant factor Δ  = Δ −1 .We chose the ratio as  ≈ 1.03 to ensure that all the meshes are approximately squared, i.e., Δ  ≈   Δ.This grid resolution was chosen after a convergence study described in Appendix C1.

The central engine model
In our simulations, the computational domain does not contain the central engine, which is considered as being embedded in the central part of the star (at  <  in ).We evolve the masses of the disk and the BH  disk and  BH and their angular momenta  disk and  BH as follows: where  fall,BH and  fall,disk are the rates of the mass accretion that directly infalls onto the BH and onto the disk, respectively. fall,BH and  fall,disk are the momentum accretions rates respectively associated to the BH and the disk.We evaluate the fraction of the infalling matter that directly falls into the BH by considering the competition between the specific angular momentum of the infalling matter and that of the innermost stable circular orbit (ISCO)  ISCO .More precisely, if the specific angular momentum is smaller than  ISCO , the infalling mass accretes onto the BH, if instead it is larger than  ISCO , it becomes a part of the disk.To determine  ISCO , we follow the prescription of Bardeen et al. (1972).We first evaluate the BH spin parameter : we then compute the ISCO radius,  ISCO , in terms of  following: where  1 and  2 are given by: We then define the specific angular momentum at the ISCO at the first order as follow: The different accretion rates are then estimated at the inner boundary  =  in as: where  (, ) is the specific angular momentum and Θ() is the Heaviside step function.
The mass and angular momentum transfer between the disk and the BH,  acc and  acc , are estimated as: with  acc the accretion time scale which is a free parameter in our models.Similarly, the contribution of the disk wind is evaluated as follows: where  w is the wind time scale and  disk =  disk / disk the specific angular momentum of the disk.Viscosity-driven mechanism is one of the possible mechanisms of the disk outflow.In this scenario, the magnetorotational instability results in the turbulent state in the disk, which acts as the effective viscosity (Balbus & Hawley 1991;Balbus & Hawley 1998).The viscous heating in the disk then drives the wind.There may be another origin of the viscosity: in the surface region of the disk, there is a velocity shear between disk matter and infalling envelope.This may induce the Kelvin-Helmholtz instability, which enhances the magnetic fields leading to the development of turbulence and dissipating the kinetic energy of the infalling matter (Obergaulinger et al. 2010).The wind time scale may be different for the different origin of viscosity.There may also be other mechanisms for launching the wind from the disk.For example, the magnetocentrifugal force by large scale magnetic fields can also work to launch the outflow (Blandford & Payne 1982).We, therefore, set  w as a free parameter not to specify the mechanism for the wind and to investigate more general central engines.

Inner boundary condition
Until the accretion disk forms, we apply an outflow condition at the inner boundary, thus allowing the material to inflow toward the central engine for  <  in .Once the disk is formed, i.e.  disk > 0, the wind outflow is injected from the inner boundary with the rate of equation ( 24) within an half opening angle  w .
In Fig. 1 we present the geometry used in our simulations.The outflow half opening angle at the inner boundary and its orientation can be freely chosen in our code.For this work we set  w = /4 and we direct it along the equatorial plane.Within the opening angle of 2 w , the wind density is set such that the following relation is satisfied: where  * 1 and  * 2 are the angles of the edges of 2 w , and  w is the radial velocity of the wind. 1 In this work we decided to describe the outflow density  w using a parabolic density profile defined as: with  0 being derived from the integration of equation ( 26).In equation ( 27), we set the parameter  = −1/cos 2 (/2 −  w ) so that the density is zero at the edges of the opening angle (i.e. at  1 = /4 and  2 = 3/4), and reaching maximum value  0 at  = /2, i.e. along the equatorial direction.We define the energy of the disk wind at the inner boundary having as a fraction of the energy related to the disk escape velocity  esc : where  esc is: with the disk radius  disk defined as: In equation ( 28) the internal energy of the wind is given by  int,w / w = (1/2)  therm  2 w . therm is a free parameter in our simulations, measuring the fraction of the wind kinetic energy assumed to be corresponding to its internal energy. is a fudge factor used to represent the uncertainties coming from the lack of knowledge of the precise disk structure (Hayakawa & Maeda 2018), and it is a free parameter in our simulations.The pressure of the outflow is computed using the tabulated EOS with the density  w and the internal energy  int,w as input parameters.equation ( 28) indicates that the asymptotic velocity of the injected matter, the velocity of the matter at infinite distance in the case the total specific energy (1/2) 2 +  int / + Φ is conserved, is  esc .
A part of the injected matter could fall back to the central engine, affecting the disk mass.This is the case when the ram pressure of the matter at the inner boundary is larger than that of the injected matter.
To avoid the recycling of the injected matter, in our study we set the angular momentum in the ghost cells to zero.In this way we do not allow the injected matter to fall back to the disk, but only to the BH.
For  <  * 1 or  >  * 2 , the boundary condition is set to prevent the matter from inflowing from the central engine to the computational domain.To achieve that we set zero fluxes (reflecting boundary condition) if the radial velocity in a first active cell is positive, while we allow the mass infall to the central engine if it is negative.

The equation of state
The thermodynamical properties of the star are described by a tabulated equation of state (EOS) that includes the ion, the radiation, the electrons, and  − - + pair.In this work we use an oxygen-based EOS, i.e an EOS using oxygen as the only component of the ion (i.e.  = 0.5), resulting in a 16 O mass fraction of 1.This decision is made considering the composition of our progenitor model dominated by oxygen outside the iron core (see Aguilera-Dena et al. 2020).

Diagnostics
In the following subsection, we describe the method used to calculate the properties of the ejecta and injected matter.
In our simulations, we define the ejecta mass  ej as the sum of unbound matter mass.The explosion energy  expl is the energy carried by the unbound matter.Several criteria are used to define the unbound matter in hydrodynamic simulations (citation here).For our study we chose to use the Bernoulli criterion which takes into account the thermal effect on the matter and the effects of the gravitational potential, and is defined as follows: Using the Bernoulli criterion we track the evolution of the ejecta mass2 and energy at every time step by integrating the equations: The injected mass  inj represents the matter coming from the central engine with a positive mass flux at the inner boundary  in .It is defined as: We consider the injected energy  inj as the energy carried by  inj with positive binding energy.We, then, compute the injected energy  inj applying the Bernoulli criterion as follows:

Parameters and initial condition
In this work, as mentioned above, we consider the scenario of a failed CCSN and we assume that the mass of the innermost region (1.28  ⊙ ) of our progenitor corresponds to the initial mass of the BH.In our model, the beginning of the disk formation is when the condition to launch the wind from the inner boundary is met for the first time.For our simulations, the density and velocity structure of the wind are fixed as explained in Sec.2.6.We then use four more free parameters, which are the wind time scale  w , the ratio between the accretion and wind time scales  acc / w , the ratio between the radial velocity of the outflow and the escape velocity  , and  therm which measures the fraction of the wind kinetic energy assumed to be corresponding to its internal energy (see equation ( 28)).In this work, we fix the direction of the outflow, its opening angle, and the density profile as we want to investigate the parameter space of the other quantities.Specifically, we set the wind along the equatorial plane using an half-opening angle  w of /4, and the density profile  w described in equation ( 27).
Using this setup, we investigate the parameter space for  w ,  acc / w ,  and  therm .We sample  w in a wide interval, (0.1, 1, 3.16, 10) s, using also more extreme values like 0.1 or 10 s (usually  w is few seconds as shown by Wang & Burrows 2023) to survey a parameter space as large as possible and to analyse the condition to reach the energies and the amount of 56 Ni produced in high-energy SNe. acc is set through the ratio  acc / w and the value of  w .We vary  acc / w in the interval (1,3.16,10,∞).The accretion time scale controls the accretion rate onto the BH from the disk, therefore it allows to track the dynamics of the central engine (see also Kumar et al. 2008).If  acc has small values so that it is shorter than the infalling time scale of the envelope (which is given by  fall,disk / disk ), then the accretion rate onto the BH tracks the rate at which mass is falling onto the disk.On the contrary, for longer  acc up to the extreme case of  acc = ∞ the mass infall onto the disk dominates.Varying  acc / w from 1 to ∞ allows us to investigate the effect of these two very different scenarios on the explosion and on the 56 Ni production.We assume the wind time scale and accretion time scale to be constant throughout the explosion in order to model the central engine as simple as possible.
In our simulations we use  2 = (0.1, 0.3) following the approach of Hayakawa & Maeda (2018) who used  2 = 0.1 in their work.We increase it because of our interest in the high energy explosions.
Finally, we set  therm as (0.1, 0.01) following the typical values of the wind internal energy in the literature (as in Hayakawa & Maeda 2018).In our work, we, then, test several combinations of these parameters.

Tracer particles and nucleosynthesis
To obtain thermodynamic histories of the ejecta, we use tracer particles following the method described in Fujibayashi et al. (2023b).In this method, the evolution of tracer particles is followed backward in time.Hence, they are placed every time interval Δ from the end of the simulation at radius  =  out in the range of 0 ≤  ≤ .The mass of each particle is defined as Δ =    ext 2 ΔΩΔ, where ΔΩ is the solid angle element.The time interval Δ is defined as as Δ :=  out Δ/⟨  ⟩, where Δ is the interval of the polar angle and ⟨  ⟩ is the average radial velocity of the ejecta at  =  out .This formulation of Δ ensures an optimal distribution of tracer particles in time.
This method is also utilized to judge whether a given fluid element is the injected matter from the inner boundary or the one that originates from the stellar envelope: A particle is tagged as an injected matter if it crosses the inner boundary during the back-tracing.On the other hand, if a particle stays inside the computational domain until the initial snapshot, it is an ejecta component originating from the stellar envelope.
A disadvantage of using a post-process particle tracing method is that the accuracy of the thermodynamics histories is limited by the frequency of the output (see, e.g., Sieverding et al. 2023).For most of the model, the time interval of the output is ∼ 70 ms.To check the systematic error that stems from this limitation, we performed a set of simulations using a higher time resolution of Δ ∼ 10 ms for selected models, and performed particle tracing using this time interval.This convergence study is presented in Appendix C2.
For tracer particles originating from the stellar envelope, the full thermodynamical history is available.To obtain the nucleosynthetic yield based on the density and temperature evolution along such particles, we perform nucleosynthesis calculation with the open-source code torch (Timmes et al. 2000) using 200 isotopes.The initial composition of the calculation is set by the stellar composition at the initial position of each tracer particle.Due to the lack of knowledge about the thermodynamical history of the injected matter (coming from the inner boundary), we perform nucleosynthesis calculation only for the stellar envelope.
Since we carry out the nucleosynthesis calculation without evolving the stellar composition from that of the original star (hence assuming always symmetric matter) and without taking into account neutrino interactions, the 56 Ni mass evaluated might be slightly higher than in reality.Yet, the nucleosynthesis calculation is performed only for those particles which are inside the computational domain for the whole simulation (i.e. at  >  in = 10 8 cm) where the neutrino interaction is not significant and hence not expected to strongly affect the 56 Ni production.

RESULTS
In this section, we will first describe the hydrodynamical evolution of the supernova explosion using one of our simulations as the characteristic model (subsection 3.1).Then, in the subsection 3.2, we will compare the results of some observables between our different models, using this to analyse the effect of  w ,  acc / w ,  and  therm on the explosion and 56 Ni production.

Hydrodynamical evolution
To present the general outline of the evolution of our simulations, we describe the model, M20_10_1_0.3_0.1, with parameters  w = 10 s,  acc / w = 1,  2 = 0.3 and  therm = 0.1 as example.In Fig. 2, we plot the evolution of the total mass, BH mass, disk mass, mass enclosed in our computational domain, and ejecta mass for this model.In the first ∼ 10 seconds, the infalling matter accretes only onto the BH because of the small angular momentum of infalling matter.The disk formation starts after about 10 s, which is the condition to trigger the wind injection in our simulation.A part of the disk mass accretes to the BH (equation ( 22)) and the other part contributes to the wind according to equation ( 24).The injection of the wind does not occur for the first 10-20 s due to the small ram pressure of the injected matter compared to that of the infalling matter.In this model it starts at  ≈ 20 s and it leads to the increase in the ejecta mass.The disk mass peaks at about 1  ⊙ , followed by a decrease because of the rate of the mass accretion into the BH and injection as the wind mass exceeds that of the supply to the disk from the stellar envelope.After ∼200 s, the mass components reach approximately constant values In cyan the total mass, in orange the mass enclosed in our computational domain, in green the disk mass, in red the BH mass and in purple the ejecta mass. total is plotted to show the conservation of mass throughout the simulation.and  BH becomes ∼ 11  ⊙ .The ejecta mass reaches a temporal maximum beyond 5  ⊙ after ∼ 32 s and then converges to 4.6  ⊙ .
In Fig. 3 we show the time evolution of the injected and explosion energies of the characteristic model, M20_10_1_0.3_0.1.It shows a wind injection beginning at ∼ 20 s and lasting ∼ 20 s.Within the first 40 seconds from the start of the simulation the wind has already been almost injected, reaching an energy  inj of ∼ 19.33 × 10 51 erg.In this model, we see that towards the end of the wind injection period,  expl reaches a temporal maximum, before plateauing after the wind injection is finished at  expl ∼ 12.30 × 10 51 erg.
The dynamics of the explosion of M20_10_1_0.3_0.1 can be followed in the left panel of Fig. 4 where we present the time evolution of the mass outflow rate Φ m averaged over the injection angle.The Table 1.Model description and key results.From left to right, the columns contain wind time scale, the ratio of the accretion and wind time scales, the squared ratio of the asymptotic velocity of injected matter to escape velocity of the disk, the internal to kinetic energy ratio of injected matter, cumulative injected energy, ejecta mass, explosion energy, average ejecta velocity, the mass of ejecta component that is originated from the computational domain and experienced temperature higher than 5 GK, the mass of the 56 Ni synthesized, the mass of ejecta component originated from the injected matter.angular-averaged mass outflow rate is given by: where  1 and  2 are the edges of the angle within which we average the mass outflow rate, in this case they limit the injection angle between /4 and 3/4, corresponding to  * 1 and  * 2 introduced in Sec.2.6.We focus here on the mass outflow rate averaged over the injection angle because it results to be almost equivalent to that averaged over the entire computational region (i.e., over [0 − ]).This indicates that the explosion is quasi spherical.The right panel of Fig. 4 displays the time evolution of the ratio between the averaged ram pressure of the expanding matter, Pexp ram , and that of the infalling matter, Pinfall ram .These are defined as: Pexp ram () = The left panel highlights the formation of a strong mass outflow at about the onset of the injection starting from the inner boundary and reaching the outer boundary.The positive mass rate dominates over the infalling envelope in the outer layers for the first ∼ 200 s of the simulation which corresponds to the time within which the ejecta mass converges.This region of the plot perfectly corresponds to that with the highest value of Pexp ram / Pinfall ram , in the right panel.This comparison suggests that the explosion in this model is driven by the injected wind which has a ram pressure larger than that of the infalling matter at the onset of the injection.The competition between the ram pressure of the wind and that of the infalling envelope can be analysed more in detail in the upper panel of Fig. 5, where we compare the ram pressure of the expanding (solid lines) and infalling matter (dotted lines) along the equatorial plane in this simulation at  = 19 s, soon after the injection has started.In the figure the expanding and infalling matter close to the inner boundary represents the wind component and the infalling envelope respectively.The upper panel confirms that after the injection starts, for M20_10_1_0.3_0.1, the ram pressure of the injected matter dominates.As a result, most of the injected matter can expand outward without falling back, leading to a highly-energetic explosion with ∼ 10 52 erg.
If the wind injection is weak, the injected matter is not able to efficiently push forward the stellar envelope determining a subenergetic explosion with < 10 51 erg.The dynamics of such explosion is presented in Fig. 6 using the model M20_1_1_0.3_0.1.This model has the same parameters as M20_10_1_0.3_0.1 ( acc / w = 1,  2 = 0.3 and  therm = 0.1) apart from the wind timescale which is  w = 1 s (and hence  acc = 1 s).In this simulation we measure lower injected and explosion energy, i.e.  inj = 3.08 × 10 51 erg and  expl = 0.088 × 10 51 erg.In the upper and bottom panels of Fig. 6 we show the time evolution of Φ m (left panels) and Pexp ram / Pinfall ram (right panels) averaged over the injection angle and outside that, respectively.The dynamics of the explosion of M20_1_1_0.3_0.1 looks different from that of M20_10_1_0.3_0.1 (see Fig. 4).In this case there is no positive Φ m dominating at all radii from the inner boundary to the outer boundary within the injection angle (see left panels of Fig. 6).This means that most of the injected matter cannot expand outwards without falling back.Indeed the upper-left panel of Fig. 6 shows that a negative averaged mass rate always dominates the innermost region around 2 × 10 8 cm, within the injection angle.This infalling mass stops the expansion of the injected matter.However a positive Φ m is present at larger radii.The bottom-left panel of Fig. 6 shows that the expanding component of the mass flux seems to dominate not only at  ≳ 5 × 10 8 cm, but also in the innermost region outside the injection angle, between ∼ 180 s and ∼ 700 s.Nonetheless, the expanding matter is blocked by the infalling matter at around 2 × 10 8 cm here as well.These regions of dominating positive (negative) Φ m in the left panels correspond to regions in which the ram pressure of the expanding (infalling) matter dominates in the right panels.This confirms that the dominating component of the ram pressure determines the direction of the muss flux, as suggested for M20_10_1_0.3_0.1.Therefore, since in M20_1_1_0.3_0.1 the ram pressure of the injected material close to the inner boundary seems to be on average always weaker than that of the infalling, contrary to what happens in M20_10_1_0.3_0.1, the explosion in this model is unlikely to be driven by the wind.The competition between the ram pressure of the injected wind and that of the infalling matter of this model is further investigated in the lower panel of Fig. 5.In this case even though we find some matter with positive velocity, the injected matter has a ram pressure smaller than that of the infalling envelope and some of it is found also at smaller radii, i.e. ∼ 10 8 cm, among the wind, limiting the explosion energy.
These results leads to the conclusion that if the ram pressure of the injected matter is smaller than that of the infalling matter, the explosion is determined by another mechanism: the infalling envelope bounces on the wind, launching the shock wave propagating outward.The outer layer of the star is then swept by the shock wave, being unbound.In this case, the energy source of the explosion is not the energy injected, but the released gravitational binding energy of bounced matter.Such a case is shown in the lower panel of Fig. 5 for the model M20_1_1_0.3_0.1.The plot shows that the ram pressure of the infalling envelope (dashed line) dominates over that of the injected matter (solid line).However the amount of the injected matter is sufficient to act as a "wall" causing the bounce of the infalling envelope, which occurs at around  ≲ 10 9 cm.In this figure the shock wave launched is also visible at larger radii, i.e. at  = 2 × 10 9 cm, propagating outwards and then sweeping the outer layer of the star.Since this unbound mass is located only at radii of tens of thousands of kilometer, where the density is small, also the amount of the unbound mass is small.At the same time the shock propagation decreases the infalling matter velocity and its ram pressure so that if the latter becomes sufficiently small, then the injected matter can move over it and go outwards becoming another ejecta components.However this happens after hundreds of seconds, when the energy budget -which is determined by the mass supply of the infalling envelope to the disk -is low determining also a very low  inj .
The comparison of the  ram competition between the two models shows that whether the ram pressure of the disk wind can overcome the ram pressure of the accretion flow or not determines a distinction between highly-energetic explosions and sub-energetic explosions.The efficiency of the explosion mechanism can also be measured by the ratio  expl / inj which indicates the fraction of the injected energy transferred to the ejecta.For M20_10_1_0.3_0.1, the explosion energy is the ∼ 64% of the injected energy, while in the case of M20_1_1_0.3_0.1,  expl / inj ≈ 0.029.
We find that some models with sub-energetic explosions have  expl >  inj .For example, model M20_1_1_0.1_0.01 has  expl ≈ 5 × 10 49 erg, while  inj < 10 49 erg (see Table .1).This can happen because the energy source of such sub-energetic explosion is different from that of the injected matter, it is the gravitational binding energy released by the bouncing infalling envelope (see Sec. 3.1).

Bimodality of 𝐸 expl
We plot the injected and explosion energies of all models with different parameters against respectively the ejecta mass (left panel of Fig. 7) and the average ejecta velocity  ej = √︃ 2 expl / ej (right panel of Fig. 7).The figures show that the models fall into two categories having different explosion energies, while such a distinction is not seen for the injection energy, which shows a rather continuous distribution.The first category is made of highly-energetic explosions, characterized by the explosion energy of about 10 52 erg.The second category is made of sub-energetic explosions with an energy of approximately 10 50 erg.This bimodal distribution seems to be mainly controlled by the wind time scale.Models with shorter  w (i.e. w ∼ 0.1 − 1 s) are located in the left, low-energy side of the plots, while those with longer  w (∼ 3.16 − 10 s) belong to the high explosion energy group, on the right of our plots.
Another feature distinguishing the highly-energetic and subenergetic explosions is the fraction of the injected energy converted in explosion energy which can be measured by the ratio  expl / inj .Both panels in Fig. 7 shows that the gap between  inj and  expl decreases with increasing  ej and  ej , indicating a correlation between  expl / inj and  expl .
To explain the bimodality of the explosion, we compare the evolution of two models belonging to the highly-energetic and sub-energetic categories respectively.Since M20_10_1_0.3_0.1, described in Section 3.1, lays in the right-hand side of both panels of Fig. 7, we compare it to the model M20_1_1_0.3_0.1 having the same parameters ( acc / w = 1,  2 = 0.3 and  therm = 0.1), but a different  w = 1 s (and hence  acc = 1 s) which is instead located on the left, low-energy side.In Fig. 8, we compare the time evolution of injected energy (upper panel) and explosion energy (lower panel) of the two models in the first 70 seconds of the simulation.Compared to the injected and explosion energies of M20_10_1_0.3_0.1 (see Sec. 3.1), M20_1_1_0.3_0.1 reaches only  inj = 3.08 × 10 51 erg and  expl = 0.088 × 10 51 erg (these values are also shown in Table 1).In the model with  w = 10 s the injected energy grows faster with a steeper slope after the onset of the injection.
It is clear from equation ( 24) that, for a given disk mass, the shorter wind time scale leads to the higher mass injection rate.The upper panel of Fig. 9 shows that the difference in the disk mass is always less than an order of magnitude in the first 70 s of the evolution.Since the wind time scales are different by a factor of ten for those two models, the mass injection rate in model M20_1_1_0.3_0.1 may be always larger than that in M20_10_1_0.3_0.1.Indeed, the lower panel of Fig. 9 shows that the mass injection begins slightly earlier in model M20_1_1_0.3_0.1 than the other model, which indicates that the wind power is stronger in the model in the earlier phase.However, the mass injection rate is somehow smaller in the later phase ( ≳ 20 s).This stems from the smaller escape velocity of the disk,  esc , as found in Fig. 10, especially in the later phase.The escape velocity plays an important role in determining the mass injection rate and affects it more than  w because the flux  wind defined in equation ( 24) is actually computed according to equation ( 26) solving the Riemann problem at the inner boundary. 3The smaller escape velocity leads to the smaller specific energy of injected matter, which is proportional to  2  2 esc .Hence, the wind is injected less efficiently in model M20_1_1_0.3_0.1 in the later phase.Therefore, the injected mass and energy saturate earlier.
As found in equation ( 29), which can be rewritten as  esc = √ 2  BH /  disk with equation ( 30), the larger specific angular momentum of the disk leads to smaller escape velocity.Using Eqs. ( 9), (11), and ( 22)-( 25), the evolution equation of the disk specific angular momentum is written as: where  fall :=  fall,disk /  fall,disk is the specific angular momentum of the infalling matter.The two terms of the equation are the contribution of the mass accretion onto the BH and the contribution of the infalling envelope respectively.Due to the contribution of the first term, if it dominates over the second, the disk specific angular momentum  disk =  disk / disk always increases because  disk >  ISCO .The second term does not always add a positive contribution to  disk since it can be also negative when  disk >  fall .If the absolute value of these negative contributions is smaller that the first term, then the accretion onto the BH dominates and  disk keeps increasing.In our simulations we find that the second term becomes also negative, i.e.  disk >  fall , however its absolute value is on average smaller that the first term.In this case, it is evident from the equation that the time scale of the increase in  disk is determined by  acc .Therefore, the decrease in  esc and thus the energy injection efficiency drop faster in model M20_1_1_0.3_0.1 than those in model M20_10_1_0.3_0.1.

Parameter dependence of 𝐸 expl and 𝑀 ej
In this section we analyse the effects that the parameters of our model have on the explosion.To do that we show in Fig. 11 the distribution of our models in the  expl - ej plane.We display our results (filled markers) together with the observational data for broad-lined type Ic SNe taken from Taddia et al. (2019b) and for stripped-envelope SNe and superluminous SNe from Gomez et al. (2022) (empty markers).We also show the results obtained by Fujibayashi et al. (2023a) who did a first-principled general relativistic neutrino-radiation viscoushydrodynamics simulation using the same progenitor model.The values they measured for the ejecta mass of  ej = 2.2  ⊙ and the explosion energy of  expl = 2.2 × 10 51 erg can be considered as the lower limits for the self-consistent simulation since they are still growing at the end of their simulation.The wind time scale of our models is indicated by the color of the marker, while the shape Both models have the same parameters apart from  w (see Table 1).Model M20_10_1_0.3_0.1 is represented using red lines, while model M20_1_1_0.3_0.1 is represented using blue lines.Lower panel: evolution of the explosion energy in the first 70 seconds of the simulation for the same two models.
distinguishes models with different  2 and  therm .The lines connect models with the same  w and  acc , but different  2 or  therm .
The first parameter we study is the wind time scale,  w , sampled in the interval (0.1, 1, 3.16, 10) s.Fig. 11 shows that most of the models with  w ≥ 3.16 s lay on the right side of the panel having  expl ≳ 1 × 10 51 erg and  ej ≳ 2.5  ⊙ , while all models with  w = 0.1 s occupy the lower-left side of the plot presenting  expl ≲ 10 51 erg and  ej ≲ 0.07  ⊙ .Longer wind time scales lead to higher explosion energy because they keep a larger escape velocity for longer time, as explained above and shown in Fig. 10.Only in models with  w = 1 s this parameter seems not to be predominant with respect to the  2 and  therm in determining whether an explosion in highly or subenergetic.This case and the difference with models having other  w will be discussed in a following paragraph, after analysing the general effects of the two other parameters  and  therm .
Considering  2 , the figure shows that a higher value of the parameter increases the explosion energy and the mass ejecta.This effect is evident by following the gray solid lines in Fig. 11 from the circular to the squared markers which represent models with all the same parameter but  2 = 0.1 and  2 = 0.3 respectively.The lines show that in all simulation an increase in  2 makes the point move towards the upper-right side of the plot, i.e. towards higher ejecta mass and explosion energy.This is consistent with our model of the wind (see equation ( 28)) in which we use  2 to set the asymptotic velocity of the injected matter as  esc .Hence, a higher  corresponds to a larger kinetic energy of the wind and enhances the energy budget for the explosion.Similarly,  expl and  ej increase by decreasing  therm .According to equation ( 28),  therm determines the ratio between the internal to the kinetic energies of the wind.Since the sum of these energies is provided by (1/2) 2  2 esc , reducing  therm increases the kinetic energy by decreasing the fraction of the thermal energy.
Observing the distribution of our numerical results, we note that, despite the variation of the parameters in wide ranges, the explosion of the  ZAMS = 20  ⊙ tend to remain located along a line and it does not spread to cover the space in the  ej −  expl as the observational data do.In other words, the correlation between the explosion energy and the ejecta mass is tighter in our simulations than that observationally measured by Taddia et al. (2019b) and Gomez et al. (2022).We also note that the results from our simulations presented here are in good agreement with that obtained in Fujibayashi et al. (2023a), especially considering that their values of  expl and  ej are still growing at the end of their simulation, so higher values were expected if they ran the simulation longer, which will be comparable to our  2022).The magenta plussign denotes the result obtained in a general relativistic neutrino-radiation viscous-hydrodynamics simulation with the same progenitor (Fujibayashi et al. 2023a).
results.This indicates that our model contains the case studied in Fujibayashi et al. (2023a) and explores various possibilities with a variety of wind properties.

56 Ni production
In this section, we then present the results of the 56 Ni production in our models and its dependence on the free parameters of the simulations to understand if it is possible to reproduce the observational data like those presented by Taddia et al. (2019b) and Gomez et al. (2022), especially for the high-energy SN ( expl > 10 52 erg).
The results of our nucleosynthesis calculations are presented in Table 1, where we list: the 56 Ni mass synthesized in the stellar component of the ejecta (the part of the ejecta originated from the computational domain),  stellar ej,Ni and the mass of ejecta component originated from the injected matter (i.e. from the disk),  inj ej .In Table 1 we also show the mass of the stellar component of the ejecta reaching temperature higher than 5 GK(= 5 × 10 9 K),  stellar ej,>5GK , as the 56 Ni production primarily occurs for  ≥ 5GK.
Although the injected matter explains 0.9 − 40% of the total ejecta mass, it is hard to accurately estimate the 56 Ni production for this component since the complete thermodynamical history is not available.We, therefore, estimate the upper limit of the mass of 56 Ni in the ejecta,  ej,Ni , considering that the injected matter entirely becomes 56 Ni:  ej,Ni =  stellar ej,Ni +  inj ej .It is found that it ranges from ∼ 0.02 to ∼ 2.09  ⊙ , which corresponds to ∼ 2.2 − 47% of the total ejecta mass.
Fig. 12 shows  ej,Ni and  stellar ej,Ni , the 56 Ni ejecta mass that originates from the stellar component, as a function of the explosion energy (left panel) and average ejecta velocity (right panel).Together with the results of our simulations, represented by the up and down-pointing filled triangles, we also display the observational data for broad-lined type Ic SNe taken from Taddia et al. (2019b) and for The yellow stars denote those obtained in a general relativistic neutrino-radiation viscous-hydrodynamics simulation with the same progenitor (Fujibayashi et al. 2023a).The open markers display the observational data for stripped-envelope SNe, some of which are broad-lined type Ic SNe, taken from Taddia et al. (2019b) and Gomez et al. (2022).
The simulation results highlight that, for most models,  inj ej is likely to dominate the total 56 Ni mass produced.They also show that the difference between  ej,Ni and  stellar ej,Ni is larger for higher explosion energies and, equivalently, for higher averaged ejecta velocities.This means that the value  ej,Ni we estimated is affected by the uncertainty of the thermodynamical history of  inj ej , and this is particularly true for the highly-energetic models.
Looking at Fig. 12 it is found that our numerical results for the 56 Ni mass reproduce the relation between  Ni and  expl or  Ni and  ej for very-high energy SNe with  expl > 2×10 51 erg and sub-energetic SNe with  expl < 0.1×10 51 erg of the observational data, suggesting that these SNe may be driven by a wind-driven explosion modeled as in this work.Moreover, Fig. 12 shows that also considering the 56 Ni mass produced, the models fall into two categories which correspond to those of highly-energetic and sub-energetic explosions observed in Fig. 7. Sub-energetic explosions produce a smaller amount of 56 Ni (< 0.2  ⊙ ) while highly-energetic explosions produce 0.2 − 2.1  ⊙ of 56 Ni.
Our estimate of the 56 Ni ejecta mass is in good agreement also with that obtained in Fujibayashi et al. (2023a), as it was for  expl and  ej (see Fig. 11), especially considering that their values of  expl and  Ni were still growing at the end of their simulation.
ej,>5GK can be a first indicator of the 56 Ni mass produced.It is, then, informative to compare  stellar ej,>5GK and  stellar ej,Ni .We limit this analysis to the stellar component of the ejecta as it is the matter for which we know the complete thermodynamical history.
This comparison is shown in Fig. 13 where the two quantities are displayed as a function of the explosion energy.In most of the models  stellar ej,>5GK is a good approximation of  stellar ej,Ni .For high energy ( expl ≳ 1.5 × 10 51 erg) the amount of matter experiencing T>5GK is larger than the 56 Ni mass, however, is not always the case, especially at lower energy ( expl ≲ 0.25 × 10 51 erg) when for many explosions we measure  stellar ej,>5GK <  stellar ej,Ni , which means that 56 Ni starts being produced already at  < 5GK.
The results of the 56 Ni mass are expected to be different if the wind is composed of lower electron fraction matter (i.e.,   ≪ 0.5).In this case the nucleosynthesis result is not supposed to peak at 56 Ni, but in heavier nuclei (Siegel et al. 2019) and the amount of 56 Ni produced in the injected component is not expected to be significant.Therefore the 56 Ni mass is not supposed to be dominated by  inj ej .

Variety of disk wind-driven explosions
The progenitor model used in this work is, as mentioned in Sec. 2, we a rapidly rotating, low metallicity, rotationally mixed quasichemically homogeneously evolving star with the zero-age mainsequence mass  ZAMS = 20  ⊙ taken from Aguilera-Dena et al. ( 2020), which has a very compact core, and it is supposed to fail the explosion (Ertl et al. 2016 andMüller et al. 2016).Our results show a large variety of explosion energies with  expl ranging from ∼ 0.049 × 10 51 erg to ∼ 34 × 10 51 erg (see also Table 1).Moreover, the distribution of these points in the plots can be divided into two categories: the sub-energetic explosions and the highly-energetic ones.In the plots in Fig. 7, every single point is the result obtained for the same structure of the progenitor in which we varied the parameters controlling the wind injection.It means that, in reality, a given structure of the progenitor results in only a single point in these plots.This means that, in reality, the 20 ⊙ star we employed results in either a highly-or sub-energetic explosion depending on the power of the wind injection and the competition between the ram pressure of the wind and the infalling envelope.It can be, then, interesting to study this competition for other massive progenitors with different structures to further investigate the variety of explosion properties and to verify if stars with other structures also present sub-and highly-energetic explosion branches depending on the power of the wind injection.The first parameter worthy of consideration is the progenitor mass.For instance, generally speaking, more massive stars tend to present a more compact CO-layer, and thus, have higher mass-infall rates leading to a larger amount of matter available to form the disk and to a higher thermal energy budget for the explosion energy.Therefore, employing such progenitors could allow us to model even more energetic SNe.
Another parameter that affects the fate of the core collapse and plays a role in the explosion after the formation of a BH is the rotation.The initial angular momentum profile of the progenitor has an impact on the explosion energy and the ejecta mass as shown by Fujibayashi et al. (2023a).Specifically, a star with fast rotation is expected to determine a more energetic explosion and to enhance the mass ejection.Then, we can also assume more 56 Ni will be produced.Therefore it can be informative also to perform simulations employing the same progenitor star but using different rotational profiles.
In order to investigate the variety of explosion properties depending on the progenitor structure and wind injection, it can be also interesting to use our central engine model for a failed CCSN for those progenitors that are more likely to undergo a neutrino-driven explosion according to Ertl et al. (2016) or Müller et al. (2016).Comparing the results obtained using our code for these two categories of progenitors (those supposed to fail the explosion and those expected to succeed in a neutrino-driven explosion) can be useful also to investigate the different dependencies of the explosion.Indeed Fujibayashi et al. (2023a) showed that the explosion in failed CCSN is sensitive to the mass-infall rate at the disk formation, i.e., a higher mass-infall rate (usually from the carbon-oxygen layer of the star) enhances the viscous and shock heating rates around the inner region of the disk, which result in the larger explosion energy, while Ertl et al. (2016) found that the neutrino-driven explosion is more likely to be sensitive to the compactness of inner domain, i.e. the iron-silicon layer.Hence, using the same general central engine model can allow us to verify these effects.
As mentioned in Sec.3.3, another interesting result of our study is that the correlation between  expl and  ej found in our models is stronger than that of observations (Fig. 11).This illustrates the incapability of our simulations to fully explain the variety of observational data by utilizing a single progenitor model.However, changing the mass and the rotational profile of the progenitor may fill the areas in Fig. 11 not covered by the present results.Performing explosion simulations using progenitors with different  ZAMS and rotational profiles is, thus, demanded to further investigate the correlation between the explosion energy and the ejecta mass.

The model of the wind injection
In this work, we adopt a simple prescription to set the wind by evolving the disk and the BH using Eqs.( 9)-( 12).In this model, the wind time scale and the accretion time scale are free parameters, which remain constant throughout the simulation.The flexibility in setting the parameters allows us to set engine models with diverse characteristics, unrestricted by specific scenarios.In Sec.3.3 and 3.4, we showed that the results of our models with  w ∼  acc ≳ 3.16 s reproduce  ej and  expl obtained by a general relativistic neutrinoradiation viscous-hydrodynamics simulation by Fujibayashi et al. (2023a) (see Fig. 11).Our models can also account for the 56 Ni mass measured in their work within the uncertainty of the mass fraction of 56 Ni of the injected component.These similar results can be explained by the similar time scales in the two simulations at the onset of the injection (i.e.,  ≲ 20 s).In the model used by Fujibayashi et al. (2023a)  w and  acc later evolve in time since they are evaluated from the viscous time scale that depends on the disk radius.However, this happens after the wind injection and they found that at early time in their 20  ⊙ model both  w and  acc are of several seconds, hence they are comparable to those used in our simulations with  w ∼  acc ≳ 3.16 s.Fujibayashi et al. (2023a) showed that in the viscosity-driven wind scenario, the viscous heating plays a role in determining the injected energy budget.After the efficiency of neutrino cooling drops, the outflow from the disk is driven by viscosity determining the evolution of the viscous time scale and hence of both  w and  acc .Therefore, focusing on this scenario, a future implementation of our model can be taking into account the effects of the time evolution of  w and  acc .This will allow us to provide more quantitative models for more sophisticated, although specific, scenarios in which we can describe the accretion flow during both the neutrino dominated accretion flow phase (NDAF) (Narayan & Yi 1994;Popham et al. 1999) and the advection dominated accretion flow phase (ADAF) (Narayan & Yi 1994, Kohri et al. 2005and Hayakawa & Maeda 2018).

The effect of GRB ejecta
As mentioned in Sec. 2, we exclude the central engine from the computational domain, and it is considered as being embedded in the central part of the star and characterized by the presence of a BH and a disk evolved according to Eqs. ( 9)-( 12).
If the dimensionless spin of the BH is large, in the presence of electromagnetic fields, the Blandford-Znajek effect could play an important role (see Blandford & Znajek 1977), i.e., it could launch an energetic jet or outflow along the spin axis of the BH.If a relativistic jet is produced, a gamma-ray burst will be also launched (see Izzard et al. 2004 andGottlieb et al. 2022 for simulation works).In the presence of the jet, more energy can be injected into the stellar matter, and hence, the energy budget available for the explosion and the 56 Ni production increases.Therefore performing relativistichydrodynamic simulations including the injection of relativistic jets will be one of our follow-up works.

SUMMARY AND CONCLUSIONS
We studied the hydrodynamics and nucleosynthesis for the explosion of a massive star to explore the properties of ejecta and the 56 Ni production in the collapsar scenario.Our main goal was to investigate the explosion mechanism of Type Ic SNe in the collapsar scenario.
We implemented a new feature that solves the cylindrically symmetric gravitational potential Φ to the open-source multidimensional hydrodynamics code Athena++.We used it to simulate the explosion of a rapidly rotating, rotationally mixed, quasichemically homogeneous  ZAMS ∼ 20  ⊙ star employing the progenitor model from Aguilera-Dena et al. (2020).For this work we also built a semi-analytical model for the central engine by taking into account the BH and disk evolution connected through matter and angular momentum transfer, to which we added the contribution of the disk wind following Hayakawa & Maeda (2018).
We tested different models by varying the parameters that control the properties of mass and energy injection to thoroughly investigate their influence on the final ejecta.The parameter setups of the simulations and the main results are listed in Table 1.
In all of our models, we found that the energy and mass injection occurs roughly between 10 and 20 s after the disk formation.After the wind injection, the competition between the ram pressure of the injected and infalling matter leads the disk wind-driven explosion to be sub-or highly-energetic with an explosion energy < 0.1×10 51 erg and > 10 52 erg respectively.This distinction originates from whether the  ram of the injected matter can overcome  ram of the infalling envelope and efficiently push the stellar envelope outwards.When the first one is larger than the second one, most of the matter can expand outwards without falling back, leading to a highly-energetic explosion with ∼ 1 × 10 52 erg.In the case of the sub-energetic explosions ( expl < 1 × 10 51 erg), instead, the shock wave transfers the energy from the infalling matter bouncing on the wind.Propagating then outwards, the bounce shock causes the ejection of the matter outside the star.
Studying the impact of the parameters on the final ejecta, we found that the wind timescale strongly affects  expl .In particular, models with longer wind time scales tend to reach higher explosion energies due to a high mass-infall rate.We also noticed that models with higher  2 (i.e. 2 = 0.3), and hence, higher wind kinetic energy, or smaller  therm (that is  therm = 0.01) represent highly-energetic explosions.In contrast, smaller values of  2 or larger  therm lead the explosions to be sub-energetic.
We found that sub-energetic explosions produce smaller amounts of 56 Ni (< 0.2  ⊙ ) while highly-energetic explosions produce 0.2 − 2.1  ⊙ of 56 Ni.The 56 Ni mass was evaluated separately taking into account the stellar and the injected components of the ejecta since the whole thermodynamical history of the particle tracer is only available for the former component.Our results show that  inj ej dominates and the difference between  stellar ej,Ni and  stellar ej,Ni +  inj ej is larger for highly-energetic explosions.
We also compared our numerical results with the observational data for stripped-envelope SNe, some of which are broad-lined type Ic SNe, taken from Taddia et al. (2019b) and Gomez et al. (2022).We found that the distribution of the 56 Ni mass of our models reproduces the relation between  Ni and  expl or  Ni and  ej for very-high energy SNe with  expl > 2 × 10 51 erg and sub-energetic SNe with  expl < 0.1 × 10 51 erg of the observational data.Moreover, we measured a tighter correlation between the explosion energy and the ejecta mass in our simulations than that observationally measured by Taddia et al. (2019b) and Gomez et al. (2022).
To better investigate the variety of explosion properties and to verify whether stars with different structures present sub-and highlyenergetic explosion branches, we plan to perform numerical simulations by varying the mass and the rotational profile of the progenitor in our follow-up work.More massive stars can, for instance, have a larger amount of matter to form the disk and hence a higher energy budget for the explosion energy.Moreover, Fujibayashi et al. (2023b) showed that the initial angular momentum profile of the progenitor affects the explosion energy and the ejecta mass.Therefore we expect that varying the progenitor mass and rotational profile may explain the variety of observational data.
Finally, in this work, we found that our models with  w ∼  acc of several seconds reproduce the results for  ej ,  expl and the 56 Ni mass obtained by Fujibayashi et al. (2023a) in a general relativistic neutrino-radiation viscous-hydrodynamics simulation that utilizes the same progenitor model.This is due to the similar time scales in our and their works at early times (i.e.,  < 20s), at the onset of the injection.In their model for the viscosity-driven wind scenario, however,  w and  acc later evolve in time with the viscous time scale.However this happens after the wind injection which is the moment that we model in this work and we are focusing on.At this early time  w and  acc in our work and that of Fujibayashi et al. (2023a) are similar.Considering the model used by Fujibayashi et al. (2023a), it may be interesting implement the evolution of these time scales in our model and use it to investigate the NDAF and ADAF phases in this specific scenario.
it is not part of the computational domain.Therefore we assume a spherically matter distribution in that central region (with  <  Since we are interested in computing the potential inside the computational domain in which  >  in , then  0 ( ′ , ) is always given by  0 ( ′ , ) =  ′2 /.Therefore, equation (A11) becomes:

APPENDIX B: MODEL OF THE DISK WIND
If we consider the evolution if the specific angular momentum of the disk  disk as described in equation ( 39), we can notice that if the contribution of infalling matter is small, the average specific angular momentum of the disk increases with time.As a consequence the disk radius, which is defined as  disk  2 disk /(  BH ), can become larger than the inner boundary radius.Indeed if we assume that the wind carries the average specific angular momentum of the disk and we require the angular momentum conservation, considering a uniform angular velocity,  wind is evaluated as: Since our model does not describe the injection from a disk with  disk >  in , we cannot describe such a system consistently.So we describe the disk outflow according to equation (28) considering that the injected matter does not carry angular momentum.

APPENDIX C: CONVERGENCE STUDIES C1 Resolution study
The resolution of the computational domain we use in this work has been chosen considering its effect on the duration of the simulations and the convergence of the results.In order to verify that (  ,   ) = (128, 220) is a sufficient discretization of our domain, for M20_10_1_0.3_0.1 we also performed a simulation increasing the resolution to (  ,   ) = (256, 460).Fig. C1 shows the disk mass and the injected mass measured in the two simulations.The difference of these quantities in the two runs ranges from ∼ 0.01% to ∼ 2.9%, confirming that using resolution of (  ,   ) = (128, 220) we obtain converged results.The good agreement of the two simulations is also visible comparing the evolution of the explosion energy, as done in Fig. C2.Table C1.Model description and key results for models an increased frequency (every 10 ms) of output.From left to right, the columns contain wind time scale, the ratio of the accretion and wind time scales, the squared ratio of the asymptotic velocity of injected matter to escape velocity of the disk, the internal to kinetic energy ratio of injected matter, cumulative injected energy, ejecta mass, explosion energy, average ejecta velocity, the mass of ejecta component that is originated from the computational domain and experienced temperature higher than 5 GK, the mass of the 56 Ni synthesized, the mass of ejecta component originated from the injected matter.

C2 Dependence of particle tracing on time interval of outputs
As mentioned in Section 2.10, the accuracy of the thermodynamical history evaluated in the post-process particle tracing is strongly affected by the frequency of the output.For most of the model we use a time interval of the output of 70 ms.We check the systematic performing the particle tracing with an output frequency increased to every 10 ms for selected model.The results obtained for  stellar ej,>5GK and  stellar ej,Ni with the increased frequency are displayed in parenthesis Table C1.By comparing them with the results obtained for the same models but with lower output frequency listed in Table 1, they show that the uncertainty on the accuracy of the particle tracing is of 10% level.

Figure 1 .
Figure 1.Schematic picture of the explosion in the collapsar scenario.2 w represents the angle for which we allow the wind outflow.Outside of this angle, the matter is only allowed to infall towards the central engine.The figure also shows the rotation axis and the equatorial plane.

Figure 2 .
Figure 2. Masses evolution for the characteristic model , M20_10_1_0.3_0.1.In cyan the total mass, in orange the mass enclosed in our computational domain, in green the disk mass, in red the BH mass and in purple the ejecta mass. total is plotted to show the conservation of mass throughout the simulation.

Figure 4 .PPFigure 5 .
Figure 4. Left panel: time evolution of the mass outflow rate Φ m averaged over the injection angle (i.e., in [ /4 − 3/4  ]).Right panel: space-time diagram of the ratio between the ram pressure averaged over the injection angle of the injected matter, Pexp ram , and of the infalling envelope, Pinfall ram .These plots are obtained for the model M20_10_1_0.3_0.1.

Figure 6 .
Figure 6.Upper-left panel: time evolution of the mass rate Φ m averaged over the injection angle.Upper-right panel: time evolution of the ratio between the ram pressure averaged over the injection angle of the injected matter, Pexp ram , and of the infalling envelope, Pinfall ram .Lower panels: they show the same time evolution diagrams as the upper panels but here the mass flux and the ram pressure are angled-averaged outside the injection angle, i.e., in [0 − /4] + [3/4  −  ].These plots are obtained for the model M20_1_1_0.3_0.1 (lower panels).

Figure 7 .Figure 8 .
Figure 7.The  expl (filled markers) and the  inj (open markers) against the ejecta mass  ej (left panel) and against the ejecta velocity  ej (right panel) for the models studied in this work.Results for models with different wind time scales  w = 0.1s, 1s, 3.16s and 10s are distinguished by different colours, respectively red, blue, cyan and green.

Figure 9 .Figure 10 .
Figure 9. Upper panel: evolution of the disk mass in the first 70 seconds for the simulation of the same models as in Fig. 8. Lower panel: evolution of the injected mass in the first 70 seconds of the simulation for the same two models.

Figure 11 .
Figure11.Parameter dependence with respect to the observable pair of ejecta mass  ej and explosion energy  expl .The color distinguishes the wind time scale  w .The lines connect models with all the other parameters fixed but different  2 or  therm .The open markers display the observational data for stripped-envelope SNe, some of which are broad-lined type Ic SNe, taken fromTaddia et al. (2019b) andGomez et al. (2022).The magenta plussign denotes the result obtained in a general relativistic neutrino-radiation viscous-hydrodynamics simulation with the same progenitor(Fujibayashi et al. 2023a).

Figure 12 .
Figure 12.Relations between the explosion energy and the 56 Ni mass (left) and average velocity of the ejecta and the 56 Ni mass (right).Each grey line connects  stellar ej,Ni (triangles) and  stellar ej,Ni +  inj ej (down-pointing triangles) of the same model to show the possible range of 56 Ni mass.The yellow stars denote those obtained in a general relativistic neutrino-radiation viscous-hydrodynamics simulation with the same progenitor(Fujibayashi et al. 2023a).The open markers display the observational data for stripped-envelope SNe, some of which are broad-lined type Ic SNe, taken fromTaddia et al. (2019b) andGomez et al. (2022).

Figure 13 .
Figure 13.Comparison between the  stellar ej,Ni (blue triangles) and  stellar ej,>5GK(orange down-pointing triangles) displayed against the explosion energy.Each grey line connects  stellar ej,Ni and  stellar ej,>5GK of the same model.

Figure C2 .
Figure C2.Evolution of the explosion energy of M20_10_1_0.3_0.1 at two different resolutions: The red solid line shows the results for the resolution used for all the simulations in this work (  ,   ) = (128, 220) and the dashed, blue lines displays the quantities for the higher resolution of (  ,   ) = (256, 460).
The results are for an output frequency of 70 ms.