A Study of Primordial Very Massive Star Evolution. II. Stellar Rotation and Gamma-Ray Burst Progenitors

We calculate new evolutionary models of rotating primordial very massive stars, with initial mass from 100 M ⊙ to 200 M ⊙, for two values of the initial metallicity Z = 0 and Z = 0.0002. For the first time in this mass range, we consider stellar rotation and pulsation-driven mass loss, along with radiative winds. The models evolve from the zero-age main sequence until the onset of pair-instability. We discuss the main properties of the models during their evolution and then focus on the final fate and the possible progenitors of jet-driven events. All tracks that undergo pulsational-pair instability produce successful gamma-ray bursts (GRB) in the collapsar framework, while those that collapse directly to black holes (BH) produce jet-driven supernova events. In these latter cases, the expected black hole mass changes due to the jet propagation inside the progenitor, resulting in different models that should produce BH within the pair-instability black hole mass gap. Successful GRBs predicted here from zero metallicity, and very metal-poor progenitors, may be bright enough to be detected even up to redshift ∼20 using current telescopes such as the Swift-BAT X-ray detector and the JWST.


Introduction
The first stars (Population III) that lit up in our Universe might have been far more massive than those forming nowadays, with an initial mass function peaking at ; 100 M e (Bromm et al. 1999;Abel et al. 2002).The reason is the absence of metals, which are the most efficient coolants within molecular clouds during star formation.In a previous work (Volpato et al. 2023, hereafter, Paper I, to which we refer for an introduction on Population III and very metal-poor stars), we investigated the effect of pulsation-driven mass loss (see the prescription in Nakauchi et al. 2020) on the evolution and final fate of primordial very massive stars (VMSs).We found that pulsation-driven and radiation-driven mass loss dominates during different evolutionary phases.In Paper I, the most massive stars can eject several solar masses of material already during the main sequence (MS) due to radial pulsations.This is in contrast with the modest mass loss expected from radiationdriven winds at these metallicities.We find that almost all models experience pair-creation instability during the last phases of the evolution after the ignition of carbon, neon, or oxygen, depending on the mass of the star.Models with M i = 100 M e and Z = 0 are somewhat uncertain and may avoid pair-instability.Our stars with 300 M i /M e 1000 should directly collapse into a black hole (DBH, Bond et al. 1984;Farmer et al. 2020), while models with M i = 150, 200 M e should produce pair-instability supernovae (PISN), leaving no remnant (Fowler & Hoyle 1964;Barkat et al. 1967;Rakavy & Shaviv 1967;Fraley 1968;Heger & Woosley 2002;Heger et al. 2003;Takahashi et al. 2018;Takahashi 2018;Woosley & Heger 2021;Farag et al. 2022;Costa et al. 2023).Models with M i = 100 M e and Z = 0.0002 should enter into the pulsational pair-instability supernova (PPISN) regime (Woosley et al. 2002;Chen et al. 2014;Yoshida et al. 2016;Woosley 2017;Farmer et al. 2019;Woosley & Heger 2021;Farag et al. 2022).Instead, those with the same mass and Z = 0 could end their evolution either with a failed core-collapse supernova (CCSN) or a PPISN.In the former case, these models could provide a new formation pathway for black holes (BHs), potentially helping to alleviate the black hole mass gap puzzle (see also Croon et al. 2020;Farmer et al. 2020;Sakstein et al. 2020;Costa et al. 2021;Farrell et al. 2021;Tanikawa et al. 2021;Vink et al. 2021;Farag et al. 2022, for different formation scenarios).
In Paper I, we did not consider stellar rotation, which is one of the most influential phenomena in the evolution of massive and very massive stars (Heger et al. 2000;Meynet & Maeder 2000;Brott et al. 2011;Ekström et al. 2012;Paxton et al. 2013;Yusof et al. 2013;Limongi & Chieffi 2018;Goswami et al. 2021Goswami et al. , 2022;;Higgins et al. 2022;Martinet et al. 2023).Rotation changes the gravity of the star and thus modifies the stellar geometry from the usual spherical symmetry also affecting the surface temperature (e.g., Costa et al. 2019a, 2019b, andreferences therein).The fact that the star is rotating implies different kinds of turbulent mixing, e.g., meridional circulation and shear instability.These processes increase the mixing of chemical elements within the star, affecting the lifetimes of the main nuclear burning phases as well as the chemical composition at the star's surface (Maeder 2009, and references therein).Turbulence caused by rotation transports angular momentum from the stellar core to the envelope (Heger et al. 2000), where it is removed by stellar winds.Different authors (e.g., Heger et al. 2000;Georgy et al. 2011, and references therein) have suggested that mass loss is enhanced by rotation due to lower effective gravity, which is caused by centrifugal forces.However, the interplay between rotation and mass loss is still under debate, with different results reported in the literature.For instance, more sophisticated multidimensional radiation-driven wind models have not confirmed the enhancement of mass loss due to rotation (Müller & Vink 2014).
While the effect of rotation on mass loss in VMSs has been the subject of several studies (e.g., Meynet & Maeder 2000;Ekström et al. 2012;Yoon et al. 2012Yoon et al. , 2015;;Murphy et al. 2021;Martinet et al. 2023), its impact in combination with pulsation-induced mass loss has never been examined.In light of the results of Paper I, this topic clearly deserves to be investigated.We do so using the PARSEC code (Bressan et al. 2012;Costa et al. 2019aCosta et al. , 2019b;;Nguyen et al. 2022) to follow the evolution of a set of VMSs until the occurrence of pairinstability, in line with our previous work.
Rotation can also heavily affect the final fate of massive and very massive stars, and it is a necessary condition to produce successful GRB events (Woosley 1993;Woosley & Heger 2006;Yoon et al. 2006Yoon et al. , 2012Yoon et al. , 2015;;Woosley & Heger 2012).The collapsar model (Woosley 1993) is the most widely accepted theory for the formation of long GRBs, in which an accretion disk forms during the collapse of a massive star, powering the jets that produce the GRB event.The accretion disk can form only if the infalling material has enough angular momentum to avoid direct accretion onto the BH (Woosley 1993).The propagation of the jet through the progenitor is a major issue for massive stars with extended envelopes, preventing the production of a successful gammaray burst (GRB) event.The reason is that the crossing timescale for the jet to reach the surface of the star can be longer than the accretion timescale by more than 3 orders of magnitude.In this case, the jet is not powerful enough to break out from the star and produce a successful GRB (Yoon et al. 2006(Yoon et al. , 2012(Yoon et al. , 2015;;Woosley & Heger 2012).This is why, when considering Population III stars, Yoon et al. (2012) proposed the chemically homogeneous evolution as the main channel for the star to retain enough angular momentum within its core and at the same time to avoid the redward evolution in the H-R diagram.In this way, the star does not retain a very extended envelope and this facilitates the jet propagation.In this work, we investigate models of VMS that should experience pairinstability (PI) but with lower masses compared to Yoon et al. (2015) to have less extended envelopes.We find another possible pathway for the evolution of successful GRB progenitors.The absence of the stellar envelope due to pulsational pair-instability mass loss eases the jet propagation through the star.This decreases enormously the crossing timescale for the jet to reach the stellar surface, thus producing successful GRB events.
The paper is organized as follows.Section 2 summarizes the main ingredients of the PARSEC code.In Section 3 we present stellar models computed with and without rotation accounting only for radiation-driven winds and for both radiation and pulsation-driven mass loss.We discuss the main aspects of the models, especially how dredge-up episodes affect the internal structure, the surface abundances, and the wind ejecta.Then, we focus on the final fates, the compact remnants, and the possible jet-driven events for different models' structures and angular momentum configurations.In Section 4 we summarize our results and conclude.

Stellar Evolution Calculations
We compute stellar evolution calculations with the PARSEC code (Bressan et al. 2012) in its version 2.0 (Costa et al. 2019a(Costa et al. , 2019b(Costa et al. , 2021;;Nguyen et al. 2022, and reference therein).The input physics and the code parameters are as described in Paper I, whereas the details of the implementation of rotation in the PARSEC code can be found in Costa et al. (2019aCosta et al. ( , 2019b)).We use scaled-solar abundances (Caffau et al. 2011), where the initial solar metallicity follows the calibration in Bressan et al. (2012), Z initial,e = 0.01774.Here, we briefly summarize the main aspects of the PARSEC code that concern the evolution of massive stars and, in particular, the mass loss prescriptions adopted in this work.For the standard mass loss prescription (M rdw  for consistency with Paper I), we adopt the formulation by Vink et al. (2000Vink et al. ( , 2001) ) and the mass loss rates predicted by de Jager et al. (1988) for stars with an effective temperature higher and lower than 10,000 K, respectively.The metallicity dependence is µ , where we use the initial metallicity as a proxy for the iron content6 (which remains constant along the evolution).We also consider the enhancement of mass loss due to the proximity of the Eddington factor to 1 (Gräfener & Hamann 2008;Vink et al. 2011).In this case, we use the same scheme as in Chen et al. (2015), with the following details µ  , where α = 2.45 − 2.4Γ e , with Γ e the Eddington factor and α between 0 and 0.85.In each evolutionary stage, we take the maximum mass loss rate between Vink et al. (2000Vink et al. ( , 2001) ) or de Jager et al. (1988) and Vink et al. (2011), which accounts for the Eddington factor dependence.For Wolf-Rayet stars with X < 0.3 and > ( ) T log 4 eff , we use the mass loss recipe from Sander et al. (2019) with the metallicity dependence proposed by Costa et al. (2021).
For pulsation-driven mass loss (M pdw  ), we use the analytical expressions provided by Nakauchi et al. (2020) where α 1 , α 2 , β 1 , β 2 , and γ are coefficients that depend on the initial metallicity of the model (see Nakauchi et al. 2020, for more details).Finally, the last mass loss prescription is M max  , which takes the maximum between M rdw  and M pdw  at each time step during the evolution of the models.
Concerning the interplay between mass loss and stellar rotation, we use the prescription by Heger et al. (2000), which reads is the mass loss rate in the nonrotating case.Then, v is the surface tangential stellar velocity, while v crit is the break-up velocity defined as We take into account also the mechanical mass loss when the star reaches the critical rotation (Georgy et al. 2013).It is worth noticing that there are different prescriptions for the treatment of stellar rotation (e.g., Maeder & Meynet 2000), and this subject is still under investigation.Major uncertainties in stellar rotation are tied to different aspects.For instance, stellar rotation is inextricably linked to magnetic field generation (e.g., Braithwaite & Spruit 2017;Brun & Browning 2017).The interaction between rotation and magnetic fields is not completely understood despite its influence on stellar activity, which affects processes such as starspots, flares, and the stellar wind.Furthermore, rotation influences mass and angular momentum loss from stars, thus affecting their evolution.Predicting the rate and pattern of mass loss due to rotationinduced instabilities is still a source of uncertainties in evolutionary models.
We adopt the mixing-length theory by Böhm-Vitense (1958), with a mixing-length parameter α mlt = 1.74.We use the Schwarzschild criterion to define the border of the convective regions, while for core overshooting, we adopt the ballistic approximation (Bressan et al. 1981).In this latter, the core overshooting parameter (λ ov = 0.4) times the pressure scale height corresponds to the mean free path of the eddies across the border of the convective region.We also account for overshooting at the base of the convective envelope below the formal Schwarzschild border (Alongi et al. 1991;Bressan et al. 2012;Nguyen et al. 2022).Recent calibration of the red giant branch bump luminosity in a large sample of Globular Clusters with metallicity −2 < [M/H] < 0 using updated α-enhanced PARSEC models (Fu et al. 2018) suggests using an envelope overshooting parameter Λ env between 0.5 and 0.7, with the latter value being more appropriate in more metal-poor systems.Similar, if not greater, values are required to reproduce blue loops in star clusters and in low metallicity dwarf irregular galaxies (Alongi et al. 1991;Bressan et al. 2012;Tang et al. 2014).In the present calculations, we use Λ env = 0.7.
To inhibit density inversion in the inefficient convective regions of the stellar envelope, we follow the temperature gradient limitation described in Chen et al. (2015).
, convection becomes more efficient, preventing the numerical instabilities caused by density inversion.
As in Paper I, we consider two different initial chemical compositions (Z = 0, Y = 0.2485) and (Z = 0.0002, Y = 0.24885), with Z and Y the initial abundances in mass fraction of metals and helium, respectively.We compute stellar evolution models with initial mass M i = 100, 150, 200 M e and initial rotation rate ω = 0.0, 0.2, 0.3, 0.4, 0.5; where ω = Ω i /Ω crit with Ω i the initial angular velocity and Ω crit the critical angular velocity.For each combination of (M i , Z, ω), we compute two stellar models adopting two different mass loss recipes, namely M rdw  and M max  .We found that models encounter progressively greater numerical difficulties in the computation toward the highest values of mass, rotation velocity, and metallicity in the explored range.In particular, two models out of 60 could not be brought to convergence (the models with M i = 150 M e and M i = 200 M e , having Z = 0.0002, ω i = 0.5, and computed with the M rdw  prescription), so we exclude them from the following discussion.
We note that the calculations of Nakauchi et al. (2020) are based on nonrotating models, while in the present study, the interplay between rotation and stellar pulsation should be taken into account.However, we found the timescale of pulsation to be always much shorter than the rotation period in all our models, indicating that the interaction between the two processes can be safely neglected (see Appendix C for further details).
We follow the evolution of our models from the zero-age main sequence (ZAMS) until the onset of pair-creation instability (see Section 2 in Paper I).This occurs after the ignition of carbon, neon, and oxygen in the stellar core, depending on the initial mass, metallicity, rotation, and massloss prescription adopted for the models (see Tables 1 and 2).

General Properties of the Stellar Evolutionary Tracks
In Table 1 and 2 we summarize the main properties that highlight the evolution and final outcome of the models.
Figure 1 presents the evolution of all models in the H-R diagram (HRD), where we can see that there is a positive correlation between the models' luminosity and their initial rotation velocity, ω.This is most evident for models computed with Z = 0.0002, while in the case of Z = 0 the evolutionary tracks run almost superimposed except for the nonrotating ones.This behavior can be explained by the helium enrichment at the surface during MS, which forces both luminosity and effective temperature to increase.In turn, this helium enrichment is caused by two factors.First, the enhancement of the convective core between rotating and nonrotating models increases with metallicity (Groh et al. 2019).Second, at Z = 0.0002, the region above the convective core experiences greater rotational mixing as the rotation rate increases, compared to the Z = 0 models at the same rotation rates.
Another effect of rotation is to reduce or even quench the blue loops during the core helium burning (cHeB) phase.We can see that in all panels of Figure 1, there is at least one nonrotating star that evolves toward higher effective temperatures.With the addition of rotation, this is not the case anymore.Models with Z = 0 and ω > 0 evolve after the MS toward the red part of the HRD; while for Z = 0.0002 there are two models with ω = 0.2 that perform a blue loop, but in the case of ω > 0.2 no model evolves back to higher effective temperatures.
After the MS phase, due to rotational mixing and the occurrence of dredge-up (DUP) episodes, 16 stellar models reach a surface hydrogen abundance X < 0.3 (starred symbols in Figure 1).When massive stars evolve into red supergiants, the convective envelope inflates and cools, while at low densities, opacity is dominated by electron scattering.These lead to increasing atmospheric opacity and favor the development of convection in progressively deeper layers of the star, causing a DUP episode.It is worth noticing that the efficiency    of the DUP depends on the envelope undershooting parameter, as shown in Costa et al. (2021).According to our definition (see Section 2.1 in Paper I), these 16 stars exhibit a surface chemical composition similar to Wolf-Rayet (WR) stars.Still, they are not hot and spend most of their evolution with a low effective temperature in the red part of the HRD.For this reason, we refer to this kind of star as WR-manqué (WRm; note that the mass loss recipe adopted for these stars is that by de    We plot the value of initial mass (in M e ) for all tracks.Jager et al. 1988).This drastic change in the surface chemical composition occurs too late to have a great impact on the models' effective temperature evolution.If the hydrogen surface abundance were to decrease below ∼20% during MS (due to high rotation mixing), then the models would have followed the so-called chemically homogeneous evolution.In this case, they would have evolved toward higher effective temperatures, completely avoiding the red-supergiant phase (Woosley & Heger 2006;Yoon et al. 2012).We do not find any occurrence of chemically homogeneous evolution in our models.Among these 16 tracks, six stellar models reach a total amount of hydrogen between 0.26% and 0.1% of their total mass.Section 3.3 examines the impact of dredge-up and rotation on surface abundances and chemical ejecta, especially for these extreme WRm stellar models.

Internal Structure
Figure 2 shows the Kippenhahn diagrams of four different stellar models.They have the same initial mass M i = 200 M e , but different metallicities, rotational velocities, and mass loss prescriptions.At the beginning of the MS, all these stars have a convective core equal to almost 90% of their total mass.This is because these models have a very large initial mass, and rotation increases the extent of their convective cores.For example, comparing the models in panels (a) and (c), we can notice that the star in this latter panel has a slightly bigger convective core due to its higher initial rotational velocity.Then, the model in panel (a) is evolving almost at constant mass.The reason is that for this model, we adopted the massloss recipe M rdw  with Z = 0. On the other hand, for models in panels (b), (c), and (d) the total mass of the star decreases during the evolution.This is most evident for the stars in panels (b) and (d) during the first part of the MS, when the pulsationdriven mass loss is operating.
After the depletion of hydrogen in the stellar core and the start of the cHeB phase, these four models experience DUP episodes.In Figure 2 there are different degrees of envelope penetration, which imply different internal mixing.DUP episodes affect the structure of these models by reducing the mass of the helium core, which can be crucial for the final fate of the stars (see Section 3.4).They also impact the chemical composition of these four models.Depending on the efficiency of the DUP, these stars can remain classical red supergiants, or they can become WRm stars.Despite the DUP, the model in panel (a) retains enough hydrogen at the surface such that it does not become a WRm star, while the opposite occurs for models in panels (b), (c), and (d) (see also Section 3.3).The red vertical line shows when X < 0.3 at the surface, and we see that this does not occur for the model in panel (a).Moreover, the ) model in panel (b) becomes an extreme WRm star, with only 0.11% of hydrogen out of its total final mass.Therefore, this star could be considered a pure helium star with a very small carbon-oxygen core (green line in Figure 2(b)) due to the very deep penetration of the stellar envelope.
The surface abundance evolution and wind ejecta masses of all extreme WRm models, along with other selected tracks, are discussed in the Section 3.3.

Surface Chemical Abundances and Mass of the Ejecta
In this section, we consider the surface abundances and the wind ejecta mass of some selected models.Both of these are affected by rotation-induced mixing, dredge-up episodes, and mass-loss history.
Figure 3 presents the surface abundance evolution of H, He, C, N, O, and Z eff = 1 − X − Y.Each panel shows the results for the same model computed with two different mass-loss prescriptions, M rdw  (dotted line) and M max  (solid line).In Figure 3 there are all six extreme WRm models.
During MS, we see the effect of rotation-induced mixing, especially in the four Z = 0 models (panels (a), (b), (c), and (d)).That is, for example, the gradual increase in nitrogen surface abundance, along with a slower enhancement in carbon and oxygen.Approximately ∼10 5 yr before the end of computations, all six models start to experience DUP episodes.During cHeB, the hydrogen surface abundance lowers below 0.3 (vertical red line in each panel), while the effective metallicity increases very steeply.At the end of computations, all extreme WRm models in Figure 3 have Z eff > 0.6, reaching ∼0.77 in panel (d) for the highest initial rotational velocity.The interesting thing to notice is that at Z = 0, the mass-loss prescription accounting also for pulsation-driven mass loss favors the metal enrichment during the evolution, while at Z = 0.0002 the extreme WRm models are computed with M rdw  .
In panels (c), (d), and (f), both models can be considered WRm stars (two red vertical lines), but only one in each plot is almost completely hydrogen depleted.In these three cases, the extreme WRm models are those that meet the condition X < 0.3 earlier on during their evolution.At the end of calculations, these six extreme WRm stars have a surface composition mainly of helium, nitrogen, and oxygen.
Figure 4 shows the wind ejecta mass for He, C, N, O, Ne, and Mg, while tables of wind ejecta for all rotating models can be found on Zenodo. 7In the six panels of Figure 4 there are all the models presented in Figure 3 along with stellar tracks computed with their same mass and metallicity, but different initial rotational velocity and mass-loss prescription.In panels (b), (c), and (e) are the extreme WRm models discussed in Sections 3.1 and 3.3.In these six panels, we can see the effect of different initial rotational velocities on wind ejecta masses, as well as the impact of the DUP episodes.Instead, panels (a)-(b), (c)-(d), and (e)-(f) show the incidence of the two mass-loss prescriptions M rdw  and M max  .The general trend is that, also accounting for pulsation-driven mass loss, the ejecta mass of the considered elements increases.This is expected since the definition of the two mass-loss recipes was adopted.Considering the initial metallicity of the models, at higher Z i we have higher ejecta masses due to higher mass-loss rates.However, the most impactful process for example in panels (a), (b), (c), and (e) is the occurrence of DUP episodes.
The extreme WRm models eject much more metals compared to the others, and this is because of a deeper penetration of the stellar envelope (DUP) coupled to the mass loss history of the model.In panel (a) there are no extreme WRm models, but still, we see major differences between stars computed with ω = 0.2, 0.4, and 0.5 with respect to those with ω = 0.0, and 0.3.The reason for this is the metal enrichment of the models due to DUP episodes (see also Figure 3(a)).Finally, we can say that the differences in the ejecta are due to differences in the mass loss history, except when dealing with models with a very high metal enrichment due to DUP mixing episodes.In these latter cases, the relative ejecta mass of the elements considered is much higher, meaning that mass loss alone can not explain these differences and the main driving process is the mixing due to DUP episodes.

Final Fate
Figure 5 shows the final helium core mass (M He ) at the end of calculations for all tracks.We use M He as a proxy for the final fate of the star, following Woosley (2017), Farmer et al. (2019).M i and M He are positively related; however, both panels present stellar models that do not follow this trend.This is caused by DUP episodes that reduce the helium core of the models and in some cases even affect their final fate.The helium core mass is defined by the chemical composition of the envelope, and for this reason, for the 16 WRm stars (see Section 3.1) we consider two extreme possible cases for the final fate.In the first one, as usual, we use the M He core to determine the final fate.In the second case, we use the total final mass of the star at the end of the computations (M f ) to derive the final fate.The combined effects of rotation and DUP episodes affect the stellar tracks in different ways.In panel (a), for example, we can see that rotation increases the helium core mass with respect to the nonrotating case for M i = 100 M e (see also Table 2).At the same time, DUP episodes and higher mass loss rates prevent a steady growth of the helium core with an increasing initial rotation rate.
Stars with M i = 150, 200 M e have He cores that are massive enough to cause a PISN, which disrupts the whole star.There are some exceptions, though.In panel (a), seven tracks do not follow the main pattern.Four of these models were computed with M max  , and only one has M i = 150 M e .This latter is computed with w = M , 0.5 rdw  and has three possible fates.If we consider the strict definition of its helium core, this latter has a mass of 64.3 M e .With this core, the model is inside the The blue regions correspond to the star's convective core; the pink area represents the convective envelope, semiconvective zones at the boundary of the helium convective core, and convective shells.The yellow, cyan, and purple hatch regions show the hydrogen, helium, and carbon-burning core/shells, respectively.The continuous black line indicates the total mass of the star, the orange one represents the helium core, and the green one corresponds to the carbon-oxygen core.The red arrow marks the time when the star enters the unstable region with 〈Γ 1 〉 = 4/3 + 0.01; while the red vertical line shows when the star becomes a WRm with X < 0.3.Panels (a), (c): models computed with the standard PARSEC mass-loss prescription for radiation-driven winds.Panels (b), (d): models that account also for pulsationdriven winds.
uncertainty strip around the PI-PPI boundary.The lower limit of this strip is set to 63.91 M e , while the upper one is 65.24 M e .They are the He core masses of the T135D and T140D models from Woosley (2017), respectively.These two models are computed with Z = Z e /10 and = M 0  , simulating the evolution of a zero-metallicity star.Even though these two models have a He core mass very close, and in one case exceeding the 64 M e limit, Woosley (2017) found that they produced a stellar-mass ) star could imply that the helium core corresponds to the total final mass of the star, M He = 150 M e .In this case, the final fate of the model should be the direct collapse to a BH (DBH).This latter scenario also applies to the ( ) and the ) models if we consider M He = M f (see Tables 1 and 2 for the different final fates considered for each model).
The other tracks in panel (a) with M i > 100 M e that do not become PISN are the models ) with ω = 0.2, 0.3, 0.4 and 0.5.These four models are extreme WRm stars due to their almost H-free envelopes (see Sections 3.1 and 3.3).By considering their helium core masses, these models should undergo PPISN or failed corecollapse supernova (fCCSN, see Figure 5 and Tables 1 and 2).On the other hand, in the case where M He = M f , their helium cores are above the 135 M e upper threshold for PISN; thus, these stars should directly collapse, forming a massive BH.
In , ω = 0.3, 0.4), ( ) models can either collapse directly to a BH or get totally disrupted in a PISN, depending on the helium core mass considered.On the other hand, the ( ) star is in the PPISN or PISN regime in the case M He < 64 M e or M He = M f , respectively (see also Table 1).
In both panels of  (2022), and references therein).In DBH cases, we take the final mass of the star as a first approximation for M BH , always accounting for mass loss due to neutrino emission.In Tables 1  and 2, we give two (or even three) possible outcomes and the corresponding BH masses for all models with an uncertain fate based on the boundaries in Figure 5. Notice that for the model, we assume a 1% error on the 135 M e threshold for the upper limit for PISN in the fit formula by Mapelli et al. (2020).It is worth mentioning that these limits are based on nonrotating models, as long as the fit formula from Mapelli et al. (2020).There are different studies on the pair-instability limits and the effects of rotation, e.g., Glatzel et al. (1985), Woosley (2017), Marchant & Moriya (2020), Woosley & Heger (2021).For example, Woosley & Heger (2021) mention that high rotation rates could shift the lower limit for pair-instability from 64 M e to ∼70 M e .For the case where the He core mass is close to the 64 M e limit ( w ), we can not calculate the BH mass with the fit formula from Mapelli et al. (2020).For this reason, we use a linear interpolation between M He and M BH of the models T135D and T140D.In this way, we can give an estimate of the BH mass for the case of PPISN (see also Table 2).
Tables 1 and 2 and Figure 6 show the results obtained with different mass-loss prescriptions and initial rotation rates.In Figure 6, whenever there are multiple symbols for the same stellar model, it indicates that for that particular star, there is more than one possible outcome.For Z = 0, the most massive BHs reach ∼180 M e , while for Z = 0.0002 they reach ∼172 M e .The complex interplay between DUP and rotation sets the final mass of the BHs we expect from our models.It is worth noticing that rotation favors the occurrence of DUP episodes, shifting, in some cases, the mass of the possible remnant from zero to hundreds of solar masses.Within Section 3.5, the remnants' mass will be rediscussed and adjusted according to results from the analysis on possible jetdriven events (see also Tables 3 and 4).

Progenitors of Jet-driven Events
Within the collapsar scenario (Woosley 1993), the most important characteristics for being a gamma-ray burst (GRB) progenitor are a massive core to produce a BH, the lack of an extended hydrogen envelope to facilitate the jet outward propagation and high core specific angular momentum for the formation of an accretion disk.There are, however, different studies that suggest the possibility of jet propagation even through the very massive envelope of Population III stars (e.g., Ohkubo et al. 2009;Suwa & Ioka 2011;Toma et al. 2011;Nagakura et al. 2012;Wei & Liu 2022).Moreover, there are also multiple criteria for the minimum core specific angular momentum and the mass coordinate to consider when evaluating this latter (e.g., Woosley & Heger 2006;Yoon et al. 2006;Meynet & Maeder 2007;Yoon et al. 2012;Aguilera-Dena et al. 2018;Aloy & Obergaulinger 2021;Obergaulinger & Aloy 2022).
First of all, we have to consider all the shells that have a sufficient specific angular momentum to form a disk instead of falling directly into the newly formed BH (see Appendix B for a different approach in the calculations of the shell inertial moment).Figure 7 shows four example models that summarize all possible cases within our current work.This choice is based on different angular momentum configurations within the progenitors and also the different final fates for the jet-driven events.In each panel, we consider the two limiting cases for the minimum specific angular momentum needed to support mass at the last stable orbit (LSO) of a Schwarzschild and a maximally rotating Kerr BH.Then, the general case considers a BH with mass and angular momentum within the specific mass coordinate in the stellar model (see Bardeen et al. 1972, for the exact expressions).We can see that the specific angular momentum profile is very different between the four models, and therefore the expression j > j crit where we expect disk formation follows different patterns in these four progenitors.We assume that the inner 3 M e of material would form the BH, and therefore do not contribute to the accretion disk.For this reason, we exclude the inner 3 M e from subsequent calculations throughout this work for all our models (see Appendix A, where we consider the case of M BH = 5 M e for the only model where this implies a sizeable difference).
Having enough angular momentum and forming a disk are necessary but not sufficient conditions for a successful GRB.On top of that, the lack of an extended envelope is needed to allow the jet outlet.To analyze the jet propagation through the progenitor in further detail, we proceed as in Yoon et al. (2015) by computing the accretion rate from the disk.et al. (2015), in this work such accretion rates should be considered as lower limits.This is because our models do not reach the precollapse stage, and hence the density should increase with respect to the considered one.
Figure 8 shows the accretion rate, the freefall timescale, and the crossing timescale for those models presented in Figure 7. Also in these four panels we highlight where the shells have enough specific angular momentum to support a disk (see Figure 7 and related discussion).Panels (b) and (c) do not show the entire progenitor.This is because in the PPISN cases, we plot only the model up to the mass coordinate calculated with the fit formula by Mapelli et al. (2020).In this way, we consider the model after the loss of the envelope due to pulsational pair-instability (Ekström et al. 2008), assuming that the pair-induced pulsations do not affect the angular momentum of the stellar core (Aguilera-Dena et al. 2018).As a first approximation for the velocity of the ejected material due to PPI, we consider the escape velocity at the mass coordinate given by the fit formula from Mapelli et al. (2020).We find that v esc ∼ 10 4 km s −1 and the timescale between the pulsational pair-instability and the collapse of the star should be of the order of ∼100 yr.We can safely stop considering the ejected material in our calculations since it should be at distances of ∼1 pc when the jet forms, thus reducing the density of the ejected envelope by a factor (1 pc) 3 .Instead, in panels (a) and (d) we show the entire stellar model since we consider the  1 and 2. Panel (a): models calculated with Z = 0. Panel (b): models calculated with Z = 0.0002.
M He = M f case and therefore the stars do not undergo PPI (see Section 3.4 and Tables 1, 2).
The models in Figure 8 represent the four possible cases for a jet-driven event among all our rotating stellar tracks, except those that undergo PISN.For this reason, we assign the letters a, b, c, or d in the final fate column of Tables 3 and 4, according to the structure and final fate of the models following one of these four possible cases.Tables 3 and 4 also summarize different physical properties of possible jet-driven events for all rotating models considered in this work.
The model in Figure 8(a) does not experience pulsational pair-instability and it is within the DBH scenario (taking the   8; (9) remnant mass(II)the bottom of the envelope here is defined as the first point where Y < 10 −3 ; (III) assuming an error of 1% on upper limit for PISN in fit formula by Mapelli et al. (2020); a Considering M He = M f .
upper limit M He = M f , note in Table 3 that the final fate should not change even considering the lower value for M He ).Therefore, we plot the entire stellar structure where the final mass of the model is 143 M e .Panel (a) shows that there could be three distinct accretion episodes due to the distribution of specific angular momentum within the star from matter at 3 M e -14 M e , 16 M e -38 M e and 111 M e -143 M e , respectively.
Combining the first two accretion episodes, the accretion timescale for matter within the stellar core is 5.66 s.However, this is much lower compared to the crossing timescale for jet propagation through the envelope, which is ∼10 3 s (see Table 3).This shows that the jets powered by core accretion are not expected to produce a successful GRB, due to a much difficult jet propagation.Moreover, the crossing timescale is more than two orders of magnitude smaller than the freefall timescale (τ cross = τ ff ), and thus the jets would not remain collimated during propagation within the star (Yoon et al. 2015).Following Suwa & Ioka (2011), to calculate the energy of a jet, we assume the following expression for the luminosity: where M  is the mass accretion rate from Equation (5) and the accretion-to-jet conversion efficiency is η = η 0 a 2 = 10 −3 a 2 , where a is the dimensionless spin parameter (a )) of the central BH with M BH = 3M e and J the corresponding angular momentum.The first term η 0 comes from Suwa & Ioka (2011), while the dependence on the BH spin is from Blandford & Znajek (1977).In those models where there are two accretion episodes within the stellar core, to calculate the luminosity of the second jet, we assume M BH = 3 M e + M acc , where M acc is the total mass accreted by the BH before the second accretion episode.Thus, the accretion-to-jet efficiency changes too, because the dimensionless spin parameter is not calculated for a BH of 3 M e as for the other jets.The energy that the jets pump into the stellar envelope is ∼3.3 • 10 52 erg.Because these jets can not break out from the progenitor and they are going to spread out instead of remaining collimated, we sum their energy to compare it with the binding energy of the envelope.We define the bottom of the stellar envelope as the first point where X < 10 −3 (note the only exception in Table 3 to estimate more accurately M BH ), hence the total envelope binding energy is ∼1.6 • 10 50 erg.Because of this, the whole envelope would be ejected by the jets' energy, therefore hindering the third accretion episode.In this case, the final outcome of the model should be a jet-driven supernova (SN).
In panel (b), the model has enough angular momentum in its inner core to support the formation of an accretion disk.The first accretion episode lasts 0.65 s, while the crossing timescale for the jet to reach the surface of the star is 0.36 s.This means that in this case, there could be a successful GRB event powered by the accretion of matter between 3 M e and 15 M e .The crossing timescale is more than four times smaller compared to the accretion timescale for matter between 22 M e and 44 M e .This implies that the second accretion episode can not form another jet before the first one breaks out from the star.Here we assume as an upper limit that the first jet does not blow away any mass of the outer stellar core, though this mass could be very small (Zhang et al. 2004).Thus, in this case, we have two accretion episodes that produce successful GRBs."Case b" corresponds to progenitors of successful GRB events that are powered by two accretion episodes, where in Tables 3 and 4 there are the results for the two jets in the τ acc , L acc , and E acc columns.
Case c is very much similar to the latter.The only difference is that the entire model has enough angular momentum to form a disk: there is a single accretion episode that involves the whole star.Even in this case, the accretion timescale is higher concerning the crossing one, 2.83 and 0.39 s, respectively.Therefore, this model should also produce a successful GRB event.Here we assume that the whole mass of the star is going to be accreted, and this gives an upper limit for the accretion timescale and hence the total accretion energy of GRB progenitors that follow the kind of structure presented in this case c.
The model in Figure 8(d) is very similar to that in panel (a) (taking the upper limit M He = M f ).As between cases b and c, the most important difference is that in Figure 8(d) there should be only one accretion episode from matter within the stellar core.Even in this case, the accretion timescale is much shorter compared to the crossing timescale (see Table 4).Hence, also in this model, the jet powered by core accretion should produce a jet-driven SN, to which we refer as final fate d in Tables 3 and  4. On the other hand, if we consider the lower value of M He for this specific model, the jet-driven event changes deeply.Because of the lack of an extended envelope, the jet could reach the stellar surface and produce a successful GRB (see Table 4 for the numerical details).Tables 3 and 4 also present the expected mass of the remnant after the jet-driven event, where also here we account for neutrino emission and take 90%.8 of the considered final mass for the expected M BH Fryer et al. (2012), Rahman et al. (2022), and references therein.For jet-driven SNe in cases a and d, we consider only the mass of the core because the envelope is ejected by the energy of the jet(s).Instead, in cases b and c, we take into account the total mass of the model after pulsational pair-instability mass loss.All values in the last column in Tables 3 and 4 are to be considered as upper limits for the BH mass.This is because there could be more mass ejected by the jet, even though it might not be that large (Zhang et al. 2004).
Tables 3 and 4 show the luminosity due to accretionpowered jets for all rotating models considered (except those that undergo PISN).The two brightest GRB events have ( ) L log 52.5.The main driver of the jet-powered events is the accretion rate shown in Figure 8, which depends on the density of the models.Now we consider different values for the accretion-to-jet conversion efficiency parameter.With η 0 = 10 −2 ; 10 −4 , the overall results do not change.Of course, we obtain respectively higher and lower luminosities and accretion energies, but the cases highlighted in Tables 3 and 4 remain unchanged.There is only one exception, which is the ( ) model in the case of the DBH scenario (see Table 4).With η 0 = 10 −4 the accretion energy results lower than the binding energy of the stellar envelope.Thus in this case the expected remnant mass from this progenitor is higher, 78 M e .It is worth noting that this latter is almost 20 M e above the value expected when the total stellar envelope is ejected by the energy of the jet.Considering this lower efficiency, we should expect a failed jet-driven SN from this progenitor, which greatly impacts the final mass of the BH.

Concluding Remarks
In this work, we study the evolution of rotating zero-metallicity and metal-poor very massive stars with initial masses between 100 M e and 200 M e .We investigate the resulting effect arising from different initial rotation rates ω = 0.0, 0.2, 0.3, 0.4, 0.5, and pulsation-driven winds (following the prescription of Nakauchi et al. 2020), accounting also for radiation-driven winds throughout the entire HRD.These are the first very massive models that consider both stellar rotation and pulsation-driven mass loss, extending the parameter space covered by the PARSEC evolutionary tracks.We follow the evolution until the occurrence of electron-positron pair-instability after carbon, neon, or oxygen ignition depending on the initial mass, metallicity, rotation, and mass-loss prescription adopted for the models (see Tables 1 and 2).For all models, we checked that rotation and radial pulsations should not influence each other, with the former dominating over the latter (see discussion in Appendix C).
As in Paper I we discuss the final fate of these stars (Section 3.4), but then accounting also for stellar rotation we analyze the possibility of jet-driven phenomena from these rotating progenitors (Section 3.5).We consider the accretionto-jet efficiency parameter η = η 0 a 2 , where η 0 = 10 −2 , 10 −3 , 10 −4 , while a is the dimensionless spin parameter of the central BH with M BH = 3 M e (see Appendix A where we explore the M BH = 5 M e case).Results in Tables 3 and 4 assume the standard value of η 0 = 10 −3 from Suwa & Ioka (2011).
Following the analysis of Yoon et al. (2015), we identify four possible cases among our rotating models with different structures and final fates.
For stars that do not lose any mass due to pair-instability and still form a BH (DBH scenario, cases a and d), the expected outcome is a jet-driven SN.In these cases, the extended stellar envelope prevents the jet breakout, but the energy of this latter is higher than the envelope binding energy, and thus the jet unbinds the outer layers of the star.The progenitors of these successful jet-driven SNe have M i = 200 M e , except for four models with M i = 150 M e (see Tables 3 and 4).Within case d when considering η 0 = 10 −4 , there is only one exception to the fact that the accretion energy of the jet is bigger than the envelope binding energy, i.e., ( ) model when considering M He = M f .This implies an increase of the expected BH mass of almost 20 M e .
All models that undergo pulsational pair-instability lose their envelope, hence fostering the propagation of the jet as opposed to the models in Yoon et al. (2015).) model, which has three possible fates (see Table 4).The major difference between the b and c scenarios is that in the former the star does not have enough angular momentum throughout the whole structure, and thus there are two distinct accretion episodes.Instead, in the latter the whole mass of the star is accreted through the disk in one single episode, which increases the accretion timescale.This difference also occurs between models in cases a and d.
A GRB event can be detected by the BAT X-ray detector up to redshift z ∼ 20 if it has L  10 52 erg s −1 s (Komissarov & Barkov 2010;Yoon et al. 2015).All successful GRB events listed in Tables 3 and 4 would be luminous enough to be detected.This changes as we consider different values for the parameter η 0 .When adopting a lower value for efficiency, all the models get below the observability threshold above.On the other hand, increasing the accretion-to-jet efficiency boosts the possibility of detecting this kind of event.Also, the afterglow of these GRBs could be of paramount importance.The reason is that with a bigger energy budget, and therefore a radio flux peaked at late times at gigahertz frequencies (Toma et al. 2011), the radio afterglow of GRBs from Population III progenitors should be much brighter than that of Population II/Population I stars (Salvaterra 2015;Burlon et al. 2016).Hence, this could be key for distinguishing GRB events from different progenitor populations.With current instruments like the Australian Square Kilometer Array Pathfinder telescope and James Webb Space Telescope with both Near-InfraRed Camera and Near-InfraRed Spectrograph, the observation of the afterglow from Population III GRB events should be within reach (Macpherson & Coward 2017).
Jet-driven events also have a deep impact on the expected M BH for some of these progenitors.Comparing the results for the expected BH mass in Tables 1 and 2 with those in Tables 3 and 4, there are 15 models with a different remnant mass.These stars are all within the a and d cases (see discussion in Section 3.5), while for models following cases b and c, the expected mass of the BH does not change.The differences are due to jet-driven SNe that unbind the stellar envelope during the accretion-powered jet's propagation.Hence, the BH mass can be reduced by more than 130 M e in the most extreme cases.
Out of 15 models with a reduced M BH , 13 stars have an expected remnant mass that falls within the pair-instability black hole mass gap (40-65 M e to 120 M e , see also Croon et al. (2020), Farmer et al. (2020), Sakstein et al. (2020), Costa et al. (2021), Farrell et al. (2021), Tanikawa et al. (2021), Vink et al. (2021), Farag et al. (2022), for different formation scenarios).Several stellar models can produce BHs with a mass close to one of the primary and the secondary BHs in the gravitational wave event GW190521 (Abbott et al. 2020) (4) and (5) fractions of MS and cHeB lifetimes in which the star is unstable to radial pulsation; (6) and (7) occurrence of blue loop and dredge-up episode; (8) final He core mass; (9) final C-O core mass; (10) final mass of the star at the onset of dynamical instability; (11) central fuel abundance of ongoing nuclear burning at the onset of dynamical instability; (12) neutrino luminosity to radiative luminosity ratio when T c = 10 9 K; (13) and (14) final fate and associated outcome (BH or complete disruption), and (15) BH mass.a Assuming an error of 1% on the upper limit for PISN in the fit formula by Mapelli et al. (2020).b Considering M He (2019)  we set the lower limit of M He for PPISN at 45 M e. b Following Woosley (2017) we set the lower limit of M He for PPISN at 34 M e. c Following Woosley (2017) we set the lower limit of M He for PISN at 64 M e. d We set the lower limit of M He for PISN at 65.24 M e , which is the M He of the T140D model inWoosley (2017).

Figure 1 .
Figure 1.H-R diagrams of the 20 sets of tracks computed in this work.Different initial rotation rates are color coded as labeled.Panels (a), (b) and (c), (d) refer to Z = 0 and Z = 0.0002, respectively.Panels (a), (c): tracks calculated with the standard PARSEC mass-loss prescription for radiation-driven winds.Panels (b), (d): tracks computed considering also the pulsation-driven mass loss.Cyan starred symbols indicate where the stellar models become Wolf-Rayet manqué stars (X < 0.3).We plot the value of initial mass (in M e ) for all tracks.

Figure 2 .
Figure2.Kippenhahn diagrams of four selected models.The horizontal axis shows the age of the models as a logarithm of time (in yr) until the end of computations.The blue regions correspond to the star's convective core; the pink area represents the convective envelope, semiconvective zones at the boundary of the helium convective core, and convective shells.The yellow, cyan, and purple hatch regions show the hydrogen, helium, and carbon-burning core/shells, respectively.The continuous black line indicates the total mass of the star, the orange one represents the helium core, and the green one corresponds to the carbon-oxygen core.The red arrow marks the time when the star enters the unstable region with 〈Γ 1 〉 = 4/3 + 0.01; while the red vertical line shows when the star becomes a WRm with X < 0.3.Panels (a), (c): models computed with the standard PARSEC mass-loss prescription for radiation-driven winds.Panels (b), (d): models that account also for pulsationdriven winds.

Figure 3 .
Figure 3. Surface chemical abundances' evolution of six selected pairs of models, from the ZAMS to 10 yr before the end of computations.In each panel, the abundances of hydrogen, helium, carbon, nitrogen, and oxygen are color coded as shown in the legend.The black line represents the effective metallicity, Z eff = 1 − X − Y.The red vertical line marks when X < 0.3.The results are presented for two different mass-loss prescriptions.The horizontal axis is as in Figure 2.

Figure 4 .
Figure 4. Wind ejecta mass of models with initial mass M i /M e = 150, 200.In each panel, there are the ejecta masses of helium, carbon, nitrogen, oxygen, neon, and magnesium for two mass-loss recipes, M rdw  and M max  .Panels (a), (b): models calculated with Z = 0. Panels (c), (d), (e), (f): models calculated with Z = 0.0002.
panel (b) there are nine WR-manqué stars, five of which are computed with M max  .Both the M i = 150 M e and M i = 200 M e models with M rdw  and ω = 0.2 have M He > 135 M e , in the case of M He = M f .Therefore, these two models could undergo DBH.Besides these two, seven stars in this panel have a double final fate.The (

Figure 5 ,
stars with M i = 100 M e are in the pulsational pair-instability regime, but outside the uncertainty region.The only exceptions are the two models computed with ω = 0.0 (see discussion in Paper I) and ( above.As in Paper I, we use the fit formula proposed by Mapelli et al. (2020, see Appendix D for details) to compute the BH mass in the PPISN scenario, while also accounting for mass loss due to neutrino emission (10% of the baryonic mass of the proto-compact object, see Fryer et al. (2012), Rahman et al.

Figure 5 .
Figure 5.Each panel shows the helium core mass, M He , as a function of the initial mass for all models of five different initial rotation rates with a fixed initial metallicity.Helium core mass is shown at the end of the calculations.Horizontal lines limit the pulsational pair-instability (PPI), pair-instability (PI) explosion, and direct collapse to a BH (DBH) regime from bottom to top (Woosley 2017; Farmer et al. 2019, 2020).The lower red strip indicates the uncertainty range of the lower limit for PPI.Lower and upper boundaries are 34 M e (Woosley 2017) and 45 M e (Farmer et al. 2019), respectively.The black line is the average.Panel (a): similar uncertainty strip between PPI and PI with boundaries 63.91 M e and 65.24 M e .These values are the helium core masses of the T135D and T140D models from Woosley (2017), respectively (see text for more details).
density within each shell.As in Yoon

Figure 6 .
Figure 6.Expected remnant mass for all stellar models presented in this work as a function of M i .The extra symbols indicate the predicted BH mass for those stars with multiple possible fates, see also Tables1 and 2. Panel (a): models calculated with Z = 0. Panel (b): models calculated with Z = 0.0002.

Figure 7 .
Figure7.Specific angular momentum profile as a function of the mass coordinate.In each panel the black line shows the specific angular momentum for the model considered, the blue line refers to the minimum angular momentum needed to support mass at the LSO of a Schwarzschild BH, while the green line denotes the minimum angular momentum for a maximally rotating Kerr BH.The red line shows the specific angular momentum at the LSO for a BH with mass and angular momentum inside the considered mass coordinate in the computed stellar model.All panels show in lilac the regions within the stars where j > j crit .In these regions, we expect the formation of an accretion disk.The two innermost vertical lines in each panel refer to the mass coordinate 3 M e and 5 M e , respectively, while the outer one corresponds to the carbon-oxygen core of the model considered.

Figure 8 .
Figure 8. Accretion rate, freefall timescale, and crossing timescale as a function of the mass coordinate.Lilac regions show the shells with sufficient specific angular momentum to support a disk, where j > j crit .The two timescales refer to the y-axis on the right, while the accretion rate is shown on the left.
Case b and c correspond only to models with M i = 100 M e , except the (

Table 1
Most Relevant Properties of Models Computed with Z = 0.0002, M rdw  and M max 
a Failed CCSN.Following Farmer et al.

Table 3
Most Relevant Properties of Possible GRB Events for Models Computed with Z = 0.0002, M rdw  and M max  The table entries are as follows: (1) star's initial mass; (2) crossing timescale; (3) freefall timescale; (4) accretion timescale, ∑ i τ acc,i ; (5) accretion luminosity, ∑ i L acc,i ; (6) accretion energy, ∑ i τ acc,i • L acc,i ; (7) envelope binding energy; (8) case for the possible fate of the GRB progenitor according to four cases outlined in Figure

Table 4
Most Relevant Properties of Possible GRB Events for Models Computed with Z = 0, M rdw  and M max  Notes.Table entries as in Table 3. a Concerning the expected BH mass, the model follows the c case since it should experience PPI.b Following Woosley (2017) we set the lower limit of M He for PISN at 64 M e. c We set the lower limit of M He for PISN at 65.24 M e , which is the M He of the T140D model in Woosley (2017).d Following Farmer et al. (2019) we set the lower limit of M He for PPISN at 45 M e. e Following Woosley (2017) we set the lower limit of M He for PPISN at 34 M e .
, respectively.Therefore, rotating primordial very massive stars could provide a new pathway for the formation of BHs within the pair-instability black hole mass gap.angular momentum able to form an accretion disk.Hence, the whole stellar core should collapse into a BH.