Metal-enriched Pair-instability Supernovae: Effects of Rotation

In this paper, we revisit metal-enriched rotating pair-instability supernovae (PISNe) models for metallicities consistent with the Small Magellanic Cloud, the Large Magellanic Cloud (LMC), and 0.1Z ⊙. By calculating multiple models, we intend to clarify mass ranges and the ejected 56Ni masses from the PISNe, and mass-loss histories for progenitors. We find that the choice of the Wolf–Rayet (WR) mass-loss rates is important, and we adopt the recently proposed rate of Sander & Vink, which covers the mass ranges for PISNe progenitors. We show that slow rotation lowers the PISN range because the core mass increases by rotational mixing. On the other hand, when we assume a typical rotation speed for observed OB stars, the mass-loss increase becomes more significant, and the final stellar masses are lower than in nonrotating models. As a result, the typical mass range for bright supernovae (SNe), with a 56Ni mass higher than 10M ⊙ for these fast-rotating models is more than 400 and 350M ⊙ for LMC and 0.1Z ⊙ metallicities, respectively. It is interesting that unlike in previous works, we find oxygen-rich progenitors for most cases. This O-rich progenitor may be consistent with the recently identified PISN candidate SN2018ibb. He-rich progenitors are seen only for relatively dim and metal-poor (Z ≲ 0.1Z ⊙) PISNe. We also discuss the black hole mass gap for metal-enriched PISNe, and we show that the upper bound for the gap is lower than in the Population III case.


INTRODUCTION
A pair-instability supernova (PISN) is a thermonuclear explosion of a massive Oxygen core induced by the creation of electron-positron pairs.After this mechanism was initially proposed (Barkat et al. 1967;Rakavy & Shaviv 1967), detailed calculations confirmed the viability of the explosion (e.g.Umeda & Nomoto 2002;Heger & Woosley 2002;Takahashi et al. 2018).These studies mostly focused on metal free Pop III stars because the mass range of PISNe is large (about 140 − 300  ⊙ for metal free stars) and stars massive enough to experience the pair instability are thought to be rare in the present universe because of large wind mass loss for metal enriched stars.Although there are some calculations for metal enriched stars (e.g.Langer et al. 2007;Yoshida & Umeda 2011;Yusof et al. 2013;Kozyreva et al. 2014;Whalen et al. 2014), all but two of these works (discussed below) assume that the main properties of PISNe, i.e., the explosion energy and the amount of 56 Ni production, depend solely on CO-core mass and not on metallicity.Kozyreva et al. (2014) and Whalen et al. (2014) calculated hydrodynamics and nucleosynthesis for metal-enriched models ( ≳ 0.001) during the exploding phase, but only two mass models are shown for each metal-enriched case.
Observationally it is not clear if there has been any evidence for the existence of PISNe.Since the nucleosynthetic pattern is quite different from core collapse SNe, if an extremely metal poor star was formed from the gas of PISN ejecta, it should not be difficult to identify (e.g.Umeda & Nomoto 2005).Recently, a single star with the signature of PISN enrichment was reported (Xing et al. 2023), but no other metal poor stars with evidence of PISN enrichment have been found, although there is a suggestion that PISN abundance patterns might be hidden in very metal poor stars (Aoki et al. 2014).
A PISN can also be observed as a nearby supernova since in principle PISNe are possible for  ≲  ⊙ /3 (Langer et al. 2007).The first suggestion for such a discovery was reported as a Type I super luminous supernova (SLSN-I), SN2007bi (Gal-Yam et al. 2009).However, spectroscopic models of PISNe were incompatible with the observations both in the photospheric (Dessart et al. 2012) and nebular (Jerkstrand et al. 2017) phases.Since then, magnetar models have more frequently been considered to explain SLSNe-I, and PISNe models are rarely taken seriously.Nevertheless there remain some SLSNe-I, such as PS1-14bj, PTF10nmn, OGLE14, and SN2020wnt, which cannot be well explained by the magnetar models and these could well be PISNe.If these candidates are PISNe, the ejected 56 Ni mass falls in the range of ∼ 0.5−10 ⊙ (Gal-Yam 2019).Very recently, there was a report that SN2018ibb might be the best PISN candidate discovered thus far (Schulze et al. 2023).SN2018ibb is a Hydrogen poor super-luminous SN at  = 0.166, and if it is a PISN, the expected 56 Ni mass is more than 30 ⊙ .Schulze et al. (2023) also mention that there is a signature of interaction with oxygen-rich CSM.
In Yoshida & Umeda (2011) we showed that if a nearby PISN with metallicity  = 0.004 ejects more than a few solar masses 56 Ni, the initial stellar mass should be more than 500 ⊙ for standard mass loss rates in the non-rotating models.On the other hand, Yusof et al. (2013) found that rotating models can be such an PISN in the mass range 150-180  ⊙ for SMC metallicity.In general, rotation increases core mass which may reduce the PISN mass range significantly.However there are some differences in the assumptions of these works and consequences of these differences is unclear.In this paper, we revisit metal-enriched rotating PISNs models to clarify model dependence of these PISNe.This is motivated by the paucity of models in the existing literature.
The remainder of the paper is organized as follows.In Section 2, we summarize our numerical methods, including initial rotational profiles (2.1), mass loss rates (2.2), angular momentum transfer (2.3), and finally the hydrodynamical calculations of the PISNe (2.4).Section 3 shows mass loss histories (3.2), stellar evolution of rotating models (3.3), hydrodynamics (3.5), nucleosynthetic post processing (3.6) and some examples of peculiar results (3.7).We then discuss our results in the context of previous work (Sec.4) before concluding in Sec. 5.

METHODS
In this paper we calculate evolution of sub-solar metallicity stars in the PISN mass range for three metallicity cases: LMC metallicity ( =  ⊙ /3), SMC metallicity ( =  ⊙ /5) and  = 0.1 ⊙ .Here,  ⊙ is the solar metallicity and we adopt  ⊙ = 0.0141.We note that the SMC and LMC metallicities mentioned above roughly corresponds to the lower metallicity bounds for these galaxies.
For this study, we use the HOSHI (HOngo Stellar Hydrodynamics Investigator) code described e.g. in Takahashi et al. (2018) and Luo et al. (2022) for the stellar evolution.We follow the nuclear burning using a nuclear reaction network of 153 species of nuclei (Takahashi et al. 2018).Nuclear reaction rates are taken from the JINA reaclib database v1 (Cyburt et al. 2010), except for the 12  (, ) 16  rate which is taken to be 1.5 times the value given in Caughlan & Fowler (1988).In this paper we adopt the overshooting parameter  ov =0.01 as same as in (Takahashi et al. 2018), which is called the 'M' model in Luo et al. (2022).Then we use the hydrodynamical code described in Takahashi et al. (2016) and Nagele et al. (2022) for the collapse and PISN explosions after the model approaches the electron-positron pair-creation instability region.
Calculation methods for stellar evolution are mostly the same as in our previous work, but there are differences in the setting of initial models, and the treatments in the mass loss rates as explained in Sec 2.2.Previously we started stellar evolution from a zero age main sequence (ZAMS) model and set the initial rotation speed for that ZAMS model.However, the definition of ZAMS is somewhat uncertain and it is thus difficult to compare with other works which may adopt a slightly different definition for ZAMS.Usually this small difference does not matter, but in this work initial rotation speed is critical, so we adopt a more concrete definition.

Initial rotation
We now describe the determination of initial rotation in this work.First we construct an expanded pre-mainsequence model which has central temperature below log  c < 7.0 (K).Then we begin the stellar evolution calculation, and when the central temperature reaches log  c = 7.0 (K), we set the initial rotational velocity.We assume rigid rotation initially and the ratio of the surface velocity,  i , to the Keplerian velocity,  K ≡ √︁  /, is taken as a parameter.As shown below  i / K = 0.1 corresponds to relatively slow rotation, and 0.2 to fast-rotation.
Figure 1 shows typical evolution for the ratio of the surface velocity to the Keplerian velocity,  rot / k .These lines are for 0.1  ⊙ ,  = 200 and 250  ⊙ with  i / k = 0.1 and 0.2, and SMC metallicity,  = 200  ⊙ with  i / k = 0.1.As shown in this figure,  rot / k increases with the initial stellar shrinkage, and has a maximum around ZAMS, where log   (K) ∼ 7.7.These maximal values are shown in Table 1 as ( rot / k ) max .We note that the SMC model has lower ( rot / k ) max than the 0.1  ⊙ model with the same  and  i / k due to slightly larger mass loss before ZAMS.
This figure also shows that ( rot / k ) max is about 0.2 for  i / k = 0.1 and 0.3-0.4 for  i / k = 0.2.Since the typical  rot / k of massive main sequence stars is about 0.3-0.4(e.g.Georgy et al. 2012), we may say that  i / k = 0.2 models rotate rapidly, and  i / k = 0.1 models are relatively slowly rotating, though we do not have any observational evidences that very massive stars considered in this paper rotate as fast as observed massive stars.

Mass loss rates
We will find that the mass loss rates, especially the rate for Hydrogen poor stars, are important for metal-enriched PISNe.Therefore, in this subsection we describe in detail the mass loss rates we use.For Hydrogen-rich stars with the surface Hydrogen mass fraction, X(H), higher than 0.5, we use a rate by Vink et al. (2001) if the effective temperature is higher than log  eff = 4.05.We use a different rate for Hydrogen-rich (X(H)> 0.3) Yellow and Red supergiants with log  eff < 3.7 as in Sylvester et al. (1998);van Loon et al. (1999); Crowther (2000); Ekström et al. (2012).If log  eff > 4.05 and X(H) < 0.3, we use the rate for Wolf-Rayet (WR) stars.If the parameter space is outside the region mentioned above, we use de Jager et al. (1988).We note that as shown in Yusof et al. (2013), the very massive stars considered in this paper typically have large log  eff , say > 4, throughout the evolution, so we are most often using the Vink et al. (2001) rate for Hydrogen-rich stars.
Previously, e.g. in Yoshida et al. (2019) and later work using the HOSHI code, we have used Nugis & Lamers (2000) for hydrogen poor stars, and used de Jager et al. (1988) outside the region of Nugis & Lamers (2000).However, we have found that the mass loss rate sometime becomes small in this approach and we had sometimes significantly larger final mass than other groups' work.Therefore, we decided to replace the Wolf-Rayet mass loss rate with more commonly used ones from recent similar studies.Specifically, we tried two types.One type is the rate used in Aguilera-Dena et al. (2018, 2020) and Yoon (2017) (hereafter Y2017 rate).The other is a theoretical rate given in Sander & Vink (2020) (SV2020).After implementing these rates, we had roughly consistent results with other works.For example, the final mass of solar metallicity ∼ 50 ⊙ models are about ten solar masses.
There are some differences in the two WR mass loss rates.As shown in SV2020, their rate is typically smaller than the Y2017 rate.There is also a notable difference in the metallicity dependence.The rate by SV2020 decreases rapidly toward lower metallicity compared with other rates in the literature.Because of these differences we found that the PISN mass range is significantly larger if Y2017 rate is used.Since the Y2017 rate is determined mainly to explain observed proper-ties of WR stars, it is not clear if we can extrapolate it to the PISN region.On the other hand, SV2020 rate can be applicable up to 500  ⊙ and metallicity down to 0.02  ⊙ .Therefore, we use the SV2020 rate in this work, though we emphasize that the mass loss rate is quite uncertain.More specifically, for the WR stars, we use the expression in the equation ( 28) of SV2020 and parameters given in their equations from ( 23) to (32).
When a star is rotating rapidly, mass loss is enhanced due to the nearly-critical rotation at the surface (the Ω − Γ limit, Langer 1998;Maeder & Meynet 2000).In order to include this effect, we consider the following expression for the mass loss rate according to Yoon et al. (2012): where  rot ,  crit ≡ √︁   (1 − / Edd) /,  KH ,  edd are the surface rotation velocity, the critical rotation velocity, the Kelvin-Helmholtz timescale, and the Eddington luminosity, respectively.There is an ambiguity in the definition of the Kelvin-Helmholtz timescale ∼   2 /, so we introduce a parameter of  (1),  KH , such that  KH =  KH   2 /.In this paper we adopt  KH = 0.5 (Kippenhahn & Weigert 1990).
As shown above, we set an upper limit on the magnitude of the mass loss rate and refer to it as  up ≡ 0.3/ KH .We refer to the magnitude of the non-critical rate as In general, the rotation induced mass loss rate increases when a star is shrinking, because both the  rot and  tend to increase.The reader might think that in these stages,  rot increases gradually and when it reaches  up , we replace it with the upper limit.However, this is not the typical realization.Since we cannot adopt very short time steps to follow the star over evolutionary timescales, in a realistic calculation, the ratio  rot / crit often approaches and then exceeds 1.We use  up if  rot / crit exceeds 1.In the usual calculations,  rot ≪  up even if  rot / crit is very close to 1. Therefore, in practice we use  rot for  rot / crit < 1 and use  up for  rot / crit ≥ 1.This inclusion of the upper limit for mass-loss rate is different from recent related work, Aguilera-Dena et al. (2018, 2020), where they instead set the upper limit to the value  rot / crit to be 0.98.This difference may affect the results for rapidly rotating cases; however, the difference is not important in this paper, since it turns out that all the models shown in this paper end up as slow rotators because of large mass loss during the late He-burning and WR stages.

Angular momentum transfer
In the HOSHI code, the diffusion approximation is applied for the transportation of angular momentum, similar to the codes in Yoon & Langer (2005); Woosley & Heger (2006) which were used to calculate progenitor models of long gamma-ray bursts (GRBs).As in these works, we assume a magnetic model and the Tayler-Spruit dynamo (TS dynamo, Spruit 2002), is applied.With the TS dynamo, angular momentum is transferred efficiently when differential rotation exists.For the inclusion of the TS dynamo effect, we follow the method in Heger et al. (2005) as described in Takahashi et al. (2014).Because of this effect, the stellar cores are relatively slowly rotating after the He burning stages.It is then difficult to produce the fast rotating cores that are necessary for GRB progenitors (Heger et al. 2000) unless chemically homogeneous evolution (CHE) occurs (Yoon & Langer (2005)).Since we have not applied our code to CHE models and also adopt a different WR mass loss rate from previous works, we tested whether we can produce CHE stars.We confirmed that CHE stars with a fast rotating Fe core can be obtained with our code for metal poor, i.e., below 0.1 Z ⊙ , ∼ 30 − 50 ⊙ stars.In this paper, however, we focus on how rotation affects the PISN explosion mass range.

Hydrodynamical and nucleosynthesis calculation
In this paper, we aim to calculate explosion and nucleosynthesis for pair-instabillity supernovae.Such calculations were done in Takahashi et al. (2018) for Pop III PISNe and we use basically the same method.We calculate stellar evolution with the HOSHI code, until the central temperature reaches around Log Tc =9.2 and then we switch to the hydrodynamical code.The HOSHI code includes the acceleration term in the equation of motion and in principle could calculate hydrodynamics.However, since this is a Henyey type stellar evolution code, energy conservation is not as good as in a hydrodynamical code and not suitable for PISN simulations.This is the reason we switch to a hydrodynamical code.
The hydrodynamics code is a one dimensional Lagrangian general relativistic hydrodynamics code (Yamada 1997;Sumiyoshi et al. 2005;Takahashi et al. 2016).The code includes energy changes due to nuclear reactions and neutrino cooling, and we use the same 153 isotope network as in HOSHI.In this paper, we use 255 radial meshes and the HOSHI models are mapped to this grid using the same procedure as in previous works (e.g.Takahashi et al. 2018) while the numerical settings are identical to those used in Nagele et al. (2023).After shock breakout, the timesteps become large and we terminate the simulation.For one LMC model (100  ⊙ ), a pulsation occurs instead of an explosion, and for this model we terminate the calculation after shock breakout, as we are not overly interested in pulsations in this paper.
After the hydrodynamics calculation is completed, we post process the nucleosynthesis with a network consisting of 300 nuclei, as in Takahashi et al. (2018).The 56 Ni masses reported in this paper are from the post processing, and are slightly lower than the values obtained from hydrodynamics.After post processing the hydrodynamics, additional post processing is carried out with fixed temperature and density until the composition is fully decayed.The reported elemental yields use the abundances at the end of this additional post processing.

PISN Mass Range
In Table 1, we show the initial stellar mass range for the models which become PISNe and the (most) abundant elements at the surface of the models in our results.For example, He/H → O ( ≥ 180) means that the surface abundance is He and H-rich for lower masses and O-rich above  = 180 ⊙ .In these models, the mass fractions of He and H are about 0.8 and 0.2, respectively.
In general, the PISN range is higher for larger metallicity since the mass loss rate is larger.The effects of stellar rotation on the PISN range can be summarized as follows.Slow rotation lowers the range slightly since the He-core becomes larger due to rotational mixing.However, in this parameter range, faster rotation increases the mass range since the stars become WR-stars earlier because of rotation induced mass loss, meaning that the final mass becomes smaller due to the efficient mass-loss during the WR stages.For much lower metallicity where the mass loss of WR stars is small, we may find CHE stars becoming a PISN.However, we have not found such a case in this metallicity range.There are a few previous studies for metal-enriched PISNe.In Whalen et al. (2014), using the GENEVA stellar evolution code (Hirschi et al. 2004;Eggenberger et al. 2008), they showed that 150 and 200  ⊙ with 0.1  ⊙ and 500  ⊙ with 0.3  ⊙ become PISNe.Our results are broadly similar though there are some specific differences.Their models assume relatively fast rotation, which corresponds to  i / k ≳ 0.2 in our definition, and our models yield O-rich PISNe within a higher mass range.This mainly originates from the difference in the mass loss rate, especially in the WR stages.
Yusof et al. ( 2013) also provided a PISN mass range using the GENEVA code.Their non-rotaing LMC model became a PISN for  ≳ 300 ⊙ .This is roughly consistent with our results, as their rotating models also have a higher mass range.Their results for rotating SMC models are quite different from ours.They show that 150 and 200  ⊙ models become PISNe with large amounts of 56 Ni production, while for us these models do not become PISNe.Our models require lower metallicity and slower rotation to make such energetic PISNe.
Yoshida & Umeda (2011) also showed the range for metallicity  = 0.004 non-rotating models.Compared with our non-rotating SMC models, their lower bound is similar but their upper bound is much larger (more than 500 ⊙ ).In addition they found only He-rich models.These differences also originate from the chosen WR mass-loss rates.In the followings we describe in more detail how we got these results.

Mass loss histories
Using some 0.1  ⊙ models, i.e., 170  ⊙  i / k = 0 and 250  ⊙  i / k = 0, 0.1 and 0.2 models, we now describe typical mass loss histories.All of these models end as PISNe though their final masses are different.Among them, the 170  ⊙ model ends as a He-rich PISN while all the 250  ⊙ models become O-rich PISNe.Total mass evolution as a function of central temperature for these models is shown in Figure 2. We find from this figure that mass loss is large at the end of H-and He-burning stages, which correspond to log  c ∼ 7.8 (K) and 8.3 (K), respectively.During the end of the H-burning stage, larger  i / k models tend to have larger mass loss rate.This is related to the surface composition shown in Figure 3 and the size of core mass.Note that we define the CO core as the region inside the mass coordinate with  ( 4 He) > 0.1 inside the Helium layer.Figure 4 shows the mass loss rates in these periods.
With the effect of rotational matter mixing, the convective core is larger for faster rotating models.For example, at log  c ∼ 7.8 (K), where the central hydrogen mass fraction is about 0.1, convective core masses for the  i / k = 0, 0.1 and 0.2 models are about 140, 160 and 170 ⊙ , respectively.Larger convective cores mean larger He-core masses and thus a star can be a He-star more easily with mass loss.As shown in Figure 3, rotating models become He-rich at the surface in this stage, and the mass-loss rates become larger than that of the equivalent non-rotating models (Figure 4).In this particular example, the 250  ⊙  i / k = 0.2 model becomes a Wolf-Rayet star in this stage.During the early He-burning stage, the evolutionary timescale is shorter than the mass loss timescale for this metallicity, and thus the total mass is roughly constant.However, at the end of the He-burning stage, stellar luminosity increases as the stellar core shrinks and mass loss rates increase again.Because of this enhanced mass loss, all the stars shown here become WR stars, though the detailed mass loss histories are rather complicated.For example, the 250  ⊙  i / k = 0 model has a larger mass loss rate at the and of He-burning than rotating models, since the surface hydrogen mass fraction is larger.Because of this effect, the final masses of 250  ⊙  i / k = 0 and 0.2 models end up being similar.The 250  ⊙  i / k = 0.1 model has larger final mass, since it has smaller hydrogen mass fraction, but has larger He core mass than the non-rotating model.
Since the 170 ⊙ model has lower luminosity than the 250 ⊙ models, the mass loss rates are smaller in general and this model ends up as a He-rich star.In Figure 4 we find rather large variation of the mass loss rate for log  c = 8.5 ∼ 8.8.This is caused by the increase of the luminosity due to the shrinkage of stellar radius as the star becomes a WR star.Though we do not go into detail here, this variation may produce some observational signatures after the supernova explosion.
Unlike the 0.1  ⊙ models, all the SMC and LMC models become O-rich stars in the PISN mass range as shown in Tables 2 and 3, because radiative mass loss is large enough to make these stars WR stars before the end of the H-burning stage.

Evolution of Rotation and Angular momentum distribution
As mentioned above, stellar rotation affects mass loss histories and the mass range for the PISN explosion.In this subsection we describe in more detail the evolution of surface rotation velocity and internal angular momentum distribution.
One of the rotational effects is the enhancement of mass loss rates according to equation (1).To see the effects of this, Figure 5 shows the evolution of  rot / crit as well as  rot / k for rotating 0.1  ⊙ , 250  ⊙ models.The  i / k = 0.1 and 0.2 models shown in this figure were discussed in the previous subsection.As we can see,  rot / crit is always far below 1.0 for these models and the mass loss enhancement is insignificant.
One may wonder, however, if we assume faster initial rotation, will  rot be larger.We show in the bottom panel just such a case for  i / k = 0.3.Indeed  rot / crit reaches 1.0 before the main sequence stage.However, after some period of the enhanced mass loss stage, the star loses angular momentum and  rot / crit never reaches 1.0 again after the main-sequence stages.The mass loss rate of this model is persistently larger than slower rotating ones, and the final mass roughly corresponds to the lowest mass PISN, which ejects only a small amount of 56 Ni.
In Figure 6, we show the distribution and time evolution of specific angular momentum  for the same models as in Figure 5. Three solid lines correspond to  at the epochs when the central hydrogen mass fraction becomes  = 0.5,  = 0.1 and the central helium mass fraction becomes  = 0.01.These three epochs corresponds to the middle and end of the hydrogen burning stage, and the end of the helium burning stage, respectively.Each panel in this figure also show a dashed curve labeled  LSO,Kerr .This is the specific angular momentum needed to get into the last stable orbit around a maximally rotating Kerr-black hole of rest-mass equal to the mass coordinate (Bardeen et al. 1972).This curve gives a rough measure of the amount of angular momentum necessary to be a GRB progenitor, and we may call the model fast rotating if  of the model exceeds  LSO,Kerr (Yoon & Langer 2005;Yoon et al. 2006).This figure show that the angular momentum loss during and after the helium burning stages is efficient and the star becomes a slow rotator even for the  i / k = 0.3 model.In fact we find that all the PISNe models shown in this paper are slow rotators, which means that CHE does not occur for this mass and metallicity range.We note that these results depend on the mass loss rate for hydrogen  poor stars, which is rather uncertain due to the paucity of observations.

Internal Abundance Evolution Before PISN Explosion
In this subsection, we explain the evolution of internal abundance distributions using the same models shown in Sec-tion 3.2, namely,  = 0.1 ⊙ , 170 ⊙ ( i = 0) and 250  ⊙ ( i / k = 0, 0.1, 0.2) models.Figures 7 to 10 show abundance distributions for these models for 4 epochs which are defined by the time when the central hydrogen mass fractions are  = 0.5 and 0.1 (Epochs 1 and 2, respectively), central helium mass fraction is  = 0.01 (Epoch 3), and at the central temperature is log  c = 9.2 (K) (Epoch 4) which corresponds to the beginning of pair-instability collapse.
First we compare two non-rotating models in Figures 7  and 8. Until Epoch 2, these models have similar abundance patterns, though the 250  ⊙ model has a larger convective core mass as we can see from the flat distributions of H and He in the core.The evolution after Epoch 3 appears different.This is because the 250  ⊙ model has a larger mass-loss rate because of larger luminosity.As a result, the star becomes a WR star during the middle of the helium burning stage.Then, the mass loss rate becomes even larger, which can be seen in Figure 4, and the model loses the entire He-rich envelope before the start of Epoch 3. The opposite applies to the 170  ⊙ model and this model maintains a thick He envelope until Epoch 4.
Next we describe rotation effects using Figures 8 to 10 for 250  ⊙ models with different rotation velocities.From these figures we find that faster rotating models have larger convective core masses because of rotational matter mixing.As a result, faster rotating models tend to have larger massloss rates, and they become WR stars earlier, which we can verify in the figures for Epoch 2. We note that as mentioned in Section 3.2, this mass-loss enhancement is not caused by the rotation induced mass loss, but simply by the increase in luminosity.
The slowly rotating model with  i / k = 0.1 in Figure 9, has a larger final core mass than the non-rotating model in Figure 8 because the increase of core mass by rotational mixing is more effective than the increase of mass-loss rate.On the other hand, for faster rotating models with  i / k ≳ 0.2, the effect of mass loss is more important and the final mass becomes smaller for faster rotation.
Though not shown in figures, with our choice of mass loss rates, low mass PISN models have He-rich envelopes only for low metallicity:  ≲ 0.1 ⊙ .As shown in Table 3, all our PISNe models with  ≥ 0.2 ⊙ have O-rich envelopes, even for the lowest mass models because of larger radiative mass loss at higher metallicities.

Hydrodynamical calculations and PISN explosions
For most models, when the central temperature reaches 10 9.2 K, we switch to the hydrodynamics code.The exceptions to this are the   /  = 0.2.After this switch, all models collapse immediately, reaching their peak temperatures within 200 to 700 seconds.After this, the velocity reverses, and a shock forms, which propagates from the cen-   ter of the star to the surface in 10-100 seconds.Maximum temperatures and explosion energies (Table 2) fall in roughly the same range as in the Pop III case (e.g.Takahashi et al. 2018).Note that models which fail to explode because they pulse or collapse to a black hole are denoted by "PPISN" and "Collapse", respectively, in the final three columns of Table 2.One difference to the Pop III case is that the 56 Ni mass is higher for comparable CO core masses.This may derive from seed metals not present in the Pop III case (Subsection 3.6), although it is hard to say for certain (Section 4).

Nucleosynthesis during PISN explosions
After performing the hydrodynamical simulations with 153 isotopes, we then post process the trajectories of those simulations with a larger network (300 isotopes).As the temperature increases, isotopes above the iron group begin to photodissintegrate.Note that this process does not occur in the Pop III case as there are no heavy isotopes.By Log T ≈ 9.5, the composition has moved away from stability towards the proton rich side (p side), forming a continuous distribution in the neutron-proton plane.As the temperature increases further, the composition shifts to higher mass and towards the neutron rich side, so that, at maximum temperature (Log T ≈ 9.8, in the extreme case), the composition spans nearly the entire 300 isotope network.Once the temperature and density decrease, the distribution retreats to the p side, before eventually decaying back to stability.Figure 11 shows the Nickel mass (symbols) of all models in the plane of the CO  2, 3).The left panel shows SMC and LMC models, while the right panel shows Z= 0.1 ⊙ models.core mass and the initial mass.The Nickel mass is generally thought to depend only on the CO core mass, but as we discuss below (subsection 3.7), using a finer grid and including metal enriched and rotating progenitors has revealed more complex behavior.Figure 12 shows the elemental yields for all models (grouped similarly to in Table 2) and Table 4 shows isotopic yields for the M = 475  ⊙ , Z = LMC,   = 0 model.As in the Pop III case, the main signature of the PISNe is the sawtooth pattern derived from the preference for even elements (Heger & Woosley 2002;Umeda & Nomoto 2002;Takahashi et al. 2018).This sawtooth pattern extends to heavier elements for more massive cores (Figure 12) which reach higher peak temperatures (Table 2).

Peculiar phenomena
In this subsection, we summarize some peculiar phenomena we have found during the calculations.It is usually expected that there is a concrete mass region for PISNe for given metallicity and initial rotational velocity, and larger initial mass leads to larger CO core mass and larger explosion energy (and 56 Ni mass) in the range.However, we find that this is not always the case because mass loss histories can be complicated.Especially near the border of the PISN range, there is often no simple relation between the initial mass and the final outcome.For example, in 0.1 ⊙  i = 0 models, the relation between the initial mass and explosion energy is not monotonic below 170 ⊙ .As a result, the 160 ⊙ model becomes a PPISN while the 150 ⊙ model becomes a PISN.We find similar "inversion" near the upper mass limits in the SMC and 0.1 ⊙ ,  i / k = 0.1 models.Related to this, we have core collapse for 330 and 340  ⊙ but a PISN explosion for 350  ⊙ in the SMC,  i / k = 0.1 models.Because of the abnormality of these findings, we also calculated a 331  ⊙ model which exploded (Table 2).These results suggest that it is not a simple task to find an exact border for the PISN range when the effects of mass-loss are large.
We also find that some models, such as SMC, 330 ⊙ ,  i / k = 0.1, 0.1  ⊙ , 200  ⊙ ,  i = 0, and 0.1  ⊙ , 170  ⊙ ,  i / k = 0, have core masses which are expected to result in PISNe, but they do not explode.These models emerge from the helium burning phase with low entropy, and thus they avoid the PISN region, despite their massive cores.We have verified the stability of these models in the hydrodynamical code and expect them to eventually undergo core-collapse and form black holes with masses right in the middle of the PISN mass gap.
Finally, we mention another peculiar result for LMC, 330  ⊙ ,  i = 0 (referred to as the 330  ⊙ model below, for brevity).This model is at the lower edge of the PISNe range, and indeed the LMC, 340  ⊙ ,  i = 0 becomes a PPISN.For the PPISN model, the first  c peak above log  c = 9.2 (K) is not high enough to blow up the entire star, and only the outer envelope is ejected during the first "pulse".In the 330  ⊙ model, on the other hand, the peak temperature is very high (log  c = 9.79 (K) ) but slightly below the photo-dissociation temperature for Fe, and thus the star explodes with very high energy and large 56 Ni production even though the CO-core mass is rather small ( CO = 63.4⊙ ).We note that similar evolution is seen in the 0.1  ⊙ , 220  ⊙ ,  i / k = 0.1 model, though the explosion energy is not as large as the 330  ⊙ model.We have verified the properties of this explosion with much stricter numerical settings and we believe that it is real.We emphasize, however, that it likely occurs only in a very narrow mass range at the lower bound of the PISNe region, and should thus be a rare occurrence.(Asplund et al. 2009) for the explosions in Tables 2, 3, where each panel corresponds to a single subsection of these two tables.The abundances are the results of the post processing, including decaying the ejecta for s (Subsection 2.4).The maximum 56 Ni mass produced by PISNe is interesting, because it determines the maximum brightness for observed PISNe.Among the models we show in Tables 2 and 3, the maximum 56 Ni mass is 40.2  ⊙ for the 0.1  ⊙ , 410  ⊙ ,  i / k = 0.2 model.This maximum mass is roughly the same as previous Pop III PISNe results (Heger & Woosley 2002;Takahashi et al. 2018), which is reasonable because the nuclear burning during pair-instabillity collapse should be roughly independent of metallicity.
Since we show our results only with 10 or 5  ⊙ intervals with respect to the initial mass, one may wonder if much higher 56 Ni mass could be obtained if we calculated models with finer grids.We thus tried to calculate with much finer grids for several cases.However, we did not find any models with 56 Ni higher than 40  ⊙ , this is partly because we do not necessary have larger CO-core mass for a larger initial mass, as mentioned in the previous sub-section.
In order to obtain continuous  CO size samples, we made the following artificial modification using the 0.1  ⊙ , 230  ⊙ ,  i / k = 0.1 model.The original model has  CO = 106.8and 56 Ni = 25.2  ⊙ .In order to have a larger final core, we artificially reduce the mass loss rate by a constant factor only for the period when log  c = 8.5 ∼ 8.9, so that we could obtain a homogeneous progenitor set with different  CO .The maximum possible  CO to be a PISN obtained in this way is 118.3  ⊙ and the results are shown in the Table 3 labeled as ' i = 230B'.As shown there, the maximum energy is 7.00 × 10 52 erg with 56 Ni = 37.5  ⊙ .This exercise demonstrates that a 56 Ni mass upper limit of ≃ 40  ⊙ is likely quite general.

DISCUSSION
First we compare our results with Pop III PISN calculations in Takahashi et al. (2018).For the Pop III models, PISNe occurred in the CO core mass rage,  CO ∼ 72 − 134 ⊙ , though the range depends on models.From Table 2, we find the range is lower for our metal-rich models,  CO ∼ 63 − 120 ⊙ , though in some parameter sets, the lower bound is higher,  CO ∼ 71 ⊙ and the upper bound is lower, There is also a difference in the relation between  CO and 56 Ni mass or explosion energy.Our metal-enriched PISNe models tend to have higher 56 Ni masses and explosion energies than the Pop III models.Without further investigations, we cannot say, at this moment, that these differences are metallicity effects, since some adopted parameters including the 12  (, ) 16  rate are different from Takahashi et al. (2018).Nevertheless, we can say that it is dangerous to simply assume that a metal-rich PISNe is the same as a Pop III PISN with the same  CO .
Comparisons of the PISN range with other previous models were made in Section 3.1 and we do not repeat it here.We simply note that the results depend on mass-loss rates, which are rather uncertain, and in particular the WR mass-loss rates are very important.
One of the significant differences in our results with previous work is that most our models become O-rich PISNe.We find He-rich PISNe only for metal poor ( ≲ 0.1 ⊙ ), low mass PISNe.If this is true, it would be difficult to find He-rich PISNe in nearby samples.This could be consistent with past observations that have not found any evidences for He-rich PISNe, and one promising PISN candidate, SN2018ibb, described in the next subsection appears to be O-rich.We even suggest that we might have already found such O-rich SNe already in the samples of many Type-I super luminous SNe (SLSNe).It is possible that we simply have overlooked such a possibility since O-rich PISNe models have not yet been reported.

Type Ic PISN candidate SN2018ibb
Very recently Schulze et al. (2023) reported that SN2018ibb is the best PISN candidate to date.Though they only say that it is Hydrogen poor, it is likely that it is Type Ic since they show evidence that this SN interacts with Oxygen-rich CSM.
If this was a He-rich SN Ib, it is difficult to explain by our models, because SN2018ibb is bright and the expected 56 Ni mass is about 30  ⊙ .On the other hand, if this was a O-rich SN Ic, several models may fit, such as LMC ∼ 470 ⊙ , SMC ∼ 380 ⊙ , 0.1  ⊙ ∼ 300 ⊙ for non-rotating models, and SMC ∼ 350 ⊙ and 0.1  ⊙ ∼ 230 ⊙ for  i / k = 0.1 models.
For SN2018ibb, evidence of CSM interaction is observed at around 600 days after the peak of bolometric luminosity.In the next subsection we discuss the possibilities of such an interaction in our model, though more detailed discussion about SN2018ibb will be given in a forthcoming paper (Nagele et al, in prep.).

Late time mass loss and SN-CSM interaction
It is well known that some SNe show evidence of interaction with CSM which has been ejected from the progenitors through mass loss.For a usual WR star with masses below 30 ⊙ , the mass loss rates are about 10 −4  ⊙ / yr or below (see e.g.Yoon 2017), which is low enough to expect SN-CSM interaction to not be observable.For SNe Ibn, the expected mass loss rates are ∼ 10 −2  ⊙ / yr or larger to explain the peak luminosity (Maeda & Moriya 2022).In Table 2 and 3, we show the mass loss rates at log  c = 9.0 (K) by |  9.0 |.Typically around this stage, the mass loss rate increases gradually with time toward the onset of pair-instabillity.Therefore, it is expected that the mass loss rates in our models are not large enough to explain the peak luminosity in the SN light curve.
To explain the observations in the late time interaction as in SN2018ibb, lower mass loss rates than ∼ 10 −2  ⊙ / yr are sufficient.The brightness of the interaction depends on the CSM density and thus depends also on mass loss wind velocity as well as the mass-loss rate.Since we use the WR mass-loss rate by SV2020, we can make use of their estimates for the wind velocity, which is typically 1000-2000 km/s for our models.For this velocity, we can estimate that that the mass-loss rate at around log  c = 9.0 is relevant to the interaction around several hundred days after the peak (Nagele et al, in prep for more detail).The tables show that |  9.0 | is typically 6 9 × 10 −4  ⊙ / yr for darker PISNe and 1 1.3 × 10 −3  ⊙ / yr for bright PISNe.Therefore, we expect that our PISNe models will have some signatures of CSM interaction.From these rough arguments we cannot definitively say that our models explain the CSM interaction of SN2018ibb, but we view the situation as promising.

Pair-instability gap in the black hole mass distribution
Since the detection of gravitational waves from black hole merger (Abbott et al. 2021), the black hole mass function has undergone increased scrutiny.To estimate the merger event rate, the existence of a gap in the BH mass distribution due to the occurrence of PISNe, called the pair-instability mass gap is important.Previously, the gap range was mentioned as 65 − 130 ⊙ (e.g. in Farmer et al. 2019;Mapelli et al. 2020), but since then several results have shown a larger value for the lower boundary.For example, in our previous paper Umeda et al. (2020), it is shown that the lower bound for Pop III stars depends on the overshooting parameter.It could be as large as 109  ⊙ for relatively small overshoot parameter, while this could be ∼ 66 ⊙ if a larger value is assumed.
The upper bound should also be model dependent and here we recall the values from Takahashi et al. (2018) for Pop III PISNe.For their non-rotating and magnetic-rotating models, the lowest BH masses above the PISN gap were 257 − 290 ⊙ for hydrogen stars, 139 − 150 ⊙ for helium cores, and 139 − 140 ⊙ for CO cores.Considering mass loss by binary interactions, the values for CO core are often mentioned as the upper bound for the gap, i.e., the upper bound is roughly 140 ⊙ .
Let us now focus on the results of this paper.We find that the maximum CO core mass for PISNe in LMC, SMC, and 0.1 ⊙ models are ≃ 109, 118, and 118  ⊙ , respectively (Tables 2 and 3).The BH masses above the PISN gap inferred from the Tables are 119.7,116.7  ⊙ for LMC, 111.0, 117.8  ⊙ for SMC and 119.2, 107.3, 115.5  ⊙ for 0.1 ⊙ models.From these results we may say that the upper bound for the gap is 107 to 120 ⊙ , which is significantly lower than the Pop III results.Furthermore, there are several models with CO core-masses inside the "gap", but with low entropy, which we expect to form black holes (Sec.3.7).
Combined with the results in Umeda et al. (2020) the PISN gap may not exist depending on the overshooting parameter if we consider both the Pop III and Pop II populations, although the gap may exist for a fixed metallicity.

What initial rotation speed is reasonable ?
We have shown that with our choice of mass loss rates, PISNe in principle occur for LMC metallicity and lower.Is this consistent with observations in which no clear evidence for PISNe has yet been obtained?We used the phrase 'in principle occur' because PISNe for LMC metallicity is possible only for slowly rotating models, say  i / k ≲ 0.1.Even for  rot / crit = 0.1, the minimum mass to be a PISN is 400  ⊙ and larger mass is required for faster rotation.
In Georgy et al. (2012), it is discussed that the typical velocity of OB stars are  i / k ∼ 0.4.If we need to apply this value for our models, the mass range of PISNe becomes very high and PISNe will be extremely rare for LMC and SMC metallicities.However, there is basically no observations for the rotation velocities for such massive stars, and we currently do not know if  > 300 ⊙ stars should initially rotate rapidly.If we find no evidence for PISNe from galaxies with LMC and SMC metallicity, the initial rotation speed may explain the non-detection.
For 0.1  ⊙ metallicity, on the other hand, PISNe are possible in a relatively low mass range, say  ∼ 250 ⊙ , even for fast rotating cases ( i / k = 0.2).We note, however, that bright PISNe with 56 Ni ∼ 30 ⊙ such as SN2018ibb requires a larger initial mass, ∼ 400 ⊙ , even for this metallicity.

CONCLUDING REMARKS
In this paper we revisited rotating metal-enriched PISN models.Previously, only a few such models had been calculated in the literature and the general properties for these models were not clear.By calculating several models using the WR mass-loss rates by SV2020, we intended to clarify, the properties of metal-enriched PISNe, such as mass ranges, produced 56 Ni mass, and mass loss histories for LMC, SMC and 0.1 ⊙ metallicities.Since we would like to know the precise values for the 56 Ni mass, for example, we carefully performed hydrodynamical simulations for collapse and explosion with a sufficiently large nuclear reaction network, as well as stellar evolution calculations.
In general, the PISNe mass range is higher for higher metallicity since mass-loss is larger.For non-rotating cases, we find that the PISN ranges are 320 − 475 ⊙ , 300 − 390 ⊙ , and 150−330 ⊙ for LMC, SMC and 0.1 ⊙ metallicities, respectively.We note that the range for 0.1 ⊙ is already close to that of Pop III SNe (Takahashi et al. 2018).
The rotation effects which we found are interesting.For slow rotation, the PISNe range shifts downwards because the effect of increased core mass is more important than increased mass loss rate which shifts the range higher.For SMC and lower metallicities, our  i / k = 0.1 models show lower PISNe range than non rotating models.On the other hand, for the LMC model, mass loss effects are larger and the range shifts higher even for  i / k = 0.1 (Table 1).For faster rotation,  i / k = 0.2 which roughly corresponds to typical rotation speed for observed OB stars, the mass ranges are larger than non-rotating models.As a result, if a PISN progenitor were to rotate as fast as  i / k ≳ 0.2, the typical mass for bright PISNe with 56 Ni mass larger than 10 ⊙ , would be more than 400 ⊙ though the mass range is lower for lower metallicities.
One of the significant differences with this work from previous similar works is that we have obtained O-rich progenitors for most of our models while previous studies obtained Herich ones.This difference mainly comes from the choice of the WR mass-loss rates, where we adopt the SV2020 rate.Since the WR mass-loss rate is uncertain especially for these very massive stars, these results should confront observations.In our case, a He-rich PISN occurs only for dimmer PISNe with metallicity  ≲ 0.1 ⊙ .It is encouraging that a recently identified PISN candiate SN2018ibb appears to be O-rich, but if we were to find bright He-rich PISNe in the future, then we would have to modify the WR mass loss rates.

Figure 1 .
Figure 1.Evolution of  rot / K for some models as a function of Log  c up to the early He burning stage.

Figure 2 .
Figure 2. Total mass evolution as a function of log  c for the selected models in  = 0.1 ⊙ .

Figure 3 .
Figure 3. Surface abundance evolution for the models in Figure 2.

Figure 4 .
Figure 4. Mass loss rates for the models in Figure 2.

Figure 11 .
Figure 11.Summary figure showing the CO core mass as a function of the initial mass for the exploding models.The symbols show the amount of 56 Ni produced in the explosion, while the colored lines connect models with the same parameters besides mass (subsections of Tables2, 3).The left panel shows SMC and LMC models, while the right panel shows Z= 0.1 ⊙ models.

Figure 12 .
Figure 12.Log elemental abundances relative to Iron, relative to solar (Asplund et al. 2009) for the explosions in Tables2, 3, where each panel corresponds to a single subsection of these two tables.The abundances are the results of the post processing, including decaying the ejecta for s (Subsection 2.4).

Table 1 .
The initial stellar mass range for becoming PISNe and the (most) abundant elements at the surface of the final models.

Table 2 .
Summary table for models denoted by the initial mass ([ ⊙ ]), metallicity, and initial rotation.The remaining columns show maximum rotation, final mass, final CO core mass (in units of  ⊙ ), mass loss rate at log  c =9.0 (K) ([ ⊙ /year]), surface composition and the mass fraction of it in ( ), explosion energy, maximum central temperature during the explosion, and 56 Ni mass as determined by the 300 isotope post processing.† This model is incorrectly denoted as an explosion in the published version.This table has the correct information.

Table 4 .
Yields (Y [ ⊙ ]) of all stable isotopes up to79Be for the 300 isotope post processing of M = 475  ⊙ , Z = LMC,   = 0.The upper section contains the yields after the composition has fully decayed.The lower section contains the yields for radioactive isotopes with Y> 10 −6 at the end of the hydrodynamics calculation ( ≈ 8000 s). Yild tables for all models are available in the supplementary material.