Alternative possibility of GW190521: Gravitational waves from high-mass black hole-disk systems
Masaru Shibata, Kenta Kiuchi, Sho Fujibayashi, Yuichiro Sekiguchi
AAlternative possibility of GW190521: Gravitational waves from high-mass blackhole-disk systems
Masaru Shibata,
1, 2
Kenta Kiuchi,
1, 2
Sho Fujibayashi, and Yuichiro Sekiguchi Max Planck Institute for Gravitational Physics (Albert Einstein Institute),Am M¨uhlenberg 1, Potsdam-Golm 14476, Germany Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto, 606-8502, Japan Department of Physics, Toho University, Funabashi, Chiba 274-8510, Japan (Dated: January 15, 2021)We evolve high-mass disks of mass 15–50 M (cid:12) orbiting a 50 M (cid:12) spinning black hole in the frameworkof numerical relativity. Such high-mass systems could be an outcome during the collapse of rapidly-rotating very-massive stars. The massive disks are dynamically unstable to the so-called one-armedspiral-shape deformation with the maximum fractional density-perturbation of δρ/ρ (cid:38) .
1, andhence, high-amplitude gravitational waves are emitted. The waveforms are characterized by aninitial high-amplitude burst with the frequency of ∼ × − at the hypothetical distance of 100 Mpc and by a subsequent low-amplitude quasi-periodic oscillation. We illustrate that the waveforms in our models with a wide range of the diskmass resemble that of GW190521. We also point out that gravitational waves from rapidly-rotatingvery-massive stars can be the source for 3rd-generation gravitational-wave detectors for exploringthe formation process of rapidly-rotating high-mass black holes of mass ∼ M (cid:12) in an earlyuniverse. PACS numbers: 04.25.D-, 04.30.-w, 04.40.Dg
I. INTRODUCTION
In their third observational run, Advanced LIGO andAdvanced Virgo announced the discovery of a high-massbinary black hole, GW190521 [1]. The estimated mass oftwo black holes are m = 85 +21 − M (cid:12) and m = 66 +17 − M (cid:12) with the estimated distance to the source of 5 . +2 . − . Gpc.Stellar evolution calculations predict that black holes ofmass gap ∼ M (cid:12) should not be formed from themassive stars because of the presence of pulsational pairinstability and pair-instability explosion for the massive-star evolution with a range of the initial stellar mass [2–4]. For this reason, the announcement of GW190521 bythe LIGO-Virgo collaboration has been attracting muchattention and several possibilities for forming the blackholes in the mass gap have been suggested (see, e.g.,Refs. [5–9]).Subsequently Nitz and Capano [10] reanalyzed thedata with different priors and waveform models fromthose in Ref. [1]. They found that more probable com-ponent mass might be m = 168 +15 − M (cid:12) and m =16 +33 − M (cid:12) , and thus, the mass of them might not be inthe mass gap range of 50–120 M (cid:12) . A very-massive blackhole of mass (cid:38) M (cid:12) could be naturally interpretedas a remnant of the collapse of a very-massive star (seebelow). Another interesting point in their work is thattheir maximum-likelihood waveform indicates that afterthe burst waveform of a few cycles with the highest ampli-tude, gravitational waves might be emitted: They showthat with their maximum-likelihood waveform in whichthe post-burst waves are present, the signal-to-noise ratiois higher than in the absence of the post-burst waves. Forbinary black hole mergers, the gravitational waveformsafter a peak amplitude is reached should be character- ized by a ringdown waveform and the amplitude shouldapproach zero quickly. The possible presence of grav-itational waves after the burst-waves emission suggeststhat the source of GW190521 might not be the mergerof high-mass binary black holes.Besides the frequency of gravitational waves, the wave-forms such as that of GW190521, which are composedinitially of burst waves of a few wave cycles and subse-quently of a small-amplitude oscillation, have been of-ten reported in the results of numerical simulations forthe systems dynamically unstable to non-axisymmetricdeformation. For example, during the collapse of veryrapidly rotating stars, the non-axisymmetric deformationcan occur and the resulting gravitational waveform canbe qualitatively similar to that of GW190521 (besidesthe frequency: see, e.g., Fig. 15(b) of Ref. [11]). Thus,it would not be non-sense to explore other astrophysi-cal possibilities for interpreting the observational resultsof GW190521, supposing that post-burst gravitationalwaves indicated in Ref. [10] may indeed be real signals.Based on this motivation, we perform a numerical-relativity simulation for the system of a massive blackhole and a massive disk. We choose the black holeof initial mass M BH , = 50 M (cid:12) and the disk of restmass M disk = 15–50 M (cid:12) in this paper. Here the black-hole mass of 50 M (cid:12) is chosen in order to reproduce thetypical frequency of gravitational waves for GW190521which is ∼
50 Hz [1] (note that the frequency is ap-proximately proportional to M − , ). Such a high-massblack hole-disk system could be formed during the col-lapse of very-massive metal-poor stars of initial masslarger than ∼ M (cid:12) [12–18]. Here the criticalinitial mass depends on the angular momentum of theprogenitor stars [15, 17]. For these stars, the mass of a r X i v : . [ a s t r o - ph . H E ] J a n the carbon-oxygen core formed at the final stage of thestellar evolution is M CO (cid:38) M (cid:12) and the collapse istriggered by the pair instability. The collapse of suchheavy cores leads to the direct formation of a high-massblack hole, avoiding the pair-instability explosion, whilefor 70 M (cid:12) (cid:46) M CO (cid:46) M (cid:12) , a pair-instability supernovais the final outcome, leaving no compact objects at thecenter [15, 17].Here, it is important to emphasize that the mass of theblack hole during its formation may not be always equalto the entire mass of the carbon-oxygen core and it couldbe smaller than M CO in the presence of rapid rotation.One reason for this is that a fraction of the mass canform disks (and ejecta from the disks). Indeed, our pre-vious work in numerical relativity shows that this is thecase: Reference [19] illustrates that after the collapse of amoderately rapidly-rotating very-massive star of the ini-tial mass 320 M (cid:12) and M CO ≈ M (cid:12) , which is obtainedby a stellar evolution calculation [15], the outcome is ablack hole of mass ∼ M (cid:12) surrounded by a disk ofmass ∼ M (cid:12) . Our subsequent work further illustratesthat in the presence of the neutrino cooling effect whichwas absent in Ref. [19], the initial mass of the black holeformed during the early stage of the collapse is likelyto be substantially smaller than M CO [20] because thecollapse of the central region is significantly acceleratedin a runaway manner toward the earlier formation of ablack hole. The temporal mass of the black hole and diskformed during the collapse could also depend strongly onthe angular velocity profile of the pre-collapse carbon-oxygen core, which is not well understood by the stellarevolution calculations. Thus, in this paper, we supposethe case in which a high angular momentum is present inthe carbon-oxygen core just prior to the collapse, and asystem of a black hole and a disk, for both of which themass is several tens of solar mass, could be formed dur-ing the collapse of rapidly-rotating very-massive stellarcores.By numerically evolving high-mass black hole-disk sys-tems as a plausible outcome formed during the collapseof rapidly-rotating very-massive stars, we show that themassive disks in such systems are dynamically unsta-ble to the one-armed spiral-shape deformation as oftenfound in previous works for low-mass disks [21–25] (butsee also Ref. [26] for the possible importance of the mag-netic field effects for low-mass disks). Then, we demon-strate that gravitational waves emitted from the non-linearly perturbed disks are indeed composed of a high-amplitude initial burst and a subsequent low-amplitudequasi-periodic oscillation. For M BH , = 50 M (cid:12) , the fre-quency of the initial burst gravitational waves is ∼ ∼ G and c denote the gravitational constant and the speedof light, respectively. II. SET-UP OF NUMERICAL SIMULATIONS
The set-up of the present numerical simulations is es-sentially the same as that in our previous papers on ax-isymmetric simulations for the massive black hole-disksystems [28, 29], but in this paper, we perform three-dimensional hydrodynamics simulations. We numericallysolve Einstein’s equation and the hydrodynamics equa-tions. Einstein’s equation is solved using the original ver-sion of the Baumgarte-Shapiro-Shibata-Nakamura for-malism [30] together with the puncture formulation [31],Z4c constraint propagation prescription [32], and 6th-order Kreiss-Oliger dissipation. For hydrodynamics, wedo not consider the neutrino emission and weak inter-action throughout this paper, although our implemen-tation can take into account the neutrino effects in thesame manner as in Refs. [33, 34]. In the presence of theneutrino emission, the disk shrinks by the neutrino cool-ing and its compactness is increased. This effect makesthe dependence of the growth timescale of the disk insta-bility and resulting gravitational waveforms on the diskmorphology obscure. Thus we decided to switch off theneutrino effects in this work. In this prescription, onlythe advection part for the equations of the electron frac-tion, Y e , is taken into account.The quantities of black holes (mass and spin) are de-termined from their area and circumferential radii of ap- TABLE I. Initial conditions for the numerical simulation. De-scribed are the model name, initial black-hole mass ( M BH , ),rest mass of the disk ( M disk ), dimensionless spin of the blackhole ( χ ), the coordinate radii at the inner and outer edgesof the disk ( r in and r out ), entropy per baryon ( s/k with k the Boltzmann’s constant), the maximum density of thedisk ( ρ max ), and the orbital period at the maximum den-sity ( P orb ). The units are M (cid:12) for the mass, GM BH , /c ≈ . M BH , / M (cid:12) ) km for r in and r out , 10 g / cm for ρ max ,and GM BH , /c for P orb , respectively.Model M BH , M disk χ r in r out s/k ρ max P orb L11 50.0 15.1 0.78 2.0 21.8 11 4.35 119L12 50.0 15.1 0.79 2.0 27.9 12 2.85 126J12 50.0 20.0 0.78 2.0 25.4 12 3.87 129J13 50.0 19.9 0.78 2.0 32.8 13 2.54 135K13 50.0 25.0 0.78 2.0 29.7 13 3.32 137K14 50.0 25.0 0.78 2.0 38.6 14 2.21 144M12 50.0 30.1 0.77 2.0 21.1 12 6.62 130M13 50.0 30.0 0.77 2.0 27.0 13 4.26 139M14 50.0 29.9 0.78 2.0 35.0 14 2.79 147M15 50.0 30.1 0.78 2.0 45.6 15 1.87 153H14 50.0 40.0 0.76 2.0 28.9 14 4.39 150H15 50.0 40.0 0.76 2.0 37.4 15 2.89 157H16 50.0 40.0 0.77 2.0 49.0 16 1.90 163E15 50.0 50.1 0.76 2.0 31.0 15 4.35 158E16 50.0 50.2 0.76 2.0 40.0 16 2.86 168E17 50.0 50.0 0.77 2.0 52.6 17 1.86 175 parent horizons [35], assuming that these quantities arewritten as functions of the mass and spin in the sameformulation as in the vacuum Kerr black hole. Gravi-tational waves are derived through the extraction of theoutgoing component of the complex Weyl scalar [36].As in our previous work [28, 29], we employ a tab-ulated equation of state based on the DD2 equationof state [37] for a relatively high-density part and theTimmes (Helmholtz) equation of state for the low-densitypart [38]. In this tabulated equation of state, thermo-dynamics quantities such as ε , P , and h are written asfunctions of ρ , Y e , and T where ε , P , h (= c + ε + P/ρ ), ρ ,and T are the specific internal energy, pressure, specificenthalpy, rest-mass density, and matter temperature, re-spectively. We choose the lowest rest-mass density to be0 . / cm in the table.Axisymmetric equilibrium states for black hole-disksystems are prepared as the initial conditions [28, 29]using the method of Ref. [39]. As in Refs. [28, 29], wedetermine the angular velocity profile from the relation j ∝ Ω − n , (2.1)where j = c − hu ϕ is the specific angular momentum. Ωis the angular velocity defined by u ϕ /u t with u µ the fluidfour-velocity. n is a constant that determines the profileof the angular velocity, for which we choose n = 1 / ρ and Y e for the high-density mat-ter in the same form of ρ ( Y e ) as in Ref. [28]. For thismodel, the value of Y e is 0.07 in a high density region with ρ ≥ g / cm , 0 . − (43 / [ ρ/ (10 g cm − )]for 10 g / cm ≤ ρ ≤ g / cm , and settles to 1 / ρ ≤ g / cm , for which the effect of the electron de-generacy to the neutron richness is weak. In addition,we assume that the specific entropy, s , is initially con-stant, in order to obtain the first integral of the Eulerequation easily. For the high-mass disks, we choose it s/k = 11–17 in this paper where k is the Boltzmann’sconstant: By adjusting this value, the compactness forgiven mass is controlled (for the higher entropy, the diskbecomes less compact for a given value of the disk mass).We note that our setting for the disk profile is ideal-ized. For more realistic work, it is obviously necessary toperform a simulation started from rapidly-rotating very-massive progenitor stellar cores. However, the collapsesimulation with a wide variety of the parameters (e.g.,angular momentum profile) is computationally very ex-pensive. The purpose of the present work is to under-stand the nature of the non-axisymmetric instability ofthe massive disks orbiting a black hole and resulting grav-itational waves emitted. For this, we believe that thepresent setting is acceptable.The initial black-hole mass is fixed to be M BH , =50 M (cid:12) while a wide range of the disk mass, M disk ≈ M (cid:12) , is employed. The initial dimensionless spin ofthe black hole is set to be approximately 0.8 suppos-ing that the very-massive stellar cores, which form amassive disk during the collapse, should be rapidly ro-tating. We set the inner coordinate radius of the diskto be close to the radius of the innermost stable cir-cular orbit around the rapidly spinning black hole as r in ≈ GM BH , /c , and the outer edge is varied in therange of r out / ( GM BH , c − ) = 21–53 (see Table I). Theseare the plausible sizes of the dense region of the diskformed during the collapse of rotating very-massive stel-lar cores [19]. With these setups, the maximum density ofthe disks is of the order of 10 g / cm (see Table I). Thus,the matter in the high-density region is mildly neutron-rich.Numerical simulations are performed with a fixedmesh-refinement implementation used for the neutron-star mergers [40–42]. We prepare 10 refinement levelswith each region composed of 241 × ×
121 uniform gridpoints for x - y - z (the reflection symmetry with respect tothe equatorial plane is imposed). For the finest level, thegrid spacing is 0 . GM BH , /c , and the outer boundaryalong each axis is located at L ≈ GM BH , /c for allthe simulations. We checked that with such setting theblack hole can be accurately evolved within the fractionalchange of its mass and dimensional spin of 0 .
1% beforethe onset of the mass accretion resulting from the diskinstability.
FIG. 1. Snapshots of the rest-mass density profile of the disk in the equatorial plane for model M14. P orb denotes the initialorbital period of the disk at the density maximum (see Table I). The color bar shows log ρ with the unit of ρ being g / cm We do not add any perturbation initially, except fornumerical noises. By a small random perturbation au-tomatically introduced in numerical evolution, the dy-namically unstable deformation of the massive disks isenhanced for all the simulations.
III. EVOLUTION OF DISKS ANDGRAVITATIONAL WAVESA. Evolution of dynamically unstable disks
Irrespective of the initial conditions prepared in thispaper, the disk is dynamically unstable to the one-armedspiral-shape deformation. Figure 1 plots the snapshots ofthe rest-mass density profile of the disk in the equatorialplane for model M14 (see Table I for the initial condi-tion of this model). It is found that a small numericalnoise introduced during the simulation triggers the subse-quent exponential growth of the one-armed spiral-shapedeformation mode. We note that in the present con-text (i.e., in fully general relativistic computation), theenhanced density perturbation in the disk enforces theblack hole to slightly shift to the higher-density side ofthe disk, because the center of mass of the system shouldbe preserved. This enhances the gravitational attractionbetween the black hole and the high-density part of the disk, and accelerates the growth of the one-armed densityperturbation. Note that this effect cannot be taken intoaccount in the simulations with fixed background space-time, and thus, a fully general relativistic simulation isneeded to explore the one-armed spiral-shape instabilityfor high-mass disks.After the significant growth of the spiral-shape de-formation, the angular momentum transport inside thedisk, which is caused by the gravitational torque fromthe non-axisymmetric density perturbation as well as bythe black-hole tidal field, is enhanced. As a result, asignificant mass accretion onto the black hole proceeds.During the initial growth of the one-armed spiral-shapedeformation until its saturation, 20–30% of the total restmass of the disk falls into the black hole. Subsequently,a weakly perturbed disk is formed, and from such a disk,quasi-steady mass inflow to the black hole continues (e.g.,Refs. [24, 25]). For all the cases, the black hole spins upby the mass accretion and the dimensionless spin alwaysbecomes higher than 0.85 after a substantial fraction ofthe mass falls into the black hole. In particular for themassive and compact disk models such as M12, H14, E15,and E16, the dimensionless spin exceeds 0 . (a) t/P orb − − − − − − − δ m m = 1, M14m = 2m = 3m = 4 (b) t/P orb − − − − − − − δ m m = 1, M12m = 1, M13m = 1, M14m = 1, M15 (c) t/P orb − − − − − − − δ m m = 1, E15m = 1, E16m = 1, K13m = 1, K14m = 1, L11m = 1, L12 (d) t/P orb − − − − − − − δ m m = 1, E15m = 1, H14m = 1, M13m = 1, K13m = 1, J13m = 1, L12 FIG. 2. Mode amplitude for the density perturbation defined by | C m | /C as a function of time. (a) Evolution of the m = 1–4modes for model M14; (b) Evolution of the m = 1 modes for models M12, M13, M14, and M15; (c) Evolution of the m = 1modes for models L11, L12, K13, K14, E15, and E16; (d) Evolution of the m = 1 modes for models L12, J13, K13, M13, H14,and E15. The time is shown in units of P orb , which denotes the initial orbital period of the disks at the maximum rest-massdensity (see Table I). ulations of relatively short duration, we find that 1–2%of the total mass of the disk is ejected from the system.As a result of the mass ejection (as well as of the massaccretion into the black hole), the maximum density ofthe disk decreases (compare the first and last panels ofFig. 1). With a longer-term simulation taking into ac-count the viscous effect, we could find that a more frac-tion of the disk matter will be ejected from the system.However, exploring such a long-term mass ejection pro-cess is beyond the scope of this paper.Figure 2 displays the evolution of the mode amplitudeof the density perturbation for the selected models. Todefine the mode amplitude, we first extract each non-axisymmetric mode of the density perturbation by cal-culating C m = (cid:90) ρ ∗ e imϕ d x, (3.1)with m = 0–4 and ϕ = tan − ( y/x ). Here, the volumeintegral is performed only for the outside of the apparenthorizon. We then define the normalized mode amplitudeby δ m = | C m | C . (3.2) In the following, we refer to δ m as the mode amplitude.For all the simulations, δ m initially increases by thenumerical noise. However, after δ reaches ∼ − , theone-armed spiral instability, which is developed from theinitial seed perturbation that results from the numericalnoise, makes the dominant perturbation. Thus, in thefollowing, we focus only on the growth rate for δ ≥ − .Figure 2(a) compares the mode amplitude of m = 1–4for model M14. This clearly shows that only the m = 1mode initially grows exponentially with time (focusing onthe plot only for the stage of δ (cid:38) − ). During the earlyexponential growth of the m = 1 mode, the other modeswith m ≥ m = 1 modeamplitude exceeds 10 − , we find the exponential growthfor the mode amplitudes of m = 2 and m = 3 but thegrowth of these modes appears to result from the non-linear growth of the m = 1 mode because the growth ratefor m = 2 is approximately twice of that for m = 1. For m = 4, the mode amplitude is approximately constantwith δ ∼ − . This high value of δ is an artifact thatresults from the use of the Cartesian coordinates for thesimulation. Since this does not grow in time, no impact(and no harmful effect) is made by it for the evolution of TABLE II. The growth timescale of the dynamical instability in units of P orb , the maximum value of δ , and the signal-to-noiseratio (SNR) for gravitational waves in the detection at the hypothetical distance to the source of 100 Mpc with the mostoptimistic direction assuming the detection by a single advanced LIGO detector with the designed sensitivity. The last rowdescribes the type of the waveform (see the text for the classification).Model L11 L12 J12 J13 K13 K14 M12 M13 M14 M15 H14 H15 H16 E15 E16 E17 τ dyn /P orb δ , max the system.Figure 2(b)–(d) display the m = 1 mode amplitude, δ ,as a function of time for several models. These show that δ increases exponentially with time for the stage of δ (cid:38) − , for which we extract the growth timescale, τ dyn , byfitting δ in the form of ∝ exp( t/τ dyn ). It is found thatthe exponential growth timescale for the models chosenin this paper is 0.5–2.5 times of the initial orbital period, P orb , of the matter at the highest-density region of thedisk (cf. Table II). This implies that the instability foundhere is indeed the dynamical instability.Figure 2(b) compares the m = 1 mode amplitude as afunction of time for models M12–M15 for which the com-pactness of the disk defined by GM BH / ( c r out ) is differ-ent for a fixed value of the disk mass. In Fig. 2(c) we alsocompare the m = 1 mode amplitudes between E15 andE16, between K13 and K14, and between L11 and L12to see the dependence of the growth timescale of the in-stability on the disk compactness for a fixed value of thedisk mass. It is always found that the growth timescaleof the dynamical instability (in units of P orb ) is shorterfor the more compact disk models. We note that if wemeasure the growth timescale in units of GM BH , /c , thistrend is further remarkable. Thus, we can conclude thatmore compact disks are more subject to the dynamicalinstability. On the other hand, this result also suggeststhat a sufficiently less compact disk may be dynamicallystable. Figure 2(b) also shows that the maximum valueof δ is larger for more compact disks (see also Table II).Figure 2(d) compares the m = 1 mode amplitude formodels L12, J13, K13, M13, H14, and E15 for which c r out / ( GM BH ) is ∼
30. This shows that the growthtimescale of the dynamical instability is shorter for moremassive disks with a given disk compactness. In partic-ular for M disk /M BH , (cid:38) . P orb . On the other hand, this figure suggeststhat for a disk of mass smaller than a critical value thedisk may be dynamically stable [24]. This figure also il-lustrates that the maximum value of δ is larger for moremassive disks (see also Table II).Irrespective of the models, after δ reaches the max-imum of δ (cid:38) .
1, i.e., the density perturbation be-comes non-linear, the growth of the density perturba-tion saturates. Then, the degree of the deformation isapproximately preserved with the density perturbationof δ = 0 . δ will become smaller in the long-term evolution of t (cid:29) P orb .In Table II, we summarize the exponential growthtimescale, τ dyn , and the maximum value of δ , δ , max .Broadly speaking, for more compact or more massivedisks, τ dyn /P orb is smaller and δ , max is larger. As shownin Sec. III B, the dependence of these quantities on thecompactness and mass of the disks is reflected in thegravitational waveforms. We will find that three differenttypes of gravitational waveforms are produced by varying r out and the disk mass.Before closing this subsection, we note the followingpoint. In the presence of magnetic fields, the magnetoro-tational instability should occur [43] and subsequently, aturbulence could be excited in the disk [26, 44]. Then theturbulent viscosity is enhanced, and the angular momen-tum in the disk would be redistributed. If the viscous an-gular momentum transport efficiently works in the disk,the growth of the one-armed spiral instability studied inthis paper should be suppressed [26]. The timescale ofthe viscous angular momentum transport is τ vis := R ν , (3.3)where R and ν denote the cylindrical radius in the diskand the shear viscous coefficient, respectively. For the ge-ometrically thin disks, ν can be rewritten as ν = α ν c s H where α ν is the dimensionless alpha parameter [45], c s isthe sound velocity, and H is the scale height of the disk.Numerical simulations for accretion disks have shown α ν = O (10 − ) (e.g., Ref. [44]).Then the ratio of τ vis to P orb is τ vis P orb ≈ (2 πα ν ) − R Ω c s H ≈ (2 πα ν ) − (cid:18) RH (cid:19) , (3.4)where we used the approximate relation of the force bal-ance with respect to the vertical direction of the disk as H Ω ≈ c s . Because R/H (cid:38) τ vis /P orb >
10 as long as α ν is not unusually large, i.e., α ν (cid:46) .
05. Thus, for themassive and compact disks of τ dyn (cid:46) . P orb which we -3-2-1 0 1 2 3-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 L11 h ( a t M p c ) t ret (s) -1-0.5 0 0.5 1-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 L12 h ( a t M p c ) t ret (s) -3-2-1 0 1 2 3-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 J12 h ( a t M p c ) t ret (s) -1.5-1-0.5 0 0.5 1 1.5-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 J13 h ( a t M p c ) t ret (s) -3-2-1 0 1 2 3-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 K13 h ( a t M p c ) t ret (s) -1.5-1-0.5 0 0.5 1 1.5-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 K14 h ( a t M p c ) t ret (s) FIG. 3. Gravitational waves observed along the rotational axis at the hypothetical distance of 100 Mpc for the relativelylow-mass disk models L11, L12, J12, J13, K13, and K14. The time axis is chosen so that the time at the maximum amplitudebecomes approximately zero. The solid and dashed curves denote the plus and cross modes, respectively. We note that foreach panel the scale of the vertical axis is different. study in this paper, the timescale of the viscous angu-lar momentum transport is much longer than the growthtimescale of the one-armed spiral-shape instability. Thisis in particular the case for the models that result in thetype I waveforms which we are most interested in (seeSec III B). Therefore, the magnetic/viscous effects couldbe safely neglected for the models employed in this paper.However, as shown in Ref. [26], the magnetic/viscous ef-fects are likely to play an important role for low-mass oror less compact disks for which τ dyn (cid:38) P orb . B. Gravitational waves
Figures 3 and 4 display gravitational waves observedalong the rotation axis (i.e., perpendicular to the orbital plane) for all the models employed in this paper. Figure 3is for relatively less massive disks with 15 M (cid:12) ≤ M disk ≤ M (cid:12) and Fig. 4 is for 30 M (cid:12) ≤ M disk ≤ M (cid:12) .For many of the models studied in this paper (themodels except for L12, J13, K14, and M15), the wave-forms are characterized by initial burst waves of a highamplitude and subsequent quasi-periodic waves of a lowamplitude. We refer to this waveform type as type I.The initial burst waves are most efficiently emitted whenthe value of δ becomes the maximum. The numberof the wave cycle of the burst waves is 1–3 for thistype, and is smaller for more compact or more mas-sive disks for a given value of c r out / ( GM BH , ). Themaximum amplitude of the burst gravitational waves is10 − –10 − at the hypothetical distance to the sourceof D = 100 Mpc and is higher for more compact disks -10-5 0 5 10-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 M12 h ( a t M p c ) t ret (s) -6-4-2 0 2 4 6-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 M13 h ( a t M p c ) t ret (s) -4-3-2-1 0 1 2 3 4-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 M14 h ( a t M p c ) t ret (s) -1.5-1-0.5 0 0.5 1 1.5-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 M15 h ( a t M p c ) t ret (s) -8-6-4-2 0 2 4 6 8-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 H14 h ( a t M p c ) t ret (s) -4-2 0 2 4-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 H15 h ( a t M p c ) t ret (s) -3-2-1 0 1 2 3-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 H16 h ( a t M p c ) t ret (s) -10-5 0 5 10-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 E15 h ( a t M p c ) t ret (s) -6-4-2 0 2 4 6-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 E16 h ( a t M p c ) t ret (s) -4-3-2-1 0 1 2 3 4-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 × -22 E17 h ( a t M p c ) t ret (s) FIG. 4. The same as Fig. 3 but for the high-mass disk models M12, M13, M14, M15, H14, H15, H16, E15, E16, and E17.
10 15 20 25 30 35 40 45 50 55 20 25 30 35 40 45 50 55
Type Is Type IwType II M d i s k / M s un c r out / GM BH,0
FIG. 5. The types of the gravitational waveforms are sum-marized in the plane of c r out / ( GM BH , ) and M disk /M (cid:12) .The upper and lower dashed lines are M disk /M (cid:12) = 3 + c r out / ( GM BH , ) and M disk /M (cid:12) = −
10 + c r out / ( GM BH , ),respectively. with smaller values of c r out / ( GM BH ) or for more mas-sive disks. On the other hand, the maximum amplitudeof the quasi-periodic waves does not depend on the valueof c r out / ( GM BH ) as strongly as that of the initial burstwaves and is always (cid:46) − for D = 100 Mpc.For the waveforms of models L12, J13, K14 and M15,the amplitude of the initial burst gravitational waves isonly by a factor of 2–3 higher than the subsequent quasi-periodic waves in contrast to the models of the type Iwaveform. We refer to this waveform type as type II.For the type II waveform, the maximum amplitude is (cid:46) × − for D = 100 Mpc irrespective of the disk mass.In terms of δ , max , the type I waveform is obtained for δ , max (cid:38) .
15 while the type II waveform is for δ , max (cid:46) .
15 (see Table II).We can consider that there are two types in the type Iwaveform. For the models of the very compact or massivedisks, the burst waves are characterized by only one waveof a very high amplitude (see the waveforms for modelsM12, H14, E15, and E16). We refer to this waveformas type Is. The short duration of the highest amplitudephase is due to the fact that the peak values of δ for thesecases are relatively large and hence by the gravitationaltorque exerted from the non-axisymmetric structure, alarge amount of the disk matter is ejected from the high-density region of the disk settling to a less compact diskin a short timescale. For other type I models (L11, J12,K13, M13, M14, H15, H16, E17), the gravitational torqueat the moment that δ is highest is not large enough toexpel a large fraction of matter from the high-density re-gion of the disk. As a result, the burst waves are emittedfor a couple of wave cycles. We refer to the waveformtype as type Iw for this case. In terms of δ , max , the typeIw waveform is obtained only for δ , max (cid:46) .
25. We notethat the waveforms for models M13 and H15 may be clas-sified as type Is. Thus, we consider that these models are located near the boarder for distinguishing types Iw andIs (note that the waveform should change continuouslywith the change of the compactness and mass of the disk,and hence, it would not be possible to definitely classifythe waveform into a particular type).The numerical results presented here illustrate thatthe compactness and mass of the disk (which are de-termined by the angular momentum profile of the pro-genitor star prior to the collapse) are the key parame-ters for determining the type of the waveforms. Specifi-cally, for the models employed in this paper, the type IIwaveforms are obtained broadly for M disk /M (cid:12) < −
10 + c r out / ( GM BH , ), and otherwise, the type I waveform isthe result (see Fig. 5). Also, broadly for M disk /M (cid:12) > c r out / ( GM BH , ), the waveform becomes the type Is.We note that this would be quantitatively valid only for c r out / ( GM BH , ) ≥
20 and for the initial angular mo-mentum profile employed in this paper. However, weinfer that qualitatively similar relations would be satis-fied irrespective of the disk configuration for the massivedisks.For the type I waveform, the amplitude of the burst-wave part is appreciably larger than the quasi-periodicwaves subsequently emitted. Thus, if gravitational wavesof this type are detected by gravitational-wave detectorswith a low signal-to-noise ratio, it would not be very easyto confirm the detection of the quasi-periodic wave part,and essentially, the attention would be paid only to theburst-waves part. For the case of GW190521, the detec-tion with a sufficiently high signal-to-noise ratio appearsto be done only for the first ∼ ∼
1, which is not asmany as that of GW190521.Figure 6 displays the dimensionless effective amplitudeof gravitational waves for the selected models: The leftpanel compares the spectra for given values of the diskmass with different compactness and the right panel com-pares the spectra for approximately the same value of c r out / ( GM BH , )( ∼
30) with different disk mass. Here,0 (a) f [Hz]10 − − − − h e ff ( D = M p c ) E15E16K13K14L11L12Advanced LIGO (b) f [Hz]10 − − − − h e ff ( D = M p c ) E15H14M13K13J13L12Advanced LIGO
FIG. 6. Effective amplitude of gravitational waves as a function of the frequency (a) for models L11, L12, K13, K14, E15, andE16; (b) for models L12, J13, K13, M13, H14, and E15 for which c r out / ( GM BH , ) ∼
30. The designed sensitivity of advancedLIGO is taken from https://dcc.ligo.org/cgi-bin/DocDB/ShowDocument?docid=2974. the effective amplitude is defined from the Fourier trans-form of gravitational waves, h ( f ), as h eff := | f h ( f ) | ,where f denotes the frequency of gravitational waves.The hypothetical distance to the source is 100 Mpc inthis figure and we assume the observation from the ro-tation axis (i.e., from the most optimistic direction) asin Figs. 3 and 4. For the setting with M BH , = 50 M (cid:12) ,irrespective of the disk compactness, the frequency at thepeak value of h eff becomes f peak ≈ ∼
40 Hz, for the models withmore massive disks such as H14–H16 and E15–E17. Ourinterpretation for this is that the total mass inside theregion for which the disk has the highest density, M tot ,is appreciably larger than M BH , , and thus, the angularvelocity of the one-armed spiral mode is slower (in in-versely proportional to M tot ) for a give compactness ofthe disk.The peak amplitude of h eff is higher for more compactdisks (i.e., for smaller values of c r out / ( GM BH , )) with agiven value of the disk mass and for more massive disks,as expected from Figs. 3 and 4. The noteworthy featuresfound from Fig. 6 are as follows: Irrespective of the wave-form type, (i) the effective amplitude increases with thefrequency for f ≤ f peak and (ii) the effective amplitudesteeply drops for the frequency beyond f peak . Althoughthe feature (ii) is similar to those of binary black-holemergers, the feature (i) is totally different from them.Thus, in the detection of gravitational waves with a highsignal-to-noise ratio, it would be easy to distinguish thewaveforms of the massive black hole-disk systems fromthose of binary black hole mergers. However, for the de-tection with a low signal-to-noise ratio (which is in partic-ular the case for the low-frequency region of (cid:46)
30 Hz withthe gravitational-wave detectors operated in the third-observational run [46]), such difference is unlikely to beconfirmed.Figure 7 plots the spectrogram (dominant frequencyas a function of time) of gravitational waves for models L11, J12, K13, M13, M14, M15, H16, and E17. As ex-pected from Figs. 3 and 4, one prominent spot sharplyappears at f ∼ f ∼
40 Hz while for othermodels, it is ∼
50 Hz. On the other hand, for model M15of the type II waveform, not one prominent spot but anextend bright region appears. This feature is universallyfound for models L12, J13, and K14, i.e., for the low-mass or less compact models, and is different from thatof GW190521.The signal-to-noise ratio of gravitational waves withrespect to the designed sensitivity of advanced LIGOat the hypothetical distance of 100 Mpc is listed in Ta-ble II. Here, the signal-to-noise ratio is calculated for thedetection by a single detector. Note that the numeri-cal simulations were stopped at 0.3–0.5 s after the peakgravitational-wave amplitude is reached. Since long-lasting quasi-periodic waves with the duration longerthan 0.5 s is likely to be present, the signal-to-noise ratioshown here should be considered as the minimum value.Taking into account that the sensitivity of advancedLIGO in the third observational run was by a factor of ∼ ∼ FIG. 7. Spectrograms (relation between the dominant frequency as a function of time) are displayed for models L11, J12, K13,M13, M14, M15, H16, and E17. tively late-time universe. However, currently, we do nothave rich observational information on this, although thisquestion may be answered by the future observation [47].It is also worthy to note that in the future detectors currently proposed, such as Einstein Telescope [48] andCosmic Explorer [49], the signal-to-noise ratio for theevents studied in this paper will be more than 10 timesas high as the detectors currently in operation. For such2future detectors, the horizon distance will be deeper than1 Gpc and the formation rate of the very-massive stars islikely to be higher than in the late-time universe. Con-sidering that the typical frequency of gravitational wavesis approximately proportional to the inverse of the black-hole mass and that the future detectors will have a highsensitivity down to (cid:46)
20 Hz [48, 49], the detectability willbe also high for black holes of higher mass of ∼ M (cid:12) .Thus, in the future, gravitational waves not only frombinary black hole mergers but also from high-mass blackhole-disk systems could be an interesting target for ex-ploring the formation process and formation history ofmassive black holes of mass ∼ M (cid:12) . IV. DISCUSSION
Just prior to the formation of the systems composed ofa massive black hole and a massive disk which we consid-ered in this paper, a stellar core collapse of very-massivestar, triggered by the pair instability, should occur. Sincethe black hole is likely to be formed immediately afterthe collapse (e.g., Refs. [19, 20]), electromagnetic coun-terparts associated with the usual core-collapse super-novae are unlikely to accompany with this. However,after the formation of a black hole, the material could beejected from the massive disk surrounding the black holefor the long-term evolution, and such material injects theoutgoing momentum and energy into the extended en-velop of radius ∼ cm surrounding the carbon-oxygencore [19]. For a sufficiently high energy injection, theenvelop will explode as in the core-collapse supernovae,and the ejecta may be the source for an electromagneticcounterpart in the optical-infrared bands with a high lu-minosity. The brightness of the supernova-like explosionis likely to depend strongly on the kinetic energy, E kin ,of the mass outflow from the disk [19]. Our previouswork [29] shows that 10–20% of the disk material couldbe ejected by a viscous process with the typical veloc-ity of 0 . c from massive disks surrounding central blackholes. Assuming these values, we have E kin ≈ erg (cid:18) M eje M (cid:12) (cid:19) (cid:16) v eje . c (cid:17) , (4.1)where M eje and v eje are the rest mass and velocity of theejecta. Thus for the ejecta mass of 1–10 M (cid:12) , the injectionenergy is broadly considered to be 10 –10 erg. We notethat the ejecta energy may be increased by subsequentnucleosynthesis.Using the Arnett’s model [50], Ref. [19] estimatedthat for E kin = 10 –10 erg and the envelop mass of ∼ M (cid:12) , the absolute luminosity of an explosion wouldbe ∼ –10 erg/s with the duration of O (1 yr). Thus,for a hypothetical distance of 100 Mpc, the apparent mag-nitude is ∼ ∼ M (cid:12) would be rare andhence the typical distance to this collapse event wouldbe large (cid:29)
10 Mpc. Thus, the possible electromagneticcounterparts of this kind of the event are not likely to bevery luminous, and it might not be easy to detect it in theabsence of the information of the sky localization. As westudied in this paper, gravitational waveforms from thecollapse of rapidly-rotating very-massive stars could besimilar to that of the merging binary black holes witha value of high chirp mass (cid:38) M (cid:12) . Thus an alert forsuch an event could be issued as a candidate for a bi-nary black hole merger as in the case of GW190521 inthe future. Although a careful follow-up observation byelectromagnetic telescopes is usually absent for the can-didates announced as binary black holes, in the future,it will be interesting to perform follow-up observationsfor the candidates of very high-mass binary black holes.The problem is that quite unfortunately the LIGO-Virgocollaboration never announce the chirp mass (i.e., typicalfrequency) of the candidate events of gravitational waves.Thus it will not be possible to distinguish the interestingevents among a huge amount of the “binary black hole”candidates unless the policy of the LIGO-Virgo collabo-ration is changed in the future.In this paper, we focus only on the excitation of anunstable mode in massive disks surrounding a high-massblack hole and associated gravitational waves emitted.Another interesting aspect of this system is that a largeamount of matter could be ejected from the massive diskduring the enhancement of the one-armed spiral densitywave and probably through the subsequent long-term vis-cous evolution. For a reliable estimation of E kin men-tioned above, it is important to perform a long-term sim-ulation for exploring the mass ejection process [28, 29]and to quantitatively explore the value of E kin . In theejecta from the disk, nucleosynthesis could also occur andthe radio-active energy could be released for the addi-tional source of the electromagnetic counterparts. If anamount of Ni with mass larger than ∼ M (cid:12) is synthesizedin the ejecta, the event may be similar to the luminoussupernovae (e.g., Ref. [52]). The detailed exploration ofthe mass ejection and nucleosynthesis is an interestingtopic for the subsequent work. It is also interesting toexplore how much mass can be ejected from the system.If a significant fraction of the matter forms a massivedisk and is subsequently ejected from the system, the fi-nal mass of the formed black hole can be smaller thanthe progenitor mass, and may be smaller than 120 M (cid:12) ,i.e., may come into the mass-gap range.There is another possible phenomenon for emittinggravitational waves similar to GW190521. In this pa-per, we suppose that a black hole is formed during thecollapse. If the centrifugal force in the central region of3the collapsing very-massive stars is strong enough, not ablack hole but a very massive spheroidal or toroidal ob-ject can be formed (see, e.g., Refs. [11, 22] for other con-texts). If such an outcome is very compact and dynam-ically unstable to one-armed spiral mode or bar-modedeformation, gravitational waves with a high amplitude,which would have the same order of magnitude as thatfrom massive black hole-disk systems, can be emitted.Assuming the formation of a compact object, the esti-mated frequency of gravitational waves is f = 1 π (cid:115) GM obj R ≈
58 Hz (cid:18) M obj M (cid:12) (cid:19) − (cid:18) R obj GM obj c − (cid:19) − / , (4.2)where M obj is the mass of the spheroidal/toroidal objectand R obj is the typical radius for the highest-density partof such an object at the onset of the dynamical insta-bility, respectively. Thus, the formation of a compactspheroid/toroid from a rapidly rotating very-massivestellar core can be another scenario for GW190521.Finally, we note that for forming a rapidly-rotatingvery-massive star, a metal-poor star would be necessary.Such a star is unlikely to be formed frequently in thelate-time universe, and thus, it would be rare that thecollapse of the rapidly-rotating very-massive star is ob-served at the distance of several hundred Mpc. How-ever, for given uncertainty on the formation rate of suchstars in metal-poor galaxies [47], it is important to im-pose the constraint using the observational results. Fromthis perspective, the gravitational-wave observation to-gether with the electromagnetic-counterpart search canprovide us a valuable opportunity.To summarize, in this paper, we have suggested thatburst gravitational waves of frequency ∼
50 Hz, which areobserved with a small signal-to-noise ratio and are com-posed only of a few cycles of high-amplitude waves, can be interpreted not only by those of binary black holesbut by those by other scenarios. We suspect that thesource of GW190521 might not be a merger of binaryblack holes but a stellar collapse of a very massive starleading temporarily to a black hole of mass ∼ M (cid:12) anda massive disk of several tens of solar mass that is dy-namically unstable to the one-armed spiral-shape defor-mation. The point of this scenario is that the progenitorshould be heavy enough to form a high-mass black holeof mass ∼ M (cid:12) and a high-mass disk, which can emitgravitational waves of low frequency ∼
50 Hz, althoughwe do not need the fine-tuning of the disk mass. Oneshould keep in mind that it is not easy to definitely con-clude the source of this-type of burst events only from agravitational-wave observation with a low signal-to-noiseratio. In the future, electromagnetic counterpart searcheswill play a key role for pining down the most probablescenario and an alert from gravitational-wave observa-tion community suitable for such searches is obviouslynecessary.In the present paper, we prepare an initial condition ofa massive black hole and a massive disk supposing thatit would be the outcome formed during the collapse ofa rapidly-rotating very-massive stellar core. Obviously,this setting is idealized. For more realistic work, it isnecessary to perform a simulation started from a rapidly-rotating very-massive progenitor star. We plan to per-form such simulations in the subsequent work.
ACKNOWLEDGMENTS
We thank Koh Takahashi, Masaomi Tanaka, andTakami Kuroda for helpful conversations. This workwas in part supported by Grant-in-Aid for ScientificResearch (Grant Nos. JP16H02183, JP18H01213, andJP20H00158) of Japanese MEXT/JSPS. Numerical com-putations were performed on Sakura and Cobra clustersat Max Planck Computing and Data Facility. [1] R. Abbott et al., Phys. Rev. Lett. , 101102 (2020).[2] S. E. Woosley, Astrophys. J. , 244 (2017).[3] S. E. Woosley, Astrophys. J. , 49 (2019).[4] T. Yoshida, H. Umeda, K. Maeda, and T. Ishii, Mon.Not. R. Astron. Soc. , 351 (2016).[5] G. Fragione, A. Loeb, and F. A. Rasio, Astrophys. J. , L26 (2020).[6] E. J. Farrell et al., arXiv:2009.06585.[7] T. Kinugawa, T. Nakamura, and H. Nakano, Mon. Not.R. astron. Soc, , L49 (2021).[8] M. Safarzadeh and Z. Haiman, Astrophys. J. , L21(2020).[9] B. Lie and V. Bromm, Astrophys. J. , L40 (2020).[10] A. H. Nitz and C. D. Capano, arXiv: 2010.12558.[11] M. Shibata and Y. Sekiguchi, Phys. Rev. D , 024014.[12] H. Umeda and K. Nomoto, Astrophys. J. , 385(2002). [13] A. Heger and S. E. Woosley, Astrophys. J. , 532(2002).[14] K. Takahashi, T. Yoshida, H. Umeda, K. Sumiyoshi, andS. Yamada, Mon. Not. R. Astron. Soc. , 1320 (2016).[15] K. Takahashi, T. Yoshida, and H. Umeda, Astrophys. J. , 111 (2018).[16] C. L. Fryer, S. E. Woosley, and A. Heger, Astrophys. J. , 372 (2001).[17] S.-C. Yoon, A. Dierks, and N. Langer, Astron. Astrophys. , A113 (2012).[18] S.-C. Yoon, J. Kang, and A. Kozyreva, Astrophys. J. , 16 (2015).[19] H. Uchida, M. Shibata, K. Takahashi, and T. Yoshida,Astrophys. J. , 98 (2019).[20] H. Uchida, M. Shibata, K. Takahashi, and T. Yoshida,Phys. Rev. D , 041402 (2019).[21] J. F. Hawley, Astrophys. J. , 496 (1991). [22] B. Zink, N. Stergioulas, I. Hawke, C. D. Ott, E. Schnet-ter, and E. M¨uller, Phys. Rev. D , 024019 (2007).[23] O. Korobkin, E. B. Abdikamalov, E. Schnetter, N. Ster-gioulas, and B. Zink, Phys. Rev. D , 043007 (2011).[24] K. Kiuchi, M. Shibata, P. J. Montero, and J. A. Font,Phys. Rev. Lett. , 251102 (2011).[25] E. Wessel, V. Paschalidis, A. Taokaros, M. Ruiz, and S.L. Shapiro, arXiv: 2011.04077.[26] M. Bugli, J. Guilet, E. M¨uller, L. Del Zanna, N. Buc-ciantini, and P. J. Montero, Mon. Not. R. Astron. Soc. , 108 (2018).[27] C. Reisswig, C. D. Ott, E. Abdikamalov, R. Haas, P.M¨oesta, and E. Schnetter, Phys. Rev. Lett. , 151101(2013).[28] S. Fujibayashi, M. Shibata, S. Wanajo, K. Kiuchi, K.Kyutoku, and Y. Sekiguchi, Phys. Rev. D , 083029(2020).[29] S. Fujibayashi, M. Shibata, S. Wanajo, K. Kiuchi, K.Kyutoku, and Y. Sekiguchi, Phys. Rev. D , 123014(2020).[30] M. Shibata and T. Nakamura, Phys. Rev. D ,5428(1995): T. W. Baumgarte and S. L. Shapiro, Phys.Rev. D , 024007(1998).[31] 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).[32] D. Hilditch, S. Bernuzzi, M. Thierfelder, Z. Cao, W.Tichy, and B. Br¨ugmann, Phys. Rev. D , 084057(2013).[33] S. Fujibayashi, K. Kiuchi, N. Nishimura, Y. Sekiguchi,and M. Shibata, Astrophys. J , 64 (2018). [34] S. Fujibayashi, S. Wanajo, K. Kiuchi, K. Kyutoku,Y. Sekiguchi, and M. Shibata, Astrophys. J. , 122(2020).[35] M. Shibata and K. Taniguchi, Phys. Rev. D , 084015(2008); K. Kyutoku, M. Shibata, and K. Taniguchi, Phys.Rev. D , 044049 (2010).[36] T. Yamamoto, M. Shibata, and K. Taniguchi, Phys. Rev.D , 064054 (2008).[37] S. Banik, M. Hempel, and D. Bandyophadyay, Astro-phys. J. Suppl. Ser. , 22 (2014).[38] F. X. Timmes and F. D. Swesty, Astrophys. J. Suppl. , 501 (2000).[39] M. Shibata, Phys. Rev. D , 064035 (2007).[40] Y. Sekiguchi, K. Kiuchi, K. Kyutoku, and M. Shibata,Phys. Rev. D , 064059 (2015).[41] Y. Sekiguchi, K. Kiuchi, K. Kyutoku, M. Shibata, andK. Taniguchi, Phys. Rev. D , 124046 (2016).[42] K. Kyutoku, K. Kiuchi, Y. Sekiguchi, M. Shibata, andK. Taniguchi, Phys. Rev. D , 023009 (2018).[43] S. A. Balbus and J. F. Hawley, Rev. Mod. Phys. , 1(1998).[44] J. F. Hawley, X. Guan, and J. H. Krolik, Astrophys. ,84 (2011).[45] N. I. Shakura and R. A. Sunyaev, Astron. Astrophys. ,337 (1973).[46] R. Abbott et el., arXiv:2010.14527.[47] T. Kojima et al. arXiv: 2006.03831.[48] M. Punturo et al., Class. Quantum. Grav. , 194002(2010).[49] B. P. Abbott et al. Class. Quantum. Grav. , 044011(2017).[50] W. D. Arnett, Astrophys. J. , 541 (1980).[51] M. J. Graham et al., Phys. Rev. Lett. , 251102 (2020).[52] K. Maeda and K. Nomoto, Astrophys. J.598