Merger of multiple accreting black holes concordant with gravitational wave events

Recently, advanced Laser Interferometer Gravitational-Wave Observatory (aLIGO) has detected black hole (BH) merger events, most of which are sourced by BHs more massive than $30~M_\odot$. Especially, the observation of GW170104 suggests dynamically assembled binaries favoring a distribution of misaligned spins. It has been argued that mergers of unassociated BHs can be engendered through a"chance meeting"in a multiple BH system under gas-rich environments. In this paper, we consider the merger of unassociated BHs, concordant with the massive BH merger events. To that end, we simulate a multiple BH system with a post-Newtonian $N$-body code incorporating gas accretion and general relativistic effects. As a result, we find that gas dynamical friction effectively promotes three-body interaction of BHs in dense gas of $n_\mathrm{gas}\gtrsim 10^6\,\mathrm{cm}^{-3}$, so that BH mergers can take place within $30$ Myr. This scenario predicts an isotropic distribution of spin tilts. In the concordant models with GW150914, the masses of seed BHs are required to be $\gtrsim 25~M_\odot$. The potential sites of such"chance meeting"BH mergers are active galactic nucleus (AGN) disks and dense interstellar clouds. Assuming the LIGO O1, we roughly estimate the event rates for PopI BHs and PopIII BHs in AGN disks to be respectively $\simeq 1-2\,\mathrm{yr}^{-1}$ and $\simeq 1\,\mathrm{yr}^{-1}$. Multiple episodes of AGNs may enhance the rates by roughly an order of magnitude. For massive PopI BHs in dense interstellar clouds, the rate is $\simeq 0.02\,\mathrm{yr}^{-1}$. Hence, high-density AGN disks are a more plausible site for mergers of chance meeting BHs.

Excepting the GW151226 and GW170608 events, the BH pair in each event includes a BH more massive than 30 M ⊙ . Abbott et al. (2016c) have argued that such massive BHs are unlikely to originate in metal-rich stars owing to mass loss by stellar wind. As models for the BH merger events, several binary evolution scenarios have been proposed. They include a binary of metal-free or low-metallicity stars accompanied by mass transfer or common envelope ejection (e.g. Kinugawa et al. 2014;Belczynski et al. 2016), binary evolution in a tidally distorted field (e.g. de Mink et al. 2016), binary evolution driven by fallback accretion (Tagawa et al. 2018) and dynamical interaction in dense stellar clusters (e.g. O'Leary et al. 2009;Samsing et al. 2014;Rodriguez et al. 2016). Also, BH binaries may be hardened within gas-rich environments (Escala et al. 2004;Chapon et al. 2013), especially in active galactic nucleus (AGN) disks (Kocsis et al. 2011;McKernan et al. 2012McKernan et al. , 2014Bartos et al. 2017;Stone 2017;McKernan et al. 2017). McKernan et al. (2012McKernan et al. ( , 2014 predicted the occurrence of intermediate BH mergers originating in AGN disks. McKernan et al. (2017) also considered binary formation of unassociated BHs through angular momentum exchange. Baruteau et al. (2011) investigated inward migration of massive stellar binaries hardened whitin a dense gaseous disk in the Galactic center. GW observations can provide information about component spins through measurements of an effective inspiral spin parameter χ eff , which can potentially be used to distinguish different formation channels. Isolated binary evolution does not result in a significant spin misalignment, since mass transfer and tides are to align spins with the orbital angular momentum. The GW170104 event exhibits χ eff = −0.12 +0.21 −0.30 , which disfavors spin configurations with both component spins positively aligned with the orbital angular momentum (Abbott et al. 2017a), although the less massive BH merger in GW151226 has a preference for spins with positive projections along the orbital angular momentum (Abbott et al. 2016b). The observation of GW170104 hints towards dynamically assembled binaries favoring a distribution of misaligned spins rather than near orbitaligned spins. Recently, Tagawa et al. (2015Tagawa et al. ( , 2016 have proposed mergers of unassociated BHs through a "chance meeting" in gas-rich environments, without making a priori assumption of a BH binary. They have demonstrated that a multiple stellar-mass BH system embedded in dense gas can engender mergers of BHs through gas dynamical friction and three-body interaction, which predicts an isotropic distribution of spin tilts. In this paper, we consider BH mergers by chance meetings in a multiple BH system, especially focusing on the massive BH merger events (GW150914, GW170104, and GW170814). Favorable gas-rich environments for BH mergers are provided in nuclear regions of galaxies which have density of n gas 10 7 cm −3 at 1 pc (Goodman 2003;Namekata & Umemura 2016). Another possible site is dense interstellar cloud cores of n gas = 10 5−7 cm −3 (Bergin et al. 1996;Stahler 2010), or interstellar clouds of n gas < 10 5 cm −3 (Spitzer 1978). Simulations are performed with a highly accurate post-Newtonian N -body code, where such general relativistic effects as the pericenter shift and GW emission are taken into consideration. In these simulations, the effects of gas dynamical friction and Hoyle-Lyttleton mass accretion by ambient gas are incorporated. Changing initial masses of BHs, ambient gas density, and distributions of BHs, we derive the range of BH mass that is concordant with the GW events, and thereby assess the mass of accreting gas before mergers. Also, we roughly estimate the event rates of such BH mergers both in galactic centers and in dense interstellar clouds.
The detailed description of numerical schemes is given in Tagawa et al. (2016).
The equations of motion are integrated using a fourth-order Hermite scheme (Makino & Aarseth 1992). Our simulations incorporate the effects of gas dynamical friction and gas accretion onto BHs. The general relativistic effects are dealt with post-Newtonian prescription up to a 2.5PN term (Kupi et al. 2006), where 1PN and 2PN terms correspond to pericentre shift, and a 2.5PN term to GW emission.

Setup of Simulations
The key parameters in our simulations are initial BH mass (m 0 ), initial typical extension of BH distributions (r typ ), ambient gas number density (n gas ), and accretion efficiency (ǫ). We set gas accretion rate onto each BH to be the accretion efficiency (ǫ ≤ 1) times the Hoyle-Lyttleton accretion rate (ṁ HL ), i.e., where v i is velocity of i-th BH, c s is sound speed, G is the Gravitational constant, and m H is the hydrogen mass.
The effect of radiation pressure on Hoyle-Lyttleton accretion (Watarai et al. 2000;Hanamoto et al. 2001) is incorporated as Tagawa et al. (2016). This gas accretion prescription allows super-Eddington accretion, which is verified in spherical symmetric systems (e.g. Inayoshi et al. 2016). We consider multiple BHs that are embedded in high-density gas, e.g., in galactic nuclear regions of 1 pc or in dense interstellar clouds. Then, typical extensions of BH distributions at an initial epoch (r typ ) are assumed to be from 0.01 to 1 pc. Additionally, to scrutinize dependence on ambient gas density, we consider a relatively wide range of gas density from 10 2 cm −3 to 10 10 cm −3 . We initially set up five BHs with equal mass of 20, 25, or 30 M ⊙ . Because of the uncertainty concerning the actual mass accretion rate, we vary the gas accretion efficiency ǫ in a range of 10 −3 to 1. We set BHs in a uniform gas sphere whose mass is 10 5 M ⊙ . Therefore, according to the choice of gas density, the radius of gas sphere, R gas , is changed. The gas temperature is assumed to be 1000 K (therefore c s = 3.709 km s −1 ) as Tagawa et al. (2016). Initial positions of BHs are set randomly in a x − y plane within r typ which is smaller than R gas . Initial velocity of each BH is given as the sum of a circular component and a random component. Circular velocity is given so that the centrifugal force should balance the gravity by gas in the x − y plane. In addition, we impose random velocities in the xyz space, according to the probability of a Gaussian distribution with the same dispersion as the circular velocity.
We adjudicate that two BHs merge, when their separation is less than 100 times the sum of their Schwarzschild radii. The evolution is pursued for 10 Gyr, since we consider mergers within the cosmic time. We terminate the simulation, when the first BH merger occurs.

MODELS CONCORDANT WITH GW EVENTS
Changing the set of parameters, we have simulated 264 models, of which 135 produce a binary BH merger within 10 Gyr. We have found sixteen models to match GW events, where the final BH masses fall within the estimated mass range in the observations. In Table 1, they are listed with the assumed sets of parameters. The columns are the model number, initial mass of BHs (m 0 ), ambient gas number density (n gas ), accretion ef- Table 1. Sets of parameters in which BHs merged at the masses of the GW events.   ficiency (ǫ), initial extension of BH spatial distributions (r typ ), radius of a gaseous sphere (R gas ), final masses of merged BHs (m 1 , m 2 , m 1 > m 2 ), merger time (t merge ), and the merger type in each run. Tagawa et al. (2016) scrutinized merger mechanisms in gas-rich environments. They found that gas dynamical friction is indispensable for BH mergers. First, the BH orbits contract due to gas dynamical friction, and then a subsequent merger is promoted through different mechanisms, which are classified into four types: a gas drag-driven merger (type A), an interplay-driven merger (type B), a three-body driven merger (type C), and an accretion-driven merger (type D). Figure 1 demonstrates the evolution until the first merger in Model 3 (see Table 1 for the simulation parameters), where (a) accretion rate, (b) mass and (c) velocity of a heavier BH in merged BHs, and also (d) separation of the closest pair within all BHs are shown as a function of time. Panels (c) and (d) demonstrate that the velocity decays and the separation of BHs shrinks owing to gas dynamical friction within 2 Myr. In this stage, the BH velocity oscillates between subsonic and supersonic one (the sound speed being c s = 3.709 km s −1 ), and the accretion rate intermittently reaches a super-Eddington accretion rate. In this phase, a binary forms due to energy loss by gas dynamical friction. The binary is hardened by kicking another BH through three-body interaction at around 2 Myr, which is represented by discontinuous change of the separation. Then, the BH velocity becomes highly supersonic and therefore the accretion rate is reduced to a level much lower than an Eddington accretion rate. A component in the BH binary is sometimes replaced by another one as a result of three-body interaction. Actually, such exchange occurs at 3.7 and 5.1 Myr. Three-body interaction is repeated until 13Myr, and eventually the binary merges into a massive BH due to GW radiation. Prior to the BH merger, the masses of merged BHs are enhanced by about ten M ⊙ . As shown in panel (b), most of gas accretes in an early three-body interaction phase with subsonic velocity. The mass accretion rate shortly before the merger emitting GW is reduced to less than 10 −5 Eddington accretion rate, owing to the high circular velocity of the BH binary. In practice, the final accretion rate is dependent on ambient gas density. In the other models listed in Table 1, the final accretion rate is 10 −4 Eddington accretion rate.
In Figure 2, we plot the masses of two BHs shortly before first mergers in the 135 models out of simulated 264 models, and compare them to the estimated mass range in the GW150914, GW170104, and GW170814 events. We find that the masses of two BHs are consonant to the GW150914 event in twelve models, and to the GW170104 event in two models, and to the GW170814 event in three models. Especially, Model 14 matches the GW170104  Table 1. Panels (a), (b), and (c) represent mass accretion rate in units of the Eddington mass accretion rate (ṁ E = L E /cη,η = 0.1), mass, and velocity for a heavier BH in merged BHs, respectively. Panel (d) shows the separation of the closest pair within all BHs. and GW170814 events, simultaneously. It worth noting that their merger types are type C (three-body driven mergers), except for Model 16 assuming extremely high density gas. Also, it has turned out that, in these twelve models, gas of several M ⊙ can accrete onto BHs in early three-body interaction phases. Abbott et al. (2016c) have argued that if a strong stellar wind is assumed, a BH more massive than 25M ⊙ should originate in a metal-free (PopIII) or ultra-low metal star. Even for a weak wind model, the progenitors should be of sub-solar metal abundance. Hence, the present results imply that metal poor stars are preferred as the progenitors of the GW150914 BHs. Figure 3 shows the accumulated mass on each BH before the merger. Since three-body interaction is a chaotic process, accreting mass in type C changes in a cataclysmic fashion. Supposing Bondi accretion, there must be uncertainties of ∼ 10M ⊙ in accreting mass, since the merger time can fluctuate within a factor of two according to the adopted seed random number (Tagawa et al. 2015). Taking into consideration the fact that the mass uncertainties in the observations are ∼ 7 M ⊙ , about a half of the models that match the masses in the GW150914 event may be missed. As shown in Table 1, the accretion efficiency (ǫ) in the concordant models are 0.01, except for the models assuming extremely high or low density gas (Model 1, 12, 13 and 16). Hoyle-Lyttleton-type accretion is a nonlinear function of mass, and therefore the accreting mass is a steep function of ǫ and n gas . The value of accretion efficiency is roughly determined by the balance between accretion timescale and merger timescale (Tagawa et al. 2016). In other words, the accumulated mass is regulated by these timescales. Actually, the timescales accord when the accretion efficiency is around 0.01.

Merger sites
We consider preferable sites for the present merger scenario. The first possibility is AGN disks, where the density is as high as 10 7 cm −3 and the size is as compact as 1 pc (Sirko & Goodman 2003;Burtscher et al. 2013). For a gas disk surrounding a central supermassive BH (SMBH), the Toomre Q value is estimated to be for disk temperature of 10 3 K, where M SMBH and M disk is the masses of a SMBH and an AGN disk, respectively. Hence, if M disk is lower than 10 5 M ⊙ , the disk is stabilized by the SMBH. However, a more massive disk should be stabilized by additional heating sources such as massive stars formed within the disk (Sirko & Goodman 2003). The viscous timescale of a disk is assessed by t vis ≃ 10 8 yr r 1 pc where α is the standard viscosity parameter (e.g. Umemura, Fukue, & Mineshige 1997), although the mass accretion may be flickering in ∼ 0.1 Myr (King & Nixon 2015). The AGN lifetime can be estimated by the duty cycle, P duty = N AGN t AGN /t H (z), where t AGN is the duration of a single AGN episode, N AGN is the number of AGN episodes, and t H (z) is the Hubble time at redshift z. Shankar, Weinberg, & Miralda-Escudé (2009) have derived P duty as a function of redshift and BH mass. For M SMBH = 10 7 M ⊙ , P duty ≃ 0.03 at z = 0.3 and P duty ≃ 3 × 10 −3 at z = 0. This can be translated into N AGN t AGN = 300 Myr at z = 0.3 and 41 Myr at z = 0, while N AGN t AGN = 10 Myr at z = 0.3 and 1 Myr at z = 0 for M SMBH = 10 9 M ⊙ .
BHs whose orbits are originally misaligned with AGN disks tend to be aligned due to gas dynamical friction.
The alignment timescale is estimated to be where v z is the z-component of BH velocity, h disk is the aspect ratio of an AGN disk (Goodman 2003), and h ini is the aspect ratio of an initial BH orbit against a AGN mid-plane. Since t align should be shorter than t AGN , M SMBH is constrained to be 10 7 M ⊙ . In the process of alignment, the velocity relative to the disk rotation leads to the epicyclic motion of a BH. The relative velocity is decaying due to dynamical friction, and simultaneously the circular orbit shrinks in the disk. When multiple BHs having residual reciprocal velocity interact with each other in the disk, the dynamics similar to the present simulations is expected. Also, the situation is analogous to the formation of protoplanets from planetesimals in a protoplanetary disk (Kokubo & Ida 2000).

1
Another possibility for the merger site is giant molecular clouds (GMCs).
Taking into account these timescales in the two possible sites, BH mergers should occur within 30 − 100 Myr. From Table 1, this condition requires n gas 10 6 cm −3 . Therefore, dense galactic nuclear disks and dense GMCs are potential sites for the mergers concordant with the GW events. Besides, dynamically assembled BH binaries in the present simulations predict an isotropic distribution of spin tilts without alignment with the orbital angular momentum, which is preferred to account for the misaligned spins in the GW170104 event.

Event rate in AGNs
We estimate the event rates for mergers of massive stellar-mass BHs in the first advanced LIGO observation run (LIGO O1). The horizon distance of massive BH mergers is D h ≈ 3 Gpc (z ≈ 0.3) (Belczynski et al. 2016), corresponding to a comoving volume V c ≈ 50 Gpc 3 . Here, we assess the event rates for massive BH mergers in AGN gas disks.
First, we consider remnant BHs of massive population I stars formed in a galaxy, say, PopI BHs. Due to inward migration of BHs by stellar dynamical friction, about N BH ∼ 2 × 10 4 BHs may exist within 1 pc from a SMBH in a Milky Way (MW)-sized Galaxy (Miralda-Escude & Gould 2000;Antonini 2014). To produce massive BHs with 25 M ⊙ , the initial progenitor mass is required to be 70 M ⊙ (Belczynski et al. 2016). Supposing the Salpeter initial mass function with an upper mass limit of 100M ⊙ , ∼ 20% of produced BHs are expected to be massive (f massive ∼ 0.2). Hence, the fraction of massive BH pairs is f 2 massive ∼ 0.04. Since the aspect ratio (h ini ) represents the fraction of BHs which can align to the AGN disk, the fraction of aligned BHs in t AGN is given by f align ≃ 0.2(t AGN /100 Myr) 1/4 from equation (4). In order for a merger to take place, the condition of t merge ≤ t AGN should be satisfied, where is t merge ∼ 10 Myr for n ∼ 10 7 cm −3 from the present simulations. Thus, it is required that t AGN ≥ 10 Myr and therefore N AGN ≤ P duty t H (z)/10 Myr. In the range of M SMBH ≤ 10 7 M ⊙ and 0 z 0.3, we have 3 × 10 −3 P duty 3×10 −2 (Shankar, Weinberg, & Miralda-Escudé 2009).
Then, N AGN ≤ 4 at z ∼ 0 and N AGN ≤ 31 at z ∼ 0.3. Using these assessments, the merger rate per Milky-sized galaxy is estimated to beṄ merge/gal ∼ P duty f align f 2 massive N BH /t AGN = f align f 2 massive N BH N AGN /t H (z) ≃ 10 − 20 Gyr −1 for N AGN = 1, and ≃ 30 − 300 Gyr −1 for the maximum of N AGN . From the Schechter function fit of local galaxies, the number density of MW-sized galaxies is n gal ∼ 2 × 10 6 Gpc −3 (Marzke et al. 1998). Using these values, the number of MW-sized galaxies involved in an observable volume is N gal ∼ V c n gal ∼ 1 × 10 8 . Under these assumptions, the event rate for mergers of massive PopI BHs in AGN disks in the first observing run of aLIGO is estimated to be R O1,AGN,PopI ∼Ṅ merge/gal N gal ≃ 1 − 2 yr −1 for N AGN = 1, and ≃ 3 − 30 yr −1 for the maximum of N AGN . The volumetric event rate is R vol,AGN,PopI ∼ N merge/gal n gal ≃ (2−4)×10 −2 Gpc −3 yr −1 for N AGN = 1, and ≃ 0.06 − 0.6 Gpc −3 yr −1 for the maximum of N AGN .
Next, we consider remnant BHs of population III stars (PopIII BHs). Although there are large uncertainties, roughly ten PopIII BHs are possibly born in a minihalo of 10 5 − 10 6 M ⊙ (Susa et al. 2014;Valiante et al. 2016). In this case, ∼ 10 6 PopIII BHs are expected to exist in a MW-sized galaxy (Ishiyama et al. 2016). Then, if the ratio of PopIII BHs to PopI BHs is assumed to be constant in a whole galaxy, the number of BHs in central subparsec regions is N BH ∼ 2 × 10 2 . Besides, if taking into consideration the possibility that BHs within ∼ 10 pc can migrate into subparsec regions, the number of BHs at 1 pc can increase by about one order of magnitude (Miralda-Escude & Gould 2000). So, we suppose N BH ∼ 2 × 10 3 PopIII BHs exist in an AGN disk in a MW-sized galaxy. We assess the fraction of massive ones in all PopIII BHs to be f massive ∼ 0.5 (Heger & Woosley 2002;Susa et al. 2014). Then,Ṅ merge/gal ∼ P duty f align f 2 massive N BH /t AGN ≃ 6 − 10 Gyr −1 for N AGN = 1, and ≃ 20 − 200 Gyr −1 for the maximum of N AGN . Under these assumptions, we estimate the event rate for mergers of massive PopIII BHs in AGN disks to be R O1,AGN,PopI ∼Ṅ merge/gal N gal ≃ 1 yr −1 for N AGN = 1, and ≃ 2−20 yr −1 for the maximum of N AGN . The volumetric event rate is R vol,AGN,PopI ∼ N merge/gal n gal ≃ 2 × 10 −2 Gpc −3 yr −1 for N AGN = 1, and ≃ 0.04 − 0.4 Gpc −3 yr −1 for the maximum of N AGN .

Event rate in GMCs
Here, we estimate the event rates for BH mergers in GMCs. A MW-sized galaxy contains ∼ 10 8 BHs in the volume of ∼ 100 kpc 3 (Remillard & McClintock 2006), where the fraction of massive BHs is 0.2 as discussed in previous section. There are ∼ 1000 GMCs in a galaxy and they occupy the volume 10 −3 kpc 3 (Ruffle 2006). Hence, we can assume that ∼ 200 massive BHs reside in GMCs. Considering the fact that velocity dispersion of PopI massive stars is ∼ 20 km/s (Binney & Merrifield 1998;Nordstrom et al. 2004) and the escape velocity of GMCs of ∼ 10 km/s (Dobbs et al. 2011;Dale et al. 2012), ∼ 40 percent of PopI BHs can be captured by GMCs. According to probability distributions, about 3 GMCs possess more than two massive BHs. Also, the volume filling factor of dense cores in GMCs is f core ∼ 0.05 (Bergin et al. 1996). Besides, stars which leave BHs more massive than 25 M ⊙ should be metal poor (≤ 0.3 solar metallicity) and they should be low velocity dispersion (Nordstrom et al. 2004). Most of such stars exist in outer galaxies of 10 kpc (Martinez-Medina et al. 2017), where the stellar mass is ∼ 0.1 of the total galactic stellar mass. Since BHs are redistributed in the dynamical time of a galaxy t dyn ∼ 100 Myr, the merger rate in a MW-sized galaxy isṄ merge/gal ∼ 3 × 0.1f core /t dyn ≃ 0.2 Gyr −1 . Under these assumptions, the event rate for mergers of PopI BHs in GMCs is estimated to be R O1,GMC,PopI ∼Ṅ merge/gal N gal ≃ 0.02 yr −1 and R vol,GMC,PopI ∼Ṅ merge/gal n gal ≃ 3 × 10 −4 Gpc −3 yr −1 .

CONCLUSIONS
In this paper, we have considered the mergers of unassociated BHs through a chance meeting in gas-rich environments. To elucidate the merger condition concordant with the recently detected gravitational wave events, we have conducted highly accurate post-Newtonian N -body simulations on a multiple BHs system embedded in dense gas, incorporating dynamical friction, Hoyle-Lyttleton mass accretion, and general relativistic effects such as pericentre shift and gravitational wave emission. Consequently, we have found the following: (1) Gas dynamical friction works effectively to promote three-body interaction of BHs in dense gas of n gas 10 6 cm −3 . Eventually, mergers are caused within 30 Myr. This scenario predicts an isotropic distribution of spin tilts, which is compatible with the spin misalignment seen in the GW170104 event.
(2) Before BH mergers, gas of several M ⊙ accretes onto each BH. However, gas accretion takes place predominantly during early three-body interaction phases, and the final mass accretion rates shortly before GW emmision are 10 −4 Eddington accretion rate. Thus, the electromagnetic counterparts of GW events might not be so luminous.
(3) We have found sets of model parameters concordant with the massive BHs detected in the GW events. In the concordant models, the initial extension of BH distributions is smaller than 1 pc. To account for the GW150914 event, the masses of seed BHs are required to be 25M ⊙ . Hence, metal poor stars are preferred as the progenitors of the GW150914 BHs.
(4) We have roughly estimated the event rates by the first observing run of LIGO advanced detectors. The event rates for massive PopI BHs and PopIII BHs in AGN disks are assessed to be ≃ 1 − 2 yr −1 and ≃ 1 yr −1 , respectively. If multiple episodes of AGNs are taken into consideration, the rates can be enhanced by roughly an order of magnitude. For massive PopI BHs in dense interstellar clouds, the rate is ≃ 0.02 yr −1 . Hence, high-density AGN disks are a more plausible site for mergers of chance meeting BHs.
In the present simulations, we have assumed a fairy simple configuration of matter. However, taking realistic situations into consideration, we should construct a more concrete model of gas distributions in a dense cloud/disk and gravitational potential, including stellar distributions and a central supermassive black hole. Also, the back reaction due to gas dynamical friction may alter the BH dynamics. These effects will be explored in the future analysis.
We thank the anonymous referee for useful comments. Numerical computations and analyses were carried out on Cray XC30 and computers at Center for Computational Astrophysics, National Astronomical Observatory of Japan, respectively. This research is also supported in part by the European Research Council under the European Unionfs Horizon 2020 Programme, ERC-2014-STG grant GalNUC 638435 and Interdisciplinary Computational Science Program in Center for Computational Sciences, University of Tsukuba, and Grant-in-Aid for Scientific Research (B) by JSPS (15H03638).