"Super-Kilonovae"from Massive Collapsars as Signatures of Black-Hole Birth in the Pair-instability Mass Gap

The core collapse of rapidly rotating massive ~10 Msun stars ("collapsars"), and resulting formation of hyper-accreting black holes, are a leading model for the central engines of long-duration gamma-ray bursts (GRB) and promising sources of r-process nucleosynthesis. Here, we explore the signatures of collapsars from progenitors with extremely massive helium cores>130 Msun above the pair-instability mass gap. While rapid collapse to a black hole likely precludes a prompt explosion in these systems, we demonstrate that disk outflows can generate a large quantity (up to>50 Msun) of ejecta, comprised of>5-10 Msun in r-process elements and ~0.1-1 Msun of $^{56}$Ni, expanding at velocities ~0.1c. Radioactive heating of the disk-wind ejecta powers an optical/infrared transient, with a characteristic luminosity $\sim 10^{42}$ erg s$^{-1}$ and spectral peak in the near-infrared (due to the high optical/UV opacities of lanthanide elements) similar to kilonovae from neutron star mergers, but with longer durations $\gtrsim$ 1 month. These"super-kilonovae"(superKNe) herald the birth of massive black holes>60 Msun, which, as a result of disk wind mass-loss, can populate the pair-instability mass gap 'from above' and could potentially create the binary components of GW190521. SuperKNe could be discovered via wide-field surveys such as those planned with the Roman Space Telescope or via late-time infrared follow-up observations of extremely energetic GRBs. Gravitational waves of frequency ~0.1-50 Hz from non-axisymmetric instabilities in self-gravitating massive collapsar disks are potentially detectable by proposed third-generation intermediate and high-frequency observatories at distances up to hundreds of Mpc; in contrast to the"chirp"from binary mergers, the collapsar gravitational-wave signal decreases in frequency as the disk radius grows ("sad trombone").


INTRODUCTION
The astrophysical locations which give rise to the synthesis of heavy nuclei via the rapid capture of neutrons onto lighter seed nuclei (the r-process; Burbidge et al. 1957;Cameron 1957) remains a topic of active debate (see Horowitz et al. 2019;Cowan et al. 2021;Siegel 2021 for recent reviews). Several lines of evidence, ranging from measurements of radioactive isotopes on the sea floor (e.g., Wallner et al. 2015;Hotokezaka et al. 2015) to the abundances of metal-poor stars formed in the smallest dwarf galaxies (e.g., Ji et al. 2016;Tsujimoto et al. 2017), suggest that the dominant site of the r-process is much rarer than ordinary core collapse supernovae (SNe), both in the early history of our Galaxy and to-day. The most promising contenders are the mergers of neutron star binaries (e.g., Lattimer & Schramm 1974;Symbalisty & Schramm 1982) or rare channels of core collapse SNe, such as those which give birth to a rapidly spinning magnetar (Thompson et al. 2004;Metzger et al. 2007;Winteler et al. 2012;Nishimura et al. 2015) or a hyper-accreting black hole ("collapsar"; e.g., Surman et al. 2008;Siegel et al. 2019; see also Grichener & Soker 2019). Perhaps not coincidentally, the same two types of events−neutron star mergers and collapsars−are the leading models for the central engines of gamma-ray bursts (GRB) of the short-and long-duration classes, respectively (e.g., Woosley & Bloom 2006;Berger 2014).
The radioactive decay of r-process elements in the ejecta of a neutron star mergers power a short-lived optical/infrared transient known as a kilonova (Li & Paczyński 1998;Metzger et al. 2010; . However, the large quantity of r-process ejecta 0.02 − 0.06M inferred from the kilonova to accompany GW170817, as well as the relatively low inferred outflow velocity ∼ 0.1 c of the bulk of this material (e.g., Cowperthwaite et al. 2017;Drout et al. 2017;Villar et al. 2017), do not agree with predictions from numerical relativity for the mass ejected during the early dynamical phase of the merger (see Metzger 2019; Siegel 2019; Margutti & Chornock 2020; Nakar 2020 for reviews). Instead, the dominant ejecta source in GW170817, and likely in the majority of neutron star mergers, are delayed outflows from the accretion disk which forms around the black hole (BH) or neutron star remnant (e.g., Fernández & Metzger 2013;Just et al. 2015;Siegel & Metzger 2017Fujibayashi et al. 2018). General relativistic magnetohydrodynamical (GRMHD) simulations of the long-term evolution of post-merger disk find that up to ∼ 30% of its original mass is unbound in outflows with average velocities ∼ 0.1 c (Siegel & Metzger 2017;Fujibayashi et al. 2018;Fernández et al. 2019;Christie et al. 2019;Fujibayashi et al. 2020), broadly consistent with the kilonova observed from GW170817.
As emphasized by Siegel, Barnes, & Metzger (2019), similar accretion disk outflows to those generated in neutron star mergers occur also in collapsars (see also Mac-Fadyen & Woosley 1999;Janiuk et al. 2004;Surman et al. 2006;Miller et al. 2020;Just et al. 2021). Unlike in the merger case, the collapsing stellar material feeding the disk is composed of roughly equal numbers of protons and neutrons (electron fraction Y e 0.5). However, for mass accretion rates above a critical threshold value ( 10 −1 − 10 −3 M s −1 , which depends on the effective viscosity and BH mass; Chen & Beloborodov 2007;, the inner regions of the disk are electron degenerate and act to "self-neutronize" via electron captures on protons (e.g., Beloborodov 2003), thus maintaining a low electron fraction Y e ≈ 0.1 in a regulated process (Siegel & Metzger 2017). As a result, the collapsar disk outflows, which feed on this neutronrich reservoir, can themselves possess a sufficiently high neutron concentration to enable an r-process, throughout much of the epoch over which the GRB jet is being powered. However, the details of the synthesized composition−particularly the partitioning between light and heavy r-process elements−are sensitive to the impact of neutrino absorption processes on the electron fraction of the outflowing material (Surman et al. 2006;Miller et al. 2020;Li & Siegel 2021).
In comparison to neutron star mergers, collapsars hold several complementary advantages as r-process sources . Firstly, as a result of being generated promptly from very massive stars and empirically found to occur in small dwarf galaxies at low metallicity (e.g., Fruchter et al. 2006), collapsars naturally explain the r-process enrichment in ultra-faint dwarf galaxies such as Reticulum II (e.g., Ji et al. 2016) and metal-poor stars in the Galactic halo (e.g., Brauer et al. 2021). Furthermore, if the gamma-ray luminosity of GRBs scales with BH accretion rate in the same way in mergers as in collapsars, then from the relative rate and gammaray fluence distributions of long-versus short-duration GRBs, one is led to conclude that the total mass accreted through collapsar disks over cosmic time (and hence the integrated amount of disk wind ejecta) could exceed that in neutron star mergers . Arguments based on the chemical evolution history of the early Milky Way galaxy have been made both in favor of rare SNe/collapsars Siegel et al. 2019;van de Voort et al. 2020;Yamazaki et al. 2021;Brauer et al. 2021) and mergers (Shen et al. 2015;Duggan et al. 2018; Macias & Ramirez-Ruiz 2019; Bartos & Márka 2019;Holmbeck et al. 2019;Tarumi et al. 2021) as sources of early r-process.
While the kilonova from GW170817 provided ample evidence that neutron star mergers can execute an rprocess, the same signature has not yet been seen from the SNe observed in coincidence with long GRBs. However, this fact is not necessarily constraining yet, insofar that r-process material is easier to hide in the collapsar case. In particular, a prompt and powerful supernova explosion may be required to explain the large masses of 56 Ni inferred from GRB supernova light curves (e.g., Cano 2016; Barnes et al. 2018; however, see Zenati et al. 2020). By contrast, the r-process-generating disk outflows occur over longer times, up to tens of seconds or more after collapse, commensurate with the observed duration of the long GRB. Unless efficiently mixed to the highest velocities, the r-process elements (and any associated photometric or spectroscopic signatures) are therefore buried behind several solar masses of "ordinary" supernova ejecta (dominated by α-elements such as oxygen).
Nevertheless, if present in the inner ejecta layers, rprocess elements could manifest as a late-time infrared signal ) arising from the high opacity of heavy r-process nuclei Tanaka & Hotokezaka 2013). This signal is challenging to detect given the typically large distances to GRB SNe and the late times required (at which point the emission is faint). The detection prospects will improve with the advent of the James Webb Space Telescope (JWST), particularly if the nebular spectra of lanthanide-rich material also peaks in the infrared .
Collapsars of the type observed so far as SNe may also only represent a subset of accretion-powered core collapse events. The progenitor stars which give rise to long GRBs are typically believed to possess ZAMS masses 40M with helium cores at death of 10M (e.g., Woosley & Heger 2006). Upon collapse of their iron cores, these events first go through a rapidly rotating proto-neutron star phase (Dessart et al. 2008), in which a millisecond magnetar is formed (Thompson & Duncan 1993;Raynaud et al. 2020). The strong and collimated outflow from such a magnetar during the first seconds after its birth (Thompson et al. 2004;Metzger et al. 2007), before it accretes sufficient matter to collapse into a BH, may play an important role in shock heating and in unbinding much of the outer layers of the star and generating the required large 56 Ni masses (e.g., Shankar et al. 2021).
On the other hand, the number of long GRBs with detected SNe only number around a dozen, and the majority of these are associated with the volumetrically more common but physically distinct "low luminosity" class of GRBs (e.g., Liang et al. 2007). It thus remains unclear whether the more energetic, classical long GRBs always occur in coincidence with 56 Ni-powered SNe. Indeed, luminous SNe have been ruled out for a few nominally long GRBs (e.g., Fynbo et al. 2006;Gehrels et al. 2006), though the nature of these events (e.g., whether they are actually short GRBs masquerading as collapsars) remains unclear (e.g., Zhang et al. 2007).
Within this context, we consider in this paper the fate of initially much more massive stars, those with ZAMS masses M ZAMS 260M which are predicted to evolve helium cores by the time of core collapse above the pairinstability (PI) gap 130M (e.g., Woosley et al. 2002, Woosley 2017, Renzo et al. 2020b, Farmer et al. 2020, Woosley & Heger 2021. If the initial mass function (IMF) is an indication, such stars are potentially much rarer than the ordinarily considered collapsar progenitors with M ZAMS 40M . On the other hand, if such stars are rapidly spinning (e.g., Marchant et al. 2019;Marchant & Moriya 2020)-possibly because of continuous gas accretion throughout their life (e.g. Jermyn et al. 2021;Dittmann et al. 2021)-and if these form collapsar-like disks upon collapse in proportion to their (much higher) helium core masses, their resulting yield of r-process ejecta in disk winds could be substantially greater.
Another key difference is that a prompt explosion (e.g., as attributed to a proto-magnetar above, or fallback accretion; Powell et al. 2021) is more challenging to obtain for these very massive stars. This is because (1) the nominal timescale for BH formation is much faster, within 0.1 s, due to large masses exceeding 2M and high compactness of their iron cores (e.g., Renzo et al. 2020b); (2) their larger 10 53 erg gravitational binding energies relative to those of lower mass helium cores 10 52 erg exceed the rotational energies of even maximally spinning neutron stars. As a result of the assuredly failed initial explosion of post-PI cores, these systems are unlikely to eject a large quantity of prompt, shock-synthesized 56 Ni and unprocessed stellar material (however, see Fujibayashi et al. 2021). Instead, the bulk of the ejecta will arise over longer timescales from disk outflows, which-scaling up from low-mass collapsarscould amount to 10M of r-process and Fe-group elements (including 56 Ni).
Rather than the usual picture of GRB SNe, the type of collapse transient above the PI-gap we envision is in some ways more akin to a scaled-up neutron star merger. At risk of committing etymological heresy, we therefore refer to these massive collapsar transient events as "super-kilonovae" (superKNe). As we shall discuss, if superKNe exist, their long durations and red colors may render them potentially identifiable through either follow-up infrared observations of long GRBs (e.g., with JWST), or blindly in surveys with the Vera Rubin Observatory (Tyson 2002) or the Nancy Grace Roman Space Telescope (Roman; Spergel et al. 2015).
The gravitational wave observatory LIGO/Virgo detected a binary BH merger, GW190521, for which both binary components of masses ∼ 85M and ∼ 66M , respectively (Abbott et al. 2020), were inside the nominal PI mass gap. 1 Tentative evidence suggests an effective high BH spin of the progenitor binary, albeit with the spin axis misaligned with the orbital momentum axis (however, see Mandel & Fragos 2020;Nitz & Capano 2021). These unusual properties have motivated a number of theoretical studies proposing new ways to populate the PI mass gap, such as through dynamical stellar mergers (e.g. Di Carlo et al. 2019Carlo et al. , 2020Renzo et al. 2020a), hierarchical black hole mergers in dense environments (e.g., Antonini & Rasio 2016;Yang et al. 2019;Tagawa et al. 2021;Gerosa & Fishbach 2021), modifying stellar physics at low metallicity (e.g. Farrell et al. 2021;Vink et al. 2021), or through external gas accre-tion (e.g., Safarzadeh & Haiman 2020). As we shall describe, if both the BHs acquired their low expected masses and high spin as a result of inefficient disk accretion, superKN events of the type envisioned here provide a novel single-star channel for filling the PI mass gap "from above". This paper is organized as follows. Sec. 2 presents a semi-analytic model for the collapse of rotating massive stars, their accretion disks and disk wind ejecta, and resulting heavy element nucleosynthesis, which builds on earlier work in Siegel et al. (2019). Calibrating the model such that collapsars generate BH accretion events consistent with the observed properties of long GRB jets, we then apply the model to more massive 130M progenitors above the PI mass gap. Using our results for the disk wind ejecta, in Sec. 3 we calculate the light curves and spectra of their superKN emission by means of Monte Carlo radiative transfer simulations. Sec. 4 explores the prospects for discovering superKNe by future optical/infrared surveys or in follow-up observations of long GRBs. Section 5 discusses several implications of our findings, including gravitational-wave emission from self-gravitating phases of the collapsar disk evolution; the astrophysical origin of GW190521; and the luminous radio and optical emission that results from the superKN ejecta interacting with surrounding gas. Sec. 6 summarizes our results.

Stellar models
To model the pre-collapse structure of the superKN progenitors, we employ the MESA stellar evolution models of Renzo et al. (2020b), publicly available at https: //zenodo.org/record/3406357. The simulations start from naked helium cores of metallicity Z = 0.001 which are then self-consistently evolved from helium core ignition, through possible (pulsational) pair-instability (PPI), to the onset of core-collapse (defined as when the radial in-fall velocity exceeds 1000 km s −1 ). We label the input stellar models according to their initial helium core mass, e.g., model 200.25 corresponds to M He,init = 200.25M , and focus on models "above" the PI gap, which do not experience pair-instability driven pulses.
The models are computed using a 22-isotope nuclear reaction network, which is sufficient to capture the bulk of the energy generation throughout the stellar evolution, but cannot accurately capture the weakinteractions in the innermost core (e.g., Farmer et al. 2016). However, the deepest layers of the core promptly fall into the newly formed BH (see Sec. 2.2) and hence do not contribute to the accretion disk and its outflows.
These models were evolved without rotation, which we instead artificially impose at the point of core collapse (see Sec. 2.2). The main effect of rotation during the precore-collapse evolution is mixing at the core-envelope interface, which leads to more massive helium cores for a given initial mass. In the extreme case of chemically homogeneous evolution (Maeder & Meynet 2000), the entire star may become a helium core. This will impact how many stars develop core masses reaching into the PI/pulsational PI regime or beyond, and thus the predicted population statistics. However, because Renzo et al. (2020b) only simulate the helium core, this does not affect our present study. Rotation can also enhance the wind mass-loss rate (e.g., Langer 1998), and increase the radius in the outer layers at the rotational equator by up to 50%, which is neglected in the progenitors we use. Finally, by adding centrifugal support to the core, rotation can modestly increase the PI/pulsational PI mass range (e.g., Glatzel et al. 1985). 2 For sufficiently large initial core masses M He,init 200 M , the final mass at collapse would nominally produce a BH above the PI mass gap (neglecting subsequent mass-loss in accretion disk outflows, as explored in the present study). This arises because the gravitational energy released by the PI-driven collapse acts to photo-disintegrate nuclei produced in the thermonuclear explosion, instead of generating outwards bulk motion (e.g., Bond et al. 1984). Since these models do not experience pulses of mass loss, their pre-collapse total mass is determined by the assumed wind mass-loss prescription: the minimum final helium core mass above the PI mass gap for our MESA models is M He,fin 125 M . Renzo et al. (2020b) estimated the corresponding final BH mass (again, neglecting post-collapse disk outflows) as the total baryonic mass with binding energy < 10 48 ergs, which effectively corresponds to the total final mass within a few 0.1 M (e.g., Farmer et al. 2019;Renzo et al. 2020b).
A key ingredient in modeling fallback accretion is the radial density profile of the star at collapse (Sec. 2.2). Despite their large masses, helium stars above the PI gap remain compact throughout their lives, never expanding 2 For example, using a setup similar to Renzo et al. (2020b), Marchant & Moriya (2020) study the impact of an initial rotation frequency ω/ω crit = 0.9, where ω crit ≡ (1 − L /L Edd )GM /R 3 and L /L Edd is the stellar luminosity in units of the Eddington luminosity. They found a ∼ 4% (∼ 15%) increase in the maximum BH mass below the PI mass gap assuming angular momentum is transported by a Spruit-Tayler dynamo (assuming no angular momentum transport). The stronger angular momentum coupling found by  would likely result in an even more modest effect. as a result of PI pulses. Their typical radii R ≈ 10 R are similar to the normally considered Wolf-Rayet progenitors of GRBs (e.g., Woosley & Heger 2006; Fig. 1). If stars in this mass range reach core collapse with their hydrogen envelope intact (e.g., for sufficiently low metallicity, as in population III stars) their radii could be considerably larger; however, no red supergiants of this mass have yet been observed.
We also employ the 15 M and 20 M single (hydrogen-rich) star models from Heger et al. (2000) to test our collapsar model on more canonical long-GRB progenitors (see Appendix D for a discussion of results).
These are computed starting from a surface equatorial velocity of 200 km s −1 at ZAMS and assume that mean molecular weight gradients do not impede rotational mixing (f µ = 0, "weak molecular weight barriers"). They are labeled E15 and E20 respectively, and are publicly available at https://2sn.org/stellarevolution/rotation/.

Collapsar Model
The masses and composition of the superKNe ejecta are computed by modeling the collapse of a progenitor star. Depending on the stellar angular momentum profile, the collapse and fallback of envelope material leads to the formation of an accretion disk, which gives rise to massive neutron-rich disk outflows Miller et al. 2020;Just et al. 2021). Although rotation profiles of massive stars at the time of core collapse, in particular of those above the PI mass gap considered here, are highly uncertain (e.g., Heger et al. 2000;Marchant & Moriya 2020), the specific angular momentum j z generally increases with stellar radius. In-falling stellar material thus circularizes at increasingly larger radii from the BH with time.
We endow the stellar models with mass M and radius R = R He,fin at the time of core collapse (Sec. 2.1) with an angular momentum profile that assumes rigid rotation on spherical shells, with angular velocity Ω(r, θ) = Ω(r). This results in where r, θ are the radial and polar angle coordinates, respectively. We adopt a general parametrized angular momentum profile of the form where r b , p, and f K are free parameters. This corresponds to a low-density 'envelope', composed primarily of helium in the models considered here, rotating at a fraction f K < 1 of the local Keplerian angular momentum j K = GM enc (r)r, where M enc (r) is the mass enclosed interior to radius r, and an inner 'core', in which rotation is suppressed by a power-law with index p relative to the fraction of local break-up rotation adopted for the envelope. Although the parameter values (r b , p, f K ) are uncertain, as we discuss below, they can be "calibrated" to produce the timescales and energetics of the disk accretion consistent with the observed properties of long GRB jets. Figure 1 illustrates the parametrized rotation profile for model 250.25.
Assuming an axisymmetric rotating star, we discretize the progenitor stellar model into (n r , n θ ) mass elements, logarithmically spaced in stellar radius r, and uniformly spaced in cos θ. The angular resolution is chosen sufficiently high (typically n θ = 1001) that the accuracy in numerically computing global quantities by integration (total mass, total fall-back mass, etc.) is dominated by the finite radial resolution of the stellar progenitor models. Defining t = 0 as the onset of core-collapse, a given stellar layer at radius r will start to collapse onto the centre upon the sound travel time t s (r) = r 0 c −1 s (r)dr from the centre to r. Due to its finite angular momentum, a given fluid element will do so on an eccentric trajectory and circularize on the equatorial plane at time (cf. Kumar et al. 2008) t circ (r, θ) = t s (r) (3)
The innermost parts of the stellar core may not possess sufficient angular momentum to circularize in an accretion disk and, instead, directly collapse into a BH. We define the initial BH as a 'seed BH' formed by the innermost stellar layers up to radius r •,0 with enclosed mass M •,0 = 0.5 M , a safe assumption for all stellar models considered here. This seed BH has dimensionless spin parameter where J enc (r) is the enclosed angular momentum, G is the gravitational constant, and c is the speed of light. The corresponding innermost stable circular orbit (ISCO) is given by (Bardeen et al. 1972) Upon initial BH formation, we follow the collapse of the outer stellar layers according to Eqs. (3) and (4) and distinguish between mass elements that circularize outside the BH to form a disk (r circ (r(t), θ) > r ISCO (t)), giving rise to a 'disk feeding rate'ṁ fb,disk , and those that directly fall into the BH without accreting through a disk (r circ (r(t), θ) ≤ r ISCO (t)), giving rise to a direct fallback rate onto the BHṁ fb,• . Here, r(t) refers to the radius of a stellar element at polar angle θ which circularizes at time t in the equatorial plane. We denote the associated rates of angular momentum supplied to the disk and the BH byJ fb,disk andJ fb,• , respectively.
We follow the evolution of the BH, disk, and ejecta properties solving the following equations: Here, r 2 ISCO − a • r g (r ISCO r g /2) 1/2 + a 2 • r 2 g /4 r ISCO r 2 ISCO − 3r ISCO r g /2 + a • r g (r ISCO r g /2) 1/2 1/2 is the specific angular momentum of a fluid element at the ISCO of the BH with mass M • , spin a • = cJ • /GM 2 • , and gravitational radius r g = 2GM • /c 2 ( Bardeen et al. 1972). Mass is accreted onto the BH at a ratė where is the viscous timescale of the disk, with α being the standard dimensionless disk viscosity (Shakura & Sunyaev 1973), is the Keplerian angular velocity of the disk, and h z,disk its scale height (we take h z,disk ≈ 0.5 as a fiducial value). The disk radius r disk (t) is defined by the current disk mass and angular momentum, The disk accretion flow gives rise to powerful outflows with mass-loss at a ratė and associated angular momentum loss ratė Neutrinos cool the disk effectively above the critical "ignition" accretion rate for weak interactions (Chen & Beloborodov 2007;Siegel et al. 2019;De & Siegel 2020), which is approximately given by (see Motivated by the findings of GRMHD simulations of neutrino-cooled accretion flows (Siegel & Metzger 2018;Fernández et al. 2019;Siegel et al. 2019;De & Siegel 2020), we assume that for high accretion rates >Ṁ ign a fraction 1 − f acc ≈ 0.3 of the disk mass is unbound in outflows. This fraction is assumed to increase to ≈ 0.6 belowṀ ign , under the assumption that inefficient cooling will result in excess heating and outflow production (e.g., Blandford & Begelman 1999;De & Siegel 2020). Similarly, we assume that enhanced outflow production occurs also at very high accretion rates, for which neutrinos become effectively trapped in the optically thick accretion disk and are advected into the BH before radiating. This threshold "trapping" accretion rate is given by (see Appendix A.3) Insofar asṀ ν,trap scales in the same way with the (growing) BH mass asṀ ign , we find this trapped regime is of little practical importance in our models. In summary, the accretion efficiency is given by Eqs. (9)-(14) allow a calculation of the total ejecta mass M ejecta obtained from a particular collapsar model.
We evolve this set of coupled differential equations numerically until all stellar progenitor material has collapsed and has either been accreted onto the BH or been ejected into outflows. Note that these equations explicitly conserve mass and angular momentum. Time stepping is equidistant in log t and chosen sufficiently high, such that i) the accuracy of the total fallback mass is dominated by the radial resolution of the provided stellar model (see Appendix C) and that ii) conservation of total mass and angular momentum in Eqs. (9)-(14) is achieved to better than 10 −14 relative accuracy for all model runs.
Once a disk forms around the BH and its accretion rateṁ acc exceeds 10 −4 M s −1 , we assume that a relativistic jet emerges, powerful enough to drill through the remaining outer layers in the polar region. This threshold is motivated by typical GRB luminosities L γ ∼ 2 × 10 50 erg s −1 (Goldstein et al. 2016), which, if accretion powered, require an accretion rate of at leasṫ M ∼ L γ /c 2 ∼ 1.1 × 10 −4 M s −1 . If this threshold is surpassed, we ignore any remaining material in the polar regions θ < θ jet and θ > 180 • − θ jet for the subsequent fallback process. This material has little effect on the total quantity of material accreted through the disk as it predominantly falls into the BH directly due to the low angular momentum in these regions. However, it has a slight indirect effect on nucleosynthesis by modifying the BH mass (see below). As a fiducial value, we take θ jet = 30 • . We further justify the existence of such a successful jet a posteriori by the fact that our models reach the regimeṀ > 10 −4 M s −1 favorable for powering typical observed long GRBs including the time necessary for the jet to drill through the stellar envelope (see Sec. 2.3).
The fallback process may in some cases give rise to massive, gravitationally unstable accretion disks. In this limit, the disk mass becomes comparable to the BH mass, and our assumption of a Kerr metric would not be justified anymore. We estimate this instability region by monitoring the ratio of self-gravity to external gravitational acceleration by the BH potential (Paczynski 1978;Gammie 2001), where Σ is the disk's surface density. If Q < 1, we remove excess disk mass by enhancing accretion and wind production such as to restore Q = 1. This is motivated by the fact that gravitationally unstable disks tend to self-regulate by increased angular momentum transport via gravitationally driven turbulence, thereby increas-ing the accretion rate and reducing the disk mass until Q > 1 (e.g., Gammie 2001). The composition of the disk wind ejecta at a given time depends most sensitively on the instantaneous accretion rate . Following Siegel et al. (2019), Miller et al. (2020), and Li & Siegel (2021), we define the following accretion regimes: Here,Ṁ ν,r−p represents a threshold between production of lanthanides and first-to-second peak r-process elements only, accounting for the fact that increased neutrino irradiation at high accretion rates tends to raise the electron fraction above ≈ 0.25 required for lanthanide production (e.g., Lippuner & Roberts 2015). We assume this threshold scales with the accretion rate above which the inner disk becomes optically thick to neutrinos, which we estimate as (see Appendix A.3) This expression has been normalized using numerical results by Siegel et al. (2019) and Miller et al. (2020) for M • ≈ 3M . Additionally including the effects of neutrino fast flavor conversions may increaseṀ ν,r−p significantly (Li & Siegel 2021), possibly up to ≈ 1M s −1 or higher for such light BHs. We therefore treat the normalization as a free parameter and explore different scenarios in which the value is scaled up by a factor of ten. Below the ignition rateṀ ign , r-process production ceases abruptly and nucleosynthesis in the outflows with roughly equal numbers of neutrons and protons (Y e 0.5) only proceeds up to iron-peak elements . A large fraction of the outflowing material in this epoch remains, however, as 4 He instead of forming heavier isotopes. This is due to the slow rate of the triple-α reaction needed to create seed nuclei when Y e ≈ 0.5, relative to the much faster neutron-catalyzed reaction 4 He(αn, γ) 9 Be(α, n) 12 C that operates when Y e 0.5 (Woosley & Hoffman 1992). Here, we employ a simple model to estimate the yield of 56 Ni in such Y e 0.5 disk outflows, similar to Siegel et al. (2019) A requisite for the synthesis of 56 Ni in disk outflows is that nuclei from stellar fallback material are dissociated into individual nucleons once entering the inner part of the accretion disk. At late times during the accretion process, the disk densities and temperatures may not be high enough to ensure full dissociation. We estimate the transition time t diss to this state by evaluating the conditions under which only 50% of α particles are dissociated in the disk (see Appendix B). For t > t diss we ignore any potential further nucleosynthesis in disk outflows.

Collapsar Model Results
We start in Sec. 2.3.1 by walking through the evolution of the collapse and mass-ejection process for a representative model corresponding to a star above the nominal PI mass gap. Appendix D presents the results of our model when applied to "ordinary" low-mass collapsars (with BH masses below the PI mass gap), demonstrating that for the fiducial range of parameters considered in this work, we obtain properties in agreement with observed GRBs and previously predicted r-process ejecta. Using the same parameters (now "calibrated" to reproduce the properties of ordinary collapsars) we present in Sec. 2.3.2 a parameter exploration of ejecta masses and nuclear compositions for massive collapsars above the PI mass gap. Figure 2 illustrates the collapse evolution of model 250.25 with representative rotation parameters of p = 4.5, f K = 0.3, and r b = 1.5 × 10 9 cm. Upon seed BH formation, the BH grows rapidly in mass and spin through low-angular momentum material at small radii directly falling into the seed BH before circularizing. Only after a few seconds and accretion of ≈ 10 M of the inner layers, first material starts to circularize outside the BH horizon to form an accretion disk (cf. Fig. 2, top and bottom panel), and initiate accretion onto the BH through a disk in addition to direct infall.

Basic Model Evolution
Direct fallback onto the BH subsides after the accretion of about 20 M in this model (cf. Fig. 2, top panel) when a significant fraction of low angular momentum material residing in the polar region of the progenitor model has fallen into the BH. Further BH growth then proceeds almost entirely through disk accretion. This initial direct fallback episode partially clears up the polar regions for a relativistic jet to propagate through the outer stellar layers to eventually break out of the star and generate a long gamma-ray burst. Around the same time, a significant fallback rate onto the disk sets in (cf. Fig. 2, top panel) to establish a heavy ∼ 15 M accretion disk on a timescale of a few seconds (cf. Fig. 2, bottom panel). The disk accretion rate onto the BH, m acc , quickly exceedsṀ ign and we assume that a relativistic jet forms. This removes the remaining lowangular momentum material in the polar regions and thus results in suppression of direct fallback onto the BH, which becomes negligible compared to disk fallback (cf. rect consequences on nucleosynthesis in the disk winds due to its effect on the BH mass (cf. Eq. (26)). For somewhat larger values of r b , the situation changes and direct fallback onto the BH may extend to late times even in the presence of a jet, due to the overall lower angular momentum budget of the progenitor star outside the polar cone with opening angle θ jet . For more extreme scenarios, fallback onto the disk may become close to non-existent. As soon as the disk forms, most angular momentum resides in the disk rather than the BH in this model (cf. Fig. 2, center panel). The majority of this is being blown off in the ejecta, while a subdominant amount is transferred to the BH as disk matter gradually accretes through the ISCO onto the BH. For significantly larger values of r b this trend reverses, and most angular momentum is transferred to the BH rather than the ejecta as less material accretes through a disk.
The top panel of Fig. 3 shows the history of ejecta production in the model discussed above. Shown is the instantaneous accretion state M disk /t visc of the disk as a function of the cumulative ejected wind mass, together with the nucleosynthesis regimes defined in Eq. (26). This evolution shows a 'sweep' through most nucleosynthesis regimes, typical of the models considered here. Nucleosynthesis regimes change during the evolution as a result of the BH mass growth and can be more dramatic in some cases than illustrated here. Outflows are first created in the regime of a main r-process with lanthanide production, during which the bulk of the wind ejecta is produced. The remaining ≈ 10 M of ejecta originate in a regime that mostly ejects α-particles and ∼ 0.1 M of 56 Ni. The bottom panel of Fig. 3 illustrates 56 Ni production in this regime. Shown are the mass fraction of 56 Ni produced in disk outflows according to Eq. (B14) as well as the mass fraction of α-particles in the accretion disk according to Eq. (B18). The vertical dashed line indicates the dissociation time t diss after which < 50% of α-particles are dissociated into individual nucleons in the accretion disk (Sec. 2.2, Appendix B). As a conservative estimate, for t > t diss , we ignore any further production of 56 Ni according to Eq. (B14) as the required free nucleons become unavailable. However, this represents only a slight correction in most cases, as by far the dominant amount of 56 Ni is typically produced before t = t diss .

Parameter Study of Massive Collapsars
Before systematically applying our model across the parameter space of massive collapsars, we first apply it to 'ordinary' collapsars of stars well below the PI mass gap, the results of which we describe in Appendix D. We use the progenitor models of Heger et al. (2000) as representative of typical stellar progenitors of canonical long GRBs (MacFadyen & Woosley 1999). Our results for the nucleosynthesis yields of the disk outflows as a function of the parameters {f K , r b } which enter the progenitor angular momentum profile (Fig. 16), broadly agree with those previously presented in Siegel et al. (2019), though some quantitative differences arise due to our more detailed treatment of different regimes of BH accretion (see Appendix D for a discussion). Our low-mass collapsar models also exhibit BH accretion timescales and energetics of putative jet activity in agreement with long GRB observations. We can therefore claim a rough "calibration" of our model across the adopted parameter space of progenitor angular momentum properties, allowing for more confidence when extrapolating to the regime of more massive collapsars described below.   Figs. 4 and 18 summarize our results for the ejecta and GRB properties for model 250.25 as a representative example of a stellar model above the PI mass gap, in the parameter space {f K , r b }. The top panel of Fig. 4 shows that, even for a progenitor mass M = 150 M at the onset of collapse (that is, well above the PI mass gap), the final BH remnant can populate the entire mass gap between ∼ 55 M − 130 M (for typical parameter values), depending on the rotation profile at the onset of collapse. Labelled contours indicate the inferred primary mass of GW190521, together with its 90% confidence limits. We focus on this region of the parameter space in what follows, insofar that superKNe generated from such events probe BHs formed in the PI mass gap.
As in case of the low-mass collapsars (Appendix D), our results are not sensitive to the precise value of the power-law coefficient p, which we thus ignore in what follows. We find ubiquitous r-process production throughout the parameter space, ranging between ∼ 0.1−2.3 M of heavy (A > 136) r-process material including lanthanides and ∼ 1 − 29 M of light (A < 136) r-process elements. Additionally, between ∼ 0.05−1 M of 56 Ni are synthesized in the ejecta.
Interestingly, the region of highest r-process production is well aligned with intermediate final BH masses in a range similar to the GW190521 confidence region (Sec. 5.3). For large r b the outer stellar layers possess too little angular momentum to form massive accretion disks that give rise to copious r-process ejecta, as most material directly falls into the BH. On the other hand, for small values of r b and high values of f K massive disks form; however, high angular momentum leads to large disk radii r disk and associated viscous timescales, such that the accretion rate drops below the required thresholds for r-process production for most of the accretion process. This occurs despite the presence of spiral modes in this regime, which tend to increase the accretion rate (Sec. 2.2). Most r-process material (both light and heavy) is synthesized for small values of both r b and f K , which represents the optimal compromise between high angular momentum and sufficient compactness of the accretion disk. We discuss the possible contribution of massive collapsars to the long GRB population in Sec. 4.3.
For use in our subsequent light curve models (Sec. 3), we decompose the ejecta content of the collapsar models further into mass fractions of several constituents of interest. Assuming full mixing of all ejecta content (see also Sec. 3.1), we calculate the mass fraction X La of lanthanides (atomic mass number 136 A 176) based on the amount of main r-process material, assuming the solar r-process abundance pattern (Arnould et al. 2007) motivated by the results of Siegel et al. (2019). A mass fraction X lrp for light r-process elements is based on the combined mass fraction of light r-process material only plus the fraction of main r-process ejecta with A < 136 when applying the solar r-process abundance pattern. Finally, we also compute the mass fraction X Ni of 56 Ni. Results are depicted in Fig. 5. For concreteness, we select several models along iso-mass contours for the final BH mass within the GW190521 confidence region and report the corresponding ejecta parameters in Tab. 1.

SUPER-KILONOVA EMISSION
As the disk outflows expand away from the BH, the ejecta shell they form eventually gives rise to optical/infrared emission powered by radioactive decay (the "superKN").

Analytic Estimates
We begin with analytic estimates of the superKN properties. The total ejecta mass M ej is comprised of up to three main components: (1) radioactive rprocess nuclei, mass fraction X rp ; (2) radioactive 56 Ni, X Ni ; (3) non-radioactive 4 He, X He = 1 − X rp − X Ni (also a placeholder for other non-radioactive elements). Typical values for our fiducial models (Sec. As described in the previous section, the total r-process mass fraction can be further subdivided into that of light r-process nuclei X lrb and of lanthanides X La . For simplicity, throughout this section we assume the ejecta are mixed homogeneously into a single approximately spherical shell. Physically, such mixing could result from hydrodynamic instabilities that develop between different components of the radial and temporally-dependent disk winds and or due to its interaction with the GRB jet (e.g., Gottlieb et al. 2021).
The light curve will peak roughly when the expansion timescale equals the photon diffusion timescale (e.g., Arnett 1982), where v ej is the average ejecta velocity. The effective gray opacity κ varies in kilonovae from 1 cm 2 g −1 for ejecta dominated by light r-process species, to ∼ 10 − 30 cm 2 g −1 for ejecta containing a sizable quantity of lanthanide atoms and ions (e.g., Kasen et al. 2013;Tanaka et al. 2020). However, κ will be smaller than 9.0 9.2 9.4 9.6 9. these estimates in the superKN case due to the large mass fraction of light elements, X He ∼ 0.5 − 0.9, which contribute negligibly to the opacity. For our analytical estimates below, we linearly interpolate κ between 0.03 cm 2 g −1 (at X La = 10 −4 ) and 3 cm 2 g −1 (at X La ≥ 0.2), which we find results in reasonable agreement with the detailed radiation transport calculations present in Sec. 3.2. The peak luminosity and effective temperature can also be estimated using analytic formulae (e.g., Metzger 2019), where we have used the radioactive heating rate of rprocess nuclei from Metzger et al. (2010) with an assumed thermalization efficiency of 50%. Near peak light t ∼ t pk ∼ 100 d, the specific radioactive heating rate of 56 Ni is ∼ 10 − 30 times higher than that of rprocess elements (e.g., Siegel et al. 2019). Given values X Ni /X rp ∼ 0.01 − 0.05 for most of our disk outflow models, L pk is moderately underestimated by Eq. (30), which neglects 56 Ni heating. Fig. 6 shows the predicted peak timescale, luminosity, and effective temperature of the superKN emission in the parameter space {f K , r b } for the fiducial model 250.25. For the same parameters which generate remnant BHs with masses in the PI gap, we predict peak luminosities L pk ∼ 10 42 erg s −1 and characteristic durations of months. Though similar to other types of SNe in duration, superKNe are characterized by significantly cooler emission (T eff ≈ 1000 K), as confirmed by radiative transfer calculations presented in the next section.

Model Selection and Parameters
To elaborate on the estimates of §3.1, we carried out detailed radiation transport simulations for five ejecta models whose properties (M ej , M Ni , M lrp , and M La ) span the space defined by the subset of simulations that produced BHs within the mass gap (60 M • /M 130), i.e., models that fall between the two cyan contours of Fig. 5. (See also Woosley 2017;Farmer et al. 2019;Renzo et al. 2020b;Farmer et al. 2020;Costa et al. 2021;Mehta et al. 2021). The parameters of the mass gap models are largely confined to a plane in M ej -M Ni -M lrp -M La space, making it straightforward to select a handful of characteristic parameters from the full set. We used the KMeans routine of sklearn (Pedregosa et al. 2011) to divide our models into four clusters, and adopted the positions of the cluster centers as four representative super-kilonova models. However, a small fraction of the mass-gap models occupy a distinct region of the parameter space, having large M ej , but little to no nucleosynthetic products heavier than He. Since these models were not captured by our clusters, we added  a fifth model to explore the edge case of a high-mass, nickel-free outflow. Our five models are listed in Tab. 2. We performed for the models of Tab. 2 onedimensional radiation transport calculations carried out with Monte Carlo radiation transport code Sedona (Kasen et al. 2006;Kasen et al. in prep.). We adopted for each model a density profile such that the mass external to the velocity coordinate v follows a power-law, Above, the minimum ejecta velocity v min is determined by the characteristic velocity v ej = (2E kin /M ej ) 1/2 (with E kin the ejecta kinetic energy), and the choice of power- We take α = 2.5 and v ej = 0.1c for all models, consistent with predictions of accretion disk outflow velocities (e.g. Fernández et al. 2015;Siegel et al. 2019). The opacity of the outflowing gas, and therefore the nature of the transients' electromagnetic emission, is sensitive to the abundance pattern in the ejecta. Specifically, lanthanides and actinides, and to a lesser extent elements in the d -block of the periodic table, contribute a high opacity, while the opacities of s-and p-block elements is significantly lower Tanaka et al. 2020).
In this work, we predict the synthesis of helium, 56 Ni, and light and heavy r -process material, but do not carry out detailed nucleosynthesis calculations, e.g. by postprocessing fluid trajectories. The composition of each model is then solely a function of its M ej , M Ni , M lrp , and X La . We assume that heavy (A > 136) r -process material is 41% lanthanides and actinides by mass, equal to the solar value of M La /M A>136 . The remainder is split between d -block and s/p-block elements (54% and 5% by mass, respectively). For light r -process material, X La = 0. We estimated it comprises 95% (5%) d -block (s-/p-block) elements by mass.
The composition adopted for our radiation transport models is limited by both our imperfect knowledge of the details of nucleosynthesis and incomplete atomic data of the sort necessary to calculate photon opacities in the ejecta. Lanthanide and actinide mass (M La ) is divided among lanthanide elements following the solar pattern, with one adjustment: because the required atomic data is not available for atomic number Z=71, we redistribute the solar mass fraction of Z=71 to Z=70.
Atomic data is also unavailable for most of the dblock elements produced by r -process (whether heavy or light). We thus distribute d -block mass evenly among elements with Z = 21 − 28 (excluding Z = 23 for lack of data), artificially increasing the mass numbers to A ∼ 90 to avoid overestimating the ion number density. All rprocess s-and p-block material is modeled by the lowopacity filler Ca (Z = 20). 4 He and 56 Ni (as well as its daughter products 56 Co and 56 Fe) are straightforward to incorporate into the composition.
Our radiation transport simulations include radioactivity from both the 56 Ni decay chain and from the rprocess. We explicitly track energy loss by γ-rays from 56 Ni and 56 Co, and assume that positrons from 56 Co decay thermalize immediately upon production. We model r -process radioactivity using the results of Lippuner & Roberts (2015) for an outflow with (Y e , s B , τ exp ) = (0.13, 32k B , 0.84 ms), with s B the initial entropy per baryon and τ exp the expansion timescale. To account for thermalization, we adjust the absolute radioactive heating rate following the analytic prescription of Barnes et al. (2016).

Radiation Transport Results
The bolometric light curves of models A through E are presented in Fig. 7. For comparison, we also show the light curves of typical SNe of various subtypes: Type Ia SN 2011fe (Tsvetkov et al. 2013), Type Ic-bl SN 2002ap (Tomita et al. 2006), Type IIp SN 2013ab (Bose et al. 2015), and the electron-capture SN 2018dz (Hiramatsu et al. 2021) . The superKN light curves exhibit considerable diversity, which is not surprising given the large ranges of ejecta and radioactive masses these systems may produce. As would be expected from simple Arnett-style (Arnett 1982) arguments, higher masses are generally associated with longer light-curve durations. This can be seen in the progression from model A to model D.
However, as model E demonstrates, the shape of the light curve also depends on the presence of 56 Ni in the ejecta. While the mass of r -process material burned in superKN outflows greatly exceeds that of 56 Ni, the energy generated by the 56 Ni decay chain, per unit mass, exceeds that of r -process decay by orders of magnitude (e.g., Metzger et al. 2010;Siegel et al. 2019). When 56 Ni is present, it can be the main source of radiation energy for the transient. As a result of the long half-life of the 56 Ni daughter 56 Co (τ Co 1/2 ≈ 77 days), the energy generation rate for 56 Ni-producing systems is declining slowly just around the time the light curves reach their maxima. The effect is a more extended light curve (see Barnes et al. 2021 for more detailed discussions).
Model E, which produces no 56 Ni, has a relatively short (∼month) duration, despite its high mass (M ej = 60M ), owing to the steep decline of the r -process radioactivity that is its only source of energy. The qualitative difference between models that burn even small amounts of 56 Ni and models that burn none points to the importance of a careful treatment of nucleosynthesis in disk outflows.
As is apparent from Fig. 7, the diversity of superKN light curves allows them to mimic other types of SNe. While superKNe do not produce sufficient 56 Ni to approach the luminosity of SNe Ia, they can, at various epochs, mimic the bolometric light curves of SNe Ic-bl, SNe IIp as well as electron-capture SNe. However, the high opacity of the r -process-enriched ejecta pushes the superKN emission to redder wavelengths than what is observed for other classes of SNe. This is illustrated in Fig. 8, which shows the normalized spectra for models A through E at bolometric peak.
Unlike other types of SNe, most of the superKN flux emerges at near-and even mid-infrared wavelengths. This is likely due to a combination of lower radioactive heating per unit ejecta mass, as well as the high opacity from r -process elements (particularly lanthanides and actinides) and the high M ej , which work in concert to increase the optical depth across the ejecta and push the photosphere out to the exterior where temperatures are cooler.
A second distinguishing feature of superKNe is their broad absorption features. These are a product of our assumed ejecta velocities (v ej = 0.1c), which are higher than what is inferred for all supernova other than the hyper-energetic SNe Ic-bl. And while SNe Ic-bl produce spectra with similarly wide absorption features, in the case of Ic-bl these features are found at much bluer (4000Å λ 8000Å) wavelengths. Thus, despite their bolometric similarities, superKNe are spectroscopically unique among SNe.
The peak photospheric temperatures of superKNe ∼ 1000 K are also similar to those required for solid condensation, suggesting the possibility of dust formation in the ejecta (e.g., Takami et al. 2014;Gall et al. 2017). Insofar as the optical/NIR opacity of ∼µm sized dust is roughly comparable to that of lanthanide-enriched ejecta, dust formation would not qualitatively impact the appearance of the transient. However, this does imply potential degeneracy between the photometric signatures of superKNe and other dust-enshrouded explosions unrelated to r-process production, including stellar mergers (e.g., Kasliwal et al. 2017). This degeneracy with dusty transients can generally be broken by the predicted broad spectral features of superKNe (v ej ∼ 0.1c).

10000
20000 30000 wavelength [Å] normalized flux + offset A t pk = 14.5 days B t pk = 29.5 days C t pk = 53.5 days D t pk = 58.5 days E t pk = 16.5 days Figure 8. The flux per unit wavelength at bolometric peak for each of the five models defined in Tab. 2. All spectra have broad absorption features consistent with a high-velocity outflow, and a low-temperature, pseudo-black body spectrum, consistent with a high-opacity composition. These spectra distinguish superKNe from other classes of SNe, which are much bluer, and from other dust-enshrouded explosions, in which broad absorption features are absent.

DISCOVERY PROSPECTS
In this section we explore the discovery prospects of superKNe with future optical/infrared transient surveys and via late-time infrared follow-up observations of energetic long GRBs. We then discuss how superKN emission could be enhanced by circumstellar interaction for collapsars embedded in AGN disks.

Volumetric Rates
We begin by estimating the volumetric rate of su-perKNe. One approach is to scale from the observed rates of ordinary collapsars. The local (redshift z 0) volumetric rate of classical long GRBs is ≈ 0.6 − 2 Gpc −3 yr −1 (Wanderman & Piran 2010), which for an assumed gamma-ray beaming fraction f b = 0.006 (Goldstein et al. 2016), corresponds to a total collapsar rate of ≈ 100 − 300 Gpc −3 yr −1 . Under the assumption that ordinary collapsars originate from stars of initial mass M ZAMS 40M , then the more massive stars M ZAMS 250M which generate helium core masses above the PI mass gap (M BH 130M ) will be less common by at least a factor ∼ (40/250) α−1 ∼ 0.1 − 0.3 for an initial-mass function (IMF) dN /dM ∝ M −α , where we consider values for the power-law index between α = 2.35 for a Salpeter IMF and a shallower value α ≈ 1.8 (Schneider et al. 2018). This optimistically assumes that (i) stars that massive exist (e.g., de Koter et al. 1997;Crowther et al. 2016), and that (ii) these can form helium cores such that M He M ZAMS , for instance because of rotational mixing (e.g., Maeder & Meynet 2000;Marchant et al. 2016; or continuous accretion of gas (e.g., Jermyn et al. 2021;Dittmann et al. 2021). Various processes act to remove mass from a massive star during its evolution, and generally the more massive the star, the larger its mass loss rate. Some of these mechanisms (e.g., continuumdriven stellar winds and eruptive mass loss phenomena, see also Renzo et al. 2020a) might occur even at low metallicity.
With the above estimate and caveats, we obtain an optimistic maximum local rate of superKNe from massive collapsars of ∼ 10 − 100 Gpc −3 yr −1 . On the other hand, the long GRB rate increases with redshift in rough proportion to the cosmic star-formation rate (SFR ∝ (1 + z) 3.4 for z 1; e.g., Yüksel et al. 2008) and hence the maximum rate of superKNe is larger at redshift z 1 by a factor ∼ 10 than at z 0, corresponding to a maximum superKN rate of ∼ 100 − 1000 Gpc −3 yr −1 at z 1.
The superKN rate question can be approached from another perspective: What is the minimum birth-rate of BHs in the PI mass gap to explain GW190521-like merger events (Sec. 5.3) via the massive collapsar channel? The rate of GW190521-like mergers at z 0 was estimated by LIGO/Virgo to be ∼ 0.5 − 1 Gpc −3 yr −1 (Abbott et al. 2020). This rate is smaller than the maximum superKN rate estimated above, consistent with only a small fraction of BHs formed through this channel ending up in tight binaries that merge due to gravitational waves at z ≈ 0.

Discovery with Optical/Infrared Surveys
We now evaluate the prospects for discovering su-perKNe with impending wide-field optical/infrared surveys.
First, we explore the expected observable rates within the Legacy Survey of Space and Time (LSST) conducted with the Vera Rubin Observatory. LSST is currently set to commence in early 2024 and will explore the southern sky in optical wavelengths to a 5σ stacked nightly visit depth of ∼ 24 mag. We inject the set of SEDONA light curves of models described in Tab. 2 into the publicly available LSST operations simulator, OpSim (Delgado & Reuter 2016). We use the most recent baseline scheduler (baseline v1.7) to calculate LSST pointings, limiting magnitudes, and expected sky noise across a full simulated 10 year survey in ugrizY bands. We additionally apply dust reddening following the dust maps of Schlegel et al. (1998). For each model, we inject a superKN randomly 300 times within the full LSST simulation (including both the wide-fast-deep survey and deep-drilling fields) at redshift bins of 0.01.
We find that superKNe discovered with LSST are confined to the local universe, with z < 0.1. Assuming that the superKN rate traces star-formation with a local rate of 10 Gpc −3 yr −1 , we expect LSST to discover ∼ 0 − 0.2 superKNe annually, resulting in up to 2 events over its 10-year nominal duration. We note that the larger the ejecta masses (i.e., Models b and e) the most likely the detection with LSST.
Given the expected red colors of the superKN emission (Fig. 8), we additionally explore the possibility of discovering superKN with the Nancy Grace Roman Space Telescope, expected to launch in the mid 2020s. Although not fully defined, Roman expects to conduct a ∼ 5 year, 10 deg 2 SN survey, primarily targeted at Type Ia SNe for cosmological distance measurements. We assume a survey cadence of 30 days and single-visit, stacked 5σ depth of 27th magnitude, corresponding to roughly an hour of integration time (in F158 band). We inject the same set of models using the Roman F062, F158 and F184 filters, corresponding to central wavelengths of 0.62, 1.58 and 1.83 µm, respectively. We assume observations are taken in each filter at the same epoch, and consider superKNe with three or more 3σ detections to be detectable. Assuming the Roman widefield survey footprint is chosen to minimize galactic dust, we do not account for any galactic reddening.
We find that Roman is most sensitive to models with the largest Lanthanide fractions. Assuming that the su-perKN rate traces star-formation with a local rate of 10 Gpc −3 yr −1 , we expect a 5 year Roman survey as described would find roughly 1-20 superKNe, most favoring the Lanthanide-rich Model B. These superKNe will be observable out to a redshift of z ∼ 0.9. We note that longer cadences significantly decrease the number of su-perKN detections possible with Roman, at least within the 3-detection discovery criterion we have adopted.

Energetic Long GRB Accompanied by SuperKNe
SuperKNe could also be detected following a subset of long GRBs. Figure 18 summarizes the GRB properties for our fiducial massive collapsar model 250.25. We find accretion timescales comparable to those of ordinary collapsars from lower mass progenitor stars (Appendix D). These mass gap collapsars are therefore candidates for contributing to the observed population of long GRBs, except that they may be a factor of ∼ 10 times more luminous and energetic than typical GRBs if the gammaray luminosity tracks the BH accreted mass. Furthermore, if the fraction of massive stars above the PI mass gap which form or evolve into collapsar progenitors is greater at lower metallicity, this could imprint itself on the redshift evolution of the long GRB luminosity function (for which there is claimed evidence; Petrosian et al. 2015;Sun et al. 2015;Pescalli et al. 2016).
In the local universe, long GRBs with supernovae are commonly accompanied with the luminous hyperenergetic Type Ic SNe with broad lines (Ic-BL; e.g., Woosley & Bloom 2006;Japelj et al. 2018;Modjaz et al. 2020). The superKN transients we predict from the birth of more massive BHs are of comparable or moderately lower peak luminosities than ordinary collapsar SNe (e.g., Cano 2016) but significantly redder (Figs. 7,8). Luminous optical SNe have been ruled out to accompany a few nominally long duration GRBs (Fynbo et al. 2006;Gehrels et al. 2006). One of these events, GRB 060614, was found to exhibit a red excess which Jin et al. (2015) interpreted as a kilonova. However, the luminosity and timescale of the excess could also be consistent with superKN emission from a massive collapsar of the type described here. We encourage future deep infrared follow-up observations of energetic long GRB with Roman or JWST on timescales of weeks to months after the burst to search for infrared superKN emission.

SuperKNe Embedded in AGN Disks
The optical emission from superKNe could be significantly enhanced by circumstellar interaction if they are embedded in a gas-rich environment. Graham et al. (2020) reported a candidate optical wavelength counterpart to GW190521 in the form of a flare from an active galactic nucleus (AGN). The flare reached a peak luminosity L pk ∼ 10 45 erg s −1 in excess of the nominal level of AGN emission and lasted a timescale t pk ∼ 50 days, over which it radiated a total energy of E rad ∼ 10 51 erg. Shibata et al. (2021) propose a scenario for GW190521 as a massive stellar core collapse generating a single BH and a massive accretion disk 10 − 50M rather than a binary BH merger. Although our results in Secs. 5.2 and 5.3 challenge this interpretation, our present work shows that a prediction of this scenario is a superKN counterpart with M ej ∼ 3 − 20M and v ej ∼ 0.1 c. Though the predicted peak timescale, t pk ∼ 50 days (Eq. 29), of the superKNe emission roughly agrees with that observed by Graham et al. (2020), the luminosity powered by radioactivity ∼ 10 42 erg s −1 (Fig. 7) is too small to explain the observations by several orders of magnitude.
This problem could be alleviated if the collapsing star is embedded in a dense gaseous AGN disk (e.g., Jermyn et al. 2021;Dittmann et al. 2021). If the density of the AGN disk at the star location is sufficiently high, ρ 10 −15 g cm −3 , runaway accretion of mass might help building up very massive and fast rotating helium cores. The mass accretion might be interrupted as the AGN turns off (on a few Myr timescale), and depending on the balance between mass loss processes and the previous accretion phase one might obtain a superKN progenitor. At its collapse, the shock-mediated collision between the superKN ejecta and the surrounding disk material could power a more luminous optical signal than from radioactive decay alone, akin to interactionpowered super-luminous SNe (e.g., Smith et al. 2007).
Given the large kinetic energy of the superKN ejecta, E kin ∼ M ej v 2 ej /2 ∼ 1 − 5 × 10 53 erg, the Graham et al. (2020) transient could be powered by tapping into only ∼ 1% of E kin by shock deceleration. Insofar as such luminous shocks are radiative and momentumconserving, the swept-up gaseous mass in the AGN disk required to dissipate E rad ∼ 10 51 erg is only M sw ∼ (E rad /E kin )M ej ∼ 0.1 − 1M . Treating the swept-up material as being approximately spherical and expanding at ∼ v ej , the optical diffusion time through M sw is (Eq. 29), where κ is now normalized to a value more appropriate to AGN disk material. Insofar that t pk,min is significantly shorter than the observed ∼ 50 d rise time of the Graham et al. (2020) counterpart, this implies the rise of the putative counterpart would instead need to be limited by photon diffusion through the unshocked external AGN disk material (e.g., Graham et al. 2020;Perna et al. 2021).
A bigger challenge for this scenario is the typically much closer source distance for GW190521 that would be predicted if this resulted from a self-gravitating collapsar disk instead of a binary BH merger (redshift z 0.05; Sec. 5.2), compared to that of the AGN identified by Graham et al. (2020) at redshift z = 0.438.

Luminous Slow Radio Transients
In addition to their prompt optical/IR signal, su-perKNe produce synchrotron radio emission as the ejecta decelerates by driving a shock into the circumburst medium (e.g., Nakar & Piran 2011;Metzger & Bower 2014). This emission can be particularly luminous because the kinetic energy of the superKN ejecta E kin ≈ 1 − 5 × 10 53 erg can be one to two orders of magnitude higher than those of ordinary collapsar SNe.
The radio transient rises on the timescale required for the ejecta to sweep up a mass comparable to their own, where n is the particle density of the external medium. The peak luminosity at a frequency ν = 1 GHz can be estimated as (e.g., Nakar & Piran 2011) where the fraction of the shock energy placed into relativistic electrons e and magnetic fields B are normalized to characteristic values, respectively, and we have assumed a power-law index p = 2.3 for the energy distribution of the shock accelerated electrons, dN/dE ∝ E −p . For characteristic circumstellar densities n ∼ 0.1 − 10 cm −3 the peak radio luminosity is comparable to that of rare energetic transients, such as those from binary neutron star mergers that generate stable magnetar remnants (e.g., Metzger & Bower 2014;Schroeder et al. 2020). However, the predicted timescale of the radio evolution of decades to centuries is much longer in the superKN case due to the large ejecta mass. This slow evolution makes it challenging to uniquely associate the radio source with a known GRB or gravitational wave event, or to even identify it as a transient in radio timedomain surveys (e.g., Metzger et al. 2015). We note that luminous radio point sources are in fact common in the types of dwarf galaxies which host collapsars (e.g., Eftekhari et al. 2020). Ofek (2017) place an upper limit on the local volumetric density of persistent radio sources in dwarf galaxies of luminosity 3 × 10 38 erg s −1 of N 5 × 10 4 Gpc −3 . Assuming the superKN radio emission remains above this luminosity threshold for a time t det ∼ 10t radio ∼ 10 3 yr, this constrains the local rate of superKNe to obey 10 − 100 Gpc −3 yr −1 , consistent with the estimates given in Sec. 4.1.

Gravitational Wave Emission
The accretion disks formed in superKN collapsars may become susceptible to gravitational instabilities if their disk mass approaches an order-unity fraction of the BH mass during the fallback evolution process (Sec. 2.2). As shown in Fig. 9, only progenitor cores with high angular momentum (small r b and/or high f K ) lead to fallback accretion that result in gravitational instabilities. Lowangular momentum cores instead form heavier BHs with relatively smaller disk masses.
The onset time of the instability of typically a few seconds (Fig. 9), representative of all superKN progenitor models investigated here, is determined by the progenitor structure, its rotation profile, and the free-fall timescale. Once triggered, subsequent fallback material collapsing onto the disk continues to excite these instabilities in the collapsar disk for a timescale of sec-onds to hundreds of seconds (Fig. 9), until viscous draining of the disk becomes fast compared to the free-fall timescale of the remaining outer layers of the progenitor star (roughly ∼ 10 s for our fiducial model in Fig. 2).
The onset of the instability manifests itself as the exponential growth of a non-axisymmetric one-arm (m = 1) density mode in the disk with growth time on the order of the orbital period of the disk, typically followed by exponential growth of an m = 2 mode (e.g., Kiuchi, K. et al. 2011;Shibata et al. 2021;Wessel et al. 2021). These non-axisymmetric density perturbations give rise to gravitational-wave emission with dominant frequency at the orbital and twice the orbital frequency, respectively (e.g., Wessel et al. 2021).
As long as further fallback keeps the disk in the instability regime defined by Eq. (25), we assume that the dominant gravitational-wave frequencies of these modes are determined by the evolving angular frequency Ω K,disk of the disk (Eq. 18) with radius r disk (t) (Sec. 2.2). Since r disk (t) monotonically increases with time as the black hole grows and material with larger specific angular momentum enters the disk, the gravitationalwave frequency decreases, sweeping down with a rate and amplitude that depends on the density and angular momentum structure of the progenitor star envelope. The gravitational-wave signal thus exhibits a "sad-trombone" pattern in the time-frequency spectrogram, as opposed to a "chirp" signal generally associated with gravitational waves from compact binary mergers. Examples of the frequency evolution of the disk for different mass models and for high and low specific angular momentum of the progenitor envelope are shown in Fig. 10. Over a large range of the parameter space and progenitor models explored here, superKN collapsars are strong emitters of quasi-monochromatic gravitational waves of duration ∼ 1 − 100 s with a decreasing frequency trend (between ∼ 0.1−40 Hz for the l = m = 2 and ∼ few × 10 −2 − 25 Hz for the l = 2, m = 1 mode) characteristic of their progenitor stellar structure (see Figs. 9 and 19 for a representative example). If detected, such gravitational-wave signals could reveal information about the rotation profiles of and angular momentum transport in evolved massive stars. The "sadtrombone" feature simultaneously followed by typically two dominant modes separated in frequency space by the instantaneous characteristic disk rotation frequency may prove useful in searching and detecting such sources with gravitational-wave detectors.
We calculate the gravitational wave strain of emitted gravitational waves as described in Appendix E. Figures  11-14 Figure 10. Disk frequency evolution (Eq. 18) for three progenitor models (250.25, 220.25, 200.25) with rotation parameters p = 4.5, r b = 1.5 × 10 9 cm, and overall small (fK = 0.3; top) or large (fK = 0.6; bottom) Keplerian angular momentum parameter. Plotted is twice the orbital angular frequency, which corresponds to the gravitational-wave frequency of the m = 2 density mode of the gravitationally unstable disk. The frequency evolution is largely controlled by fK, with all models reflecting the 'sad-trombone' nature of the gravitational-wave signal.
perKN events are expected to occur once every ∼ 3 years for our fiducial local superKN rate of 10 Gpc −3 yr −1 (Sec. 4.1). Figure 11 shows the time evolution of the plus and cross polarization strain calculated for the fiducial progenitor model (Fig. 2) assuming a face-on orientation of the collapsar disk (ι = 0). The maximum characteristic strain h c (typically h c ∼ 10 −24 − 10 −22 ) and the frequency range of the gravitational wave emission vary considerably across the {f K , r b } parameter space (Figs. 19 and 20,Appendix E).
SuperKN collapsars are multi-band gravitationalwave sources. Figures 12 and 14 compare the gravitational-wave signal in frequency space to the sensitivity of advanced LIGO (aLIGO), Cosmic Explorer   Amplitude spectral density (ASD) of gravitational-wave emission from the collapsar disk, shown for three progenitor models (250.25, 220.25, 200.25) and stellar rotation parameters p = 4.5, fK = 0.3, r b = 1.5 × 10 9 cm at an assumed source distance of 200 Mpc. The shaded region for each curve shows the unphysical frequency regime above the maximum disk frequency as plotted in Fig. 10, which is ignored in the SNR calculations. Shown for comparison are the measured or predicted noise curves for aLIGO, CE, ET, DECIGO, and BBO with sensitivity curve data from https://dcc.ligo.org/LIGO-T1500293/public and Yagi & Seto (2011). deciherz regime of DECIGO and BBO as the disk expands. The relative strain amplitude in these two different bands encodes information about the total mass and mass profile of the progenitors (Fig. 12). Lighter progenitors typically give rise to louder gravitational-wave signals over a narrower frequency band for the same rotation profile.
The overall magnitude of the amplitude spectral density is largely determined by the progenitor angular momentum as illustrated in Fig. 14. In the limit of high angular momentum (large value of the parameter f K ) for fixed r b , the instability and gravitational-wave emission are triggered earlier than for smaller values of f K (cf. Fig. 9). This is because matter deposition in the disk at early times is enhanced (rather than direct fallback onto the black hole). Under these conditions, the gravitational-wave signal is relatively weak due to the small disk and BH mass. Owing to enhanced viscosity and enhanced accretion during the instability epoch, disks that become unstable early on tend to stay relatively light; the gravitational-wave signal thus remains relatively weak throughout the fallback process. As a result, these signals tend to peak late and thus in the decihertz regime, which may only render them detectable there for Mpc distances. A non-detection in the highfrequency band may thus be indicative of the angular momentum budget of the progenitor star. In the other limit of low angular momentum (small value of the parameter f K and large r b ), the accretion disk may never become susceptible to the instability and gravitationalwave emission may be negligible (cf. Fig. 9). Hence, there exists an intermediate regime of progenitor angular momentum (intermediate values of f K ) in which the gravitational wave strain becomes maximal. For the given parameters of our fiducial progenitor model, this optimum is reached for f K ≈ 0.2, which is also reflected by the detection horizons (Figs. 13, 21).
We calculate a detection horizon for these events assuming an optimal matched filter and an SNR of 8 (see Appendix E for details). We find a detection horizon of ∼5 Mpc (aLIGO), ∼300 Mpc (ET), ∼250 Mpc (CE), and ∼425 Mpc (DECIGO) for our fiducial model with mass 250.25M , f K = 0.3, and r b = 1.5 × 10 9 cm. A parameter space study of the detection horizons is presented in Fig. 13, showing that third-generation detectors (ET, CE) as well as DECIGO are able to detect gravitational waves from superKN collapsars at distances of typically a few hundred Mpc up to a few Gpc. Detection horizons for aLIGO are typically limited to 100 Mpc (Fig. 21, Appendix E). BBO will be particularly sensitive to the lowest-frequency sources with low angular momentum in the progenitor 'core' (medium to large values of r b ) and typically reach several hundred Mpc to several Gpc. (Fig. 21, Appendix E).

GW190521
Our work has several potential implications for the gravitational wave event GW190521 (Abbott et al. 2020). Firstly, as already discussed, in the standard in- . Detection horizons of gravitational waves from our fiducial model shown in Fig. 2 with p = 4.5, fK = 0.3 and r b = 1.5 × 10 9 cm for Cosmic Explorer (top), the Einstein Telescope (center), and DECIGO (bottom), assuming optimal matched filtering and a signal-to-noise ratio of 8. For progenitors with medium to low rotation, these detectors may be able to detect gravitational waves from superKN collapsars at distances of typically a few hundred Mpc up to a few Gpc. These estimates are based on the corresponding physical frequency regime as indicated in Fig. 10. The sharp decrease in detection horizon for CE at log r b 9.4 is due to low-frequency emission below 10 Hz (below CE's sensitivity band  Figure 14. Same as Fig. 12 but for the progenitor model 250.25 with p = 4.5, r b = 1.5 × 10 9 cm and different values of the Keplerian fraction, fK = 0.2, 0.4 and 0.6. Relatively low to medium angular momentum models (here fK ≈ 0.2) generate a stronger signal in all detectors (with detection horizon ∼ 400 Mpc with aLIGO at an SNR of 8) compared to cases with higher angular momentum (large values of fK).
The fK = 0.6 model would only be detectable by BBO with a distance up to 100 Mpc at an SNR of 8 (cf. Fig. 21).
terpretation of GW190521 as a binary BH coalescence, mass loss associated with the birth of one or both of the constituent BHs can place them in the nominal PI mass gap "from above", even if they would have been above the PI gap if all of the star's mass were accreted at the time of core collapse (Fig. 4, top panel). To generate a BH with a mass consistent with the more massive member of GW190521 of ∼ 88M (Abbott et al. 2020) from a star with a helium core nominally above the gap, would require the ejection of 50M of ejecta (most of it r-process enriched; Sec. 2.2). In a direct sense, su-perKNe probe one channel for forming BHs in the PI mass gap.
Our scenario requires a fast rotating pre-collapse star and predicts that the magnitude of the BH spin would be nearly maximal (a BH,fin ∼ 1; Fig. 18, bottom panel). Although a low orbit-aligned spin χ eff 0.35 (90% confidence) was measured for GW190521, there is some evidence for a large spin component in the binary plane (Abbott et al. 2020). However, assuming that the progenitor stars can retain large rotation rates (see however Spruit 2002; and the progenitor of GW190521 formed from an isolated stellar binary through common envelope evolution (e.g., Belczynski et al. 2016), stable mass transfer (e.g., van den Heuvel et al. 2017;van Son et al. 2021), or via chemically homogeneous evolution driven by tidal interactions (e.g., Maeder & Meynet 2000;Marchant et al. 2016), one would expect the stellar angular momentum vector-and hence that of the BHs formed from the collapse-to be aligned with the orbital angular momentum .
In the case of rapidly rotating progenitors, we speculate that misaligned spins could arise from a kick imparted to the BH by mass loss in the disk winds. Our calculations in Sec. 5.2 indicate that the formed disks can become self-gravitating and hence will be subject to bar-mode like instabilities, generating non-axisymmetric spiral density waves. The latter could impart a nonaxisymmetric component to the wind mass loss, which would endow the BH with an effective kick. To significantly misalign the spins without breaking the binary, the natal kick must be comparable to the pre-collapse orbital velocity of the system, v kick ∼ 300 km s −1 (e.g., Kalogera 1996;Callister et al. 2021). Given the characteristic wind ejecta speed v ej ∼ 0.1 c, from momentum conservation an asymmetry in the disk mass-loss rate or velocity at the level of v kick /v ej ∼ 10 −2 would be sufficient to impart significant spin-orbital misalignment. Although self-gravitating instabilities result in non-axisymmetric disk density fluctuations at the level δρ/ρ 0.1, quantifying the extent to which these impart non-axisymmetric mass-loss will require additional GRMHD simulations of the disk outflows in the regime of massive, self-gravitating disks.
In an alternative approach, Shibata et al. (2021) interpret GW190521 as gravitational waves from a nonaxisymmetric instability similar to those discussed in Sec. 5.2 in a massive BH-disk system, thought to originate from the collapse of a massive star. In contrast to the BH-torus systems of fixed mass numerically evolved by Shibata et al. (2021), in an astrophysical setting mass is continuously fed to the disk at a rateṁ fb,disk due to fallback of the progenitor envelope (Sec. 2.2). If the accretion disk approaches the gravitationally unstable regime its mass evolution M disk (t) is dominated by the addition of fallback material throughṁ fb,disk . The associated timescale over which the fallback rate changes is τṁ fb,disk ∝ t α , where α ≈ 1 (Siegel et al. 2019 and Sec. 2.2). This is reflected by the fact that disks for the progenitor models considered here become gravitationally unstable on a timescale of a few seconds after BH formation, and this instability phase then lasts for a few seconds to tens of seconds (Fig. 9, Sec. 5.2). Generating a gravitational-wave signal of only a few cycles and duration t GW ∼ 0.1 s as required by GW190521 (Abbott et al. 2020) with comparatively negligible amplitude thereafter thus requires a fallback rate given by the total mass of the BH-disk system divided by the duration of the signal ofṁ fb,disk > (M • + M disk )/t GW ∼ 650 − 1000 M s −1 for the configurations considered by Shibata et al. (2021). While fallback rates of the order of up to ∼ 30 − 40 M s −1 may be reached realistically (cf., e.g., Fig. 2), fallback rates that are larger by one or two orders of magnitude seem implausible even with the most compact progenitor models possible (Sec. 2.1). These fallback rates also set limits on the compactness and thus on the gravitational-wave frequency of possible BH-disk systems-the frequency of gravitationalwave emission may be in tension with GW190521 as well. While a frequency around ∼ 60 Hz of GW190521 may not be impossible per se, our unstable BH-disk systems are typically not compact enough to reach such high frequencies even for an l = m = 2 mode and instead strongly prefer maximum frequencies below 20-30 Hz (Fig. 10). Another consequence of the weaker gravitational wave signals that we predict from massive collapsars are much closer detection horizons, which for Advanced LIGO at O3 sensitivity amount to 100 Mpc (Fig. 21, Appendix E); it is unlikely that a superKN transient from GW190521 would have gone undetected by existing wide-field optical surveys at such close distances.

Implications for Galactic r-Process Enrichment
Although the massive collapsars studied here are probably less common by a factor 10 − 30 than the bulk of ordinary collapsars (Sec. 4.1), they can in principle generate ∼ 10 times more r-process ejecta mass for similar progenitor angular momentum structure. The contribution of superKNe to total r-process production in the Universe could therefore be non-negligible relative to that of ordinary collapsars.
SuperKNe are probably too rare to explain the occurrence of individual r-process pollution events in small dwarf galaxies (e.g., Ji et al. 2016), but this does not exclude them from contributing to larger stellar systems. The total mass of r-process elements in the Milky Way is only ∼ 10 4 M , so given an r-process yield of 10M per superKNe, the number of contributing events must be 1000 and probably 100 (accounting for the dominant contribution likely coming from other channels such as lower mass collapsars and neutron star mergers).
Depending on the efficiency of gas mixing and retention in the environments of superKNe, subsequent generations of star formation could produce a modest number of extremely r-process-enriched stars. The dilution mass of the interstellar medium into which the superKN ejecta is mixed, can be estimated as (e.g., Macias & Ramirez-Ruiz 2019) M dil ≈ 2 × 10 7 M E kin 5 × 10 53 erg 0.97 where we have assumed 10 km s −1 for the sound speed of the interstellar medium. The value of M dil in su-perKNe is larger than in ordinary SNe (E kin ∼ 10 51 erg) or ordinary collapsars (E kin ∼ 10 52 erg) by a factor of 10 − 100.
The total amount of r-process material generated by superKNe is larger than ordinary collapsars by a similar factor ∼ 10, while the production ∼ 0.1 − 0.5M of 56 Ni and hence 56 Fe (Fig. 5) are similar to ordinary collapsars. If a superKN were to occur in otherwise pristine material at very low metallicity (perhaps an questionable idealization given that SNe tend to be spatially and temporally clustered), the next generation of stars which form from this material could possess a metallicity as low as [Fe/H] ∼ −5 and a Europium abundance as high as [Eu/Fe] ∼ 5, much higher than the current record holder (Reichert et al. 2021). This abundance combination would also contrast strongly with the chemical signatures of PI SNe (e.g., Woosley et al. 2002;Aoki et al. 2014), for which a large quantity of iron group elements but no r-process elements are produced.

CONCLUSIONS
We have explored the collapse of rotating very massive 130M helium stars and predicted their nucleosynthetic, electromagnetic, and gravitational waves signatures. Our conclusions can be summarized as follows.
• Building on Siegel, Barnes, & Metzger (2019), we present a semi-analytic model for the BH accretion disk in collapsars and its associated outflows which predict the quantity and composition of the disk wind ejecta, as well as the final BH mass and spin, given an assumed angular momentum structure of the progenitor star. The accretion regimes are calibrated based on the results of numerical GRMHD simulations and analytic scaling relations (Appendix A). Although the radial angular momentum structure of the progenitor star at collapse is uncertain theoretically, our approach allows us to cover a wide portion of the physically allowed parameter space. Applied to "ordinary" low-mass collapsars, the model predicts accretion luminosities and durations broadly consistent with long GRB observations.
• Our main application is to massive collapsars, originating from progenitor stars with final helium core masses ∼ 125 − 150M , which avoid pair instability SNe and nominally (in the case of zero mass ejection) would create BHs above the PI mass-gap. Analogous to lower-mass collapsars, as the fall-back accretion rate declines in time, the composition of the disk outflows systematically evolve from heavier to lighter elements (Fig. 3). Across a wide parameter space of progenitor rotational properties, we find total wind ejecta masses ∼ 10 − 50M , of which ∼ 10 − 60% is composed of r-process nuclei, including a sizable quantity of lanthanide elements associated with heavy rprocess production. The remaining ejecta is primarily unprocessed material (assumed to be 4 He in our models) and a modest quantity ∼ 0.1−1M of 56 Ni, formed from the brief hot, proton-rich phases of the disk evolution.
• The radioactive decay of r-process nuclei and 56 Ni in the ejecta of massive collapsars powers a months-long transient with a peak luminosity ∼ 10 42 erg s −1 (Fig. 7), which we refer to as a "su-perKN". The spectral energy distribution of su-perKNe near maximum light peaks at several microns due to the large opacity of the lanthanide elements (Fig. 8), similar to lanthanide-rich kilonovae from neutron star mergers. Although the bolometric light curves of superKNe are broadly similar to common types of core-collapse SNe, their combination of extremely red colors and high-velocity spectral features (v ej ∼ 0.1 c) should render superKNe distinguishable from other transient classes.
Our radiative transfer calculations have assumed a homogeneous ejecta structure; if the ejecta instead exhibits significant radial stratification, particularly a low lanthanide abundance in the highest velocity outermost layers, then the early light curve could be substantially brighter and bluer than our baseline predictions.
• Even for a progenitor stars well above the PI mass gap at collapse, the final BH remnant can populate the entire mass gap between ∼ 55 − 130M due to disk wind mass-loss (e.g., Fig. 4; Tab. 1). Su-perKNe therefore probe one mechanism for populating the PI mass gap "from above". The BHs formed through this channel are predicted to be rapidly spinning due to the large accretion of angular momentum, with final Kerr parameter a BH,fin ∼ 1. If the BH is formed in a binary, we speculate that its spin angular momentum axis could become misaligned with that of the binary angular momentum due to non-asymmetric mass-outflows associated with the gravitationallyunstable phases of the accretion (Sec. 5.2). Future numerical simulation work is necessary to explore this possibility quantitatively.
• One avenue to discover SuperKNe is via wide-field optical/infrared surveys. A 5-year survey with the Roman Space Telescope similar to that planned for Type Ia SNe, could potentially detect ∼ 1 − 20 superKNe out to redshift z 1, for an assumed z = 0 superKN rate of ∼ 10 Gpc −3 yr −1 . Su-perKNe could also be discovered by LSST, but the detection rate is lower because the predicted emission peaks at redder wavelengths than covered by the LSST bands. Measurements or limits on the occurrence rate of superKNe would constrain the birth rate of PI mass gap BHs via this channel (for comparison, the local rate of GW190521-like mergers is ∼ 1 Gpc −3 yr −1 ; Abbott et al. 2020). Su-perKNe may also be detectable following (particularly energetic) GRBs with JWST after the GRB afterglow has faded.
• The large kinetic energies of the superKN ejecta 10 53 erg results in a bright, long-lived synchrotron radio transient as the ejecta decelerates via shock interaction with the circumstellar medium (Sec. 5.1). However, the slow evolution of the radio emission for typical circumstellar densities will render these radio sources challenging to identify as radio transients (they may appear as luminous persistent sources in star-forming dwarf galaxies, for example; Eftekhari et al. 2020).
If superKNe occur inside gaseous AGN disks, shock interaction with the dense disk material could substantially enhance the optical luminosity of the event relative to that powered by radioactivity alone. This offers a speculative explanation for the claimed optical counterpart of GW190521 (Graham et al. 2020), provided it represents a gravitational wave burst from a core collapse event  instead of a black hole merger (however, see Sec. 5.3).
• The massive accretion disks from massive collapsars can become gravitationally unstable, generating gravitational wave emission as a result of non-axisymmetric density fluctuations. The predicted duration of the gravitational waves is several seconds or longer (calling into question the core-collapse origin for GW190521 proposed by Shibata et al. 2021, Sec. 5.3), while the frequency range overlaps the sensitivity window of ground-based (e.g., LIGO/CE/ET) and spacebased intermediate-frequency gravitational-wave detectors (e.g., DECIGO, BBO). Unlike the gravitational wave signal of compact binary mergers, which increase in frequency and amplitude with time ("chirp"), the gravitational wave signals of collapsar disks decreases in frequency as the disk radius grows ("sad-trombone"). Our simple estimates suggest that gravitational waves from massive collapsar disks are detectable by CE/ET/DECIGO to distances of up to several hundred Mpc (Figs. 12,13,14), interior to which the event rate could be as high as once every few years.
• SuperKNe are unlikely to contribute dominantly to the total production of r-process elements in the Universe, compared to neutron star mergers or ordinary low-mass collapsars, because the progenitors are disfavored by the initial mass function of stars. However, their extremely r-process-rich but iron-poor ejecta could in principle seed the creation of a small fraction of stars with abundance ratios more extreme than currently known metalpoor r-process-enhanced stars (e.g., [Eu/Fe] ∼ 5).
Scaling r to the radius of the innermost stable circular orbit, r ISCO ∝ M • , we have r ∝ M • and Ω ∝ M −1 • , thus givingṀ A derivation of this scaling in the general-relativistic context is given in De & Siegel (2020). The scalinġ M ign ∝ α 5/3 has been verified in 1D models (Chen & Beloborodov 2007) for BH masses M • ≈ 3M . However, the scaling should hold for the higher values of M • of interest in this paper because the assumptions which enter the above derivation (i.e., optically thin cooling, radiation pressure dominating over gas pressure at the advective transition) only strengthen with increasing BH mass. In particular, the ratio of radiation to gas pressure, increases with M • . Other forms of neutrino cooling, particularly arising from the annihilation of electronpositron pairs into neutrino-antineutrino pairs, will dominate over the pair-capture rates entering the radiative cooling rate at sufficiently high T 3 /ρ (e.g., Qian & Woosley 1996), changing the above scalings. However, the weak M • -dependence in Eq. (A7) suggests that this is not likely to occur for even the highest BH masses M • ∼ 100M of interest in this paper.

A.2. Neutrino Opaque and Trapping Thresholds
We now consider the accretion rateṀ ν,r−p , above which the inner disk is optically thick to neutrinos. The vertical optical depth through the disk obeys where σ ∝ T 2 is the energy-dependent absorption cross section with which electron neutrinos or antineutrinos are absorbed by free neutrons or protons, respectively. Using Eqs. (A3) and (A4) for T and ρ, and again taking H/r ∼ const., we find: Evaluating the marginally thick condition τ ν = 1 at r ≈ r ISCO , we obtainṀ ν,r−p ∝ αM

A.3. Neutrino Trapping Threshold
Finally, consider the "trapping" accretion rateṀ tr , above which the thermal energy released by the accretion flow is advected into the BH faster than it can be radiated through neutrinos. Equating the neutrino diffusion timescale out of the disk midplane t diff ∼ τ ν (H/c) with the inwards flow time ∼ r/v, we find that neutrinos are trapped for τ ν > (c/v)(r/H). Using v ∼ ν/r for the inflow velocity of a steady disk and Eq. (A8), and again taking H/r ∼ const., the trapping condition can be written Finally, recalling that νΣ ∝Ṁ for steady accretion, then using Eq. (A3) and evaluating Eq. (A11) at r = r ISCO , we findṀ B. PRODUCTION OF 56 NI IN DISK WINDS AND DISSOCIATION THRESHOLD Assuming the dominant seed particle formation process in the disk wind outflow is 4 He(2α, γ) 12 C, rather than the neutron-catalyzed reaction 4 He(αn, γ) 9 Be(α, n) 12 C in neutron-rich environments Y e 0.5 (Woosley & Hoffman 1992), the destruction of α-particles in disk winds proceeds as (Roberts et al. 2010) where Y α is the abundance of α-particles, λ 3α (T ) is the temperature-dependent triple-alpha rate coefficient, and the factor of 14 is due to assuming α-captures cease at 56 Ni. Furthermore, dτ = −(τ d /3T )dT , where τ d is the expansion time of the outflow around the point of αparticle formation. We take τ d ≈ 5 − 30 ms typical of expansion time scales of such accretion disk winds (Siegel & Metzger 2018;Siegel et al. 2019). Integration of Eq. (B13) yields the resulting abundance of seed particles (Roberts et al. 2010) Here, s f is the final entropy in k B per baryon, which we estimate by s f ≈ s m = s m,rad + s m,N , where is the entropy of radiation in the disk midplane and s m,N ≈ 13.3 + ln (T m /1MeV) 3/2 ρ m /10 7 g cm −3 (B17) is the entropy in the non-relativistic nucleons. The disk midplane density ρ m and temperature T m are calculated using an α-disk model as in Appendix A (see also Siegel et al. 2019). We translate Eq. (B14) into a mass fraction X56 Ni and X He ≈ 1 − X seed by assuming that seed particles are mostly within the iron peak (charge numbers 24 ≤ Z ≤ 28), using the nucleosynthesis results of Siegel et al. (2019). At late times during the fallback process (low accretion rates), the accretion disk must be hot and dense enough to dissociate α-particles as well as heavier nuclei of the infalling stellar material into individual nucleons, a requisite for the synthesis of 56 Ni in disk outflows. We estimate the time of this transition to a state in which dissociation into individual nucleons becomes suppressed by considering a fluid of neutrons, protons, and α-particles with a proton fraction Y e = 0.5 and determining the state in which half of the α-particles are dissociated. Assuming that nuclear statistical equilibrium holds (a good assumption as our disk midplane temperatures typically stay above 0.5 MeV), this can be estimated from the Saha equation (Shapiro & Teukolsky 1983) X 4 nuc = 1.57 × 10 4 X α ρ −3 m,10 T 9/2 m,10 exp − 32.81 T m,10 , where X nuc = 2X p = 2X n is the mass fraction of nucleons, and ρ m,10 and T m,10 denote the midplane density and temperature in 10 10 g cm −3 and 10 10 K, respectively. Expressing X nuc = 1 − X α and setting X α = 0.5, Eq. (B18) can be solved numerically for the transition time t diss upon inserting ρ m and T m as obtained from the numerical evolution of Eqs. (9)-(14). The timedependence of the functions ρ m and T m is determined by the evolution of M • , M disk , and r disk . At t > t diss we set Y seed = 0 in Eq. (B14), discarding any potential 56 Ni production after this time.

C. RESOLUTION STUDIES
We have performed numerical convergence tests to check convergence of nucleosynthesis results from our collapsar model (Sec. 2.3) and to determine optimal resolution for our numerical collapsar evolution calculations (Sec. 2.2). Figure 15 shows results of two convergence tests to determine optimal discretization for the polar coordinate and for time integration. The top panel illustrates that at our fiducial resolution in the polar angle of n θ = 1001 the relative error in computing the total mass of the star by numerical integration is dominated by the radial resolution of the progenitor model (relative error of 0.0018). The bottom panel indicates that for sufficiently high angular resolution (n θ ∼ 1001) and sufficiently large number of time steps of several hundred to 10 3 , the relative error in computing the stellar mass by numerical integration of Eqs. (9)-(14) is again dominated by the radial discretization of the progenitor model. The relative error of 0.002 in this case is consistent with the error budget obtained for the corresponding convergence test in the polar angle. This shows that our results are converged with roughly 10 3 grid points both in θ an in time, which we employ for all model runs.  Figure 16 presents results for model E20 of Heger et al. (2000), one representative case for ordinary collapsars below the PI BH mass gap. We vary the parameters of the adopted progenitor rotation profile (cf. Eq. (2)) within ranges motivated by the structure of the stellar evolution models (see Sec. 2.2), making nearly identical assumptions regarding the rotation profile as for massgap collapsars (Sec. 2.3.2). Our results are almost insensitive to the exact value of the power-law index p, which we thus fix to p = 4.5 for simplicity here. For model E20, we find ≈ 0.04 − 0.7 M of r-process material, including ≈ 0.01 − 0.04 M of heavy (A > 136) r-process material and ≈ 0.03 − 0.65 M of light (A < 136) r-process material, and ≈ 0.14 − 0.26 M of 56 Ni.
In comparison to Siegel et al. (2019), the updated model presented here tends to predict moderately less heavy r-process material, more light r-process material, and more 56 Ni. This is the result of i) a more detailed treatment of the disk accretion rate onto the BH, ii) an additional nucleosynthesis regime of light rprocess material only at high accretion rates >Ṁ ν,r−p (cf. Eq. (26)), and iii) a detailed evolution of the nucleosynthesis regimes throughout the accretion process as a result of BH growth. Overall, however, the mass ranges of all nucleosynthesis products found here broadly agree with Siegel et al. (2019). In particular, our refined analysis still predicts a sizable amount of lanthanide-bearing r-process ejecta of ≈ 0.04−0.7 M across various progenitor models of Heger et al. (2000). These results remain consistent with Miller et al. (2020), insofar that the mass accretion and generation of disk winds occur over a wide range of accretion rates, which drift through nucleosynthesis regimes characterized by varying degrees of neutrino irradiation (cf. Eq. (26)). Interestingly, our new models result in disk-wind 56 Ni yields approaching the values required to explain the light curves of observed GRB SNe (e.g., Cano et al. 2016), without a prompt shock-heated explosion (e.g., Barnes et al. 2018) or explosive nucleosynthesis at larger radii in the disk (Zenati et al. 2020).
Our collapsar model is also in good agreement with properties of observed GRBs. We check for consistency of our collapsar model with observed GRBs in terms of GRB timescales and energies. We assume that the accreted mass onto the BH is proportional to the radiated γ-ray energy, that is, L γ ∝ ηṁ acc c 2 , where L γ is the observed gamma-ray luminosity and η is an efficiency parameter.
A necessary requirement for collapsar accretion to explain observed GRBs is that the evolution time of the accretion rate be smaller or equal to the typical time required to generate a GRB in the engine frame, i.e., τṁ acc τ GRB . The timescale τṁ acc increases with time, typically expected as a power-law τṁ acc ∝ t α , where α 1. Following Siegel et al. (2019), we define Furthermore, let t GRB denote the time relative to the onset of disk accretion at which the equality τṁ acc = τ GRB is reached. Consistency then requires The GRB duration in the engine frame is determined by (Bromberg et al. 2012;Sobacchi et al. 2017) where τ γ is the observed duration of a GRB in the engine rest frame and τ b is the time required for the jet to drill through the stellar envelope. Assuming a typical value of τ b = 57 +13 −10 s (Sobacchi et al. 2017), and a typical observed GRB duration of τ γ = T 90 /(1 + z) = 9 s, with a characteristic T 90 27 s and redshift z 2 (Narayana Bhat et al. 2016), one finds τ GRB ≈ 66 s. The top panel of Fig. 17 shows that t GRB ∼ τ GRB essentially throughout the parameter space. Similar results are found for the other models of Heger et al. (2000). Hence, consistency with observed GRB durations according to Eq. (D20) holds.
Consistency with typical observed GRB energies requires that where E γ,iso ∼ 1 × 10 53 erg is the typical isotropicequivalent gamma-ray energy of observed GRBs and f b 0.006 is the beaming fraction (Goldstein et al. 2016), L jet is the luminosity of the accretion-powered jet, and M acc,GRB is the accreted mass onto the BH through the disk during the GRB timescale t GRB . With these values, the condition translates into M acc,GRB 2.5 × 10 −3 M , which we find is satisfied throughout the parameter space where the peak accretion rate reachesṁ acc > 10 −4 M s −1 (cf. Fig. 17 . We therefore find that our ordinary collapsar models are consistent both with typical GRB duration times and energies, including drill time. Figure 18 shows a parameter-space scan for model 250.25, a typical mass gap collapsar model. The GRB properties are in good agreement with observational constraints. While the GRB durations are typically similar to ordinary collapsars, the accreted mass during the GRB phase may be up to a factor ∼ 10 higher. One may thus speculate that these models give rise to GRBs that may be a factor ∼ 10 more luminous or energetic, if the gamma-ray luminosity tracks accreted mass.

E. GRAVITATIONAL-WAVE EMISSION
We calculate the gravitational-wave strain of emitted gravitational waves by approximating 'the lump' of the unstable disk (assumed to correspond to an overdensity of δρ/ρ 0.1; Shibata et al. 2021;Wessel et al. 2021) and the central BH as two orbiting point masses. The frequencies of gravitational-wave emission can thus be directly predicted from the evolution of the disk angular velocity according to the collapsar model in Sec. 2.2. The maximum and minimum gravitational- wave frequencies vary considerably across the {f K , r b } parameter space (see Fig. 19 for our fiducial model).
According to the quadrupole formula, assuming that the orbital radius only slowly changes with respect to the orbital frequency, the plus (h + ) and cross (h × ) polarizations of the gravitational waves at distance r and inclination ι of the disk with respect to the observer can be written as h + (t) = 4G rc 4 µr 2 disk Ω 2 K,disk 1 + cos 2 ι 2 cos[Φ(t)], (E24) h × (t) = 4G rc 4 µr 2 disk Ω 2 K,disk cos ι sin[Φ(t)], where Φ(t) ≡ t t0 2Ω K,disk (t ) dt , with t = t 0 denoting the start time of the gravitational instability. These expressions apply to the l = m = 2 mode, while Ω K,disk is replaced by 0.5Ω K,disk for the m = 1 mode. Furthermore, µ = M lump M • /(M lump + M • ) is the reduced mass of the lump-BH system, and we set M lump = f lump M disk with f lump = 0.2. Uncertainties in the value of f lump can be absorbed into uncertainties of the scale height of the disk and the threshold mass fraction f disk,thr ≡ M • /M disk at which gravitational instabilities set in (cf. Eq. (25)). We neglect corrections ∝μ to Eqs. (E24) and (E25) due to a time-dependent reduced mass µ as a result of accretion onto the black hole, assuming that M • changes only weakly over the course of gravitational-wave emission.
We characterize gravitational-wave emission in the frequency domain (positive frequencies f only) by computing the characteristic strain, defined as  Fig. 2 with p = 4.5, fK = 0.3 and r b = 1.5 × 10 9 cm for advanced LIGO at design sensitivity (top) and Big Bang Observer (bottom), assuming optimal matched filtering and a signal-to-noise ratio of 8. While for aLIGO the detection horizon is typically limited to 100 Mpc, BBO will be able to detect such sources up to typically several hundred Mpc to several Gpc, with particular sensitivity for progenitors with low-angular momentum 'cores' (medium to large values of r b ). Contours delineate final BH masses as in previous figures.
For an estimate of the horizon distance we assume that the detector receives a signal from a directly overhead source and hence the optimal strain response at the detector can be written as whereh + andh × are the Fourier transforms of the respective polarization strain amplitudes that we compute employing a Tukey window function limited to the physical frequencies between the maximum and minimum frequency expected from disk evolution (Fig. 19). We compare gravitational wave signals to detector sensitivity in terms of the amplitude spectral density √ S h (Moore et al. 2015), where S h denotes the power spectral density, and calculate the signal-to-noise ratio (SNR) using an optimal filter (Moore et al. 2015), Characteristic strains of superKN collapsars range from ∼ 10 −24 − 10 −22 depending on the rotation profile of the progenitor. A representative example is shown in Fig. 20. Detection horizons for advanced LIGO and BBO assuming SNR = 8 are shown in Fig. 21, while those for CE, ET, and DECIGO are presented in Fig. 13.