Viscous evolution of a massive disk surrounding stellar-mass black holes in full general relativity

Long-term viscous neutrino-radiation hydrodynamics simulations in full general relativity are performed for a massive disk surrounding spinning stellar-mass black holes with mass $M_{\rm BH}=4$, $6$, and $10M_\odot$ and initial dimensionless spin $\chi \approx 0.8$. The initial disk is chosen to have mass $M_{\rm disk}\approx 0.1$ or $3M_\odot$ as plausible models of the remnants for the merger of black hole-neutron star binaries or the stellar core collapse from a rapidly rotating progenitor, respectively. For $M_{\rm disk} \approx 0.1M_\odot$ with the outer disk edge initially located at $r_{\rm out} \sim 200$ km, we find that $15$%-$20$% of $M_{\rm disk}$ is ejected and the average electron fraction of the ejecta is $\langle Y_e \rangle = 0.30$-$0.35$ as found in the previous study. For $M_{\rm disk} \approx 3M_\odot$, we find that $\approx 10$%-$20$% of $M_{\rm disk}$ is ejected for $r_{\rm out}\approx 200$-$1000$ km. In addition, $\langle Y_e \rangle$ of the ejecta can be enhanced to be $\gtrsim 0.4$ because the electron fraction is increased significantly during the long-term viscous expansion of the disk with high neutrino luminosity until the mass ejection sets in. Our results suggest that not heavy $r$-process elements but light trans-iron elements would be synthesized in the matter ejected from a massive torus surrounding stellar-mass black holes. We also find that the outcomes of the viscous evolution for the high-mass disk case is composed of a rapidly spinning black hole surrounded by a torus with a narrow funnel, which appears to be suitable for generating gamma-ray bursts.


I. INTRODUCTION
A stellar-mass spinning black hole surrounded by a massive disk is believed to be a frequent remnant for the merger of neutron star binaries (binary neutron stars and black hole-neutron star binaries), and for the stellar core collapse of massive and rapidly rotating progenitors. Such remnants have been speculated to be the central engines of gamma-ray bursts [1][2][3][4] and kilonovae [5][6][7][8]. This fact motivates the community to explore in detail the formation process of the black hole and disk surrounding it (see, e.g., Refs. [9,10] for a review of numerical simulations for the formation of black holes surrounded by disks) and subsequent long-term evolution of the disk by magnetohydrodynamics/viscoushydrodynamics processes with various sophisticated levels (e.g., Refs. [11][12][13][14][15][16][17][18][19][20][21][22][23][24]). In our previous paper [23], we performed a long-term viscous radiation hydrodynamics simulation for the system of a spinning black hole with low mass M BH = 3M and of a disk with mass M disk = 0.03-0.5M in full general relativity, paying particular attention to the mass ejection in the post-merger evolution of neutron-star binaries.
In this paper, we extend our previous work [23] and explore the dependence of the viscous evolution process of the disk and resulting mass ejection on the black-hole mass, M BH , and disk mass, M disk . We consider the systems of M BH = 4, 6, 10M and of M disk ≈ 0.1 and 3M with the initial dimensionless spin of the black hole χ ≈ 0.8. These are plausible remnants for the merger of black hole-neutron star binaries [25][26][27][28][29] or the stellar core collapse from a rapidly rotating progenitor to a black hole [30,31]. The purpose of this study is to explore the quantitative dependence of the disk evolution and subsequent mass ejection as well as nucleosynthesis on the mass of the black hole and disk. In particular, we pay attention to the evolution of the high-mass disk model in this paper.
The system of a spinning black hole and a high-mass disk can be a plausible model as remnants for the collapse of rapidly rotating progenitor stars to a black hole. The evolution of high-mass disks surrounding a black hole with relevant microphysics (such as realistic equation of state and neutrino cooling) has not been studied selfconsistently in previous work because not only general relativistic gravity for the black hole but also the selfgravity of the disk have to be taken into account. In the present study, we employ the framework of full general relativity, and thus, we can explore the high-mass disk model with no assumption (except for the assumption of axial symmetry). We will show that in the presence of a high-mass disk, the neutrino luminosity is enhanced, and as a result, the viscous evolution timescale of the disk is increased significantly. This effect modifies the property of the ejected matter and resulting elements produced in the nucleosynthesis. In addition, we show that the out-come of the evolution of the system is a rapidly spinning black hole with a geometrically-thick torus and a narrow funnel. This system appears to be suitable for generating gamma-ray bursts, i.e., for driving a collimated jet toward the direction of the rotation axis.
The paper is organized as follows. In Sec. II, we briefly summarize the basic equations employed in the present numerical simulation and initial conditions employed in this work. Section III presents numerical results for the simulations. We present the evolution process of the disk and black hole, the properties of the ejecta and nucleosynthesis in the matter ejected from the disks, and the final outcome of the disk evolution. Section IV is devoted to a summary. Throughout this paper, G, c, and k denote the gravitational constant, speed of light, and Boltzmann's constant, respectively.

II. SUMMARY FOR THE METHOD OF NUMERICAL COMPUTATIONS
We perform viscous neutrino-radiation hydrodynamics simulations for systems of a stellar-mass black hole and a massive disk in the framework of full general relativity with the assumption of axial symmetry using the same method and implementation as in Ref. [23]: We numerically solve Einstein's equation, the viscoushydrodynamics equations, the evolution equation for the viscous tensor, the evolution equations for the lepton fractions including the electron fraction, and neutrinoradiation transfer equations. Einstein's equation is solved using the original version of the Baumgarte-Shapiro-Shibata-Nakamura formalism [32] together with the puncture formulation [33], Z4c constraint propagation prescription [34], and 5th-order Kreiss-Oliger dissipation. The axial symmetry for the geometric variables is imposed using a cartoon method with the 4th-order accuracy in space [35,36]. The viscous-hydrodynamics equations and evolution equations for the viscous tensor are solved by the method described in Ref. [37]. For evolving the lepton fractions, we take into account electron and positron capture, electron-positron pair annihilation, nucleon-nucleon bremsstrahlung, and plasmon decay [38,39]. The quantities of black holes (mass and spin) are determined from their area and circumferential radii of apparent horizons [40], assuming that these quantities are written as functions of the mass and spin in the same formulation as in the vacuum case.
As in our previous work [23], we employ a tabulated equation of state based on the DD2 equation of state [41] for a relatively high-density part and the Timmes equation of state for the low-density part [42]. In this tabulated equation of state, thermodynamics quantities such as ε, P , and h are written as functions of ρ, Y e , and T where ε, P , h(= c 2 + ε + P/ρ), ρ, Y e and T are the specific internal energy, pressure, specific enthalpy, restmass density, electron fraction, and matter temperature, respectively. We choose the lowest rest-mass density to be 0.1 g/cm 3 in the table, and the atmosphere density for ρ * := ρu t √ −g in the hydrodynamics simulation is chosen to be 10 g/cm 3 in the central region which is smoothly decreased to 1 g/cm 3 in the outer region. Here u µ and g denote the four velocity and the determinant of the spacetime metric, respectively.
For viscous hydrodynamics, we need to input the viscous coefficient ν. Following our previous work [23], we set ν = α ν hc s H/c 2 where α ν is the dimensionless viscous coefficient (the so-called alpha parameter [43]), c s is the sound velocity, and H is a scale height. We basically employ α ν = 0.05 as in our previous work [23] supposing that the origin of the effective viscosity is the turbulence induced hypothetically by the magneto-rotational instability [44][45][46]. For comparison, we also employ α ν = 0.10 and 0.15 for a model with M BH = 10M and M disk ≈ 3M . We set H to a constant value that is approximately equal to the radius at the innermost stable circular orbit around the Kerr black hole of χ = 0.8, i.e., H = 30 km(M BH /10M ) ≈ 2GM BH /c 2 . This setting is also the same as in Ref. [23]. Note that H is a constant determined by the initial state of the black hole.
The viscous timescale (for heating and angular momentum transport) is written approximately as (2.1) where R denotes the cylindrical radius in the disk and the reference values are chosen for M BH = 10M . As we show in Sec. III, the evolution timescale for our choice of the viscous coefficient is indeed of the order of ∼ 1 s, which is much longer than the dynamical timescale (i.e., rotational period) of the disk approximately written as Thus, in the viscous hydrodynamics of the disks, the evolution should proceed in a quasi-steady manner.
Axisymmetric equilibrium states for black hole-disk systems are prepared as the initial condition [23] using the method shown in Ref. [47]. As in the previous study, we determine the angular velocity from the relation where j = c −1 hu ϕ is the specific angular momentum. Ω is the angular velocity defined by u ϕ /u t and n is a constant that determines the profile of the angular velocity (see Ref. [48] for more careful choice for deriving disks with "Keplerian profile"). In this paper, we fiducially choose n = 1/7 to align the value with that chosen in Ref. [23]. For two models for which the disk has a geometrically thin and less compact structure, we employ n = 1/5 and 1/4 (cf. Table I). These models are used to explore the dependence of the results on the initial disk compactness (cf. Secs. III C and III E). For obtaining the initial condition, we assume a relation between ρ and Y e in the same form of ρ(Y e ) as in Ref. [23]. For this model, the value of Y e is 0.07 in a high density region and approaches 0.5 for low density for which the effect of the electron degeneracy is weak. In addition, we assume that the specific entropy, s, is initially constant (in order to obtain the first integral of the Euler equation easily). We always choose s = 6k in this paper for simplicity.
The initial black-hole mass is chosen to be M BH,0 = 4, 6, and 10M with the disk mass being M disk ≈ 0.1 and 3M . For most of the cases, we set the inner edge of the disk to be r in = 2GM BH,0 /c 2 . For this case, the outer edge of the disk is located at r out ∼ 200-250 km. For two models we choose less compact disks for which r out ∼ 500 and 1000 km, M BH,0 ≈ 10M , and M disk ≈ 3M with s/k = 6, n = 1/5 or 1/4 and r in = 3.5GM BH,0 /c 2 (see Table I).
In our formalism for obtaining the equilibrium state composed of a black hole and a disk [47], we initially prepare a black-hole geometry as a seed. In this work, we set the dimensionless spin of the seed black hole to be χ = 0.8. This implies that only in the limit of M disk → 0, χ becomes 0.8. As already mentioned, the black-hole spin is measured from the area, A AH , and circumferential radii around the equatorial and meridian planes, c e and c p , for the black-hole horizon [47], assuming that A AH , c e , and c p are written as functions of the mass, M BH , and dimensionless spin, χ, of Kerr black holes as in the vacuum black-hole case. In the presence of matter outside the black hole, its geometry is modified, and for the highmass disk case with M disk ≈ 3M , the dimensionless spin becomes slightly smaller than 0.8 (see Table I).
Because the mass ratio, M disk /M BH,0 , of the initial conditions employed for the high-mass disk case is very large (3/10-3/4), the self-gravity of the disk can play an important role for its evolution. It is well-known that the self-gravitating disk is subject to non-axisymmetric deformation and fragmentation (e.g., Refs. [49][50][51][52][53]). In our axisymmetric simulation, this effect cannot be taken into account. For the high-mass disk, the so-called m = 1 mode is often non-linearly excited, and as a result, a significant fraction of the initial disk mass is likely to be swallowed by the black hole due to the angular momentum transport in the dynamical timescale of the system. The remnant after the non-axisymmetric instability proceeds will be a black hole surrounded by a disk that satisfies M BH M disk . As we show in Sec. III, in the viscous hydrodynamics, a significant fraction of the disk is also swallowed by the black hole after the viscous angular momentum transport sets in, and the resulting remnant is a black hole and a disk with M BH M disk as well. Thus, if we focus on the process only after a large fraction of the disk is swallowed by the black hole, the viscous hydrodynamics approach is acceptable (that is, the effect of the non-axisymmetric deformation is substituted by the viscous effect). However, we should keep in mind that we need a non-axisymmetric simulation to clarify the evolution of the high-mass disk system quantitatively. For the non-axisymmetric simulation, high-amplitude gravitational waves are likely to be emitted during the nonaxisymemtric instability proceeds. Thus, exploring the feature of gravitational waves will be also an interesting topic in the non-axisymmetric simulations.

III. NUMERICAL RESULTS
A. Models, setting, and diagnostics Numerical computations are performed for the black hole-disk systems summarized in the previous section (see also Table I). For most models in this paper, the simulations are performed taking into account the neutrino absorption/irradiation effect. For one model, M10H05n, we switch off this neutrino effect to examine the importance of the neutrino absorption/irradiation effect for the high-mass disk case. In the present work, we do not incorporate a heating effect by the neutrino pair annihilation [38] for simplicity (thus the mass and kinetic energy of the ejecta might be underestimated).
The viscous tensor is set to be zero initially; we prepare an equilibrium state assuming the ideal fluid with no viscosity. Because the viscous tensor changes to a non-zero profile during the very early evolution, the disk profile is modified for all the cases. For the high-mass disk models, this (artificial) modification is significant in particular for the relatively low-mass black-hole models. Specifically, the disk radially oscillates with a high amplitude for t 0.2 s for such cases, but subsequently, it relaxes to a quasi-steady state, which is evolved by a long-term viscous process. Thus in this paper, we pay particular attention to this later phase.
Following our previous work [23], we employ a nonuniform grid for the two dimensional coordinates (x, z) in the simulation: For x ≤ x uni = 0.9GM BH /c 2 , a uniform grid is used with the grid spacing ∆x = 0.015GM BH /c 2 , and for x > x uni , ∆x is increased uniformly as ∆x i+1 = 1.01∆x i where the subscript i denotes the ith grid with ∆x i := x i+1 − x i . For z, the same grid structure as for x is used. The black-hole horizon is always located in the uniform grid zone. The location of the outer boundaries along each axis, L, is 6500-11000 km; for larger values of M BH , L is larger (see Table I).
In the analysis of the simulation results, we always derive the following quantities: average cylindrical radius R mat , average specific entropy s , and average electron fraction Y e both for the matter located outside the black hole and for the ejecta (see the method for identifying the ejecta below). Here, these average quantities TABLE I. Initial equilibrium models for the numerical simulation. Described are the model name, black-hole mass, rest mass of the disk, dimensionless spin of the black hole, the coordinate radii at the inner and outer edges of the disk (rin and rout), entropy per baryon (s/k), the value of n (that determines the rotational profile), and electron fraction (Ye) for the disk, the viscous coefficient (αν ), the scale height (H), and the location of outer boundaries along each axis in units of 10 where I and M mat denote a moment of inertia and rest mass of the matter located outside the black hole defined, respectively, by out implies that the volume integral is performed for the matter located outside the black hole. The integration is practically performed for the y = 0 plane with d 3 x = 2πxdxdz. The ejecta component is identified using the same criterion as in Ref. [23]. First, we identify a matter component with |hu t | > h min c 2 as the ejecta. Here h min denotes the minimum value of the specific enthalpy in the adopted equation-of-state table, which is ≈ 0.9987c 2 . For the matter escaping from a sphere of r = r ext , we define the ejection rates of the rest mass and energy (kinetic energy plus internal energy) at a given radius and time byṀ whereê := hαu t − P/(ραu t ). The surface integral is performed on a sphere of r = r ext with dS i = δ ir r 2 ext dθdϕ for the ejecta component. r ext is chosen to be 4000-5000 km in the present work.
Here, ρ √ −gu t (= ρ * ) obeys the continuity equation for the rest mass, and thus, the time integration for it gives a conserved quantity. In the absence of gravity, ρê √ −gu t also obeys the source-free energy-conservation equation, and far from the central region, the sum of its time integration and the gravitational potential energy of the escaped component are approximately conserved. Thus, the total rest mass and energy (excluding the gravitational potential energy) of the ejecta (which escape away from a sphere of r = r ext ) are calculated by In addition, we add the rest mass for the ejecta component located inside a sphere of r = r ext , M eje,in (t), giving the total ejecta mass on each time slice, M eje = M eje,esc + M eje,in . Far from the central object, E eje,esc is approximated by where U and T kin are the values of the internal energy and kinetic energy of the ejecta at r ext → ∞, respectively. The last term of Eq. (3.10) approximately denotes the contribution of the gravitational potential energy of the matter at r = r ext [23]. Since the ratio of the internal energy to the kinetic energy of the ejecta decreases with its expansion, we may approximate U/T kin ≈ 0, and hence, E eje,esc by E eje,esc ≈ M eje,esc c 2 +T kin +GM BH M eje,esc /r ext for the region far from the central object. We then define the average velocity of the ejecta (for the component that escapes from a sphere of r = r ext ) by In this subsection, we compare the results of the lowmass disk models with M BH = 3, 4, 6, and 10M and M disk = 0.1M , to clarify the quantitative dependence of the disk evolution and mass ejection on the blackhole mass. Here, the results with M BH = 3M were already derived in Ref. [23]. For all these models, the initial condition has the same specific entropy (s = 6k), the same angular velocity profile (n = 1/7). For the higher-mass black holes, the maximum density of the disk is lower because of the larger radius of the outer edge, and thus, the larger volume of the disk.
First we summarize general processes for the viscous evolution of the disk. We note that the mass and spin of the black hole do not change much because the disk mass is much smaller than the black-hole mass, and thus, the accretion of the matter into the black hole is a minor effect for the models with M disk = 0.1M .
In the very early stage of the evolution of the system (for t 100 ms), 60%-70% of the disk in rest mass falls into the black hole from its inner region due to the viscous angular momentum transport process (and partly due to the transition process associated with the growth of the viscous tensor). The maximum rest-mass accretion rate is ≈ 3M /s and it monotonically decreases for t 10 ms leading to < 0.1M /s for t 100 ms. Subsequently, the remaining part of the disk component expands gradually due to the viscous heating and angular momentum transport processes. By these, the maximum temperature and rest-mass density decrease monotonically with time: see t > 10 ms part of Fig. 1. In the relatively early stage with t 0.5-2 s, the maximum temperature of the disk is high, i.e., kT 3 MeV and the maximum rest-mass density is also larger than ρ ∼ 10 8.5 g/cm 3 (see the top panels of Fig. 1). As a result, the neutrino emissivity is preserved to be high, L ν 10 51 erg/s, due to the electron and positron capture on nucleons (see the bottomleft panel of Fig. 1). For this stage, the thermal energy generated by the viscous heating is consumed primarily by the neutrino cooling, and hence, the viscous heating cannot be efficiently used for the outward expansion of the disk. Indeed, the viscous heating rate of the order of ∼ νM disk Ω 2 is comparable to L ν in this stage. However, for the later stage in which the typical temperature of the disk decreases below ∼ 3 MeV/k [23], the neutrino cool-ing becomes inefficient, because the neutrino emissivity depends strongly on the matter temperature T approximately as ∝ T 6 (e.g., Refs. [11,13]). Then, the viscous heating is mostly used for the outward expansion of the disk.
In particular, for the late stage for which the neutrino cooling is inefficient, the convective motion is excited and blobs of the matter viscously heated in the vicinity of the black hole are moved toward the outer region of the disk. Here, the onset of the convective motion is reflected in a short-timescale variation in the curves of T max , ρ max , and L ν . This motion activates the mass ejection of the matter from the outer part of the disk. All these features are qualitatively the same as those found in our previous paper [23]. However, the evolution process and the property of the ejecta depends quantitatively on the black-hole mass as we state in the following.
As we described above, the mass ejection sets in when the typical temperature of the disk decreases below ∼ 3 MeV/k. Figure 1 shows that for the higher-mass black holes, the temperature is lower than that in the lowermass models at given time. As already mentioned, this is due to the fact that for the higher-mass black holes, the inner-edge radius of the disk has to be larger due to the larger horizon radius, and hence, the maximum temperature and rest-mass density are smaller from the beginning of the simulation. As a consequence, the neutrino luminosity as well as the efficiency of the neutrino emission (see the bottom panels of Fig. 1) is always lower for the higher-mass black holes. Thus, the acceleration of the disk expansion due to the viscous heating in the inefficient neutrino cooling stage sets in earlier for the higher-mass black holes.
The top-left panel of Fig. 2 shows the evolution of the average cylindrical radius R mat . In the early evolution stage of the disk during which kT max 3 MeV (see Fig. 1), R mat gradually increases with time primarily due to the viscous angular momentum transport. On the other hand, for the later stage of kT max 3 MeV, the increase of R mat is accelerated. This signals the onset of the mass ejection as we described in our previous paper [23]. In this mass ejection stage, the increase of the average specific entropy due to the viscous heating is also accelerated (see the top-right panel of Fig. 2). For the higher-mass black holes, the mass ejection sets in earlier. This mass ejection is triggered primarily by the onset of the convective motion. This is found from the fact that the time from which the short-timescale variation is found in the curves of T max and ρ max agrees approximately with the time at which the steep increases in R mat and s are found.
The dependence of the evolution timescale of the disk on the black-hole mass is reflected in the composition variation of the disk. In the early stage of the disk evolution prior to the onset of the mass ejection, electrons are relativistic and in a degenerate state because the restmass density is higher than ρ max ∼ 10 8.5 -10 9 g/cm 3 . As a result, the average electron fraction Y e is relatively Here the efficiency is defined by the total neutrino luminosity, Lν , divided by the rest-mass energy accretion rate of the matter into the black hole, c 2 dM fall /dt. Note that the short-timescale oscillation in the curves of Tmax, ρmax, and neutrino emission efficiency for the late time is due to the fact that convective motion is activated, and thus, the disk is disturbed.
low. After a significant fraction of the disk falls into the black hole and the disk relaxes to a quasi-steady evolution stage, the average value of Y e increases with the disk expansion (i.e., with the decrease of the rest-mass density).
After the typical temperature of the disk decreases below ∼ 3 MeV/k, not only the neutrino luminosity drops significantly, but also the weak interaction process freezes out. As a result, the average value of Y e relaxes to constants. Our simulation clearly shows this property irrespective of the black-hole mass (see the bottom-left panel of Fig. 2). For the higher-mass black holes, the average value of Y e for given time is always higher for t 10 −2 s due to the weaker electron degeneracy. The reason for this is that the efficiency of the neutrino emission (see the bottom-right panel of Fig. 1) is lower due to the lower temperature and density of the disk, and hence, the viscous heating is more efficient for higher-mass black holes.
On the other hand, the freezeout of Y e occurs earlier for the higher-mass black holes due to the earlier drop of the temperature and density, and hence, it is not trivial whether the final average value of Y e is higher or lower for the higher-mass black holes. Our present result shows that the final relaxed average value of Y e becomes slightly higher for the higher-mass black holes. Our interpretation for this is that the electron degeneracy in the disk is weaker when the weak interaction processes freeze out, and thus, the value of Y e in the equilibrium of the electron/positron capture at the freeze out is higher for higher-mass black holes.
The average entropy is also always higher for the higher-mass black holes. The reason for this is the neutrino cooling efficiency is lower, and thus, viscous heating is more efficiently used for enhancing its entropy for the higher-mass black holes. After the freezeout of the weak interaction at t ∼ 1 s, the intermittent production of the high-energy blobs, which eventually become ejecta, also contributes to the increase of the entropy.
We note that this conclusion would hold only for the comparison among the models with the same initial radius of the disk (i.e., with the same physical values of r out ). If we compare the models with different values of r out , the results could be modified (see, e.g., Ref. [24]).
The bottom-right panel of Fig. 2 shows that ∼ 15%-20% of the initial disk mass becomes the ejecta. Here, for the higher-mass black holes, the disk is more com-  pact in the sense that r out /M BH is smaller. As a result, the ejecta mass is smaller. This agrees qualitatively with the results of Ref. [24]. The fraction of the ejecta mass agrees quantitatively with that in our previous paper [23]. We note here that all the matter located outside the black hole in the late time eventually becomes the ejecta as this figure illustrates. One non-trivial point is that the increase of the ejecta mass is slower for the higher-mass black holes. Our interpretation for this is that although the mass ejection sets in earlier, the viscous heating timescale, which should be proportional to R 2 disk /ν ∝ M BH , is always longer for the higher-mass black holes. This effect should be reflected in the curve of M eje . Figure 3 displays the evolution for the average values of Y e and velocity of the ejecta component. Reflecting the evolution of Y e in the disk during its viscous expansion, the average value of Y e for the ejecta becomes 0.30-0.35 and is higher for the higher-mass black holes. This property also agrees qualitatively with a recent report in Ref. [24]. The velocity of the ejecta is ≈ 0.05-0.07c depending only weakly on the black-hole mass.
After the viscous evolution of the disk and subsequent mass ejection, the system settles to a spinning black hole surrounded by a geometrically-thick torus and a narrow funnel in the vicinity of the rotation axis. Such an outcome appears to be suitable for generating gamma-ray bursts, i.e., for driving a collimated jet. This topic is discussed in Sec. III H.

C. Evolution of black holes for high-mass disk models
In the following three subsections, we describe the results for the high-mass disk models of M disk ≈ 3M with the initial black-hole mass, M BH,0 = 4, 6, and 10M . We first summarize the evolution of the black holes as a result of matter accretion. Figure 4 shows the evolution of the mass and dimensionless spin of the black holes (top panels) together with the rest-mass infall rate into the black holes (bottom panels) for all the models with M disk ≈ 3M . The left panels compare the results for M BH,0 = 4, 6, and 10M with α ν = 0.05 and for M BH,0 = 10M with different values of the viscous coefficient. The right ones compare the results for M BH,0 = 10M with different initial disk radii and with and without neutrino absorption/irradiation.
For these high-mass disk cases, ≈ 70%-90% of the initial disk matter falls eventually into the black hole, and as a result, the black-hole mass increases approximately to M BH,0 + ζM disk with ζ ≈ 0.7-0.9. For the initially compact disks with α ν = 0.05 (models M04H05, M06H05, M10H05, and M10H05n), ≈ 90% of the initial disk matter eventually falls into the black hole irrespective of the models (cf. the bottom-right panel of Fig. 6), and thus, the final black-hole mass approaches approximately M BH,0 +0.9M disk . The rest-mass accretion rate is quite high (∼ 50M /s for α ν = 0.05) in the early stage for these models (see the bottom panels of Fig. 4). In this high-accretion rate stage, the neutrino emission efficiency is low because of the very high mass accretion rate (cf. the bottom-right panel of Fig. 5), and hence, the advection of the matter into the black hole is the major dissipation process of the matter energy. We note that the short-term very high accretion rate in t ≤ 10 ms for the larger values of α ν (≥ 0.1) is achieved during the relaxation stage of the accretion disk in which the viscous tensor changes from zero to the quasi-steady state, and thus, it is an artifact due to the setting of the initial condition.
The rest-mass accretion rate, dM acc /dt, decreases steeply with time after the peak is reached, and then, the growth of the black hole is essentially stopped for t 5 s irrespective of the initial conditions. Here, when the matter infall rate into the black hole drops steeply, the neutrino luminosity also drops, and as a result, the mass ejection sets in. An interesting finding is that the rest-mass accretion rate is approximately and universally proportional to t −3/2 after the peak is reached. As already mentioned, before the mass ejection is activated, the dominant dissipation process of the disk is the neutrino cooling, and thus, in this stage, the disk is in a neutrino-dominate accretion flow (NDAF) state [54] (see Sec. III D for more details).
For given values of M BH,0 and r out , the fraction of the disk matter that falls into the black hole (i.e., the value of ζ) is larger than that for the low-mass disk models (cf. Sec. III B). Our interpretation for this is that the selfgravity of the disk plays a role for confining the matter in the central region of the system.
The values of ζ naturally decrease with the increase of the initial disk radius r out (compare the results for models M10H05, M10H05w, and M10H05x): For M10H05w and M10H05x, only ∼ 80% and 70% of M disk falls into the black hole, respectively (cf. the bottom-right panel of Fig. 9). It is also found that the peak rest-mass accretion rate is significantly reduced with the increase of r out , indicating that the value of M fall is decreased. Thus, the fraction of the disk mass that falls into the black hole, and hence, the final black-hole mass depend strongly on the initial disk radius. This fraction also depends on the magnitude of the viscous coefficient (cf. the bottom-right panel of Fig. 6): For the larger viscous coefficients, the values of ζ are slightly smaller, and as a result, the final black-hole mass is slightly smaller. For these models, the rest-mass accretion rate drops steeply for an earlier time of t 1 s at which the mass ejection sets in (see Sec. III D). This leads to the consequence that the ejecta mass increases with the increase of α ν .
Irrespective of the initial black-hole mass, the dimensionless spin of the black hole, χ, monotonically increases toward an asymptotic value of χ BH 0.9 for α ν = 0.05. This asymptotic value depends on the initial black-hole mass because the fraction of the accreted angular momentum to that of the initial black hole, which is proportional to M disk /M BH,0 , is higher for the smallermass black hole in our setting. However, the asymptotic values are always smaller than unity. We note that because of the presence of the disk, the dimensionless angular momentum of the entire system defined by χ tot := cJ 0 /(GM 2 0 ), where J 0 and M 0 are the total angular momentum and initial gravitational mass, respectively, exceeds unity. Specifically, χ tot ≈ 1.20, 1.15, 1.06, 1.18, and 1.28 for models M04H05, M06H05, M10H05, M10H05w, and M10H05x, respectively. Reflecting these initial conditions, for the relatively low-mass black holes (M04H05 and M06H05), the value of χ BH exceeds 0.97, but the further increase is halted for the chosen viscous coefficient. The upper left and right panels of Fig. 4, respectively, show that for the larger viscous coefficient or for the larger initial extent of the disk, the value of χ BH is smaller. Both results are reasonable: For the larger viscous coefficient, the outward angular momentum transport becomes more efficient, and as a result, the angular momentum accretion into the black hole is suppressed. For the larger initial extent of the disk, the fraction of the matter that falls into the black hole becomes smaller as the bottom-right panel of Fig. 4 shows. This also suppresses the angular momentum gain by the black hole, leading to the smaller value of χ BH .
It is interesting to note that in the absence of the neutrino absorption/irradiation effect, the final values of M BH and χ BH are slightly smaller and larger than those in the presence of this effect, respectively. The smaller value of M BH may be interpreted as a consequence of less absorbed neutrino energy in the absence of neutrino absorption by the infalling matter. Indeed, the neutrino luminosity (i.e., the fraction of neutrinos that escape to infinity) is slightly higher in the absence of the neutrino absorption/irradiation (cf. the bottom-left panel of Fig. 8). By contrast, the reason for the larger value of χ BH in the absence of the neutrino absorption/irradiation is not clear. This result indicates that in the presence of the neutrino absorption/irradiation the angular momentum is more efficiently radiated away by neutrinos. Thus, the possible interpretation for this is that in the presence of the neutrino absorption/irradiation, neutrinos are frequently reprocessed in the central region of the disk and the final emission toward infinity occurs in a relatively far region from the black hole for which the specific angular momentum of the disk matter is larger than that in the central region, and hence, the rotation effect (a relativistic effect) increases the angular momentum loss by the emitted neutrinos.

D. Evolution of disks for high-mass disk models
Next we summarize the viscous evolution of the disk for the case of M disk ≈ 3M . In this subsection, we focus on models M04H05, M06H05, M10H05, M10H10, and M10H15.
As in the low-mass disk cases, in the early stage of the evolution of the system (in t 100 ms; see the bottom ). Note again that the short-timescale variation in the curves of Tmax, ρmax, and neutrino emission efficiency for the late time is due to the fact that convective motion is activated, and thus, the disk is disturbed significantly.
panels of Figs. 4 and 6), ∼ 60% of the rest mass of the disk falls into the black hole due to the viscous angular momentum transport processes, and subsequently, the remaining part of the disk component is evolved by the viscous heating and angular momentum transport processes. 1 Quantitatively, the evolution process depends on the black-hole mass, viscous coefficient, and initial extent of the disk. This subsection focuses on the dependence on the black-hole mass and viscous coefficient. Figures 5 and 6 display the evolution of several quantities of the disk for models M04H05, M06H05, M10H05, M10H10, and M10H15, i.e., for α ν = 0.05 with M BH,0 = 4, 6, and 10M and for M BH,0 = 10M with α ν = 0.10 and 0.15. For these cases, the maximum temperature and density are by a factor of ∼ 3 higher than those for the low-mass disk cases (compare Figs. 1 and 5). Associated with this situation, the neutrino luminosity reaches 10 54 erg/s, which is more than 10 times as high as that for the low-mass disk cases. We note that the emissivity of neutrinos via the capture process of electrons and positrons by nucleons is approximately proportional to T 6 . Thus, the neutrino emissivity is not directly reflected in the luminosity (this is also observed from a lowefficiency of the neutrino luminosity in the early stage of t 10 ms: see the bottom-right panel of Fig. 5). The reason for this is that the disk is so dense and hot that the emitted neutrinos are trapped [54]. However, the total energy carried away by neutrinos is still high, ∼ 10 53 erg (≈ 0.05M c 2 ), irrespective of the black-hole mass. Thus, more than 1% of the rest-mass energy of the disk is carried away by neutrinos. All these properties are found irrespective of the black-hole mass and the viscous coefficient.
The bottom-right panel of Fig. 5 shows that the efficiency of the neutrino emission depends strongly on the black-hole mass. As we discussed in Sec. III B, the primary reason for this may be that for the lower-mass black-hole cases, the temperature and density of the disk are higher than for the higher-mass cases, and as a result, the neutrino luminosity and the efficiency of the neutrino . We note that for models M04H05 and M06H05, the density and velocity profiles are highly disturbed in the initial transition stage of t 0.2 s, and as a result, the average velocity and electron fraction oscillate due to an artifact of the initial condition. emission could be also higher. However, Fig. 5 shows that the neutrino luminosity does not depend strongly on the black-hole mass, due to the trapping of neutrinos. Our alternative interpretation for this dependence of the efficiency on the black-hole mass is that the black-hole spin is reflected. As we showed in Fig. 4, the dimensionless spin of the black holes, χ, is higher for the lower-mass black holes and close to unity for M BH = 4 and 6M . For such rapidly spinning black holes, the horizon radius steeply decreases as χ approaches unity, and thus, the rest-mass accretion rates decrease as a result of the decrease in the horizon area. Simultaneously, as χ approaches unity, the gravitational potential around the black hole becomes deeper, and hence, the efficiency for releasing the gravitational potential energy is steeply enhanced [55]. Thus, it is reasonable that the neutrino emission efficiency steeply increases with the decrease of the black-hole mass (i.e., with the increase of the resulting dimensionless spin of the black hole). The maximum neutrino emission efficiency depends very weakly on the viscous coefficient (compare the results of M10H05, M10H10, and M10H15). This agrees with the result in our previous paper [23]. The timescale to reach the maximum is approximately proportional to α −1 ν , reflecting the viscous timescale. Overall, the viscous coefficient basically changes the timescale for the evolution of the disk. The curves of the disk quantities like the neutrino luminosity, maximum temperature, and maximum density become similar if we plot these quantities as functions of α ν t.
The maximum temperature and the neutrino luminosity are maintained to be higher than 2 MeV/k and 10 51 erg/s until t ∼ 2-10 s. Here, this timescale depends on the black-hole mass and viscous coefficient. The neutrino cooling timescale is maintained to be as long as the viscous heating timescale for t ∼ 1-5 s. During this NDAF stage [54], the thermal energy generated by the viscous heating is consumed by the neutrino cooling, and after this time, the mass ejection sets in. The onset of the mass ejection is signaled again clearly by the accelerated increase of s/k (see the top-right panel of Fig. 6).
As found from the curves of the ejecta mass M eje (see the bottom-right panel of Fig. 6) as well as of s/k , the mass ejection sets in at t ∼ 1-5 s, which is later than in the low-mass disk case. This is due to the higher disk mass and resulting longer-term high neutrino-emissivity stage (compare Figs. 1 and 5). Associated with the delay of the mass ejection, the electron fraction is enhanced during the longer-term disk evolution for the higher-mass disk case. The bottom-left panel of Fig. 6 shows that irrespective of the black-hole mass, the final average value of Y e exceeds 0.4 for α ν = 0.05, which is higher than that for the low-mass disk case. Therefore, we conclude that the ejecta from the system of a stellar-mass black hole and a high-mass disk with mass M has high values of the electron fraction with its average ∼ 0.4, if the mass ejection proceeds via the viscous effects.
For the larger viscous coefficient, the average value of Y e is higher before the mass ejection sets in (t 1 s) because of the stronger viscous effect that results in the faster decrease of the electron degeneracy and higher equilibrium value of Y e by the electron/positron capture. However, its final average value is lower for the larger viscous coefficients. The reason for this is that for the larger viscous coefficients, the viscous timescale of the disk is shorter, and thus, the freezeout of the electron/positron capture occurs earlier (see also Ref. [23]). For a very high value of α ν as 0.15, the final average value of Y e can be as low as ∼ 0.35. Thus, in the presence of an efficient mass-ejection mechanism which we do not take into account in the present work, the value of Y e could be also decreased. Figure 7 displays the evolution for the average values of Y e and velocity of the ejecta component. Reflecting the evolution of Y e in the disk during its viscous expansion, the average value of Y e for the ejecta becomes higher than that for the low-mass disk case: approximately 0.35-0.45. Here, the value depends on the black-hole mass and viscous coefficient. As in the low-mass disk cases, the final average value of Y e is higher for the higher black-hole mass while it becomes lower for the larger viscous coefficients.
The final average velocity of the ejecta, V ej , is ≈ 0.05-0.07c and, for a given value of α ν , it depends very weakly on the black-hole mass as in the low-mass disk case (compare the right panels of Figs 3 and 7). The value of V ej is larger for larger viscous parameters because of the stronger viscous effect (angular-momentum transport and viscous heating effects).
As in the low-mass disk case, after the viscous evolution of the disk and subsequent mass ejection, the system relaxes to a state composed of a spinning black hole surrounded by a geometrically-thick torus and a narrow funnel in the vicinity of the rotation axis. In particular, for the high-mass disk case, the resulting black hole is rapidly spinning with the dimensionless spin 0.9. Such an outcome would be quite suitable for generating gamma-ray bursts, i.e., for driving a collimated jet accelerated by thermal energy generated from the deep gravitational potential near the spinning black hole (see the discussion in Sec. III H).

E. Dependence on the initial disk model and neutrino irradiation
This subsection briefly shows the dependence of the disk evolution on the initial disk extent and neutrino irradiation referring to the models with M disk ≈ 3M and M BH = 10M . Figures 8 and 9 compare the results for models M10H05, M10H05w, M10H05x, and M10H05n. For the larger initial disk radii (compare the results of models M10H05w and M10H05x with those of M10H05), the evolution process is obviously delayed because the viscous timescale becomes longer for the larger initial disk radii (see Eq. (2.1)). Specifically, for the larger initial disk radii, the peak maximum temperature and neutrino luminosity are reached later (at t ∼ 50-100 ms) and associated with this, the neutrino luminosity remains higher for the later stage. This fact is also inferred from the bottom-right panel of Fig. 4, which shows that the peak of the rest-mass accretion is reached at t ∼ 30 ms and 60 ms for M10H05w and M10H05x, respectively. As a result of the delay in the evolution, the disk expansion timescale becomes longer and the mass ejection sets in slightly later. In spite of the difference in the onset time of the mass ejection, the final average value of Y e depends only weakly on the initial disk radius. This would be due to the fact that the physical condition at the onset of the mass ejection such as temperature, density, and neutrino luminosity depends only weakly on the initial disk radius. One significant difference among the models is found in the total amount of the disk mass that swallowed into the black hole. For the larger initial disk radii, the fraction of the disk matter that falls into the black hole is smaller naturally, and as a result, the fraction of the ejecta can be increased. However, besides this quantitative difference, the initial disk radius does not lead to the remarkable difference in the disk evolution and the property of the ejecta; e.g., for models M10H05w and M10H05x, the final average values of Y e and velocity are as high as those for model M10H05.
By comparing the results for models M10H05 and M10H05n, we find that the effect of the neutrino irradiation modifies the disk evolution process quantitatively. For example, in the absence of the neutrino irradiation, the neutrino emissivity is slightly enhanced and as a result, the onset of the disk expansion is delayed while the final average value of Y e is decreased due to the omission of the irradiation. However, these effects do not substantially change the disk evolution.
The properties of the ejecta reflect the disk evolution. In the absence of the neutrino absorption/irradiation, the asymptotic average value of Y e is ≈ 0.41 which is lower than that in the presence the neutrino absorption/irradiation (see Fig. 9). Thus, the neutrino absorption effect plays a role for quantitatively determining the electron fraction of the ejecta. In the absence of this neutrino effect, the asymptotic ejecta velocity is slightly smaller (≈ 0.05c) than in its presence (≈ 0.06c), reflecting the absence of the neutrino absorption/irradiation (i.e., neutrino radiation pressure).

F. Nucleosynthesis
Nucleosynthetic abundances are calculated for all the models in post-processing steps by using the reactionnetwork code, rNET [58], consisting of 6300 isotopes with the range of atomic number, Z = 1-110. The adopted theoretical rates for both light-particle (nucleon and αparticle) capture (TALYS [59]) and β-decay (GT2 [60]) are based on the microscopic mass prediction of HFB-21 [61]. The experimental rates, if available, are taken from REACLIB V2.0 [62]. Neutrino-induced reactions are not included in the present nucleosynthesis calculations. The thermodynamic histories of ∼ 5000 tracer particles for each model are deduced as described in our previous work [39]. Each nucleosynthesis calculation is started with the initial mass fractions of protons and neutrons being Y e  and 1 − Y e , respectively, at which the temperature in a given tracer particle decreases to 10 10 K (=10 GK; here 1 GK= 10 9 K). At such high temperature, nuclear statistical equilibrium (NSE) immediately establishes and thus any initial compositions with the same net charge, i.e., XZ/A = Y e (here X, Z, and A indicate the mass fraction, atomic number, and atomic mass number of a given isotope, respectively) should give the same result. We note that the Y e values at 5 GK are adopted for the initial composition.
As all nucleosynthesis-relevant quantities such as Y e , s, and V ej (or expansion timescale) span wide ranges across our explored models, the trend of resulting abundance distributions cannot be interpreted based on solely a single quantity (e.g., Y e ). Before presenting the nucleosynthesis results, therefore, we analyse these relevant quantities in some detail by utilizing the criterion for the production of heavy r-process nuclei described in Ref. [64], for 0.38 < Y e < 0.5 (neutron-deficient condition) with t exp being the duration of the expansion from 5 GK (the end of NSE) to 2.5 GK (the beginning of neutron capture) in units of second. Note that the original value of the coefficient in Eq. (3.14), 5 × 10 −4 in Ref. [64], is replaced by 5.5 × 10 −4 such that f (Y e ) becomes continuous at Y e = 0.38. Here, the proton-to-nucleon ratio of the seed nuclei synthesized in nuclear quasi-statistical equilibrium (QSE) is assumed to be 0.38.
According to Ref. [64], the r-process nuclei with A = 200 are expected to be abundantly produced for the case that the condition of C r > 1 is satisfied. This criterion also can be used to inspect the productivity of nuclei beyond the seeds (A ∼ 80-90) even for C r < 1. Figure 10 shows the dependence of C r on Y e for fixed values of s and t exp . We find that a combination of higher s, lower Y e , and shorter t exp gives a larger value of C r , i.e., a condition for synthesizing heavier nuclei. Among these quantities, a change of t exp has a relatively small impact on C r as evident from Eq. (3.12) and Fig. 10. As pointed out in Ref. [64], for Y e < 0.2 (out of the range for the use of Eq. (3.12)), the heavy r-process nuclei are synthesized less dependently on s and t exp . However, such a condition is not met in our explored models except for M10H05n (without neutrino irradiation). Figures 11-14 plot the mass histograms of the ejecta as a function of Y e , s/k, t exp , and C r for all the models employed in this paper. The mass-averaged values of these quantities are also presented in Table II. Here, Y e (T = 5 GK) is the electron fraction of each tracer particle determined at the time when its temperature decreases to 5 GK. For the low-mass disk cases, the electron fraction is distributed between ∼ 0.2 and 0.5, and for the high-mass disk cases, it is between ∼ 0.3 and 0.5 (except for M10H05n, for which the heating/irradiation is absent). As we showed in Secs. III B and III D, the average values of Y e for the ejecta is between ≈ 0.30 and 0.35 for low-mass disk models and between ≈ 0.35 and 0.50 for high-mass disk models (see also Y e in Table II). Thus the distribution is consistent with the average. Note that the cutoff of the distribution at the high-Y e end in each panel of Fig. 11 is a consequence of the fact that we limit the range as Y e ≤ 0.55 in our simulations.
The distribution of the electron fraction depends quantitatively on the black-hole mass, disk mass, and the values of α ν . Figure 11 shows that the electron fraction is higher for the higher-mass black-hole models. The bottom-left panel of Fig. 11 also shows that the electron fraction is lower for the larger viscous coefficient case (i.e., for the case that the mass ejection timescale is shorter). However, even in the extremely high value of α ν = 0.15, the value of Y e is always higher than 0.25; ejecta is not extremely neutron-rich for the high-mass disk case. The bottom-right panel of Fig. 11 indicates that the changes of r in or r out only slightly modify the Y e distribution. This panel also confirms that the inclusion of neutrino Here, the value of Ye is determined at the time when the temperature of each ejecta component decreases to 5 × 10 9 K (referred to as 5 GK). The cutoff of the distribution at the high-Ye end is due to the fact that we limit the range being Ye ≤ 0.55 in our simulations.
irradiation is necessary to obtain a reliable Y e distribution. Figure 12 shows that the typical specific entropy of the ejecta is 10-30k irrespective of the black-hole mass, the disk mass, the value of α ν , and the disk size. We find, however, that the distribution extends to higher entropy as dM/ds ∝ s −3 (note that the bin spacing in Fig. 12 is proportional to log s). As a result, nonnegligible amounts of ejecta have the entropy values exceeding s/k = 100 in many of our explored models. Figure 13 shows the expansion timescales defined by t exp = t(T = 2.5 GK) − t(T = 5 GK). For the low-mass disk cases (see the top-left panel of Fig. 13), the typical expansion timescale is ∼ 100 ms irrespective of the blackhole mass. For the high-mass disk cases, this timescale is longer, typically ∼ 500 ms (see Fig. 13). A reason for this trend is that the radius, at which the temperature of the fluid becomes ∼ 5 GK, is larger for higher-mass disks due to their higher temperature. This makes the timescale for the decrease of the temperature longer, because the temperature in radiation-dominated material drops approximately as ∝ r −1 .
From the distributions of Y e , s/k, and t exp presented in Figs. 11-13, the distribution of C r is derived by using Eq. (3.12). Here, the tracer particles with electron fractions only in the range of 0.2 < Y e < 0.5, for which Eq. (3.12) is applicable, are considered. We find that the ejecta satisfying C r > 1 are subdominant or absent in our explored models. It is anticipated from Fig. 10 that the component with s/k > 100 is necessary to meet the criterion C r > 1 for 0.2 < Y e < 0.5 and 0.01 < t exp /s < 1. Therefore, the heavy r-process nuclei synthesized in our models should originate from high entropy ejecta with s/k > 100. Such ejecta cannot be the dominant component in the present cases because of the entropy distribution scaled approximately as dM/ds ∝ s −3 (Fig. 12) with the lowest value of s/k ∼ 10. This is evident also from the average values of C r , C r < 0.1, in Table II. We note that the high entropy ejecta with s/k > 100 can be only mildly neutron-rich. The viscous and neutrino heating inevitably induces weak interaction, which leads to an increasing trend of the specific entropy with Y e . As a result, the electron fractions in the ejecta with s/k > 100 become Y e 0.35.
The high-entropy component of the ejecta originates from the innermost part of the disk, for which the specific entropy of the matter is increased efficiently by the viscous heating. The material is ejected intermittently by the buoyancy force after the weak interaction in the disk freezes out. The high-entropy component with s/k 100  has a moderate electron fraction of Y e ∼ 0.35 for the lowmass disk models (e.g., M10L05), while the high-entropy component has a high electron fraction of Y e 0.45 for the high-mass disk models (e.g., M10H05). This is because the temperature of the disk is lower and the weak interaction timescale is longer for the low-mass disk models, and thus, the initial moderate-value electron fraction in the outer part of the disk is preserved more easily than that for the high-mass disk models. For higher viscosity models, M10H10 and M10H15, a part of the highentropy component has somewhat low electron fraction of Y e ∼ 0.35-0.45 because of the lower freeze-out values of Y e for the higher-viscosity cases.
The mass fraction of nucleosynthetic products as a function of atomic mass number, X(A), is shown in Fig. 15 for all the models (solid lines). The dotted lines indicate the results excluding the tracer particles with s/k > 100. The r-process residuals to the solar system abundances [57] (filled circles; hereafter referred to as the solar r-residuals) are also shown, which are vertically shifted to match the value of X(82) for M10L05 (the top-left panel of Fig. 15) or M10H05 (the other panels of Figs. 15). We find that the heavy r-process nuclei with A > 132 (beyond the neutron shell closure of N = 82) are exclusively synthesized in the ejecta with s/k > 100 as anticipated from the analysis of C r value (except for M10H05n).
For the low-mass disk models (M disk ≈ 0.10 M ), the abundance distributions follow the low-A (= 80-110) side of the solar r-residual pattern (see the top-left panel of Fig. 15). However, the heavy r-process abundances with A > 132 are deficient compared to the scaled solar rresiduals by more than three orders of magnitude. The result is consistent with the small values of C r (∼ 0.06-0.09; Table II). This is a weak r-process signature that is also found for model K8 (gray line) [23] as well as for the post-merger ejecta from massive neutron stars studied in our previous work [39]. Although the amounts are small, we find that these nuclei (A > 132) are more abundant for higher M BH models, despite a tendency of their higher values of Y e (Fig. 11 and Table II). This is due to the fact that the ejecta for higher M BH models contain larger amounts of matter with higher values of s/k (> 100; the top-left panel of Fig. 12) and thus C r ( 1; the top-left panel of Fig. 14).
For the high-mass disk models (M disk ≈ 3.0 M ), the abundance distribution is characterized by a sharp peak at A ≈ 82 (see Fig. 15), which is formed in modestly neutron-rich (Y e = 0.3-0.4) and low-entropy (s/k = 10-20) ejecta. Compared to the low-mass disk models, the overall electron fraction (Fig. 11) and expansion timescale (Fig. 13) for the high-mass disk models are higher and longer, respectively, while those of entropy are similar (Fig. 12). Such a physical condition, i.e., the modest neutron-richness, long expansion timescale, and relatively low entropy leads to nucleosynthesis in QSE rather than by neutron capture [63,64]. This also can be anticipated from the smaller values of C r (= 0.02-  Table II) in the high-mass disk models, which indicate weaker productivity of heavy nuclei beyond the seeds (A = 80-90). As a result, the abundance pattern exhibits a sharp peak at 82 Se (synthesized as 82 Ge that has the proton-to-nucleon ratio of 0.390) with negligible mass fractions of the nuclei beyond A ≈ 132.
The change of M BH modifies the abundance patterns in the range of A = 100-130 (see the top-right panel of Fig. 15) for the high-mass disk models, reflecting the difference of the entropy distribution on the high-s side (> 50k; see the top-right panel of Fig. 12). This is also evident from the top-right panel of Fig. 14, which shows the larger ejecta mass with C r > 0.1 for model M10H05 than those for other models.
Higher values of α ν (= 0.10 and 0.15) lead to the production of the nuclei with A > 132 (see the bottomleft panel of Fig. 15) because of larger amounts of highs (> 100k) components (see the bottom-left panel of Fig. 12) and thus of C r > 1 (see the bottom-left panel of Fig. 14). For model M10H10 (and M10H15), the heavy r-process nuclei with A > 140 are synthesized entirely in the high-entropy ejecta with s/k > 100. A similar trend has also been found in our previous work for models with M BH = 3 M and M disk = 0.1 M (K8h and K8s in Ref. [23], although the heavy r-process nuclei originate also from the lowest-Y e ( 0.2) components for the small disk mass). In fact, the exclusion of ejecta with s/k > 100 (dotted lines in the bottom-left panel of Fig. 15) diminishes the abundances for A > 132 (near or below the lower-bound of the plot), which confirms the origin of heavy r-process nuclei being from the high-entropy ejecta. However, these nuclei are under-abundant by about 3 orders of magnitude compared to the solar rresiduals (scaled at A = 82).
The change of the disk size slightly modifies the abundance patterns for A > 90 as found in the bottomright panel of Fig. 15 that compares models M10H05, M10H05w, and M10H05x. This is due to the overall lower entropy for the wider disk models (see the bottomright panel of Fig. 12) because of longer viscous timescale due to larger initial radius of the disk. This leads to the weaker productivity of nuclei beyond the seeds (i.e., smaller C r ; see the bottom-right panel of Fig. 14 and Table II). We note that, for the high-mass disk models, the neutrino absorption on free nucleons plays a significant role for shaping the Y e distribution by diminishing the low-Y e ( 0.3) component (M10H05 and M10H05n in the bottom-right panel of Fig. 11), which suppresses the production of nuclei with A > 132 (see the bottomright panel of Fig. 15). This is due to the high neutrino luminosity ( 10 54 erg/s at the peak; cf. Fig. 5), which are about 10 times higher than those for the low-mass disk models (cf. Fig. 1). Note that, for model M10H05n (for which the neutrino irradiation is absence), the production of nuclei with A > 132 is due to the presence of low-Y e ∼ 0.2 ejecta (the solid and dotted lines in the bottom-right panel of Fig. 15 are overlapped).

G. Implications from the Nucleosynthesis Calculation
We presume that our low-mass and high-mass disk models represent the remnants for the merger of black hole-neutron star and the core-collapse of rapidly rotating massive stars (collapsars), respectively. In this respect, the upper limits for the frequencies of such events can be obtained from the nucleosynthesis results as what follows. Here, we do not consider the contribution from the earlier ejecta (if any) of black hole-neutron star mergers or collapsar. By including this contribution, the constraint on the frequency will be even stronger. Table III (seventh and eighth columns) presents the maximum overproduction factor (p max ) and the relevant isotope for each model. The overproduction factor is defined by where X(Z, A) and X (Z, A) denote the mass fractions of the isotope with Z and A for each model and in the solar system [65], respectively. In Fig. 16, the overproduction factors for M10L05 (left) and M10H05 (right) are plotted as representatives of the low-mass and highmass disk models, respectively. The maximum Galactic fraction of black hole-neutron star mergers or collapsars represented by each model with respect to that of core-collapse supernovae (CCSNe), f max ( 1 , Table III; ninth column), can be estimated as [63] f  [65]. Here, the total ejecta mass for each model (M ej,tot ; third column in Table III) estimated as in our previous study [39] is adopted instead of the ejecta mass at the end of simulation (M ej ; second column in Table III) that is still increasing as found in Figs. 2, 6, and 9. The fraction f max can be translated to the maximum rate of black holeneutron star mergers or collapsars represented by each model in the Galaxy, R max (Table III; last column), by adopting the inferred Galactic CCSN rate ≈ 2.30 × 10 −2 yr −1 [66].
In Table II, we find f max ∼ 2 × 10 −3 and R max ∼ 40 Myr −1 for the low-mass disk models. These values are similar to what are expected for binary neutron star mergers to be the dominant r-process site in the Galaxy and in the solar system [67,68]. However, as shown in the top-left panels of Figs. 15 and 16, the low-mass disk models account for the origin of isotopes in the range A = 80-110 only. In fact, spectroscopic [69,70] and galactic chemical evolution studies [71,72] imply the presence of the production site of light trans-iron elements, such as Sr, in addition to the main r-process site. Our results indicate that the black hole-accretion disks followed by black hole-neutron star mergers can be the second sources of such elements. Spectroscopic studies have also revealed the presence of metal-poor stars that exhibit a descending trend of trans-iron abundance patterns [73,74], which is presumed to be a signature of weak r-processing [75]. Such a weak r-process-like abundance patterns may be explained by our low-mass disk models (see also a similar result for the post-merger ejecta from massive neutron star-accretion disks in our recent work [39]). It should be stressed, however, that the contribution from the early dynamical ejecta should be added to assess the full nucleosynthetic outcomes from black hole-neutron star mergers.
For the high-mass disk models, we find f max ∼ (0.02-0.9) × 10 −3 and R max ∼ 5-20 Myr −1 (see Table II). Provided that collapsars are the dominant sources of long duration gamma-ray bursts, the local volumetric rate of collapsars is estimated to be R GRB /f b ≈ 260 Gpc −3 yr −1 , where R GRB ≈ 1.3 Gpc −3 yr −1 is the local volumetric rate of long gamma-ray bursts pointing toward the Earth [85] and f b ≈ 5 × 10 −3 is the beaming factor [86]. This can be translated to the Galactic rate of long gamma-ray bursts ≈ 26 Myr −1 , i.e., about 0.1% of the Galactic CCSN rate, by using the number of Milky Way analogous galaxies ≈ 0.01 Mpc −3 yr −1 [87]. These values are in good agreement with the upper bounds for the high-mass disk models. Thus, our models, in particular M10H05 (f max = 0.94 × 10 −3 and R max = 22 Myr −1 ; Table II), may reasonably represent the black hole-accretion disks formed in collapsars.
As the sources of heavy nuclei, the black hole-accretion disks formed in collapsars can contribute only to the light trans-iron elements with A ∼ 80 according to the right panel of Fig. 16. It is noteworthy, however, that the overproduction factor of the neutron-rich isotope 48 Ca accounts for ∼ 0.1 of p max = p( 82 Se). This implies that collapsars can be the sources of 48 Ca up to ∼ 10% of its total amount in the solar system. In Table III (fourth column), the ejecta mass of 48 Ca calculated as M ej,tot X( 48 Ca) is presented for all the models. Currently, the astrophysical sources of 48 Ca are unknown, for which the suggested sites include high-density type-Ia supernovae [76,77] and electron-capture supernovae (both core-collapse [78] and thermonuclear types [79]). It should be noted that neutrino-induced reactions are not included in the present nucleosynthesis calculations. The high neutrino luminosity for the high-mass disk models (Fig. 5) may lead to a νp-process [80][81][82][83][84] that produces some proton-rich isotopes for which the astrophysical origins are uncertain (e.g., 92 Mo).
The type-Ic supernovae accompanying with long gamma-ray bursts are considered to be powered by a large amount (0.1-0.6 M [88]) of 56 Ni. Indeed, the pioneering work of collapsars [31] has predicted the mass ejection of 56 Ni as large as ∼ 1 M , assuming Y e = 0.5 in the accretion disk. However, the ejecta from our highmass disk models contain only ∼ 0.01-0.02 M of 56 Ni (calculated as M ej,tot X( 56 Ni); fifth column in Table III). This is due to the fact that the ejecta in our models have relatively small amounts of matter with Y e > 0.49 (see Fig. 11), for which 56 Ni is efficiently produced [63]. This may indicate that the bulk of 56 Ni in this type of supernovae comes from the early ejecta [89], for which the driving mechanism is unsettled ( neutrino heating or magnetic pressure). It is also suggested that the detonation of the late-time accreted material can be an additional source of 56 Ni [90]. Figure 17 displays snapshots for the profiles of the restmass density, temperature, specific entropy, and electron fraction for models M04L05, M10L05, M04H05, M10H05, M10H15, and M10H05w at late times (t ≈ 5.5-7.1 s). By these times, mass ejection has already set in and the disk has relaxed to a quasi-steady state. These figures show that the outcomes of the evolution of the system are composed of a rapidly spinning black hole surrounded by a torus and a narrow funnel in the vicinity of the rotation axis, irrespective of the black-hole mass and initial disk mass. The dimensionless spin is larger for the larger initial disk mass. The rest-mass density in the funnel region is 10 2 g/cm 3 for the initially low-mass disk models and 10 2.5 g/cm 3 for the initially high-mass disk models. The density is close to the artificial atmosphere density, and thus, in reality, the density might be even lower. The specific entropy in the funnel region is always very high, s/k 100, and this indicates that the radiation pressure dominates over the gas pressure. By contrast, the region outside the funnel, the specific entropy is of the order of 10k.

H. Relation with gamma-ray bursts
The half opening angle of the funnel region, which we define as the region with the density less than 10 3 g/cm 3 , is ∼ 0.1 rad and is narrower for the smaller-mass black holes and for the initially larger-mass disks. This angle depends weakly on the initial disk radius although it is wider for the large viscous coefficient. For all the models that we studied, the total rest mass in the funnel region is of the order of 10 −6 M in the computational region, which may be suitable for a high-energy relativistic jet passing through avoiding the baryon loading. All these facts indicate that the outcome of the viscous evolution of the system of a rapidly spinning black hole and a massive disk is suitable for generating gamma-ray bursts [3,4].
In the late-evolution stage of the disk with t > 1 s for which the matter temperature satisfies kT < 3 MeV, the neutrino emissivity is much smaller than the viscous heating. In such a stage, viscous heating is fully used not only for the disk expansion and resulting steady mass ejection but also for driving a high-velocity outflow in the vicinity of the rotation axis and torus surface. In particular, near the innermost region of the disk, the viscous heating rate is high and the high velocity outflow is generated. Since the matter density is low above such a region, the high velocity is naturally realized. The order of the viscous heating rate is estimated by where M disk,inn denotes the rest mass of the inner region of the disk which contributes to the efficient viscous heating, and we supposed that the black-hole mass would be 5-7M . Note that around rapidly spinning black holes with χ 0.95, a small orbital radius of ∼ 2GM BH /c 2 is possible and in such a case, Ω can be ≈ 10 4 rad/s for M BH = 7M .
The isotropic luminosity of typical long gamma-ray bursts is 10 51 erg/s [3,4]. If the beaming effect (by which the total luminosity can be smaller) and energy conversion efficiency to gamma-rays (by which higher luminosity is necessary) are taken into account, the viscous heating rate shown in Eq. (3.17) is a substantial fraction of the long gamma-ray burst luminosity. Unfortunately, for the viscous heating, the matter is always accompanied, and hence, it is not possible to drive the outflow of the Lorentz factor to ∼ 100 by this heating effect due to the baryon loading problem. However, an outflow by this heating effect can play an important role for cleaning up the region in the vicinity of the rotation axis to reduce the density there.
In the presence of electromagnetic fields, the Blandford-Znajek effect [91] could provide another energy injection mechanism by extracting the huge rotational kinetic energy of the remnant black hole. In the presence of a radial magnetic field on the black-hole horizon B r , the outward Poynting flux from the horizon is written as (e.g., Ref. [92]) F BZ = 2 (GM BH ) 2 c 3 (B r ) 2 ωr + (Ω H − ω) sin 2 θ, (3.18) wherer + = 1 + 1 − χ 2 , Ω H = c 3 χ/(2GM BHr+ ), and ω is a rotational frequency of the electromagnetic field on the horizon for which one often assumes ω = Ω H /2. Then, assuming that B r and ω are uniform on the horizon, we obtain the total luminosity  Here, a strong magnetic field of the order of 10 14 G would be achieved if the magnetorotational instability efficiently enhances the magnetic-field strength, B, until the equipartition is satisfied in the disk, i.e., B 2 /(4π) ∼ ρ disk c 2 s with ρ disk being the typical rest-mass density of the disk at t 1 s which is ∼ 10 8 g/cm 3 and with c s ∼ 0.1c. The luminosity, L BZ , is comparable to the isotropic luminosity of long gamma-ray bursts, and thus, the Blandford-Znajek effect could also be a reasonable energy source.
The funnel structure shown in Fig. 17 looks suitable for confining a jet launched from the central part. In the funnel region, the rest-mass density of the matter is likely to be much lower than the electromagnetic pressure B 2 /8π, and hence, a force-free magnetosphere would be established. By contrast in the geometrically thick torus next to the funnel, the rest-mass density increases steeply to 10 6 g/cm 3 for which the rest-mass energy density of the matter exceeds the magnetic pressure. Therefore, only in the narrow funnel, a clear spiral-shape magnetic field would be established and an efficient particle acceleration would be achieved, leading to generating a narrow jet in the funnel.

IV. SUMMARY
As an extension of our previous study [23], we performed viscous neutrino-radiation hydrodynamics simulations for accretion disks surrounding a spinning black hole with M BH = 4, 6, and 10M and initial dimensionless spin χ ≈ 0.8. We consider a compact disk with M disk ≈ 0.1 or 3M and with the outer edge located at r out = 200-1000 km. For M disk ≈ 0.1M , we find that ∼ 20% of M disk is ejected and the average electron fraction of the ejecta is Y e = 0.30-0.35. Here, Y e is slightly higher for the larger values of M BH . Nucleosynthesis calculation shows that only light r-process elements with mass number of 80-110 are efficiently synthesized in such ejecta, because of the absence of the low-Y e component with Y e 0.2. Thus, the results obtained are qualitatively the same as those found in Ref. [23] irrespective of the black-hole mass.
For high-mass disks with M disk ≈ 3M , we found that the luminosity of neutrinos exceeds 10 54 erg/s in the early stage of the disk evolution and neutrino cooling is the dominant cooling process for the first seconds. This timescale is several times longer than that for M disk ≈ 0.1M for given black-hole mass. We also found that (i) ∼ 10% of M disk is ejected for the compact disk models with the disk extent ∼ 200 km and (ii) the ejecta mass increases with the increase of the initial disk radius. Irrespective of the models, Y e of the ejecta is enhanced to be 0.35 because the electron fraction is increased significantly during the long-term viscous evolution of the disk until the neutrino cooling timescale becomes longer than the viscous heating timescale. Our nucleosynthesis calculation indicates that not r-process elements but trans-iron elements with atomic mass number ∼ 80 are predominately synthesized in the matter ejected from a massive torus surrounding stellar-mass black holes. This suggests that if the remnant of the collapsars is composed of a black hole and a massive disk, such remnants may not be the sites for the r-process nucleosynthesis.
By the matter accretion, the dimensionless spin of the black hole, χ, increases. For the initially high-mass disk case, the value of χ exceeds 0.9. In particular for the case that the initial ratio of the disk mass to the black-hole mass is larger than 1/2, it exceeds 0.95. This value depends on the viscous coefficient, and for the larger viscous coefficient, the resulting value of χ is slightly smaller. For the high-mass disk case, the outcome after the viscous evolution is the system composed of a rapidly spinning black hole with χ 0.9, a geometrically thick torus, and a narrow funnel in the vicinity of the rotation axis. In particular, due to the outflow activity in a late stage of the disk evolution, the rest-mass density in the rotation axis becomes quite low 10 2 g/cm 3 , and the total rest mass in the funnel region becomes as small as 10 −6 M . The outcome appears to be suitable for driving gammaray bursts. Since the viscous heating is unlikely to be its central engine, we cannot reproduce the gamma-ray bursts in our simulation. However, the final outcome that we find in this paper indicates that in the presence of magnetic fields for which the electromagnetic energy is comparable to the thermal energy of the disk together with a rapidly spinning black hole, the Blandford-Znajek mechanism is likely to provide the energy injection required for launching the gamma-ray bursts.