Mass ejection from disks surrounding a low-mass black hole: Viscous neutrino-radiation hydrodynamics simulation in full general relativity
Sho Fujibayashi, Masaru Shibata, Shinya Wanajo, Kenta Kiuchi, Koutarou Kyutoku, Yuichiro Sekiguchi
MMass ejection from disks surrounding a low-mass black hole: Viscousneutrino-radiation hydrodynamics simulation in full general relativity
Sho Fujibayashi, Masaru Shibata,
1, 2
Shinya Wanajo,
1, 3
KentaKiuchi,
1, 2
Koutarou Kyutoku,
4, 5, 3, 2 and Yuichiro Sekiguchi Max Planck Institute for Gravitational Physics (Albert Einstein Institute),Am Muhlenberg 1, Potsdam-Golm 14476, Germany Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto, 606-8502, Japan Interdisciplinary Theoretical Science (iTHES) Research Group, RIKEN, Wako, Saitama 351-0198, Japan Department of Physics, Kyoto University, Kyoto 606-8502, Japan Department of Particle and Nuclear Physics, the GraduateUniversity for Advanced Studies (Sokendai), Tsukuba 305-0801, Japan Department of Physics, Toho University, Funabashi, Chiba 274-8510, Japan (Dated: January 15, 2020)New viscous neutrino-radiation hydrodynamics simulations are performed for accretion disks sur-rounding a spinning black hole with low mass 3 M (cid:12) and dimensionless spin 0.8 or 0.6 in full generalrelativity, aiming at modeling the evolution of a merger remnant of massive binary neutron starsor low-mass black hole-neutron star binaries. We reconfirm the following results found by previousstudies of other groups: 15–30% of the disk mass is ejected from the system with the average velocityof ∼ PACS numbers: 04.25.D-, 04.30.-w, 04.40.Dg
I. INTRODUCTION
The first direct detection of gravitational waves fromthe final stage of an inspiraling binary neutron starsystem (GW170817) by advanced LIGO and advancedVIRGO [1] was accompanied with a wide variety ofthe follow-up observations of electromagnetic counter-parts [2]. This event heralded the opening of theera of the multi-messenger astronomy composed ofgravitational-wave and electromagnetic-counterpart ob-servations, and demonstrated that the observation ofelectromagnetic signals plays a key role for understand-ing the merger and subsequent mass ejection processes ofneutron-star binaries, which cannot be understood onlyfrom the gravitational-wave observation.A popular interpretation for the merger and post-merger evolution of binary neutron stars in GW170817is as follows (e.g., see Refs. [3–6]). After the merger ofa binary neutron star, a hypermassive neutron star wasformed with an accretion disk around it. Subsequently,the hypermassive neutron star survived for ∼ . . M (cid:12) , which might be the central engineof a gamma-ray burst associated with GW170817 [7, 8].At the merger and during the post-merger stage, ejec- tion of matter with mass ∼ . M (cid:12) occurred. At themerger, a neutron-rich material, which can synthesizeheavy r -process elements, was ejected in the dynamicalprocess, while in the post-merger stage, the mass ejec-tion occurred from a disk (or torus) surrounding the rem-nant hypermassive neutron star and the black hole sub-sequently formed. The observational results suggest thatthe neutron richness of the post-merger ejecta is unlikelyto be very high. A canonical interpretation for this isthat the neutrino irradiation from the hypermassive neu-tron star to the ejecta is strong enough to reduce theneutron richness. However, physically well-modeled andwell-resolved numerical simulations taking into accountfull general relativity, neutrino transport, and angularmomentum transport by magnetohydrodynamics or vis-cous hydrodynamics effects have not been performed yet.Thus it is not clear whether the canonical interpretationis really correct.Although the electromagnetic observation forGW170817 provides us rich information for theneutron-star merger, it will not be always the casethat the next events have the similar feature for theelectromagnetic counterparts as the GW190425 eventsuggests [9]. For GW170817, the total mass of twoneutron stars was not so large that the remnant couldform a massive neutron star at least temporarily. On a r X i v : . [ a s t r o - ph . H E ] J a n the other hand, for more massive binary neutron stars orfor the case of black hole-neutron star binaries, a blackhole surrounded by a disk is the expected remnant. Forthis case, the neutrino emission as strong as that fromhypermassive neutron stars is absent, and the propertiesof the post-merger ejecta can be different from those forGW170817 [10–21]. It is worthy to explore in detail theejecta properties for the case that a black hole is theimmediate remnant of neutron-star mergers.There are a lot of previous work for exploring the ejectafrom the system composed of a black hole and a disksurrounding it. However, the previous work has beenperformed in some simplifications. In the pioneer simu-lation work by Fern´andez, Metzger, and their collabora-tors [10, 11, 13], they only qualitatively took into accountthe general relativistic effects and radiation transfer ef-fects of neutrinos. In particular, the spacetime structurearound the black hole was only qualitatively considered.In Refs. [12, 14], the authors took into account the neu-trino irradiation effect carefully, but again they did not oronly qualitatively did take into account the general rela-tivistic effects and black hole spacetime. In Refs. [15–19],the authors performed a magnetohydrodynamics simu-lation taking into account the general relativistic effectwith a fixed Kerr black hole as the background space-time. However, in Refs. [15–18], the neutrino transfereffect was only approximately taken into account and/orthe equation of state employed was an approximate one(e.g., effects of heavy nuclei are not incorporated, andin some work, degenerate pressure of electrons, which isthe key in the dense disk, is not taken into account). InRef. [19], numerical simulations are performed only for ashort term, ∼ .
13 s (i.e., an optimistic magnetic fieldfor enhancing the mass ejection is initially prepared),although general relativity (background spacetime of ablack hole), a realistic equation of state, and a detailedradiation transport are incorporated together for the firsttime.In the present work, we perform a long-term fullygeneral relativistic viscous neutrino-radiation hydrody-namics simulation for black hole-disk systems approxi-mately taking into account the neutrino irradiation effect.By focusing on the viscous evolution process, this workcan provide a complementary aspect for the evolution ofblack-hole accretion disks to the magnetohydrodynamicswork. In particular, we carefully resolve the black holespacetime and the inner part of the disk, and take intoaccount the effect of the self-gravity of the disk and theblack-hole evolution by matter accretion self-consistently.Particular emphasis is put on the point that we resolvethe vicinity of the black hole including the inner regionof the disk with the resolution higher than those in theprevious simulations. This setting enables us to followaccurate viscous evolution of the disk.The paper is organized as follows. In Sec. II, we brieflysummarize the basic equations employed in the presentsimulation study, and then, we describe the method toprepare the initial condition composed of the equilibrium state of disks. Dependence of the property of the disk sur-rounding a black hole on the equation of state and veloc-ity profile is also presented. Section III presents numeri-cal results for the simulations, focusing on the propertiesof the ejecta and nucleosynthesis in the matter ejectedfrom the disks. Section IV is devoted to a summary.Throughout this paper, G , c , and k denote the gravita-tional constant, speed of light, and Boltzmann’s constant,respectively. II. METHODS OF NUMERICALCOMPUTATIONA. Basic equations
We evolve black hole-disk systems in the framework offull general relativity. For the case that the disk mass ismuch smaller than the black-hole mass, we are allowed toperform a simulation in a fixed background of black-holespacetime. However, for the case that the disk mass is asubstantial fraction of the black-hole mass, such assump-tion may break down. In the last decade, the numerical-relativity community has established methods to evolvesystems composed of black holes for a long timescale sta-bly and accurately (e.g., Ref. [22] for a review). Thus,we do not have to assume the fixed background. In thepresent work, we numerically solve both Einstein’s equa-tion and matter-field equations self-consistently.The major purpose of this paper is to clarify the vis-cous evolution of a system composed of a self-gravitatingdisk surrounding a spinning low-mass black hole, whichis the plausible outcome formed after the merger of alow-mass black hole-neutron star binary or a binary neu-tron star of high total mass. For the disk evolutionin nature, the magnetohydrodynamics or viscous heat-ing/angular momentum transport is one of the key pro-cesses. In addition, the neutrino cooling of the diskand neutrino irradiation to the matter are key processesfor the disk evolution and for determining the propertyof the matter ejected from the system. In this work,we choose basic equations as follows: Einstein’s equa-tion, the viscous-hydrodynamics equations, the evolutionequation for the viscous tensor, the evolution equationsfor the lepton fractions including the electron fraction,and neutrino-radiation transfer equations. Here, for theequations for the lepton fractions, we take into accountelectron and positron capture, electron-positron pair an-nihilation, nucleon-nucleon bremstrahlung, and plasmondecay [20]. For solving Einstein’s equation, we em-ploy the Baumgarte-Shapiro-Shibata-Nakamura formal-ism [23] together with the puncture formulation [24], Z4cconstraint propagation prescription [25], and 5th-orderKreiss-Oliger dissipation. In reality, the viscous angularmomentum transport is likely to be induced effectivelyby a magnetohydrodynamics process [26]. In the presentwork, this process is approximated by viscous hydrody-namics. All these basic equations are the same as thosein Ref. [20] in which we performed a simulation for thesystem of a massive neutron star and disk.One crucial difference of the present simulation fromthe previous ones [20, 21] is that we have to evolve aspinning black hole located at the center stably and ac-curately for a long timescale, at least ∼ χ = 0 . x , should be smaller than ∼ . GM BH /c (see Sec. III) where M BH is the initialmass of the black hole: For M BH = 3 M (cid:12) , the grid spac-ing has to be smaller than ∼
80 m. Otherwise, the areaand dimensionless spin of the black hole, respectively, in-crease and decrease spuriously and significantly. In thepresent work, we employ ∆ x = 0 . GM BH /c . We alsoperformed simulations with ∆ x = 0 . GM BH /c and0 . GM BH /c for a particular model (referred to as K8model in Table I), and show the convergence propertyin Sec. III F. We note that such a high resolution helpsaccurately resolving the inner region of the disk.As mentioned above, a sufficiently high grid reso-lution is the key for a reliable simulation. Further-more, we have to evolve the disk for the timescale of (cid:38) × GM BH /c ∼ x, y, z ) and evolve the system only in the ( x, z )plane (the y = 0 plane).As already mentioned, in addition to Einstein’s equa-tion, we solve the viscous-hydrodynamics equation, evo-lution equations for the viscous tensor and lepton frac-tions, and radiation transfer equation. These equationsare solved in the cylindrical coordinates composed of( x, z ).For the x and z directions, the following non-uniformgrid is used for the present numerical simulation: For x ≤ x = 0 . GM BH /c , a uniform grid is used, andfor x > x , the grid spacing ∆ x is increased uniformlyas ∆ x i +1 = 1 . x i where the subscript i denotes the i th grid with x = 0 at i = 0. For z , the same gridstructure as for x is used. The black-hole horizon isalways located in the uniform grid zone. The locationof the outer boundaries along each axis is chosen to be ≈ GM BH /c ≈ M BH / M (cid:12) ) km in this study. B. Initial condition
We prepare an axisymmetric equilibrium state for theblack hole-disk system in the framework of full generalrelativity as the initial condition for our numerical sim-ulation. For the gravitational field equations, we employa puncture formulation developed in Ref. [29] with thefollowing line element: ds = − α c dt + ψ [ e η ( dr + r dθ )+ r sin θ ( β ϕ dt + dϕ ) ] . (2.1)Here, we employ spherical polar coordinates and α , ψ , β ϕ , and η are functions of r and θ . Note that η and β vanish for the non-spinning black hole. In the following,we denote the dimensionless spin parameter by χ , whichis chosen to be 0.8 or 0.6.We consider the case that the fluid four velocity hasthe form of u r = 0 = u θ and u ϕ = Ω u t , and the fluid isisentropic. Here Ω is the angular velocity which is also afunction of r and θ . Then, under the assumption of theisentropy, the Euler equation is integrated to give thefirst integral in the form hu t + 1 c (cid:90) hu ϕ d Ω = C, (2.2)where h is the specific enthalpy and C is a constant. Us-ing the rest-mass density ρ , the specific internal energy ε , and the pressure P , h is written as c + ε + P/ρ . Us-ing the normalization relation for u µ , u µ u µ = − c , u t iswritten as u t = (cid:0) α − c − γ ϕϕ (Ω + β ϕ ) (cid:1) − / , (2.3)where γ ϕϕ = ψ r sin θ .In this paper, we simply assume that j = c − hu ϕ is afunction of Ω in the form of j = A n Ω − n , (2.4)where A n and n are constants (see Ref. [30] for morecareful choice). In the Newtonian limit, j ≈ u ϕ ≈ (cid:36) Ωwith (cid:36) = r sin θ , and hence, Ω ∝ (cid:36) − / ( n +1) . Thus, for n = 1 /
3, the Keplerian angular velocity is recovered inthe Newtonian limit and for n = 0, the specific angularmomentum is constant. Since the angular velocity profileof the remnant disk of neutron-star mergers is close to theKeplerian, we should try to employ the value of n as closeas 1 /
3. We find that the disk mass in the equilibriumstates with 1 / ≤ n < / (cid:46) GM BH /c (see discussion below). Thus, we employthe values of n ≤ / n and an equation of state, wehave to determine Ω, h , A n , and C . In the following,we describe the basic equations for this procedure. First,Eq. (2.4) is rewritten asΩ n (Ω + β ϕ ) hu t γ ϕϕ = A n c , (2.5)where we used u ϕ = u t (Ω + β ϕ ) γ ϕϕ . (2.6)Using Eq. (2.4), Eq. (2.2) is rewritten as hu t + A n Ω − n − n = h (cid:18) u t + u ϕ Ω c (1 − n ) (cid:19) = C. (2.7)From Eqs. (2.3) and (2.5), we also obtain A n Ω − n [ c α − γ ϕϕ (Ω + β ϕ ) ] γ ϕϕ (Ω + β ϕ ) = (cid:18) C − A n Ω − n − n (cid:19) . (2.8)This algebraic equation is used to determine Ω for givenvalues of A n and C with gravitational fields computed. A n and C are determined by choosing the inner andouter edges of the disk in the equatorial plane. Since thevalues of h are identical at such edges, Eq. (2.7) gives (cid:18) u t + u ϕ Ω c (1 − n ) (cid:19) in = (cid:18) u t + u ϕ Ω c (1 − n ) (cid:19) out , (2.9)where “in” and “out” in the subscripts indicate the quan-tities at the inner and outer edges, respectively. Here, u t and u ϕ are given by Eqs. (2.3) and (2.6), respectively.In addition, Eq. (2.4) is written to ( j Ω n ) in = ( j Ω n ) out ,which leads to ( u ϕ Ω n ) in = ( u ϕ Ω n ) out . (2.10)Equations (2.9) and (2.10) constitute simultaneous equa-tions for Ω in and Ω out , and hence, by solving these equa-tions, we first determine these angular velocities. Subse-quently, using Eqs. (2.5) and (2.7), A n and C are deter-mined for the value of h at ρ = ρ min , where the minimumvalue of h and ρ min are found for a given tabulated equa-tion of state. Then, Ω at each point is determined bysolving Eq. (2.8). Once A n , C , and Ω at each point aredetermined, h is determined from Eq. (2.7).For a given value of h , thermodynamical variables, ρ , ε , and P , are determined through the equation of stategiven. In this work, we employ a tabulated equation ofstate based on the DD2 equation of state [31] for a rel-atively high-density part and the Timmes equation ofstate for the low-density part [32]. We choose the lowestrest-mass density to be 0 . / cm in the table. In thisequation of state, ε , P , and h are functions of ρ , Y e , and T where Y e and T are the electron fraction and mattertemperature, respectively. Thus, to determine ρ , Y e , and T from h , we need to employ two conditions for relat-ing these three variables. One condition is just the as-sumption that the specific entropy s is constant becauseEq. (2.2) is derived under this condition. For the other,we adopt two relations. In one case, we assume a relationbetween Y e and ρ in the form ρ ( Y e ) (see Fig. 1) and inthe other case, we simply set Y e = 0 .
1. For the formercase, we derive an approximate relation of ρ ( Y e ) fromour numerical results for the remnant of binary neutronstar mergers [33] (see also, e.g., Refs. [34–36]), for which Y e ρ (g/cm )
10 100 k T ( M e V ) x (km) 0 0.1 0.2 0.3 0.4 0.5 Y e ρ ( g / c m ) FIG. 1. Top: Two relations for Y e employed for construct-ing equilibrium states of disks surrounding the black hole.Bottom: Density, electron fraction, and temperature ( kT ) asfunctions of the radius in the equatorial plane for a disk model(model K8 in Table I). the value of Y e in the major part of the disk is approx-imately determined by weak interaction processes. Forthis model, the value of Y e is larger for the smaller den-sity, because the effect of the electron degeneracy, com-pared to the thermal effects, becomes weaker. We basi-cally employ this as the fiducial model, and the equationof state is simply referred to by specifying the value of s/k . The equation of state with Y e = 0 . Y e = 0 . M BH = 3 M (cid:12) , χ = 0 . M disk = 0 . M (cid:12) , r in =2 . GM BH /c , r out ≈ GM BH /c , s/k = 6, and n = 1 / (cid:46) g / cm andthe maximum temperature (in terms of kT ) is ∼
20 40 60 80 100 120 0 0.05 0.1 0.15 0.2 r ou t / M B H n s =6 ks =7 ks =8 ks =7 k , Y e =0.1 s =8 k , Y e =0.1 M d i s k / M B H n s =6 ks =7 ks =8 ks =7 k , Y e =0.1 s =8 k , Y e =0.1 FIG. 2. Left: r out / ( GM BH /c ) as a function of n for the case of M disk /M BH = 1 /
30. Right: M disk /M BH as a function of n for r out = 40 GM BH /c . For both panels, r in = 2 . GM BH /c and χ = 0 .
20 40 60 80 100 120 0 0.05 0.1 0.15 0.2 r ou t / M B H n s =6 k , χ =0.8 s =7 k , χ =0.8 s =6 k , χ =0.6 s =7 k , χ =0.6 M d i s k / M B H n s =6 k , χ =0.8 s =7 k , χ =0.8 s =6 k , χ =0.6 s =7 k , χ =0.6 FIG. 3. The same as Fig. 2 but for the comparison between χ = 0 . .
6. For χ = 0 . r in = 2 . GM BH /c and2 . GM BH /c , respectively. as a result, the electron fraction is low ≤ .
1. This is atypical structure of the disk around low-mass black holeswith M disk = O (0 . M (cid:12) ).Figure 2 displays several relations for the radius of theouter edge, r out , and baryon mass, M disk , of the disk fora variety of the values of n and s/k . In this plot, we fix M BH = 3 M (cid:12) , χ = 0 .
8, and r in = 2 . GM BH /c . On theother hand, we employ a variety of the equations of statewith s/k = 6–8 and a wide range of n . The left panelshows the outer edge of the disk, r out , as a function of n for M disk = 0 . M (cid:12) . For the larger values of n toward 1 / M disk , the extent ofthe disk (i.e., r out ) needs to be increased with n . It is alsofound that for the higher value of s/k , the extent of thedisk has to be larger. The reason for this is that for thehigher value of s/k , the pressure for a given value of thedensity is larger, and hence, the overall density becomes relatively small.The right panel of Fig. 2 shows M disk /M BH as a func-tion of n for r out = 40 GM BH /c . This clarifies that forthe larger values of n , the disk mass is smaller for thegiven extent of the disk, and for the smaller values of s/k the disk mass is larger. It is also found that for the Y e = 0 . s/k . This is reasonable because for the smaller valueof Y e , the electron degeneracy pressure is less important,and hence, to enhance the pressure, a higher value of therest-mass density is needed.In Fig. 3, we compare the results for χ = 0 . χ = 0 .
6, the inner edge of the disk cannot be assmall as that for χ = 0 . χ for givenblack-hole mass. Here, we set r in = 2 . GM BH /c for χ = 0 . . GM BH /c for χ = 0 .
8. It is foundfor χ = 0 . r out becomes larger for givenparameters of s , n , and M disk than for χ = 0 .
8. It is alsofound that the disk mass is smaller for the given valuesof s , n , and r out = 40 GM BH /c than for χ = 0 .
8. Thereason for these results is that for the smaller values of χ ,the gravitational potential of the black hole is shallower,and hence, the amount of the bounded material becomessmaller. Thus, to obtain a model with the same diskmass and same value of r out / ( GM BH /c ) for χ = 0 . χ = 0 .
8, we need to prepare a smaller value of s or asmaller value of n .The simulations are performed for several models asthe initial conditions. First, we employ models of s/k = 6and n = 1 / Y e case and s/k = 8and n = 1 / Y e = 0 . χ = 0 . r in , r out ) ≈ (2 GM BH /c , GM BH /c ) and the mass ofthe black hole and the baryon rest-mass of the disk to be M BH = 3 M (cid:12) and M disk = 0 . M (cid:12) , respectively. We alsoemploy a compact disk model (D8) for which s/k = 6, n = 1 /
10, ( r in , r out ) ≈ (2 GM BH /c , GM BH /c ) withthe same mass and spin, M BH = 3 M (cid:12) , χ = 0 .
8, and M disk = 0 . M (cid:12) , as those for models K8 and Y8. We willshow that the slight difference in s and n does not dras-tically change overall dynamics of the disk and the prop-erties of the ejecta, although for initially more compactdisks, the ejecta mass becomes smaller as a natural conse-quence. In addition, we employ one model with χ = 0 . M BH = 3 M (cid:12) , and M disk = 0 . M (cid:12) in which n = 1 / s = 6 k , and ( r in , r out ) ≈ (2 . GM BH /c , GM BH /c ) toexplore the dependence of the numerical results on theblack-hole spin.We note that the outer edge of the disk shouldnot be very large for modeling the merger remnantof neutron-star binaries; numerical-relativity simulationshave shown that r out is between 100–200 km. For thisreason, we fiducially set r out ≈ GM BH /c which is ≈
180 km for M BH = 3 M (cid:12) . For the compact model(D8), r out ≈
130 km, which is also a reasonable valueIn addition to these models, we employ initial condi-tions with M disk = 0 . M (cid:12) and r out ≈ GM BH /c , andwith M disk = 0 . M (cid:12) and r out ≈ GM BH /c . For bothmodels, we employ the fiducial Y e equation of state, with χ = 0 . M BH = 3 M (cid:12) , r BH = 2 GM BH /c , n = 1 /
7, and s = 6 k , which are the same as those for model K8 (seeTable I). By performing simulations for these initial con-ditions, we examine the effect of the mass (i.e., the effectsby the density and temperature) of the disk on the evo-lution of the system. We note that the initial conditionwith M disk = 0 . M (cid:12) is a good model for neutron-starmergers, but that with 0 . M (cid:12) is not. Rather, such ahigh-mass disk may be a good model for the remnant ofthe massive stellar core collapse to a black hole.Because the mass ratio, M disk /M BH , of the initial con-ditions employed here is fairly large (1/100–1/6) andhence the disks are weakly self-gravitating, they may besubject to non-axisymmetric deformation even if the an-gular velocity profile is close to the Keplerian. The previ-ous work in general relativity (e.g., Refs. [39–41]) shows that if the self gravity of the disks is not extremely large,spiral arms are formed and contribute to angular mo-mentum transport by the gravitational torque exertedby the non-axisymmetric structure. In our simulation,such non-axisymmetric effects cannot be taken into ac-count, but the angular momentum transport is incorpo-rated through the viscous hydrodynamics. Moreover, theprevious work [39–41] shows that the density enhance-ment in the spiral arms is not very strong and the angu-lar momentum transport is much less efficient than thatin the viscous hydrodynamics with the alpha viscous pa-rameter [42], α ν = O (10 − ). Therefore we suppose thatthe non-axisymmetric deformation effects would not bevery important for the models that we employ here. C. Method for analysing ejecta
Here, we briefly summarize how we identify the mat-ter as ejecta. The unbound matter should be consid-ered as ejecta. In this work, we employ the followingcondition for identifying matter in an unbound orbit: | hu t | > h min c where h min denotes the minimum valueof the specific enthalpy h in the chosen tabulated equa-tion of state, which is ≈ . c , and u t is a negativequantity. The reason why h min is smaller than c is thatthe effect of the binding energy of nucleus is present inthe equation of state.To analyse the ejecta, we first extract the outgoingcomponent of the matter at the radius of r ext = 2000–4000 km and identify the ejecta. Here, by changing theextraction radius, we examine the convergence of theejecta mass. In addition, we analyze the matter locatedwithin a sphere of r = r ext and the component with | hu t | > h min c is identified as the ejecta. By summingup these two components, we determine the quantities ofthe ejecta.For the ejecta component escaping from a sphere of r = r ext , we define the ejection rates of the rest mass andtotal energy at a given radius by˙ M eje := c − (cid:73) ρ √− gu i dS i , (2.11)˙ E eje := c − (cid:73) ρ ˆ e √− gu i dS i , (2.12)where g denotes the determinant of the spacetime metricand ˆ e := hαu t − P/ ( ραu t ). The surface integral is per-formed at r = r ext with dS i = δ ir r dθdϕ for the ejectacomponent. Here, we note that ρ √− gu t obeys the conti-nuity equation of the rest mass and ρ ˆ e √− gu t obeys theenergy conservation equations in the absence of gravity.Hence, for the region far from the central object, the timeintegration of these quantities are conserved. Thus, byperforming the time integration, the total rest mass andenergy of the ejecta (which escape away from a sphere of TABLE I. Initial conditions for the numerical simulation. Described are the model name, black hole mass, disk mass, black-holedimensionless spin, the radii at the inner and outer edges of the disk ( r in and r out ), entropy per baryon ( s/k ), n , and electronfraction ( Y e ) for the disk, and α ν H . The units of the mass are M (cid:12) and the units of r in and r out are GM BH /c ≈ .
43 km. Thelast column shows whether the neutrino irradiation is switched on or off.Model M BH M disk χ r in r out s/k n Y e α ν H (km) neutrino irradiationK8 3.0 0.10 0.8 2.0 41 6 1/7 0.07–0.5 0.45 yesK8h 3.0 0.10 0.8 2.0 41 6 1/7 0.07–0.5 0.90 yesK8s 3.0 0.10 0.8 2.0 41 6 1/7 0.07–0.5 1.35 yesK8n 3.0 0.10 0.8 2.0 41 6 1/7 0.07–0.5 0.45 noD8 3.0 0.10 0.8 2.0 29 6 1/10 0.07–0.5 0.45 yesY8 3.0 0.10 0.8 2.0 41 8 1/8 0.1 0.45 yesK6 3.0 0.10 0.6 2.8 41 6 1/10 0.07–0.5 0.45 yesC8 3.0 0.50 0.8 2.0 58 6 1/7 0.07–0.5 0.45 yesE8 3.0 0.03 0.8 2.0 29 6 1/7 0.07–0.5 0.45 yes r = r ext ) are calculated by M eje , esc ( t ) := (cid:90) t ˙ M eje dt, (2.13) E eje , esc ( t ) := (cid:90) t ˙ E eje dt. (2.14)In addition, we add the rest mass for the ejecta compo-nent located inside a sphere of r = r ext , M eje , in ( t ), givingthe total ejecta mass of the ejecta, M eje = M eje , esc + M eje , in .We note that far from the central object, E eje , esc isapproximated by E eje , esc ≈ M eje , esc c + U + T kin + GM BH M eje , esc r ext , (2.15)where U and T kin are the values of the internal en-ergy and kinetic energy of the ejecta at r ext → ∞ , re-spectively. The last term of Eq. (2.15) approximatelydenotes the contribution of the potential binding en-ergy to E eje , esc , which cannot be neglected for r ext (cid:46) GM BH /c ≈ v eje ,is ∼ . c and ( v eje /c ) is of the order of 10 − . Sincethe ratio of the internal energy to the kinetic energyof the ejecta decreases with its expansion, we may ap-proximate U/T kin ≈
0, and hence, E eje , esc by E eje , esc ≈ M eje , esc c + T kin + GM BH M eje , esc /r ext for the region farfrom the central object. We then define the average ve-locity of the ejecta (for the component that escapes froma sphere of r = r ext ) by v eje := (cid:115) E eje , esc − M eje , esc c − GM BH M eje , esc /r ext ) M eje , esc . (2.16)We note that the correction of the gravitational poten-tial energy GM BH M eje , esc /r ext is important for r ext (cid:46) km, and just by taking into account this correction,the values of v eje become only weakly dependent on theextraction radius, r ext . III. NUMERICAL RESULTSA. Setting
Numerical computations are performed for the blackhole-disk systems summarized in the previous section (seealso Table I). For the viscous-hydrodynamics, we needto input the viscous coefficient ν . In this work, we set ν = α ν hc s H/c where α ν is the dimensionless viscous co-efficient (the so-called alpha parameter), c s is the soundvelocity, and H is a scale height. We basically employ α ν = 0 .
05 taking into account the result of recent mag-netohydrodynamics simulations of Refs. [16, 18], whichindicates that the magnitude of the effective viscous pa-rameter is high with α ν ≈ . H = 9 km ≈ GM BH /c .That is, we set it approximately equal to the radius at theinnermost stable circular orbit around the Kerr black holeof χ = 0 .
8. For the outer part of the disk, the value ofthe scale height might be larger than H ≈ c s / Ω in the standard accretiondisk theory [42]. However, for the non-stationary system,this is not likely to be the case, if we suppose that theviscosity is effectively enhanced by the turbulence causedby the magneto-rotational instability (MRI) [26], becausethe exponential growth rate of the MRI is proportionalto the local angular velocity, Ω (for the Keplerian angu-lar velocity profile, the growth rate is 3Ω / (cid:38) / (3Ω / π/ Ω ≈
310 ms for the radius of r = 1000 km with theblack-hole mass of 3 M (cid:12) assuming the Kepler motion ofthe disk, while the timescale of the mass ejection, which istriggered by the viscous heating in the inner part of thedisk, is ∼ . α ν / . − ( H/ − s (see Sec. III B).Thus, it would not be realistic to employ a large value ofit for the outer part of the disk. Employing c s / Ω as thescale height could result in an overestimated value of thescale height for a large value of c s , and the mass ejectionprocess could be induced from an outer part of the diskin an unrealistic manner for the early evolution stage ofdisks with t (cid:28) ν , we change α ν H from 0 .
45 km to 0 .
90 kmand 1 .
35 km for the model of K8 series (see Table I).We note that in our setting of ν , varying H is equiva-lent to varying α ν ; numerical results for α ν = 0 .
05 and H = 9 km are the same as those, e.g., for α ν = 0 .
01 and H = 45 km. However for simplicity, in the following, theviscous coefficients with α ν H = 0 .
90 km and 1 .
35 km arereferred to as α ν = 0 .
10 and α ν = 0 .
15 with H = 9 km,respectively.The viscous timescale (for heating and angular mo-mentum transport) is written approximately by τ vis := R ν ≈ .
37 s( hc − ) − (cid:16) α ν . (cid:17) − (cid:16) c s . c (cid:17) − × (cid:18) H (cid:19) − (cid:18) R
50 km (cid:19) , (3.1)where R denotes the cylindrical radius of the disk. Aswe show in Sec. III B, the evolution timescale for ourchoice of the viscous coefficient is indeed of the order of0.1 s. We note that the dynamical timescale of the diskis approximately τ dyn := 2 π (cid:115) R GM BH = 3 . (cid:18) R
50 km (cid:19) / (cid:18) M M (cid:12) (cid:19) − / , (3.2)and hence, τ vis is much longer than τ dyn . Thus, if the sys-tem evolves by the viscous process, the evolution shouldproceed in a quasi-stationary manner.For many models in this paper, the simulations areperformed taking into account the neutrino irradiationeffect. For one model, K8n, we switch off the neutrinoirradiation effect to examine whether it is important ornot. In this work, we do not incorporate a heating effectby the neutrino pair annihilation [20], because the neu-trino luminosity is not very high for most of the evolutionstage of the system in the models employed. Only in thevery early stage with t (cid:46)
10 ms, the neutrino luminosityis high ( (cid:38) erg / s in total), and thus, the neutrino pairannihilation heating may play a role in the evolution ofthe disk and associated mass ejection. However, in thispaper we should not consider mass ejection in such anearly stage because the system is initially in a spuriouslyvarying phase due to a rather artificial initial condition(composed of a stationary equilibrium state of the disk ofthe perfect fluid), and thus, it is not very clear whetherthe numerical results show some physical phenomena ornot for such an early stage. In one model with high diskmass (model C8: see Table I), the neutrino pair annihila-tion heating may be important, but we do not take into account also for this model because model C8 is employedjust for the comparison with low-disk mass models andthe comparison should be done in the same simulationsetting.As described in Sec. II A, numerical simulations arecarried out typically with the grid spacing of ∆ x =0 . M BH near the black hole. This grid resolution ishigher than the resolution in previous general relativisticmagnetohydrodynamics simulations [15, 16, 18, 19], andthis setting enables us to perform a well-resolved simula-tion for the inner part of the disk (i.e., the thermal andgeometrical structure of the disk can be well resolved).For model K8, we perform simulations varying the gridresolution as ∆ x = 0 . M BH and 0 . M BH to confirmthe only weak dependence of the numerical results on thegrid resolution. B. Viscous hydrodynamics of disks for M disk = 0 . M (cid:12) This subsection focuses on presenting the results forthe models with M disk = 0 . M (cid:12) , except for model D8.Because the disk is more compact for model D8 thanfor others, the fraction of the disk matter that falls intothe black hole is larger and the ejecta mass is smallerthan those for the other less-compact disk models. Be-sides this difference, the numerical results depend onlyweakly on the initial disk compactness. Thus, we onlybriefly summarize the results for model D8 in AppendixB, comparing them with those for model K8, and in thissubsection, we only show the results for other models.Figures 4 and 5 display the evolution of the profiles forthe rest-mass density, temperature, specific entropy perbaryon, and electron fraction for models K8 and K8h.Figure 6 plots (a) the evolution of the black-hole massand dimensionless spin for models K8 with three differ-ent grid resolution, K8h, and Y8, and (b) the total restmass swallowed by the black hole for all the models withthe initial disk mass M disk = 0 . M (cid:12) . The mass and di-mensionless spin of the black hole during the evolutionare approximately determined in the same method as inRef. [43]: We calculate the area, A AH , and circumferen-tial radii around the equatorial and meridian planes, c e and c p , for the apparent horizon of the black hole, andthen estimate the mass and dimensionless spin, assumingthat A AH , c e , and c p are written as functions of the mass, M BH , and dimensionless spin, χ , of Kerr black holes asin the vacuum black-hole case.The first three panels of Fig. 7 show the evolution ofthe average cylindrical radius R mat , average value of thespecific entropy (cid:104) s (cid:105) , and average value of Y e , (cid:104) Y e (cid:105) , forthe matter located outside the black hole. Here, these FIG. 4. Snapshots for the rest-mass density in units of g / cm , temperature ( kT ) in units of MeV, specific entropy per baryonin units of k , and electron fraction Y e for selected time slices for model K8. Only for the first panel ( t = 0) the plotted region is300 km ×
300 km, and for others, it is 2000 km × Y e , and temperatureof the atmosphere artificially added is ≈
10 g / cm , 0.5, and ≈ .
036 MeV /k , respectively (cf. the first panel). average quantities are defined by R mat := (cid:114) IM mat , (3.3) (cid:104) s (cid:105) := 1 M mat (cid:90) out ρ ∗ s d x, (3.4) (cid:104) Y e (cid:105) := 1 M mat (cid:90) out ρ ∗ Y e d x, (3.5) where ρ ∗ := ρ √− gu t , and I and M mat denote a momentof inertia and rest mass of the matter located outside the0 FIG. 5. The same as Fig 4 but for model K8h. black hole, defined by I := (cid:90) out ρ ∗ ( x + y ) d x, (3.6) M mat := (cid:90) out ρ ∗ d x. (3.7) (cid:82) out implies that the volume integral is performed for thematter outside the black hole. The last three panels ofFig. 7 show the total neutrino luminosity L ν , an efficiencyof the neutrino emission defined by the total neutrinoluminosity, L ν , divided by the rest-mass energy accretionrate of the matter into the black hole, c dM fall /dt , andthe total ejecta mass M eje as functions of time. Note thatthe neutrino luminosity is defined by the total neutrinoemission rate minus the neutrino absorption rate, bothof which are calculated by the volume integral.As Figs. 4–7 show, the disk evolves approximately onthe viscous timescale defined by Eq. (3.1). We note thatthe typical value of c s is 0 . c at R ∼
100 km and c s is a decreasing function of R . In the early stage of theevolution of the disk with the timescale less than 200 ms,a substantial fraction of the inner part of the disk withsmall values of R falls into the black hole. Irrespective ofthe models and grid resolutions, ∼ M disk , falls into the black hole during thisearly stage, and as a result, the black-hole mass increasesby ∼ . . M disk as found in Fig. 6(a). Subsequently,the infall rate significantly decreases (see Fig. 6(b)) andfor t (cid:38) r = 2000–4000 km), which istypically ˙ M eje ∼ − M (cid:12) /s at the peak. Figure 6(a) alsoshows that for t (cid:46)
200 ms, the dimensionless spin slightlyincreases due to the matter accretion onto the black hole.We note that the subsequent gradual decrease of χ is dueto the insufficient grid resolution (cf. Sec. III F).As Fig. 4 and Fig. 7(a) show, the outer part of thedisk expands outwards spending more than hundredsmilliseconds by the viscous effect. Since the local vis-cous timescale for the outer part of the disk with a largevalue of R is longer than the inner part, the expansiontimescale there could be much longer than the timescaleof the matter infall to the black hole. However, becausethe viscous heating/angular momentum transport in theinner part of the disk also contributes to the disk expan-sion, the outer part of the disk expands on the timescalemuch shorter than the local viscous timescale of Eq. (3.1).For α ν = 0 .
05, an outer edge of the disk in the equato-1 (a) . . . . . M B H / M B H , K8 ( α ν = 0 . α ν = 0 . Y e = 0 . . . . . . t (sec)0 . . . χ (b) . . . . . t (sec)0 . . . . . M f a ll ( M (cid:12) ) K8 ( α ν = 0 . α ν = 0 . α ν = 0 . Y e = 0 . χ = 0 . FIG. 6. Left: Evolution of the black-hole mass and dimensionless spin for models K8, K8h, and Y8. For model K8, the resultswith three different grid resolutions are shown. The gradual decrease of χ is due to the insufficient grid resolution. Right: Totalamount of the rest mass swallowed into the black hole for models K8 (with two grid resolutions), K8h, K8s, K8n, Y8, and K6. rial plane with its rest-mass density ∼ g / cm reaches ∼ km at t ∼ . α − ν as predicted from Eq. (3.1) (e.g., compareFigs. 4 and 5, and see Fig. 7).For α ν = 0 .
05, the mass ejection sets in (or strictlyspeaking, the ejecta component appears) at ∼ . ρ (cid:38) g / cm isfound for r (cid:38) km). A signal of the mass ejection isidentified by the evolution of R mat in Fig. 7(a), whichshows a steep increase of this quantity from ∼
300 kmto higher values irrespective of the models. The pri-mary driving force of the mass ejection is the viscousheating in the inner region of the disk under the neg-ligible neutrino cooling. For t (cid:38) . kT ) below ∼ t (cid:46)
100 ms during which the total neutrino luminosity is (cid:38) erg / s, the viscous heating is mostly consumed bythe neutrino emission, and thus, the viscous effect cannothave power to eject matter. Indeed, the viscous heatingefficiency, ∼ νM disk Ω , is of the same order of L ν for thetypical values of ν , M disk , and Ω of the disk.The viscous heating is always most efficient in the in-nermost region of the disk (i.e., the region closest to theblack hole). The enhancement of the specific entropythen triggers convective motion from the innermost re-gion toward outer regions, after the neutrino cooling be- comes inefficient. We evaluated the Solberg-Holland fre-quency for the convective instability [44, 45], and foundthat the innermost region of the disk at r ∼ ∼ (10 ms) − , and hence,the timescale for the convective motion is much shorterthan the viscous timescale.However, the convective activity is suppressed by thecentrifugal force, which tends to stabilize the convectivemotion. This stabilization effect is strong in the vicinityof the main body of the disk with the highest density. Asa result, the matter of the high specific entropy producedin the innermost region goes along a high-latitude regionof the disk. Nevertheless, the high-entropy convectiveblob eventually brings the thermal energy into the outerpart of the disk with a large cylindrical radius and witha high latitude. By the increase of the thermal energythere, the matter in the outer part of the disk obtainsthe energy enough to escape from the system as ejecta.Therefore the onset time of the mass ejection is deter-mined approximately by the viscous heating timescalefor the innermost region of the disk at the moment thatthe neutrino emission timescale becomes as long as theviscous heating timescale, and the major process of en-ergy transport is the convection (cf. Ref. [46]). By theconvective activity, the matter is ejected to any directionexcept along the rotation axis (see Sec. III D).We note that the viscous heating/angular momentumtransport in the high-density region of the disk also con-tribute to the expansion of the entire region of the disk.Thus not only the convective motion but also the contin-uous viscous effect and resulting disk expansion play arole for the mass ejection. This is in particular the casefor the high viscous coefficients as indicated in Sec. III D.In addition, the thermal energy stored in the innermostregion by the viscous heating drives an intermittent wind.This contributes to non-steady mass ejection in particu-lar for the late time evolution of the disk in which the2 (a) − − t (sec)10 R m a t ( k m ) K8 ( α ν = 0 . α ν = 0 . α ν = 0 . Y e = 0 . χ = 0 . (b) − − t (sec)681012141618 h s / k i K8 ( α ν = 0 . α ν = 0 . α ν = 0 . Y e = 0 . χ = 0 . (c) − − t (sec)0 . . . . . . . h Y e i K8 ( α ν = 0 . α ν = 0 . α ν = 0 . Y e = 0 . χ = 0 . (d) − − t (sec)10 L ν ( e r g/ s ) K8 ( α ν = 0 . α ν = 0 . α ν = 0 . Y e = 0 . χ = 0 . (e) − − t (sec)0 . . . . . . L ν / ( c d M f a ll / d t ) K8 ( α ν = 0 . α ν = 0 . α ν = 0 . Y e = 0 . χ = 0 . (f) − − t (sec)0 . . . . . . . M e j e ( M (cid:12) ) K8 ( α ν = 0 . α ν = 0 . α ν = 0 . Y e = 0 . χ = 0 . FIG. 7. Several quantities for models K8, K8h, K8n, Y8, and K6. Average cylindrical radius (top left), average specific entropy(top right), and average value of Y e (middle left) of the matter located outside the black hole as functions of time. The middleright panel shows the time evolution of the total neutrino luminosity, L ν . The bottom two panels show an efficiency of theneutrino emission defined by the total neutrino luminosity, L ν , divided by the rest-mass energy accretion rate of the matterinto the black hole, c dM fall /dt (left) and total mass of the ejecta component M eje (bottom right) as functions of time. Formodel K8, the results with a higher resolution run are also plotted and show their weak dependence on the grid resolution. neutrino cooling plays a negligible role and viscous heat-ing can be fully available for driving the wind.To examine the significance of the neutrino heating (ir-radiation) effect, we also performed a simulation withoutneutrino heating (model K8n). Our result shows that theneutrino heating does not play a substantial role for themass ejection; the mass ejection rate and ejecta veloc-ity are not substantially influenced by the neutrino heat-ing. The reason for this is that the neutrino luminosityis not very high during the mass ejection stage in our model (see Fig. 7(d)): Only in the early stage of the diskevolution with t (cid:46)
20 ms, the total luminosity exceeds10 erg / s, while for t >
100 ms for which the mass ejec-tion becomes active, the luminosity drops to < erg / sexponentially, resulting in the small contribution to themass ejection by the neutrino heating. This result is inbroad agreement with that of Ref. [14]. We note that the maximum energy obtained from neutrinos per t ∼ . α ν = 0 .
05 with/without neutrino irradiation), theviscous heating and resulting convection continue to playa leading role for the mass ejection, although the massejection rate gradually decreases with the decrease of thedisk density. The total ejecta mass is 15–30% of theinitial disk mass (see Fig. 7 and Table II). For a fixedviscous coefficient with α ν = 0 .
05, the ejecta mass is ≈ Y e distribution, density and velocity profiles, compactness ofthe disk, and black-hole spin). The fraction of the ejectamass agrees broadly with the results of the earlier viscoushydrodynamics work [10, 11, 14].The entropy and electron fraction in the disk increasewith its viscous and convective expansion, after the initialquick matter accretion onto the black hole for t (cid:46)
100 msceases (see Fig. 7(b) and (c)). The typical average valueof the specific entropy is 10–12 k when the mass ejection isactivated. After the initial infall stage of the disk matteronto the black hole, the average value of Y e monotoni-cally increases with the decrease of the disk density, andit is higher than ∼ . α ν = 0 .
05. In thelate time, the average value of Y e settles to constants of ∼ .
3. The mechanism for this Y e evolution is summa-rized as follows: During the disk expansion ongoing untilthe weak-interaction freezes out (i.e., the temperature is kT (cid:38) Y e is approximately determined bythe equality of the rates of electron/positron capture onnucleons: see Appendix A for more details) However, for kT (cid:46) τ β , becomes longer than the viscous timescale, τ vis [10](here τ β is approximately equal to the neutrino emissiontimescale). Thus the settled value of Y e is determinedapproximately by the condition of τ β = τ vis . The on-set time of this freeze out of the weak interaction agrees nucleon (neutron or proton) via the neutrino absorption is es-timated approximately as ∆ E ν σ ν / (4 πr ) where σ ν is the crosssection of nucleons with neutrinos and ∆ E ν is the total energyemitted by each neutrino species (here electron neutrino or antineutrino). In model K8, for each, ∆ E ν ∼ × erg / s × .
02 s ∼ erg, σ ν ∼ × − cm for neutrinos of the energy of ∼
10 MeV, and 4 πr (cid:38) cm because r (cid:38) GM BH /c ∼ ∼ GM BH /c r ) MeV pernucleon. In reality the energy gain of each nucleon is muchsmaller than this because the cooling by the neutrino emissioncould be of the same order as the heating. However, even ifwe ignore the neutrino cooling, the obtained energy gain is atmost as large as the gravitational potential energy of nucleons, ∼ GM BH m n /r , where m n c is the rest-mass energy of nucleons ∼
940 MeV. In particular, for r ≥ GM BH /c , the estimated en-ergy gain is by more than one order of magnitude smaller thanthe gravitational potential energy. approximately with the time that the mass ejection isactivated. In our numerical result, the relaxed averagevalues of Y e are not very small, i.e., (cid:38) .
25, irrespectiveof the models for α ν = 0 .
05. This has an important im-plication for the nucleosynthesis of lanthanide elementsas we discuss in Sec. III E.Figure 7(e) shows that the total neutrino luminosityfor χ = 0 . ∼ χ = 0 . χ = 0 .
8. Thereason for this is that the radius of the innermost sta-ble circular orbit is larger for the smaller value of χ , andthus, the depth of the gravitational potential in the vicin-ity of the black hole is shallower; i.e., the gravitationalpotential energy available for the dissipation is smaller.For higher viscous coefficients, the efficiency is higher,but the enhancement is not very remarkable. Thus, themaximum efficiency is basically determined by the spinof the black hole. In the late stage of the disk evolutionin which the weak interaction freezes out, the efficiencyapproaches zero because the emissivity of neutrinos expo-nentially drops. As already mentioned, in this late stage,the mass ejection is enhanced.Figure 5 displays the evolution of the profiles for modelK8h, i.e., for a larger viscosity model than model K8.For α ν = 0 .
10 and 0 .
15 (K8h and K8s models), the massejection sets in earlier than for model K8. This is clearlyfound by comparing the profiles at t = 0 . α ν = 0 .
10 and 0.15 are more enhancedthan for α ν = 0 .
05, in particular in the innermost regionof the disk.Associated with the stronger viscous heating, the spe-cific entropy and the average value of Y e in the disk in-crease more rapidly (see Fig. 7(b) and (c),early phases).The reason for this is that the temperature of the disk ishigher and the electron degeneracy is weakened at giventime for the higher viscosity models. As we mentionedabove, however, the disk expansion occurs more quicklyfor the larger viscous coefficients. As a result, the weakinteraction freezes out earlier, and also, the electron frac-tion settles to smaller values for the larger viscous coef-ficients (compare the panels in Figs. 4 and 5 at t = 2 sand also see Fig 7(c)). This modifies the Y e distributionof the ejecta as discussed in Sec. III D.4The larger viscous effect slightly suppresses the matterinfall into the black hole (see Fig. 6(b)). Associated withthis, the total ejecta mass for α ν = 0 .
10 and 0.15 becomeslarger than that for α ν = 0 .
05 by 0 . M (cid:12) and 0 . M (cid:12) ,respectively (see Fig. 7(e) and (f) as well as Table II).All these results show that the magnitude of the viscouscoefficient significantly influences the disk evolution andthe quantitative properties of the ejecta.Figure 7 clearly shows that the effect of the black-holespin is not substantial on the evolution of the system forthe astrophysically plausible values of χ as the mergerremnants (compare the results of models K8 and K6).We can find that the average value of Y e and neutrinoluminosity for model K6 are slightly smaller than thosefor model K8, because the innermost edge of the disk formodel K6 is located slightly far outside the black holethan for model K8, and as a result, the viscous heatingefficiency and resulting ejecta mass are slightly smaller.However, this spin effect is minor for modifying the evolu-tion of the system: The effect in the change of the viscouscoefficient by a factor of 2 has much stronger impact. C. Viscous hydrodynamics of disks: comparisonamong different disk mass models
The difference in the disk mass (compare models K8,C8, and E8) results in the modification of the timescaleto reach the freeze out of the weak interaction and onsettime of the mass ejection. Figure 8 displays the snapshotsof the profiles for the rest-mass density, temperature, spe-cific entropy per baryon, and electron fraction for modelsC8 at selected time slices. Figure 9 also compares theaverage cylindrical radius, average specific entropy, av-erage value of Y e for the matter outside the black hole,total neutrino luminosity, an efficiency of the neutrinoemission, L ν / ( dM fall /dt ), and ejecta mass as well as thetotal mass swallowed by the black hole among modelsK8, C8, and E8. These figures show that for the largerdisk mass, the density and temperature is always higherat given time and the freeze-out time of the weak interac-tion comes later. As a result, for the more massive mod-els, the thermal energy generated by the viscous heatingis released by the neutrino emission for a longer term, re-sulting in more luminous and long-lasting neutrino emis-sion. Then, the disk expansion timescale becomes longer(see Fig. 9(a)), the onset time of the convection becomeslater, and the mass ejection is delayed [14]. These factsare well reflected in the maximum efficiency of the neu-trino emission and the duration of the high-luminosityneutrino emission (see Fig. 9(c) and (e)). Due to the de-lay of the freeze out of the weak interaction, the averagevalue of Y e in the disk at the freeze out of the weak in-teraction becomes higher for more massive disk models(see Fig. 9(c)). Figure 8 also illustrates that besides thequantitative difference mentioned above, the evolutionprocess of the disk is qualitatively similar among differ-ent disk-mass models; e.g., the mass ejection is primarily TABLE II. Rest mass, average value of velocity, and averagevalue of Y e for ejecta. All the quantities are extracted at t = 3 . t = 2 . M ej ( M (cid:12) ) V ave /c Y e, ave X lan K8 0.020 0.06 0.31 2 . × − K8 (high res.) 0.020 0.06 0.31K8h 0.023 0.08 0.28 1 . × − K8s 0.027 0.08 0.26 8 . × − K8n 0.021 0.07 0.30D8 0.015 0.06 0.31Y8 0.023 0.06 0.32K6 0.019 0.06 0.30C8 0.092 0.06 0.32 1 . × − E8 0.006 0.07 0.31 1 . × − driven by the convective activity caused by the viscousheating in the region of the disk close to the black holefor all the models.For the remnant of neutron-star mergers, disks withmass of (cid:38) . M (cid:12) are not very likely. However, such aheavy disk surrounding a spinning black hole could beformed as a result of rotating massive stellar core col-lapses, in particular as a central engine of long gamma-ray bursts [47]. Such heavy disks are also likely to beevolved effectively by viscous-hydrodynamics processes.Our present numerical results demonstrate that for largedisk mass, the electron fraction at which the weak inter-action freezes out could be large with Y e (cid:38) . D. Properties of ejecta
Figure 10 shows the average value of Y e and averagevelocity of the ejecta as functions of time for all the mod-els (except for model D8) employed in this paper. Here,we plot the ejecta velocity determined only for the ejectacomponent that escapes from the sphere of r = 2000 km.We note that even if the extraction radius is changedfrom r = 2000 km to 3000 km and 4000 km, the asymp-totic values of Y e and the velocity are varied only within3%. Table II lists the rest mass, average velocity, andaverage value of Y e of the ejecta for all the models con-sidered in this paper. We also show the mass histogram(mass fraction) of the ejecta as a function of Y e for mod-els K8, K8h, K8s, K8n, Y8, K6, C8, and E8 in Fig. 11.Here, the value of Y e is determined at the time when thetemperature of the corresponding ejecta component de-5 FIG. 8. The same as Fig 4 but for model C8. For the last two panels at t = 4 and 6 s, the plots are made for an enlarged regionof 5000 km × creases to 5 × K (referred to as 5 GK: note that thetemperature of the ejecta component in general decreasesdown to ∼ Y e for the ejecta with α ν = 0 . ∼ . Y e .On the other hand, it decreases significantly with the increase of the viscous coefficient. As Fig. 11(a) showsthat the distribution shifts to the lower side of Y e withthe increase of α ν . This result can be expected fromFig. 7(c), which shows that for the matter outside theblack hole, the asymptotic average values of Y e is ∼ . α ν = 0 .
05, 0.10, and 0.15, respectively.These values agree approximately with the average valueof Y e for the ejecta. This agreement indicates that the6 (a) − − t (sec)10 R m a t ( k m ) C8 ( M disk = 0 . M (cid:12) )K8 ( M disk = 0 . M (cid:12) )E8 ( M disk = 0 . M (cid:12) ) (b) − − t (sec)6810121416 h s / k i C8 ( M disk = 0 . M (cid:12) )K8 ( M disk = 0 . M (cid:12) )E8 ( M disk = 0 . M (cid:12) ) (c) − − t (sec)0 . . . . . . h Y e i C8 ( M disk = 0 . M (cid:12) )K8 ( M disk = 0 . M (cid:12) )E8 ( M disk = 0 . M (cid:12) ) (d) − − t (sec)10 L ν ( e r g/ s ) C8 ( M disk = 0 . M (cid:12) )K8 ( M disk = 0 . M (cid:12) )E8 ( M disk = 0 . M (cid:12) ) (e) − − t (sec)0 . . . . . . L ν / ( c d M f a ll / d t ) C8 ( M disk = 0 . M (cid:12) )K8 ( M disk = 0 . M (cid:12) )E8 ( M disk = 0 . M (cid:12) ) (f) − − − t (sec)10 − − M e j e a nd M f a ll ( M (cid:12) ) C8 ( M disk = 0 . M (cid:12) ),ejefallK8 ( M disk = 0 . M (cid:12) )E8 ( M disk = 0 . M (cid:12) ) FIG. 9. The same as Fig. 7 but for comparison among different disk mass models C8, K8, and E8. In the bottom right panel,not only the total ejecta mass (solid curves) but also the total mass that falls into the black hole (dashed curves) are plotted. matter in the outer part of the disk that expands bythe viscous effect and the resulting convective heatingeventually becomes the ejecta component. The averagevalue of Y e for the ejecta increases with the increase ofthe disk mass (compare the results for models C8, K8,and E8). This correlation also agrees with that found inFig. 9(c).The asymptotic average velocity of the ejecta, v eje , is ∼ . c for α ν = 0 .
05 and again depends very weaklyon the initial condition for the disk profiles of density,angular velocity, and Y e , and the black-hole spin. Forhigher viscous coefficients, v eje is increased, and for α ν =0 .
10 and 0.15, it is ≈ . c . This reflects an efficientacceleration of the matter in the outer part of the disk for the higher viscous coefficients.Although the results for the mass and velocity of theejecta agree broadly with those of the earlier numericalstudies by other groups [10, 11, 14], our result for themass distribution of Y e shows a noticeable difference fromthe previous results. In the previous results, a substan-tial fraction of the ejecta has rather small values of Y e between 0 . Y e components are rare, and the peak is locatedapproximately between Y e = 0 .
25 and 0.30. We note thatthe low Y e component is present in the disk for the earlystage in which the disk density is high, (cid:38) g / cm (see,e.g., Fig. 4 at t = 0 . Y e is increased and the ejecta does7 (a) t (sec)0 . . . . . . . . Y e K8 ( α ν = 0 . α ν = 0 . α ν = 0 . Y e = 0 . χ = 0 . (b) t (sec)0 . . . . . . . v e j e / c K8 ( α ν = 0 . α ν = 0 . α ν = 0 . Y e = 0 . χ = 0 . (c) t (sec)0 . . . . . . . Y e C8 ( M disk = 0 . M (cid:12) )K8 ( M disk = 0 . M (cid:12) )E8 ( M disk = 0 . M (cid:12) ) (d) t (sec)0 . . . . . . . v e j e / c C8 ( M disk = 0 . M (cid:12) )K8 ( M disk = 0 . M (cid:12) )E8 ( M disk = 0 . M (cid:12) ) FIG. 10. Average value of Y e (left) and and average velocity (right) of the ejecta as functions of time for models with M disk = 0 . M (cid:12) (top) and for models C8, K8, and E8 (bottom). Here, the average velocity is determined only for the ejectacomponent that escapes from r = 2000 km. not have the component with Y e ≥ .
2. One possiblereason for the difference of our results from the previ-ous ones is that the equation of state and/or treatmentof the weak interaction (e.g., electron and positron cap-ture) that we employ may be different from those in theprevious studies. Another reason is that in our simula-tion we do not have any mechanism for the mass ejectionwhich is more efficient than the viscosity driven ejection.In the following, we describe these points in more detail.In our simulation, the mass ejection occurs only af-ter the viscous heating and angular momentum trans-port proceed at least for several hundreds milliseconds,which causes the matter in the disk to expand to r (cid:38) km. We do not find any other major components formass ejection like the neutrino-wind component, which isfound in Ref. [14]. As we already mentioned in Sec. III B,in this viscous evolution process, associated with the de-crease of the density and with the decrease of the degreeof electron degeneracy, the value of Y e in the disk gradu-ally increases from a low value of Y e ∼ . Y e (cid:38) . kT ∼ Y e is de-termined approximately by the condition that electronand positron capture rates are identical. After the freezeout of the weak interaction, the average of Y e approaches asymptotically ∼ . α ν = 0 . Y e among differentgroups.For the late stage with lower temperature of kT (cid:46) Y e in the disk. Since the Y e distribution of the ejecta isjust the reflection of that in the disk, the value of Y e forthe ejecta thus determined is not very small but ratherlarge as (cid:38) .
25 (for α ν = 0 . Y e . Thisindicates that the low Y e ejecta with Y e ∼ . α ν , the weak interaction freezesout earlier (much earlier than 1 s) because faster expan-sion occurs for the disk. Since the duration to increasethe value of Y e is shorter, Y e becomes smaller for larger8 (a) . . . . . . Y e ( T = 5 GK)10 − − − − − ∆ M / M tr a j K8 ( α ν = 0 . α ν = 0 . α ν = 0 . (b) . . . . . . Y e ( T = 5 GK)10 − − − − − ∆ M / M tr a j K8 ( Y e =0.07-0.5)Y8 ( Y e =0.1) (c) . . . . . . Y e ( T = 5 GK)10 − − − − − ∆ M / M tr a j K8 ( χ = 0 . χ = 0 . (d) . . . . . . Y e ( T = 5 GK)10 − − − − − ∆ M / M tr a j K8K8n (irradiation off) (e) . . . . . . Y e ( T = 5 GK)10 − − − − − ∆ M / M tr a j C8 ( M disk = 0 . M (cid:12) )K8 ( M disk = 0 . M (cid:12) )E8 ( M disk = 0 . M (cid:12) ) (f) . . . . . . Y e ( T = 5 GK)10 − − − − − ∆ M / M tr a j K8K8 (high res)
FIG. 11. Mass histogram (mass fraction) of the ejecta as a function of Y e for models (a) K8, K8h, and K8s, (b) K8 and Y8,(c) K8 and K6, (d) K8 and K8n, (e) C8, K8, and E8, and (f) K8 with two grid resolutions. Here, the value of Y e is determinedat the time when the temperature of each ejecta component decreases to 5 × K (referred to as 5 GK). values of α ν . This tendency is clearly found by compar-ing the mass histogram as a function of Y e for models K8,K8h, and K8s (see Fig. 11(a)). This suggests that in thepresence of a mass ejection process more efficient than theviscosity-driven mechanism, like a magnetohydrodynam-ics mechanism, which would be primarily not the MRIbut the magnetic winding and resulting strong Lorentzforce with a hypothetically very large poloidal magneticfield aligned with the black-hole spin direction [18], thematter with even lower Y e components may be ejected.However, our present simulations show that in the ab-sence of such an efficient mass ejection process with theejection timescale much shorter than 1 s (i.e., only withthe viscous process), the value of Y e for the ejecta can be quite large as Y e (cid:38) . Y e ejecta components with Y e ≤ . Y e is quite similar between these two models.This indicates that the Y e distribution of the ejecta de-pends only weakly on its initial condition, if we focus onlyon the viscosity-driven ejecta. We note, however, that if9 (a) π /4 π /2 θ (rad)0 . . . . . . . ∆ M / M tr a j K8 ( α ν = 0 . α ν = 0 . α ν = 0 . (b) π /4 π /2 θ (rad)0 . . . . . . . ∆ M / M tr a j C8 ( M disk = 0 . M (cid:12) )K8 ( M disk = 0 . M (cid:12) )E8 ( M disk = 0 . M (cid:12) ) FIG. 12. Mass fraction as a function of the polar direction of the mass ejection (a) for models K8, K8h, and K8s and (b) formodels C8, K8, and E8. a very efficient mechanism of mass ejection is present inthe very early stage of the disk evolution for t (cid:28)
100 ms,the resulting mass histogram as a function of Y e may re-flect the difference of the initial condition. This pointshould be explored in the simulation community.Figure 11(c) compares the mass histogram for mod-els K8 and K6, i.e., for different black-hole spin mod-els. This figure shows that for the smaller black-holespin (K6), the Y e distribution shifts slightly to a lower Y e side. Our naive interpretation for this is that for thesmaller spin, the disk mass slightly becomes smaller be-cause a larger fraction of the disk matter falls into theblack hole, resulting in slight shortening until the freezeout of the weak interaction is reached. Alternative possi-bility is that a difference in the initial profile of the diskis reflected. However, the difference in the Y e distribu-tion is quite small, and hence, we may conclude that thedifference of the spin is not very important for the prop-erties of the viscosity-driven ejecta, as far as we focus onthe astrophysically plausible values of χ .Figure 11(d) compares the mass histogram for mod-els K8 and K8n to understand the unimportance of theneutrino irradiation for the ejecta in the black hole-disksystem. This figure shows that two results are quite simi-lar each other, and hence, the neutrino irradiation indeeddoes not play an important role. This is quite reasonablein our viscous evolution models because the mass ejec-tion is activated only after the neutrino luminosity drops(i.e., after the weak interaction freezes out), whereas theneutrino irradiation could be important only for the casethat the neutrino luminosity is still high. Thus, the neu-trino irradiation would be important only for the casethat mass ejection occurs in an early stage in which theneutrino luminosity is still high (e.g., Ref. [19]).Figure 11(e) compares the mass histogram for modelsC8, K8, and E8; different initial disk mass models. It isfound that for the higher disk mass, the lowest end of Y e is larger. As already mentioned in the previous section,for the higher disk mass, the time at which the weakinteraction freezes out comes later, and the value of Y e for the disk component becomes higher. Associated withthis effect, the lowest value of Y e in the ejecta componentalso increases for the larger disk mass model. This resultsuggests that for the case that the remnant disk massis large, e.g., in the remnant of stellar core collapse to aspinning black hole [47], the ejecta from the disk may notbe very neutron rich and cannot synthesize a substantialfraction of heavy elements like lanthanide (see Sec. III E).Figure 12 displays the mass fraction as a function of thepolar angle of the mass ejection (a) for models K8, K8h,and K8s and (b) for models C8, K8, and E8. This showsthat the mass ejection occurs to a wide range of anglesexcept for the direction of the rotation axis irrespectiveof the viscous coefficient and disk mass. The absence ofthe mass ejection toward the rotational axis is naturalbecause the ejecta, which are driven in a region far fromthe black hole, have a substantial angular momentum.An interesting point is that the primary direction of themass ejection depends on the magnitude of the viscouscoefficient. For α ν = 0 .
05, the mass ejection occurs moststrongly toward the direction of θ ∼ π/ θ de-notes the polar angle. However, with the increase of α ν ,the primary angle of the mass ejection increases, and formodel K8s, the mass ejection occurs primarily to the di-rection of the equatorial plane, θ ∼ π/
2. This indicatesthat not only the convective activity but also the out-ward angular momentum transport plays an importantrole for the mass ejection with larger values of α ν . E. Nucleosynthesis in the ejecta
A nucleosynthesis calculation is performed for modelsK8, K8h, K8s, E8, and C8. For each model, a few tenthousands of tracer particles are arranged by the methoddescribed in Ref. [21], and in each tracer particle, thenucleosynthesis is explored as a post-processing step byusing the reaction network code rNET [49]. The reactionnetwork consists of 6300 species from single neutron andproton up to the isotopes with Z = 110 (see Ref. [49]0 (a) − − − − − − − − l og X A K8 ( α ν = 0 . α ν = 0 . α ν = 0 . (b) − − − − − − − − l og X A C8 ( M disk = 0 . M (cid:12) )K8 ( M disk = 0 . M (cid:12) )E8 ( M disk = 0 . M (cid:12) )solar r-residual (Goriely 1999) FIG. 13. Patterns of mass fraction obtained in the nucleosynthesis calculation for the ejecta (a) for models K8, K8h, and K8sand (b) for models C8, K8, and E8. The filled circles with the error bar denote the r -process residual to the solar systemabundance for A ≥
69 [53]. The data are shifted to match the calculated mass fraction for model K8 at A = 83. for more details). For those relevant to the r -process,the rates of both neutron capture (TALYS [51]) and β -decay (GT2 [50]) are based on the microscopic predic-tion of nuclear masses, HFB-21 [52]. In this work, weexclude fission process from the network, which is rele-vant only for Y e < .
15. Neutrino-induced reactions arenot included, either, which are unimportant as describedin previous sections. Each nucleosynthesis calculation isstarted when the temperature decreases to 10 K withthe initial compositions of free neutron and proton to be1 − Y e and Y e , respectively. Note that nuclear statisticalequilibrium (NSE) immediately establishes at such hightemperature.Figure 13 displays the mass fraction of each element asa function of the atomic mass number, A . In our ejectadata, the peak in Y e comes around Y e ∼ . . Y e (cid:46) . Y e , the heavy elements with A (cid:38)
130 are not synthesizedsubstantially. By contrast, elements with A ∼ A ∼ r -process elements [53]. This suggests that low-mass blackhole-disk systems might be a subdominant site for the r -process nucleosynthesis, assuming that the solar abun-dance pattern is universal [54] in the universe and thatthe viscosity-driven mechanism is the main channel ofthe mass ejection.For model K8 (as well as model C8), the mass fractionof lanthanide elements is quite small (see the last columnof Table II) as expected from the mass histogram shownin Fig. 11. By contrast, for the higher viscous coefficientcases (models K8h and K8s), a fraction of the lanthanideelements is synthesized, although the amount of the lan-thanide and heavier elements are still smaller than thatin the solar abundance. The reason for this dependenceis that for the higher viscous coefficients, the ejecta con-tains a fraction of neutron-rich components with Y e (cid:46) . (cid:38) ∼ times higher than that forlighter elements like irons [55–57]. Our present resultssuggest that even for the ejecta from black hole-disk sys-tem, the lanthanide fraction could be minor, in contrastto the previous understanding (e.g., Refs. [10, 11, 14]), ifthe mass ejection sets in for (cid:38) . ∼ ∼ A >
132 like lanthanide is synthesized (see also Ta-ble II). This result suggests that if the viscous process isthe dominant mechanism of mass ejection and if the vis-cous coefficient is not extremely large, the ejecta from themassive disk around black holes may not be the sourcefor the r -process nucleosynthesis of the heavy elementslike lanthanide and third peak elements (e.g., gold). Ref-erence [61] has recently illustrated that if the mass ejec-tion sets in within ∼
100 ms after the merger, in whichthe neutrino luminosity is still high, the neutron richnessof the ejecta is significantly decreased by the neutrinoirradiation. For the massive disks (Reference [61] consid-ered the disk of mass (cid:28) . M (cid:12) ), the neutrino luminosityshould be higher as our present work shows, and hence,the neutrino irradiation would be even more enhanced.Thus, if the mass ejection occurs in late time (only inthe presence of viscous mass ejection), the weak inter-action process in the disk enhances the value of Y e andif the mass ejection occurs in early time (by some pow-erful process; e.g., by a poloidal magnetic field alignedwith the black-hole spin direction), the neutrino irradi-ation process enhances the value of Y e . This suggeststhat for the massive disk around a black hole to be a sitefor the nucleosynthesis of the heavy r -process elements,a fine tuning for the timing of the mass ejection wouldbe necessary. F. Convergence on the grid resolution
Before closing Sec. III, we comment on the convergenceof our numerical results with respect to the grid resolu-tion. As found from Fig. 6(a), the accuracy for follow-ing the evolution of rapidly spinning black holes dependsstrongly on the grid resolution. Due to the truncationerror, the black-hole spin spuriously decreases and thedecrease rate can be too high for insufficient grid reso-lutions to be acceptable. For ∆ x = 0 . M BH of modelK8, the dimensionless spin decreases by ≈ .
06 in 1 s.Such a large error is not acceptable. However, the de-crease rate is suppressed significantly with improving thegrid resolution. For ∆ x = 0 . M BH and 0 . M BH ,the decrease rates are ∼ . / s and 0.01/s, respectively,so that the effect of the rapid spin for these simulations can be taken into account for t (cid:46) x = 0 . M BH is acceptable for a reliablenumerical simulation.Figure 11(f) compares the mass histogram as a func-tion of Y e for model K8 with two different grid resolu-tions. We find that two histograms agree broadly witheach other, although the mass fraction at each value of Y e does not exactly agree. IV. SUMMARY
This paper presents our first numerical results for aviscous neutrino-radiation hydrodynamics simulation ofaccretion disks surrounding a spinning black hole in fullgeneral relativity as models for the evolution of a mergerremnant of massive binary neutron stars or low-massblack hole-neutron star binaries. We reconfirm the fol-lowing results found by previous viscous hydrodynamicsstudies by other groups [10, 14]: About 15–30% of thedisk mass is ejected from the system with the average ve-locity of ∼ ∼ (cid:46) . ∼ . ACKNOWLEDGMENTS
We thank Kyohei Kawaguchi and MasaomiTanaka for helpful discussions. This work wasin part supported by Grant-in-Aid for ScientificResearch (Grant Nos. JP16H02183, JP16H06342,JP17H01131, JP18H01213, JP18H04595, JP18H05236,JP19K14720, JP16K17706, JP16H06341, JP17H06363,and JP18K03642) of Japanese MEXT/JSPS. Wealso thank the participants of a workshop entitled“Nucleosynthesis and electromagnetic counterparts ofneutron-star mergers” at Yukawa Institute for Theoret-ical Physics, Kyoto University (No. YITP-T-18-06) formany useful discussions. Numerical computations wereperformed on Sakura and Cobra clusters at Max PlanckComputing and Data Facility, on Oakforest-PACS atInformation Technology Center of the University ofTokyo, on XC50 at National Astronomical Observa-tory of Japan, and on XC40 at Yukawa Institute forTheoretical Physics.
Appendix A: Weak interaction in the disk
In this Appendix, we describe how the electron frac-tion in the disk is determined by the weak interactionprocesses. Here, we consider the electron/positron cap-ture on nucleons and nuclei as the major reaction thatdetermines the electron fraction. We suppose that theneutrino absorption process plays only a minor role forthe evolution of black hole-disk systems, in particular for3
FIG. 14. Distribution of disk matter in the ρ - T plane for model K8 at t = 0 (upper left), 0.03 s (upper right), 0.1 s (lower left),and 0.3 s (lower right). Each point shows the density and temperature at a certain grid point in the simulation result. Y e ateach grid point is displayed in color points. The black curves show the contours of Y e, eq . Behind the ρ - T distribution of thesimulation, we also plot (by shaded color) the weak interaction timescale at each bin, t eq , which is defined as the inverse ofEq. (A1) with Y e = Y e, eq . In the plane, there is a region for which Y e, eq should be higher than 0.6 [65]. We cannot calculatethe electron/positron capture rates in such a region because Y e > . Y e = 0 . the late stage of the disk evolution, because its luminos-ity is too low to significantly change the electron fractionin the disk.The reaction rates of the electron and positron captureon free nucleons and heavy nuclei are written as (see,e.g., Ref. [63]) R ec = X p λ fec + X h (cid:104) A (cid:105) h λ hec , (A1) R pc = X n λ fpc + X h (cid:104) A (cid:105) h λ hpc , (A2)where X n , X p , and X h are the mass fractions of neutron,proton, and heavy nuclei, (cid:104) A (cid:105) h is the average mass num-ber of heavy nuclei, and λ f / hec / pc is the rate of each reactionwhich has the unit of (time) − . Here, the superscripts,“f” and “h”, indicate the capture on free nucleons andheavy nuclei, respectively. The reaction rates are written by the integration with respect to the energy as [65] λ f / hec / pc = ln 2 (cid:104) f t (cid:105) f / h m e c (cid:90) ∞ dω ω ( ω − Q f / hec / pc ) × (cid:118)(cid:117)(cid:117)(cid:116) − m e c ( ω − Q f / hec / pc ) F e ( ω − Q f / hec / pc ) × Θ( ω − Q fec / pc − m e c ) , (A3)where (cid:104) f t (cid:105) is the so-called f t -value of these reactions, F e is the distribution function of electrons, which is assumedto be the Fermi-Dirac form, and Q is the Q-value of thesereactions. For the capture processes by free nucleons, (cid:104) f t (cid:105) fec / pc ≈ Q fec = ( m p − m n ) c , and Q fpc = ( m n − m p ) c . For the capture processes by heavy nuclei, wefollow the approximations in Ref. [63] for f t -values and4Q-values with Q hec = µ p − µ n ,Q hpc = µ n − µ p , (A4)and log (cid:16) (cid:104) f t (cid:105) hec / pc / s (cid:17) = . µ e < Q ec / pc )2 . µ e > Q ec / pc )2 . . /T (blocked) , (A5)where T = T / K. µ n and µ p denote the relativistic(including the mass) chemical potential of neutrons andprotons, respectively. Here, “blocked” and “unblocked”cases imply that (cid:104) N (cid:105) ≥
40 or (cid:104) Z (cid:105) ≤
20 and that (cid:104) N (cid:105) < (cid:104) Z (cid:105) >
20, respectively. (cid:104) N (cid:105) and (cid:104) Z (cid:105) denote the av-erage neutron and proton numbers of the heavy nuclei,respectively. Note that we do not consider the updatesfor the electron/positron capture on heavy nuclei in [64],which play a minor role in the present case owing to thefreeze out of weak interaction at sufficiently high temper-ature (see below).The reaction rates are functions of ρ , T , and Y e . Wecan derive the electron fraction in the equilibrium, Y e, eq ,at each density and temperature by equating the electronand positron capture rates as R fec ( ρ, T, Y e, eq ) + R hec ( ρ, T, Y e, eq )= R fpc ( ρ, T, Y e, eq ) + R hpc ( ρ, T, Y e, eq ) . (A6)Here we do not consider the blocking by neutrinos, be-cause it is important only for the case that the opticaldepth to neutrinos is large. Note that the disk materialis supposed to be optically thin to neutrinos except forthe very early epoch of the disk evolution (i.e., t (cid:46) . Y e, eq (black curves) to-gether with the distribution of Y e for the disk material(color plots) in the ρ – T plane at selected time slices, t = 0, 0.03, 0.1, and 0.3 s of the simulation for modelK8. Here, we employ the DD2 equation of state to de-termine the mass fractions, average mass number, andchemical potentials in Eqs. (A1)–(A3) as functions of ρ , T , and Y e . Due to the viscous heating/angular momen-tum transport and resulting expansion, the density andtemperature in the disk decrease. On the other hand,the specific entropy increases. Thus, the distribution ofthe disk matter moves basically to the bottom left regiongradually. Figure 14 also displays the weak interactiontimescale at each point, t eq , which is defined as the in-verse of Eq. (A1) with Y e = Y e, eq .In an early stage of the disk evolution, the condition of t > t eq is satisfied for the majority of the disk matter sothat the values of Y e approach Y e, eq . For t (cid:46) . (cid:38) g / cm ,and hence, until t = 0 . Y e < .
1. However, because of the expansion of the disk,the density goes below ∼ g / cm and then the valuesof Y e have to be larger than ∼ . t = 0 . Y e material from thedisk in the presence of the mass ejection mechanism witha short timescale of (cid:46) . t (cid:46) . Y e material may be ejected. However, for thepurely viscous mass ejection with the timescale (cid:38) . Y e (cid:46) .
2. Thisagrees totally with the results of our numerical simula-tions.The dependence of t eq on the temperature shows thatfor kT (cid:38) t eq (cid:46) . kT (cid:38) . t eq (cid:46) ∼ k , the neutrino cooling timescale, which has thesame order of magnitude as t eq in the region where thebaryon dominates the internal energy, becomes longerthan the viscous timescale, and as a result, the freezeout of the weak interaction occurs. This effect is alsoobserved in our numerical results. Appendix B: Dependence on the disk compactness
To show the dependence of the viscous evolution onthe initial compactness of the disk, we here compare theresults for models K8 and D8. Figure 15 displays theevolution of the same quantities for the matter locatedoutside the black hole as in Fig. 9 bur for models K8and D8. It is found that the values of M fall and M eje formodel D8 are larger and smaller than those for model K8,respectively. For model D8, the ejecta mass is ≈
15% ofthe initial disk mass, which is appreciably smaller thanthat for models with less compact disk like K8 and Y8 forwhich the ejecta mass is (cid:38)
20% of the initial disk mass(see Table II).Because the larger fraction of the disk mass falls intothe black hole in the early stage of the viscous evolution,the disk mass for model D8 becomes smaller than formodel K8 during the long-term viscous evolution. Dueto this reason, the timescale to reach the freeze out of theweak interaction for model D8 is slightly shorter than formodel K8 (see Sec. III C). However, besides this small dif-ference, the evolution process for two models is quite sim-ilar each other. The final average value of Y e for modelD8 is only slightly smaller than for model K8. As a re-sult, the average value of Y e and the mass histogram forthe ejecta is not very different between two models. Theaverage velocity of the ejecta is also approximately thesame for two models (see Table II).5 (a) − − t (sec)10 R m a t ( k m ) K8 ( r out = 41 M BH )D8 ( r out = 29 M BH ) (b) − − t (sec)6810121416 h s / k i K8 ( r out = 41 M BH )D8 ( r out = 29 M BH ) (c) − − t (sec)0 . . . . . . h Y e i K8 ( r out = 41 M BH )D8 ( r out = 29 M BH ) (d) − − t (sec)10 L ν ( e r g/ s ) K8 ( r out = 41 M BH )D8 ( r out = 29 M BH ) (e) − − t (sec)0 . . . . . . L ν / ( c d M f a ll / d t ) K8 ( r out = 41 M BH )D8 ( r out = 29 M BH ) (f) − − − t (sec)10 − − M e j e a nd M f a ll ( M (cid:12) ) K8 ( r out = 41 M BH ),ejefallD8 ( r out = 29 M BH ) FIG. 15. The same as Fig. 9 but for comparison between models K8 and D8.[1] B. P. Abbott et al., Phys. Rev. Lett. , 116101 (2017).[2] B. P. Abbott et al., Astrophys. J. , L12 (2017).[3] D. Kasen, B. D. Metzger, J. Barnes, E. Quataert, and E.Ramirez-Ruiz, Nature , 80 (2017).[4] B. Margalit and B. D. Metzger, Astrophys. J. , L19(2017).[5] M. Shibata et al., Phys. Rev. D , 123012 (2017).[6] A. Perego, D. Radice, and S. Bernuzzi, Astrophys. J. , L37 (2017).[7] A Goldstein et al., Astrophys. J. , L14 (2017); V.Savchenko et al., Astrophys. J. , L15 (2017).[8] K. P. Mooley et al. Nature , 355 (2018).[9] B. P. Abbott et al., arViv: 2001:01761. [10] R. Fern´andez and B. D. Metzger, Mon. Not. Royal As-tron. Soc. , 502 (2013).[11] B. D. Metzger and R. Fern´andez, Mon. Not. Royal As-tron. Soc. , 3444 (2014).[12] A. Perego, et al. on Not. R. Rstron. Aoc. , 3134(2014).[13] R. Fern´andez, D. Kasen, B. D. Metzger, and E. Quataert,Mon. Not. Royal Astron. Soc. , 750 (2015).[14] O. Just, A. Bauswein, R. A. Pulpillo, S. Goriely, and H.-Th. Janka Mon. Not. Royal Astron. Soc. , 541 (2015).[15] D. M. Siegel and B. D. Metzger, Phys. Rev. Lett. ,231102 (2017): Astrophys. J. , 52 (2018).[16] R. Fern´andez, A. Tchekhovskoy, E. Quataert, F. Fou- cart, and D. Kasen, Mon. Not. R. Astron. Soc. , 3373(2019).[17] A. Janiuk, Astrophys. J. , 2 (2019).[18] I. M. Christie et al., Mon. Not. R. Astron. Soc. , 4811(2019).[19] J. M. Miller et al., Phys. Rev. D , 023008 (2019).[20] S. Fujibayashi, K. Kiuchi, N. Nishimura, Y. Sekiguchi,and M. Shibata, Astrophys. J , 64 (2018).[21] S. Fujibayashi, K. Kiuchi, K. Kyutoku, N. Nishimura, Y.Sekiguchi, M. Shibata, and S. Wanajo, in preparation.[22] M. Shibata, Numerical Relativity (World Scientific, Sin-gapore, 2016)[23] M. Shibata and T. Nakamura, Phys. Rev. D ,5428(1995): T. W. Baumgarte and S. L. Shapiro, Phys.Rev. D , 024007(1998).[24] M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlo-chower, Phys. Rev. Lett. , 111101 (2006): J. G. Baker,J. Centrella, D.-I. Choi, M. Koppitz, and J. van Meter,Phys. Rev. Lett. , 111102 (2006).[25] D. Hilditch, S. Bernuzzi, M. Thierfelder, Z. Cao, W.Tichy, and B. Br¨ugmann, Phys. Rev. D , 084057(2013).[26] S. A. Balbus and J. F. Hawley, Rev. Mod. Phys. , 1(1998).[27] M. Alcubierre, et al., Int. J. Mod. Phys. D , 273 (2001).[28] M. Shibata, Prog. Theor. Phys. (2000), 325; M. Shi-bata, Phys. Rev. D (2003), 024033; M. Shibata andY. Sekiguchi, Prog. Theor. Phys. , 535 (2012).[29] M. Shibata, Phys. Rev. D , 064035 (2007).[30] J. Karkowski et al., Phys. Rev. D , 104017 (2018).[31] S. Banik, M. Hempel, and D. Bandyophadyay, Astro-phys. J. Suppl. Ser. , 22 (2014).[32] F. X. Timmes and F. D. Swesty, Astrophys. J. Suppl. , 501 (2000).[33] Y. Sekiguchi, K. Kiuchi, K. Kyutoku, and M. Shibata,Phys. Rev. D , 064059 (2015); Phys. Rev. D , 124046(2016).[34] F. Foucart, R. Haas, M. D. Duez, E. O’Connor, C. D.Ott, L. Roberts, L. E. Kidder, J. Lippuner, H. P. Pfeiffer,and M. A. Scheel, Phys. Rev. D , 044019 (2016).[35] D. Radice, F. Galeazzi, J. Lippuner, L. F. Roberts, C.D. Ott, and L. Rezzolla, Mon. Not. R. Soc. Astron. ,3255 (2016).[36] L. Bovard, D. Martin, F. Guercilena, A. Arcones, andL. Rezzolla, and O. Korobkin, Phys. Rev. D , 124005(2017).[37] F. Foucart, E. O’Connor, L. Roberts, M. Duez, R. Haas,L. E. Kidder, C. D. Ott, H. P. Pfeiffer, M. A. Scheel,and B. Szilagyi, Phys. Rev. D , 124021 (2015): F.Foucart, D. Desai, W. Brege, M. D. Duez, D. Kasen, D.A. Hemberger, L. E. Kidder, H. P. Pfeiffer, and M. A.Scheel, Class. Quantum Grav. , 044002 (2017). [38] K. Kyutoku, K. Kiuchi, Y. Sekiguchi, M. Shibata, andK. Taniguchi, Phys. Rev. D , 023009 (2018).[39] J. F. Hawley, Astrophys. J. , 496 (1991).[40] O. Korobkin, E. B. Abdikamalov, E. Schnetter, N. Ster-gioulas, and B. Zink, Phys. Rev. D , 043007 (2011).[41] K. Kiuchi, M. Shibata, P. J. Montero, and J. A. Font,Phys. Rev. Lett. , 251102 (2011).[42] N. I. Shakura and R. A. Sunyaev, Astron. Astrophys. ,337 (1973).[43] K. Kyutoku, M. Shibata, and K. Taniguchi, Phys. Rev.D , 044049 (2010).[44] J.-L. Tassoul, Theory of Rotating Stars (Princeton Uni-versity Press, 1978).[45] J. M. Lattimer and T. J. Mazurek, Astrophys. J. ,955 (1981).[46] W. H. Lee, E. Ramirez-Ruiz, and D. Page, Astrophys. J. , 421 (2005).[47] S. E. Woosley, Astrophys. J. , 273 (1993).[48] D. M. Siegel, J. Barnes, and B. D. Metzger, Nature ,241 (2019).[49] S. Wanajo, B. M¨uller, H-Th. Janka, and A. Heger, As-trophys. J. , 40 (2018).[50] T. Tachibana, M. Yamada, and Y. Yoshida, Prog. Theor.Phys. , 641 (1990).[51] S. Goriely, S. Hilaire, and A. J. Koning, Astron. Astro-phys. , 767 (2008).[52] S. Goriely, N. Chamel, and J. M. Persron, Phys. Rev. C , 035804 (2010).[53] S. Goriely, Astron. Astrophys. , 881 (1999).[54] J. J. Cowan, John J., C. Sneden J. E., Lawler, A. Apra-hamian, M. Wiescher, K. Langanke, G. Mart´ınez-Pinedo,and F-K. Thielemann, arXiv:1901.01410 (2019).[55] D. Kasen, N. R. Badnell, and J. Barnes, Astrophys. J. , 25 (2013): J. Barnes and D. Kasen, Astrophys. J. , 18 (2013).[56] M. Tanaka and K. Hotokezaka, Astrophys. J. , 113(2013).[57] M. Tanaka et al., Astrophys. J. , 109 (2018): arXiv:1906.08914.[58] D. Kasen, R. Fern´andez, and B. Metzger, Mon. Not. R.Astron. Soc. , 1777 (2015).[59] J. F. Hawley, X. Guan, and J. H. Krolik, Astrophys. ,84 (2011).[60] K. Kiuchi et al., Phys. Rev. D , 124039 (2018).[61] J. M. Miller et al., arXiv: 1912:03378.[62] K. Kyutoku, S. Fujibayashi, K. Hayashi, K. Kawaguchi,K. Kiuchi, M. Shibata, and M. Tanaka, in submission toAstrophys. J.[63] G. M. Fuller, W. A. Fowler, and M. J. Newman, Astro-phys. J. , 1 (1985).[64] K. Langanke and G. Mart´ınez-Pinedo, Nucl. Phys. A , 481 (2000).[65] A. M. Beloborodov, Astrophys. J.588