R-process viable outflows are suppressed in global alpha-viscosity models of collapsar disks

Collapsar disks have been proposed to be rich factories of heavy elements, but the major question of whether their outflows are neutron-rich, and could therefore represent significant sites of the rapid neutron-capture (r-) process, or dominated by iron-group elements remains unresolved. We present the first global models of collapsars that start from a stellar progenitor and self-consistently describe the evolution of the disk, its composition, and its outflows in response to the imploding stellar mantle, using energy-dependent M1 neutrino transport and an alpha-viscosity to approximate turbulent angular-momentum transport. We find that a neutron-rich, neutrino-dominated accretion flow (NDAF) is established only marginally--either for short times or relatively low viscosities--because the disk tends to disintegrate into an advective disk (ADAF) already at relatively high mass-accretion rates, launching powerful outflows but preventing it from developing a hot, dense, and therefore neutron-rich core. Viscous outflows disrupt the star within ~100s with explosion energies close to that of hypernovae. If viscosity is neglected, a stable NDAF with disk mass of about 1Msun is formed but is unable to release neutron-rich ejecta, while it produces a relatively mild explosion powered by a neutrino-driven wind blown off its surface. With ejecta electron fractions close to 0.5, all models presumably produce large amounts of Ni56. Our results suggest that collapsar models based on the alpha-viscosity are inefficient r-process sites and that genuinely magnetohydrodynamic effects may be required to generate neutron-rich outflows. A relatively weak effective viscosity generated by magnetohydrodynamic turbulence would improve the prospects for obtaining neutron-rich ejecta.


INTRODUCTION
The identification of the dominant site(s) of the rapid neutron-capture (or r-) process remains a major goal of nuclear astrophysics. Neutron-star (NS) mergers are now widely considered as the most robust r-process site (e.g. Arnould & Goriely 2020;Cowan et al. 2021), particularly after the first multi-messenger observation of a binary NS merger, GW170817 and its electromagnetic counterparts (e.g. Kasen et al. 2017;Watson et al. 2019). However, whether mergers dominate galactic r-process enrichment or other sites play comparably important roles remains poorly constrained. While current models of ordinary core-collapse supernovae face serious problems in providing r-process viable conditions (e.g. Hüdepohl et al. 2010), the situation for the rare class of magneto-rotational supernovae may be more optimistic. Neutron-rich outflows could be obtained in such events through the fast expansion of jets launched from a highly magnetized proto-NS, requiring however a very strong progenitor magnetic-field strength in order to enable the production of the heaviest r-process elements (e.g. Nishimura et al. 2015;Mösta et al. 2018;Reichert et al. 2021Reichert et al. , 2022 Another possibility, in the case of a black hole (BH) being formed (i.e. in the "collapsar" scenario) is that massive, subrelativistic winds get expelled from the hot and dense BHaccretion disk (MacFadyen & Woosley 1999;Pruet et al. 2003;Kohri et al. 2005;Surman et al. 2006;Chen & Beloborodov 2007;Metzger et al. 2008b;Nagataki et al. 2007;Siegel et al. 2019; Barnes & Metzger 2022), next to the ultrarelativistic jet that could be launched from this system and is assumed to power long gamma-ray bursts (e.g. Woosley 1993;Popham et al. 1999;Aloy et al. 2000;MacFadyen et al. 2001;Komissarov et al. 2009;Nagataki 2009;Harikae et al. 2009;Tchekhovskoy & Giannios 2015;Ito et al. 2015;Gottlieb et al. 2021).
After several studies based on 1D models and parametrized outflow trajectories (e.g. Pruet et al. 2003;Surman et al. 2006) reported rather pessimistic conditions for the activation of the r-process in the relevant regime of mass-accretion rates onto the BH,Ṁ BH , the collapsar-disk scenario has recently seen renewed interest after Siegel et al. (2019) pointed to the apparent similarity of these systems to remnant disks of NS mergers, for which theoretical models (e.g. Fernández & Metzger 2013;Just et al. 2015a;Siegel & Metzger 2018;Miller et al. 2019;Hayashi et al. 2021), tentatively supported by the observed kilonova of GW170817 (Kasen et al. 2017), suggest generically neutron-rich outflows. Using general relativistic magnetohydrodynamics (MHD) simulations, Siegel et al. (2019) assumed the ejecta produced from isolated BH disks with different initial disk masses to be representative of ejecta launched from a collapsar disk at different epochs of massinfall rates, leading them to conclude that collapsars may be even more prolific r-process sites than NS mergers. However, apart from simplifications adopted in the neutrino treatment arXiv:2205.14158v2 [astro-ph.HE] 4 Aug 2022  ), time (t NDAF ) and BH massaccretion rate (Ṁ BH (t NDAF )) of NDAF-to-ADAF transition when neutrino emission efficiency drops below 1 %, analytic estimate of NDAF-to-ADAF mass-accretion rate for isolated disks (Ṁ ign ; cf. Eq. (1)), BH mass at t disk , BH mass at t f , BH spin at t f , minimum electron fraction within unbound material at t f , mass and total (kinetic plus thermal plus gravitational) energy of gravitationally unbound material at t f , total jet energy injected into the system between t disk and t f estimated using Eq. (2). that may have led to an overestimated neutron richness (see Miller et al. 2019;Just et al. 2021, for detailed discussions of neutrino-transport effects in disks), the methodology of Siegel et al. (2019) leaves open the important question of whether a disk formed during stellar collapse produces a similar outflow signature as an isolated disk for comparable values oḟ M BH . In this study, we address exactly this question by conducting global simulations of collapsars, starting from a stellar progenitor model and self-consistently following the collapse, disk formation, as well as the evolution and disintegration of the disk under the influence of turbulent viscosity and neutrino-transport effects.
2. MODEL SETUP We adopt the strongly rotating progenitor model 16TI (Woosley & Heger 2006) having a ZAMS (collapse) mass of 16 (13.96) M , which has been used in numerous existing studies of collapsars. After following the collapse until core bounce (commencing at time t ≈ 300 ms; cf. dashed vertical line of Fig. 1(a)) we replace all material in the innermost region by a central BH (of mass M BH and dimensionless spin parameter A BH ), which we assume would have been formed several seconds later if we had retained the central NS; this collapsar scenario is supported by the self-consistent simulations of Obergaulinger & Aloy (2022). Since for this progenitor a disk is expected to form only once material from a mass coordinate of ∼ 3.8 M arrives at the innermost stable circular orbit (ISCO) (cf. Fig. 2 of Woosley & Heger 2006), we initially place, to save computational resources, the absorbing inner boundary at a radius of r in = 50 km. In order to follow the formation and evolution of the disk, we later (at about t ∼ 7 s) move r in to the arithmetic average of r ISCO and r BH (the event-horizon radius) of a Kerr BH with mass M BH and spin A BH . The parameters M BH and A BH are continuously updated using the mass and angular momentum crossing the boundary at r in , neglecting for simplicity the distinction between baryonic and gravitational mass (cf. bottom panel of Fig. 1 (a)).
All models are simulated in axisymmetry with spherical polar coordinates using the finite-volume code ALCAR-AENUS (Obergaulinger 2008;Just et al. 2015b). We assume Newtonian physics but take into account basic general relativistic effects in the gravitational potential, computed as the sum of the BH-potential by Artemova et al. (1996) and the TOVpotential for the stellar mantle (case A of Marek et al. 2006). The equation of state includes neutrons, protons, helium, and 56 Ni in nuclear statistical equilibrium, arbitrarily degenerate and arbitrarily relativistic electrons and positrons, and fully thermalized photons. The transport of ν e andν e neutrinos is treated by an energy-dependent M1 scheme (Just et al. 2015b), which takes into account the nucleonic β-processes as in Bruenn (1985) augmented by weak-magnetism corrections (Horowitz 2002). Turbulent angular-momentum transport is approximated by the α-viscosity model (Shakura & Sunyaev 1973), using the type 2 scheme of Just et al. (2015a). In order to study the impact of viscosity, we consider three viscous cases with viscosity parameters α = 0.01, 0.03, and 0.09 (with corresponding model names a01, a03, a09 and denoted here as low, medium, and high viscosity, respectively, motivated by comparisons with MHD disks in Fernández et al. 2019;Just et al. 2021), and a non-viscous case (α = 0, named a0); see Table 1. The radial domain extends from r in to r max (with r max = 10 11/12 cm for the non-viscous/viscous case) and is discretized in a logarithmic manner by ≈ 150 zones per decade, while 80 uniform zones are used in the angular direction from θ = 0 to π. Neutrino energies are sampled by 10 bins logarithmically distributed between 0 and 80 MeV.
The models are initialized by mapping the density, ρ, electron fraction, Y e , pressure, and angular velocity onto radial shells of the axisymmetric computational domain. We note that the outermost layers of progenitor model 16TI rotate super-Keplerian (Woosley & Heger 2006) and begin to expand right after starting the simulation, however, with significantly smaller speeds than encountered later in the disk-driven outflows.
Our model setup is similar to that of MacFadyen & Woosley (1999) who, however, did not consistently evolve Y e . In order to assess the validity of ignoring the proto-NS evolution, we compare in Fig. 1 (a) radial profiles of our model with those taken from a model that consistently included the proto-NS (model 16TI of Obergaulinger & Aloy 2022) at time t = 6.2 s, where we find good agreement. 3. RESULTS 3.1. Disk formation Figure 2 provides additional global properties as functions of time, while Figs. 3 and 4 show, for models a03 and a0, snapshots of the spatial distribution at selected times. Initially, material with specific angular momentum, l, lower than the corresponding Keplerian value at the ISCO, l ISCO,Kep , falls into the BH in a quasi-spherical manner. At t ∼ 13 s (corresponding to M BH ∼ 3.8 M ) gas with l ≈ l ISCO,Kep arrives at the ISCO, where it is prevented from falling into the BH by the centrifugal barrier and starts to assemble a rotationally supported disk. All models, except the one with the strongest viscosity (a09), show qualitatively the same features at the time of disk formation. Fluid elements falling in from equatorial directions accumulate at the radial outer edge of the disk, enhancing its mass (cf. Fig. 2 (a)) and leading to an immediate increase of the density ( Fig. 1 (a)), temperature ( Fig. 2 (d)), and neutrino luminosity ( Fig. 2 (b)). Since we are mainly interested in the hottest and densest regions of the flow, we define the disk here as all material within r = 500 km and with an angular momentum exceeding 50 % of the local Keplerian value. The radial deceleration due to circularization in combination with the enhanced densities and temperatures (boosting the neutrino-emission rates and leading to dissociation of nuclei into free nucleons) leaves material enough time, at least right after disk formation, to a) adapt its local Y e to the equilibrium value Y eq e dictated by the density, temperature and the local neutrino abundances (e.g. Beloborodov 2003;Arcones et al. 2010;Fujibayashi et al. 2020;Just et al. 2021), and b) radiate away large parts of the thermal energy produced in the disk. The system has now formed a hyperaccretion flow, or neutrino-dominated accretion flow (NDAF; e.g. Popham et al. 1999;Kohri et al. 2005). As is characteristic for NDAFs, the electron degeneracy (Fig. 2 (e)) close to the BH is moderately high, η e > ∼ 1, resulting in a low value of Y e ≈ Y eq e < ∼ 0.2.

Evolution of viscous models
The subsequent evolution depends sensitively on the treatment of angular-momentum transport. We first consider the viscous models. As soon as the system enters a state of nearly Keplerian, differential rotation, viscous effects become important. In the initial NDAF state, neutrino cooling is efficient enough to approximately balance thermal heating. A neces-sary condition for this balance to hold is that the neutrinocooling timescales remain shorter than or equal to the viscous heating timescales. Due to the steep temperature dependence of the former (roughly ∝ T −6 ), this requirement essentially translates into a condition on the disk temperature. Once the temperature drops below some critical value close to T ∼ 1 MeV (roughly corresponding to emission timescales of O(10 s)), viscous heating starts to dominate neutrino cooling, which triggers viscous expansion and leads to a further reduction of the temperature and the neutrino emission efficiency ( Fig. 2 (b)), resulting in the freeze-out of Y e and causing the disk to undergo a transition into a radiatively inefficient disk, a so-called advection-dominated accretion flow (ADAF; Narayan & Yi 1994). In the fiducial model with α = 0.03 this happens at about t NDAF ≈ 14.8 s (i.e. about 2 s after disk formation), whereas the low-viscosity model exhibits this transition at about t ≈ 30.6 s. In the high-viscosity model (α = 0.09) the viscosity is so strong that the disk basically never enters the NDAF phase but evolves as an ADAF right away.
The ADAF state exhibits characteristic differences to the NDAF state, in agreement with previous investigations of neutrino-cooled disks (e.g. Metzger et al. 2008a; 2)). Missing data points for disk-related quantities indicate that no material fulfills the disk-definition criteria at these times. 2021): Once an ADAF forms, the viscously disintegrating disk is characterized by markedly lower densities and temperatures, non-degenerate electrons, and (partially) recombined nuclei. These properties are crucial to the prospects for rprocess nucleosynthesis, because they imply a high equilibrium value, Y eq e ≈ 0.5 (e.g. Arcones et al. 2010). In addition, due to inefficient cooling the ADAF state is subject to strong convective activity, which is why-at least in viscous but not necessarily in MHD models (Fernández et al. 2019)outflows are launched more efficiently during the ADAF than during the NDAF phase. Finally, the NDAF-to-ADAF transition causes a steepening ofṀ BH from roughly a t −3/2 behavior to approximately t −3 (Fig. 1 (a)).
As soon as the disk has formed, an accretion shock emerges as a result of the sudden deceleration of supersonically infalling stellar material by the disk environment (Fig. 3 (c)). The shock surface expands quickly (Fig. 2 (g)), powered by the energy input from the, mainly viscously driven, disk ejecta. This process ultimately leads to the unbinding of almost the entire remaining stellar mantle (Fig. 2 (h) and Table 1). The ejecta are mainly equatorial up to the point of breakout from the star (Fig. 3 (e)) and subsequently expand in lateral directions to reach an almost spherical shape at the final simulation time of t f = 100 s. None of the expelled material has a particularly low Y e , the minimum value (achieved in the model with the lowest viscosity) being 0.453. The total energy of the viscosity-driven explosions lies within 5 − 10 × 10 51 erg, i.e. significantly higher than for ordinary core-collapse supernovae and close to that of hypernovae (see, e.g., Nagataki 2018, for a review). We note that a diskwind driven hypernova scenario was first suggested by Mac-Fadyen & Woosley (1999), but its oblate/spherical explosion geometry makes it difficult to reconcile with observations of gamma-ray bursts associated to hypernovae (Mazzali et al. 2001), which favor a prolate and therefore more likely jetdriven supernova (MacFadyen & Woosley 1999;MacFadyen et al. 2001).

Evolution of the non-viscous model
In the absence of viscous effects the disk only evolves through the effects of neutrino cooling and heating, allowing it to remain in the NDAF state for the entire simulation time of t f = 200 s, with correspondingly high densities ( Fig. 1 (a)), temperatures, and electron degeneracies, as well as low val- ues of Y e within 0.05 − 0.2 ( Fig. 2 (d)-(f)). Like in the viscous models, an accretion shock forms around the disk and keeps expanding, however, with significantly slower speeds because it is not fueled by viscous outflows. Correspondingly, the mass and energy of unbound ejecta is considerably lower than in the viscous models (though still increasing at t = t f ). The energy source powering the explosion in this model appears to be neutrino heating, which is most intense close to the BH along the surface of the disk. Exactly in this region, namely right between the main body of the disk and the polar matter inflow, a narrow outflow is launched that expands within polar angles of θ ∼ 10 − 20 • (Fig. 4 (c)). The electron fraction of this neutrino-driven wind is close to Y e = 0.5. Support for a neutrino-driven mechanism comes from the fact that the time-integrated amount of net neutrino-heating energy (taking into account only regions where neutrino heating dominates neutrino cooling) is significantly higher than in the viscous models, high enough to explain the observed explosion energy (Fig. 2 (h)), and the neutrino efficiency is the largest among all models ( Fig. 2 (b)).

DISCUSSION
In contrast to α-viscosity models of NS-merger remnant disks (Fernández & Metzger 2013;Just et al. 2015a), none of our models produces outflows with significant amounts of neutron-rich (Y e < 0.25) matter, even though they operate at comparable mass-accretion rates, 10 −3 < ∼Ṁ BH /(M s −1 ) < ∼ 1. Instead of r-process elements, the ejecta in our models, with Y e very close to 0.5 and (at least at early times) temperatures in excess of 5 GK (cf. temperature contours in Figs. 3 and 4), are likely to produce a significant amount of 56 Ni (as originally suggested by MacFadyen & Woosley 1999). One characteristic difference to previously considered models of isolated disks (i.e. disks not interacting with an imploding stellar mantle) seems to be that the NDAF state is less stable and becomes disrupted at higher mass-accretion rates, as is suggested byṀ BH (t NDAF ) (cf. Table 1) being about one order of magnitude higher than the corresponding values predicted for isolated disks (e.g. Chen & Beloborodov 2007;De & Siegel 2021;Just et al. 2021): (where we assumed A BH ≈ 0.8). An additional difference to merger disks appears to be the low efficiency by which the modeled collapsar disks release the remaining neutron-rich material right around the time of freeze-out, t NDAF . While Fig. 4.-Same as Fig. 3, but for the non-viscous model a0, with velocity arrows whose length is capped at 10 8 cm s −1 , and at the characteristic times right before (a) and after disk formation (b), during the shock expansion through the stellar mantle (c), and shortly after shock breakout (d).
the reasons for the aforementioned differences need to be explored in future, more detailed studies, we can already conclude that they must be connected to the stellar environment into which collapsar disks are embedded, because most of the remaining properties (M BH ,Ṁ BH , A BH , and m disk ) are very similar in NS-merger disks. A critical role may be played by the radial structure of the disk: While a merger disk is born close to the ISCO and evolves as a single, expanding ring of matter, a collapsar disk is continuously replenished with matter at the circularization radius corresponding to the specific angular momentum of each infalling mass shell. Hence, even though the disk mass is similar, the radial distribution of matter, and therefore its composition, may be markedly different in both cases.
Our study contains two important shortcomings (apart from the approximate treatment of gravity): First, the α-viscosity model employed here captures effects related to the magnetorotational instability (MRI) only approximately. The ability to drive neutron-rich winds during the NDAF phase is expected to be higher in MHD-disks than in α-disks (Siegel & Metzger 2018;Fernández et al. 2019). Moreover, a full MHD treatment is necessary to reliably determine the time of the NDAF-to-ADAF transition, i.e. to assess which of our choices for α is most realistic. Second, we ignore any impact of the jet that is likely to be formed by the Blandford-Znajek process (Blandford & Znajek 1977). A jet could possibly carry away low-Y e material from the inner disk during the NDAF phase (Nakamura et al. 2015) and it would probably make the explosion geometry more prolate. A rough estimate of the jet luminosity (MacFadyen et al. 2001), is provided in Fig. 2 (i), where we assume a magnetic-to-thermal pressure ratio of 10 % at the equatorial r ISCO , as well as flux conservation between r ISCO and r BH , to estimate the magnetic-field strength, B, at r BH . For α > 0.01, the (crudely) estimated total energy injected by the jet between t disk and t f (cf. Table 1) is significantly smaller than that of the viscous ejecta, and only slightly above the threshold found in Aloy et al. (2018), which reinforces the consistency of not modeling explicitly the jet generation for the most viscous models. Neglecting the jet feedback in low-viscosity models is not so well justified if Eq.
(2) is a good proxy for the jet luminosity. Despite the aforementioned shortcomings, our study suggests that the self-consistent treatment of the imploding stellar mantle is an important ingredient for obtaining reliable predictions of Y e in collapsar disks and their outflows. Our results apply to a single potentially collapsar-forming model. Other models, where the accretion disk forms earlier and at initially higher accretion rates, could exhibit an NDAF phase for a longer time and possibly produce more neutron-rich outflows. Winds from disks at very high accretion rates are, however, expected to be subject to intense neutrino irradation, which tends to drive Y e towards 0.5 (Miller et al. 2019;Just et al. 2021). Future studies need to develop a more systematic understanding of the NDAF-to-ADAF transition and explore the sensitivity to the initial rotation profile and magnetic-field distribution in the stellar progenitor in order to better constrain the composition of material released by collapsar disks.
We thank the anonymous referee for constructive remarks. OJ acknowledges support by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement No. 759253. OJ and MO acknowledge support by the DFG -Project-ID 279384907 -SFB 1245. The research leading to these results has received funding from the European Union's Horizon 2020 Programme under the AHEAD2020 project (grant agreement No. 871158). MAA and MO have been supported by the Spanish Ministry of Science, Education and Universities (PGC2018-095984-B-I00 and PID2021-127495NB-I00) and the Valencian Community (PROME-TEU/2019/071). MO acknowledges support from the European Research Council under grant EUROPIUM-667912, as well as from the Spanish Ministry of Science via the Ramón y Cajal programme (RYC2018-024938-I). SN has been supported by JSPS KAKENHI Grant Number JP19H00693 and RIKEN pioneering project "Evolution of Matter in the Universe (r-EMU)". OJ is grateful for computational support by the HOKUSAI computing facility at RIKEN and VIRGO cluster at GSI.