Impact of initial mass functions on the dynamical channel of gravitational wave sources
MMNRAS , 1–10 (2019) Preprint 26 January 2021 Compiled using MNRAS L A TEX style file v3.0
Impact of initial mass functions on the dynamical channelof gravitational wave sources
Long Wang, , (cid:63) Michiko S. Fujii, Ataru Tanikawa Department of Astronomy, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan RIKEN Center for Computational Science, 7-1-26 Minatojima-minami-machi, Chuo-ku, Kobe, Hyogo 650-0047, Japan Department of Earth Science and Astronomy, The University of Tokyo, Japan
Accepted –. Received –; in original form –
ABSTRACT
Dynamically formed black hole (BH) binaries (BBHs) are important sources of gravitational waves (GWs). Globularclusters (GCs) are one of the major environments to produce such BBHs, but the total mass of the known GCs issmall compared to that in the Galaxy, thus the fraction of BBHs formed in GCs is also small. However, this assumesthat GCs contain a canonical initial mass function like that in the field stars. This might not be true because severalstudies suggest that extreme dense and metal-poor environment can result in top-heavy IMFs, where GCs may comefrom. Although GCs with top-heavy IMFs were easily disrupted or have become dark clusters, the contribution tothe GW sources can be significant. Using a high-performance and accurate N -body code, petar , we investigate theeffect of varying IMFs by carrying out four star-by-star simulations of dense GCs with the initial mass of 5 × M (cid:12) and the half-mass radius of 2 pc. We find that the BBH merger rate does not monotonically correlate with the slopeof IMFs. Due to a rapid expansion, top-heavy IMFs lead to less efficient formation of merging BBHs. The formationrate continuously decreases as the cluster expands due to the dynamical heating caused by BHs. However, in starclusters with a top-heavier IMF, the total number of BHs is larger, and therefore the final contribution to mergingBBHs can still be more than clusters with the standard IMF, if the initial cluster mass and density is higher thanthe model we used. Key words: methods: numerical – galaxies: star clusters: general – stars: black holes
After LIGO/VIRGO detected gravitational wave (GW)events from mergers of stellar-mass black holes (BHs) andneutron stars (NSs) (Abbott et al. 2019, 2020), many stud-ies have addressed on the origins of these events: isolatedbinaries through common envelope evolution (e.g. Giacobbo& Mapelli 2018; Belczynski, et al. 2020), through chemicallyhomogeneous evolution (e.g. Marchant et al. 2016; Mandel &de Mink 2016) and stable mass transfer (e.g. Kinugawa et al.2014; Tanikawa et al. 2020), hierarchical stellar systems (e.g.Antonini et al. 2014), open clusters (e.g. Ziosi et al. 2014; Ku-mamoto et al. 2019; Di Carlo et al. 2019; Banerjee 2020a),and galactic centers (e.g. O’Leary et al. 2009). Dense stellarsystems like globular clusters (GCs) are considered as one ofthe major environment to form GW progenitors, binary blackholes (BBHs), via few-body dynamical interactions (Porte-gies Zwart & McMillan 2000; Downing et al. 2010; Tanikawa2013; Bae et al. 2014; Rodriguez et al. 2016a,b; Fujii et al.2017; Askar et al. 2017; Park et al. 2017; Samsing et al. 2018;Hong et al. 2020). Although this dynamical channel can effi- (cid:63)
E-mail:[email protected] ciently produce GW events, the total contribution seems tobe less than the events driven by the binary stellar evolution,because the total mass of GCs is a small fraction of the totalgalactic stellar mass. However, this assumes that the initialmass function (IMF) of the GCs is the same as that of thefield stars, but it might not be the true because we do nothave sufficient observational constraints on that yet.Recently, several observational evidences indicate that ex-treme dense star forming region may have top-heavy IMFs,such as the Arches cluster in the Galactic center and 30 Do-radus in the Large Magellanic Cloud, where the heavy-endof IMFs have exponent, α ≈ − . ∼ − .
9, (Lu et al. 2013;Schneider, et al. 2018; Hosek, et al. 2019). Therefor, it isnatural to expect that the old and massive GCs in the Milk-way may also contain top-heavy IMFs since its birth envi-ronment is very different from the present-day one. On thegalactic scale, Zhang et al. (2018) found that for high-redshift( z ∼ −
3) star-burst galaxies, a top-heavy (integrated) IMFis necessary to explain their star formation rate.On the other hand, there is a puzzle related to the phe-nomenon of multiple stellar populations in GCs, called the“mass budget problem”. Observations show that several GCshave a more than half of enriched stellar populations (Milone © a r X i v : . [ a s t r o - ph . H E ] J a n Long Wang et al. et al. 2017). The stellar evolution model to explain the forma-tion of element-enriched stars, especially the AGB scenario,cannot produce sufficient materials to form such a large frac-tion of young populations (Bastian & Lardo 2018, and ref-erences there in). A top-heavy IMF is a natural solution forthis (e.g. Wang et al. 2020).With top-heavy IMFs, the strong wind mass loss from mas-sive stars in the first 100 Myr significantly affects the den-sity of the systems. Later on, the massive stars leave a largenumber of BHs in the star clusters, and they will also havea strong impact on the long-term evolution of the clusters.Indeed, it has been shown that GCs with top-heavy IMFsare much easier to expand and be disruptted (Chatterjee,Rodriguez, & Rasio 2017; Giersz, et al. 2019; Wang 2020;Weatherford et al. 2021).Therefore, the dynamical evolution of GCs provides astrong constraint on the shape of IMFs. By using semi-analytic models, Marks et al. (2012) suggested that the shapeof IMFs in GCs might depend on the metallicity and initialcloud density. Their model, however, ignored the dynamicalimpact of BHs, and thus they overestimated the slopes ofIMFs for observed GCs. By comparing with scaled N -bodymodels, Baumgardt & Sollima (2017) also argued that 35 ob-served dense GCs might not have top-heavy IMFs. However,GCs with top-heavy IMFs might have existed in the pastbut have already disappeared or have become dark clusters(Banerjee & Kroupa 2011). Thus, they cannot be directly ob-served today. Considering the large number of BHs existedthere, the contribution of BBHs mergers might be significant.Chatterjee, Rodriguez, & Rasio (2017), Giersz, et al. (2019)and Weatherford et al. (2021) have performed Monte-Carlosimulations of GCs to study the effect of top-heavy IMFs onthe survival of GCs and the BBH mergers. However, Monte-Carlo method is not fully tested in the condition of top-heavyIMFs, where a large fraction of BHs exist. Rodriguez et al.(2016) compared the Monte-Carlo (the cmc code; e.g. Joshi,Rasio, & Portegies Zwart 2000) and the direct N -body meth-ods (the nbody6++gpu code; Wang et al. 2015) for million-body simulations of GCs. In this comparison, an abnormaloscillation of the core radius was found in the cmc simula-tion. Later, Rodriguez et al. (2018) showed that when BHnumber is large, the Joshi, Rasio, & Portegies Zwart (2000)code fails to produce the correct core evolution.There are also studies using direct N -body methods, whichmust be more accurate than the Monte-Carlo simulations.Wang (2020) studied the effect of top-heavy IMFs on thesurvival of star clusters, but use a simplified two-componentmodels without stellar evolution. Haghi et al. (2020) carriedout a group of N -body models, but without a realistic num-ber of stars like in the GCs. These models properly integrateindividual stellar orbits and obtained correct dynamical evo-lution of GCs, but they cannot be used to study the BBHmergers. Therefore, it is necessary to carry out accurate star-by-star N -body simulations to study the GCs with top-heavyIMFs.In this work, we carry out four models of GCs with the suf-ficient mass (5 × M (cid:12) ) and density (the half-mass radiusis 2 pc). Different slopes of IMFs are used. In Section 2.1,we describe the N -body tool, petar , used in this study. Theinitial conditions of models are shown in Section 2.2. We de-scribe the BBHs mergers and the dynamical evolution of the M ]010203040 F i n a l m a ss a t M y r [ M ] Have SN kick
Figure 1.
The initial-final mass relation of stars by applying thestellar evolution model from sse . The yellow points indicate theremnants natal kick velocity of neutron stars and black holes afterthe supernova. Due to the fallback treatment, massive BHs haveno kick velocity. models in Section 3. Finally, we discuss our results and drawconclusions in Section 4 and 5.
We use the petar code (Wang et al. 2020b) to provide the N -body models of GCs. This is a hybrid code that combinesthe particle-tree particle-particle (Oshino, Funato, & Makino2011) and slowdown algorithm regularization (Wang, Nita-dori & Makino 2020a) methods. Such combination reducesthe computing cost compared to the direct N -body methodwhile keep sufficient accuracy to deal with close encountersand few-body interactions. The code is based on the fdps framework, which can achieve a high performance with themulti-process parallel computing (Iwasawa et al. 2016; Iwa-sawa, et al. 2020; Namekata et al. 2018).The single and binary stellar evolution packages, sse/bse ,are implemented in petar (Hurley, Pols, & Tout 2000; Hur-ley, Tout, & Pols 2002). We adopt the updated version of sse/bse from Banerjee et al. (2020b). The simulations usethe semi-empirical stellar wind prescriptions from Belczynskiet al. (2010), the “rapid” supernova model for the remnantformation and material fallback from Fryer et al. (2012) to-gether with the pulsation pair-instability supernova (PPSN;Belczynski, et al. 2016). In Figure 1, we show the relationbetween zero-age main-sequence (ZAMS) masses and finalmasses of stars from 0 .
08 to 150 M (cid:12) for Z = 0 . M (cid:12) obtain the final mass of BHs fixed to 40 . M (cid:12) . Thishas a significant impact on the mass ratio distribution ofBBHs, as shown in section 3. On the other hand, the fallbackmechanism results in zero natal kick velocity for massive BHswith ZAMS mass of > M (cid:12) . They can retain in the clustersafter supernovae, while low mass BHs and most of neutronstars can gain high velocity kicks and immediately escapefrom the system. MNRAS000
08 to 150 M (cid:12) for Z = 0 . M (cid:12) obtain the final mass of BHs fixed to 40 . M (cid:12) . Thishas a significant impact on the mass ratio distribution ofBBHs, as shown in section 3. On the other hand, the fallbackmechanism results in zero natal kick velocity for massive BHswith ZAMS mass of > M (cid:12) . They can retain in the clustersafter supernovae, while low mass BHs and most of neutronstars can gain high velocity kicks and immediately escapefrom the system. MNRAS000 , 1–10 (2019) mpact of IMFs on GW sources Table 1.
The name of N -body models with corresponding α ofIMFs and the number of stars ( N )Model A1.5 A1.7 A2.0 A2.3 α -1.5 -1.7 -2.0 -2.3N 182306 312605 581582 854625 To specifically investigate how the shape of IMFs affect theformation and evolution of BBHs in GCs, we carry out fourmodels by varying the α in the Kroupa (2001) IMF with amulti-component power-law shape: ξ ( m ) ∝ m α i α = − . , . (cid:54) m/M (cid:12) < . α = − . , . (cid:54) m/M (cid:12) < . α , . (cid:54) m/M (cid:12) < . . (1)The value of α and corresponding model names are listedin Table 1. A2.3 has the canonical value of α while all theothers correspond to different degrees of top-heavy IMFs. Themaximum value of α is chosen to be 1.5. We use the full massrange of IMF, where the minimum and the maximum massesof stars are 0 .
08 and 150 M (cid:12) , respectively. For the metallicity,which determines the stellar evolution (mass loss) of stars, weadopt a typical value for GCs, Z = 0 . × M (cid:12) and the initial half-massradius, r h, , to be 2 pc. Thus, our models have a mass anddensity as high as those typical for observed GCs. In fact,the density of our model is higher than that of the previouslargest DRAGON GC model (Wang et al. 2015). Thus, oursimulations are still time-consuming even using the petar .As a start, we assume no primordial binaries in these mod-els. This simplifies the discussion since we can purely focuson the dynamical formation of BBHs. On the other hand, itis easier to perform the models. We do not apply the galacticpotential either, but we remove unbounded stars when theyreach more than 200 pc far from the cluster centre. We do not intend to create realistic GC models with primor-dial binaries and galactic tidal field, but focus on the theo-retical studies of the BBHs in GCs. Thus, we only evolve ourmodels up to 3 Gyr instead of Hubble time, using the avail-able computing resource. However, this already covers about10 initial half-mass relaxation time of the systems, thus it issufficient for analysis.
In our simulations, we do not model BBH mergers via theeffect of general relativity. We instead detect all BBHs thatare potential GW mergers in the post-process using the snap-shots of the simulations. We apply the following two steps toobtain the BBH mergers.
Table 2.
The number of BBH merger candidates ( N cand ), mergersinside GCs ( N in ) and escaped mergers ( N out ) up to the 3 Gyrevolution of GCs.Model A1.5 A1.7 A2.0 A2.3 N cand
20 10 30 43 N in N out Firstly, we estimate merger timescale for each binary basedon the formula provided by Peters (1964): t gw = 1219 c β (cid:90) e e / [1 + (121 / e ] / (1 − e ) / d ec = a (1 − e )[1 + (121 / e ] / e / β = 645 G m m ( m + m ) c . (2)By calculating t gw of BBHs in the snapshots with the timeinterval of 0 . − t merge ≡ t gw + t < t is the physical time of the cluster when t gw is calculated.However, these candidates may not actually merge becauseBBHs with high eccentricities can be perturbed by their sur-rounding stars, if they are still in the cluster. The orbits ofBBHs can also dramatically change after a strong few-bodyinteraction. If the eccentricity decreases, t gw can significantlyincreases. Thus, in the second step, we look the evolutionof each candidate and check whether they escape from thesystem or merge before their t gw increase.The numbers of candidates and confirmed mergers for eachmodel are shown in Table 2. We can notice that the numberof candidates does not monotonically depends on the α ofIMFs. More top-heavy the IMF is, less N cand tend to be,except for that of the A1.5 model (the most top-heavy one).The A1.5 model has more N cand than that of the A1.7. Wewill explain the reason in Section 3.2.Figure 2 shows the evolution of semi-major axis ( a ) of allBBHs detected in the snapshots of our simulations. The can-didates, inner and escaped BBH mergers are shown together.In this figure, we can identify two types of BBHs, the softand hard ones separated by a ≈ AU. Due to the Heggie-Hills law (Heggie 1975; Hills 1975), soft binaries are easily dis-rupted by close encounters, but hard binaries become tighter.This is seen in the Figure 2. Soft BBHs randomly appear andvanish, while hard binaries continuously evolve harder follow-ing a clear trace of decreasing a . Typically, only one or twohard BBHs appear at one time because binary-binary en-counters tend to break softer binaries. Once their a becomessmall enough, a strong encounter can eject them out of theGCs, and they disappear in the Figure. Then new hard BBHsform and repeat the same process.From the A2.3 model, we can identify a clear trend thatBBHs with large masses prefer to form first. Once they es-cape, BBHs with lower masses form and escape one by one.This feature results in a clear difference in the mass ratio dis-tribution of BBH components depending on IMFs, as shownin Figure 3. As more top-heavy IMFs have larger fractions ofstars with the ZAMS masses above 100 M (cid:12) , more equal-mass MNRAS , 1–10 (2019)
Long Wang et al.
Figure 2.
The evolution of semi-major axes of all BBHs in the fourmodels. The gray scale indicates the total masses of BBHs. Theyellow “x” represents the merger candidates. The red triangles andred crosses represents escaped mergers and mergers inside GCs.
BHs of 40 . M (cid:12) form as a result of the PPSN. These mostmassive BHs form binaries first. Thus the BBH mergers intop-heavy IMFs tend to have q closer to unity.Another major feature shown in Figure 2 is the increase ofminimum semi-major axis (maximum binding energy), espe-cially in the A1.5 model. As a result, it becomes more diffi-cult to form tight BBHs in the later evolution of GCs. Thisis reflected on the N cand as a function of time. Most of thecandidates form in the first 1000 Myr in the A1.5 model, and Figure 3.
The cumulative distribution of mass ratio q of two com-ponents in BBH candidates (lines) and confirmed (triangles) merg-ers.
50 60 70 80 m BBH
Figure 4.
The cumulative distribution of the masses of BBHscandidates (solid curves) and confirmed ones (triangles) the formation rate significantly decreases after 2000 Myr. InSection 3.3, we explain why the minimum semi-major axisincreases by analyzing the escape velocities of GCs.In Figure 4, we show the cumulative distribution of themasses of BBH merger candidates and confirmed ones. Thereare only a few confirmed mergers with m BBH < M (cid:12) inthe A2.3 model. In A1.5 models, on the other hand, most ofmergers have m BBH = 81 M (cid:12) with the mass of each compo-nent being m BBH = 40 . M (cid:12) due to the PPSN. Thus, clusterswith a more top-heavy the IMF form more massive mergers.In Figure 5, we show the cumulative distribution of t gw and t merge for BBH merger candidates and confirmed ones for allmodels. While mergers inside GCs have a short merger time( t gw < t merge < t gw > t merge =1 −
10 Gyr). Therefore, the escaped mergers from all modelswith different IMFs can be detected today.
MNRAS , 1–10 (2019) mpact of IMFs on GW sources T gw [Gyr]010203040 C u m u l a t i v e nu m b e r A1.5A1.7A2.0A2.3EscIn 10 T merge [Gyr]010203040 C u m u l a t i v e nu m b e r A1.5A1.7A2.0A2.3EscIn
Figure 5.
The cumulative number of BBH merger candidates(solid curves), mergers inside GCs (crosses) and escaped mergers(triangles) vs. t gw (upper panel) and t merge (lower panel). The stellar evolution, especially wind mass loss of massivestars in the first 100 Myr, has a significant impact on thelater-on and long-term evolution of GCs. Especially, the massloss with a top-heavy IMF is more intense due to the largerfraction of massive stars. This significantly affects the cen-tral density of GCs. As shown in Figure 6, the core radii( R c ) in the first 100 Myr become larger with more top-heavyIMFs. Once most massive stars become compact remnants,such the strong wind mass loss stops. Then, R c starts toshrink and finally core collapse happens. Figure 6 shows thatthe A1.5 model contracts the deepest during the core col-lapse, probably due to the larger masses of the BHs, and thelonger distance of sinking. This causes the densest core be-tween 350 −
450 Myr compared to those of the other models.Such high density might be the reason why the A1.5 modelhas a more number of BBH merger candidates compared tothat of the A1.7 model. Figure 2 shows that many of BBHmerger candidates appear around 400 Myr in the A1.5 model.This is consistent with the feature of core collapse.
Due to the larger masses of BHs compared to those of stars,BHs have a strong impact on the long-term evolution of starcluster (Breen & Heggie 2013). In star clusters including BHs,massive BHs sink into the cluster center due to the dynami- R c [ p c ] m B H [ M ] A1.5A1.7A2.0A2.30 100 200 300 400 500 600Time[Myr]10 n B H , c Figure 6.
The evolution of core radii ( R c ); averaged masses( (cid:104) m BH , c (cid:105) ) and numbers of BHs ( N BH , c ) inside R c . BBHs arecounted as single objects. cal friction (mass segregation) and form a dense subsystem.After that, the core collapse of the BH subsystem happensand drive the formation of BBHs, which heats the systemsvia few-body interactions. As a result of the energy-balancedevolution, the halo composed of light stars expands. Whenmore BHs exist, the expansion is faster (see also in Giersz, etal. 2019; Wang 2020). This can be identified in Figure 7, inwhich the evolution of the half-mass radii of non-BH objects( R h , noBH ) and BHs ( R h , BH ) are compared for all models.In the first 300 Myr, R h , BH decreases due to the mass seg-regation, and later it increases due to the BH heating. On theother hand, R h , noBH always increases. The stronger stellar-wind mass loss in more top-heavy models result in larger R h , BH and R h , noBH in a short time ( <
100 Myr). This is sim-ilar to the evolution of the core radii shown in Figure 6. As aresult of stronger BH heating during the long-term evolutionafter 300 Myr, R h , BH and R h , noBH increase faster in the moretop-heavy models.According to Heggie-Hill law, the boundary semi-majoraxis between soft and hard BBHs depend on the masses ofBHs and local velocity dispersion: a s / h = Gm m (cid:104) m (cid:105) σ (3)We check the density ( ρ BH ) of and 3D velocity dispersion( σ BH ) of the BHs within 10% Lagrangian radii of all BHs and R c of the BH subsystem. The results are shown in Figure 8.We find that ρ BH inside R lagr , BH , varies, as the slope ofthe IMF changes, i.e., models with more top-heavy IMFs havelower central densities. This is consistent with the behaviour MNRAS , 1–10 (2019)
Long Wang et al. R h A1.5A1.7 A2.0A2.3 BHNon-BH
Figure 7.
The evolution of half-mass radii of non-BH objects( R h , noBH ; dashed curves) and BHs ( R h , BH ; solid curves). of R h , BH (see Figure 6). In contrast, ρ BH inside R c are verysimilar among all models, different from the behaviour on alarger distance scale.We can divide the models into two groups, A1.5/A1.7 andA2.0/A2.3, based on the values of σ BH within R lagr , BH , or R c . There is a gap of ∼ σ BH between the twogroups. Models with more top-heavy IMFs have smaller σ BH ,thus their a s / h are larger due to Equation 3. This is consistentwith the hard-soft boundaries shown in Figure 2.As clusters expand, σ BH continues decreasing for all mod-els, and a s / h increases. This is identified in the A1.5 and A1.7models in Figure 2. In the A2.0 and A2.3 models, on the otherhand, we can clearly see that the masses of BBHs decreasewith time. This balances the effect of decreasing σ BH , andthus the change of a s / h is not obvious for these two models.The key quantity that controls the minimum semi-majoraxis of binaries is the central escape velocity of star clusters.As the semi-major axis of a hard binary shrinks, the binarycan gain larger kinetic energy (larger kick to the center-of-mass velocity) after a strong encounter. This process finallycauses the ejection of the binary from the center of the clusterafter a strong encounter. If the ejection velocity is below theescape velocity, the dynamical friction bring it back to thecenter. The binary can continue to encounter with othersand become harder, until the next strong encounter ejectsit again. Therefore, the escape velocity of star clusters limitsthe minimum semi-major axis of hard binaries that can bereached by few-body encounters.Based on the mechanism desribed above, we can estimateGW merger timescale ( t gw ) of an ejected BBH, which shouldbe an indicator for t gw of BBHs in a star cluster. Accordingto Equation 2 (Peters 1964), t gw ∝ a /m , where a esc is thesemi-major axis of an ejected BBH, and we assume m = m for the BBH. Since an ejected BBH has an internal velocity of ∼ v esc , a esc = Gm / (2 v ), and t gw ∝ m /v . Eventually,an ejected BBH has larger t gw when it contains heavier BHs,and when it is ejected from a star cluster with a smallerescape velocity.In Figure 9, we estimate the central escape velocity of thePlummer model, v esc = (cid:114) × . G MR h , (4) where M is the total mass of the system. Due to the larger R h in models with top-heavy IMFs, their v esc is smaller. More-over, their BHs are heavier in models with more top-heavyIMFs. Therefore, it is more difficult to form BBH mergersthere, although the number of BHs is much larger. On theother hand, v esc decreases, as the cluster expand. This ex-plains why the minimum semi-major axis increases as shownin Figure. 2 especially for the A1.5 model, because the massesof BBHs keeps a similar value up to 3 Gyr. In the case of theA2.3 model, masses of BBHs decrease, while the minimumsemi-major axis does not show a strong increasing trend.Since models with top-heavy IMFs have a stronger expan-sion and keep heavier BHs inside the clusters, the efficiencyof forming BBH mergers decrease faster. Therefore, with thesame initial total mass and size, the BBH merger rate is ac-tually lower in the GCs with top-heavy IMFs. If the clusterwith a top-heavy IMF was under a strong tidal field, the de-crease of v esc would be faster, thus it would be more difficultfor the cluster to generate merging BBHs. However, if thecluster was initially denser and more massive, it could sur-vive up to Hubble time and continue to generate mergingBBHs all the time. In such a case, the cluster with a top-heavy IMF contains a large number of BHs, and thereforethe total contribution of BBH mergers can be significantlylarger than that in a GC with the same mass and size but atop-light IMF. Such ‘dark clusters’ has been studied in theMonte-Carlo simulations (Weatherford et al. 2021). The stellar wind and BH heating drive the mass loss of starclusters in the early and later phase, respectively. Figure 10(upper panel) shows the time evolution of the total massesof BH and non-BH components in the four models. After100 Myr, the A1.5 model loses about 60% of the initial massdue to stellar-wind mass loss, while the A2.3 model loses only30%. Therefore, the impact of the stellar evolution is moresignificant in star clusters with more top-heavy IMFs. Suchstrong mass loss later drives the fast expansion of the systemas shown in Figure 6 and 7. Since our models do not havetidal field but simply remove unbound stars above 200 pc, theA1.5 model still survives and keeps about 34% initial massat 3 Gyr. In the realistic condition, where the galactic tidalfield play a role, the clusters with top-heavy IMFs tend todissolute much faster than our current models (Wang 2020).Breen & Heggie (2013) and Wang (2020) have found thatthe mass fraction of BHs in the system ( M BH /M ) evolve de-pending on the initial ratio and the tidal field in the followingmechanism. Since BHs are centrally concentrated, they arenot directly affected by the tidal field. Thus the tidal evapo-ration of BHs can be neglected. As mentioned in Section 3.3,strong close encounters between hard BBHs and intruderscan eject BHs from the cluster. This is the major scenariothat causes the mass loss of BHs.On the other hand, light stars in the halo are truncatedby the tidal field. When BH heating exists, this process isaccelerated (Breen & Heggie 2013; Giersz, et al. 2019; Wang2020). In the lower panel of Figure 10, we can identify thetwo different evolution trend of M BH /M . In the A2.3 model, M BH /M decreases with time, and finally most of BHs will beejected out from the clusters. The core collapse of light starswill happen, and a GC with a dense core will appear. MNRAS , 1–10 (2019) mpact of IMFs on GW sources B H [ M / p c ] Inside R lagr, 10% A1.5A1.7 A2.0A2.3
Inside R c Time[Myr] B H [ p c / M y r ] Time[Myr]
Figure 8.
The evolution of density of BHs ( ρ BH ) and 3D velocity dispersion of BHs ( σ BH ) within 10% Lagrangian radii ( R lagr , BH , )and R c Time[Myr] v e s c [ p c / M y r ] A1.5A1.7A2.0A2.3
Figure 9.
The evolution of the central escape velocity of the mod-els
In contrast, M BH /M in the A1.5 and A1.7 models increaseswith time. This means that light stars will finally evaporatefrom such clusters and that dark clusters form. In a strongtidal field, the mass loss of light stars is faster and such pro-cess can be accelerated. The A2.0 model is in the transitionregion.With a two-component simplified model, Wang (2020)found that the mass loss rate of BHs, M BH ( t ) /M BH (0), sim-ply depends on the mass segregation time of heavy compo-nents ( t ms ) in isolated clusters, but the dependence becomescomplex if a tidal field exists. In our models with stellar evo- M [ M ] A1.5A1.7 A2.0A2.3 BHNon-BH10 Time[Myr]0.00.20.40.6 M B H / M a ll Figure 10.
The evolution of total masses of BH and non-BH com-ponents (upper panel) and mass ratio between these two (lowerpanel). MNRAS , 1–10 (2019)
Long Wang et al. T ms ]0.60.70.80.91.0 M B H [ t ] / M B H [ ] A1.5A1.7A2.0A2.3
Figure 11.
The evolution of total mass of BHs, M BH ( t ). The massis normalized by using the value at 100 Myr. Time is in the unitof t ms (eq. 5). lution and IMFs, even without tidal field, M BH ( t ) /M BH (0)does not simply depend on t ms as shown in Figure 11. In two-component models, the definition of t ms is simple: t ms = m m t rh , , (5)where the relaxation time of light component has the form t rh , = 0 . N / R / m / G / ln Λ . (6)However, with a mass function, the precise definition of t ms is difficult to determine. Here we use the averaged mass ofnon-BH and BH components as a replacement of m and m in Equation 5 and 6. Theoretical studies have shown thatwhen multi components exist, the diffusion coefficients needto be properly averaged to obtain the correct relaxation time(Spitzer & Hart 1971; Antonini & Gieles 2019; Wang 2020).Thus, using averaged mass is not accurate to calculate t rh , .Especially, since the mass functions are different among thefour models, there probably need to introduce a correctionfactor ψ to t rh , , similar to the case in the two-componentmodels done in Wang (2020). Moreover, the stellar evolu-tion introduces another complexity. Thus, it is reasonable tosee that M BH ( t ) /M BH (0) does not simply depends on t ms .Instead, Figure 11 shows that the mass loss rate of BHs isfaster when IMF is more top-heavy. This again indicates thatstar clusters with top-heavy IMFs evolve faster. Breen & Heggie (2013) established a theory to describe thelong-term evolution of star clusters with BH subsystems. Thekey idea is based on the energy-balanced evolution of starclusters (H´enon 1961, 1975). After the BH subsystem formin the center of star clusters due to the mass segregation,BHs drive the binary heating process and provide the energyto support the whole system. The H´enon’s principle suggeststhat the energy flux from the center (BH subsystem) shouldbalance the one required by the global system. This can bedescribed by (Eq. 1 in Breen & Heggie 2013) Et rh ≈ k E BH t rh , BH , (7) where E BH and E are the total energy of the BH subsys-tem and the global system respectively; t rh and t rh , BH arethe two-body relaxation times measured at R h and R h , BH ,respectively.This relation constrains the behaviour of BHs, i.e., the den-sity of the BH subsystem, the formation rate of BBHs and theescape rate of BHs are controlled by the global system (lightstars), not the BH subsystem itself. By extending the Breen &Heggie (2013) theory to top-heavy IMFs, Wang (2020) foundthat when a large fraction of BHs exist, the energy balanceis different from the description of Equation 7. To properlymeasure the energy balance between the central BH subsys-tem and the global system, the correction factor ψ is requiredto define the relaxation time more appropriately (as discussedin Section 3.4). In Wang (2020), ψ is defined as ψ = (cid:80) k n k m /v k (cid:104) n (cid:105)(cid:104) m (cid:105) / (cid:104) v (cid:105) , (8)where n k , m k and v k are the number density, the mass of oneobject and the mean square velocity of the component k , (cid:104) n (cid:105) , (cid:104) m (cid:105) and (cid:104) v (cid:105) represent the average values of all components,respectively.Figure 1 in Wang (2020) shows that if ψ is not included, k in Equation 7 is not a constant, but depending on the in-dividual and total masses of BHs. With the correction factor ψ , the value of k becomes constant for most models with M BH /M < .
4. The N -body models in Wang (2020) are low-mass clusters with only two mass components and no stellarevolution. With more realistic models, we provide a similaranalysis by treating BHs and non-BHs as two components.Figure 12 shows the evolution of energy flux rate (mea-sured at R h and R h , BH ) with and without ψ factor, and theevolution of ψ of all objects and BHs, respectively. With thecorrection of ψ to t rh (the upper panel), the A2.0 and A2.3models show the same ratio of energy flux, while the A1.5and A1.7 models have higher ratios. Without correction (themiddle panel), on the other hand, there is no common ratioamong all four models. The lower panel shows the evolutionof ψ measured at R h and R h , BH . The value of ψ [ R h , BH ] ini-tially increases for all models but more rapidly for the top-heavy models (the A1.5 and A1.7 models). Then it slightlydecreases after 1 Gyr for the A2.0 and A2.3 models. The valueof ψ [ R h ] also increases in the beginning. For the A2.3 model,it significantly at around 100 Myr and then decreases to asimilar value to ψ [ R h ]. For the other models, it also peaks ataround 100 Myr but has a lower value compared to ψ [ R h , BH ]later on. The analysis here ignores the issue of the internal ψ factors for the BH and non-BH components discussed inSection 3.4. But the result is consistent with what is foundin Wang (2020). In this work, we did not include the general relativity effect onBH orbits during the simulation, but considered it in the post-process analysis to detect the BBH mergers. We thereforecannot detect hierarchical mergers, which can be the sourcesfor massive BHs detected by LIGO/VIRGO (Abbott et al.2019, 2020). With the same reason, in our simulations, we alsodid not find intermediate-mass black holes (IMBHs), which
MNRAS000
MNRAS000 , 1–10 (2019) mpact of IMFs on GW sources E k T r h / E k , B H B H T r h , B H A1.5A1.7 A2.0A2.310 E k T r h / E k , B H T r h , B H Time[Myr]05101520 BHAll
Figure 12.
The evolution of energy flux rate (measured at R h and R h , BH ) with (the upper panel) and without (the middle panel) thecorrection factor, ψ , respectively. ψ measured at R h and R h , BH areshown in the bottom panel. are discussed in (Portegies Zwart et al. 2004; Giersz, et al.2015; Rizzuto et al. 2020).Without primordial binaries, these models probably un-derestimate the BBH merger rates, since primordial binariescan lead to stellar-evolution-driven formation of BBHs. Theexchange of components via dynamical encounters betweenprimordial binaries and BHs can also generate BBHs. Thischannel is also missed in our current simulations.By including the tidal field, the survival timescale of A1.5and A1.7 can be much shorter. This can also affect the BBHmerger rates. In the future work, we will carry out new modelsby improving all these aspects.Since we include PPSN in our simulations, the models withtop-heavy IMFs have a large fraction of equal-mass BBHs(40 . M (cid:12) ). This can be changed by adopting the metallicitydifferent from Z = 0 .
001 or different stellar-evolution modelsfor massive stars. Thus, our finding of the mass ratio distri-bution of BBHs cannot represent all kinds of conditions. Butthe trend, i.e., top-heavy IMFs tend to lead to a higher massratio of BBHs, would be general even if the stellar evolutionmodels was changed.
In this work, we carry out four star-by-star N -body sim-ulations of GCs with different IMFs, and initially M =5 × M (cid:12) and R h = 2 pc. We find that the formation rate ofBBH mergers depends on the stellar evolution and dynami- cal process (core collapse and BH heating) in a complicatedway. There is no monotonical correlation between the slope ofIMF ( α ) and the number of (potential) BBH mergers (Fig-ure 2, Table 2). The stronger stellar-wind mass loss in thefirst 100 Myr leads to a faster expansion of GCs with moretop-heavy IMFs. As the escape velocity is lower in more top-heavy models, merging BBHs less efficiently form in theseclusters, although the number of BHs is much higher. How-ever, the A1.5 model, which shows a deeper core collapse afterexpansion, can produce a burst of BBH merger candidates atthe momentum of the core collapse (Figure 2,6).During the long-term evolution, it is more difficult to formBBH mergers in GCs with more top-heavy IMFs because theyexpand faster. This trend can be identified from the evolu-tion of the minimum semi-major axis of hard BBHs shownin Figure 2. Therefore, with the same initial mass and size,GCs with more top-heavy IMFs less efficiently produce BBHmergers within the same time interval of evolution.However, GCs with top-heavy IMFs may maintain the es-cape velocity highly enough for a long term to keep BHs insideclusters. As a result, the total number of BBH mergers canbe large in the case of the top-heavy IMFs (Weatherford etal. 2021). In other words, although the efficiency is low, high-density GCs with top-heavy IMFs can have a longer timeto produce BBH mergers. In the case of GCs with top-lightIMFs, on the ohter hand, the total amount of BBH mergersare limited to the total number of BHs, even though they aremore efficient.The comparison with the two-component models fromWang (2020) suggests that the mass loss rate of BHs dose notsimply depends on the mass segregation time ( t ms ), probablybecause it is difficult to define an accurate t ms when a massspectrum exists. However, the general trend of energy bal-ance (energy flux rate) is consistent with the result of Wang(2020). We also identify the two evolution trend of GCs. BHsescape faster in GCs with the canonical Kroupa (2001) IMF( α = − . ACKNOWLEDGMENTS
L.W. thanks the financial support from JSPS InternationalResearch Fellow (School of Science, The university of Tokyo).M.F. was supported by The University of Tokyo ExcellentYoung Researcher Program. This work was supported byJSPS KAKENHI Grant Number 17H06360 and 19H01933and Initiative on Promotion of Supercomputing for Youngor Women Researchers, Information Technology Center, TheUniversity of Tokyo, and MEXT as “Program for Promot-ing Researches on the Supercomputer Fugaku” (Toward aunified view of the universe: from large scale structures toplanets, Revealing the formation history of the universe withlarge-scale simulations and astronomical big data). Numeri-cal computations were carried out on Cray XC50 at Center
MNRAS , 1–10 (2019) Long Wang et al. for Computational Astrophysics, National Astronomical Ob-servatory of Japan.
DATA AVAILABILITY
The simulation data underlying this article are stored onCray XC50. The data were generated by the software petar ,which is available in GitHub, at https://github.com/lwang-astro/PeTar. The simulation data will be shared via the pri-vate communication with a reasonable request.
REFERENCES
Abbott B. P., et al., 2019, Physical Review X, 9, 031040Abbott B. P., et al., 2020, arXiv, arXiv:2010.14527Antonini F., Gieles M., 2019, arXiv, arXiv:1906.11855Antonini F., Murray N., Mikkola S., 2014, ApJ, 781, 45Askar A., Szkudlarek M., Gondek-Rosi´nska D., Giersz M., BulikT., 2017, MNRAS, 464, L36Bae Y.-B., Kim, C., Lee, H. M., MNRAS, 440, 2714Banerjee S., Kroupa P., 2011, ApJL, 741, L12Banerjee S., 2020, arXiv, arXiv:2011.07000Banerjee S., Belczynski K., Fryer C. L., Berczik P., Hurley J. R.,Spurzem R., Wang L., 2020, A&A, 639, A41. doi:10.1051/0004-6361/201935332Bastian, N. & Lardo, C,. 2018, ARA&A, 56, 83Baumgardt H., Sollima S., 2017, MNRAS, 472, 744Belczynski K., Bulik T., Fryer C. L., Ruiter A., Valsecchi F., VinkJ. S., Hurley J. R., 2010, ApJ, 714, 1217. doi:10.1088/0004-637X/714/2/1217Belczynski K., et al., 2016, A&A, 594, A97Belczynski K., et al., 2020, A&A, 636, A104Breen P. G., Heggie D. C., 2013, MNRAS, 432, 2779Chatterjee S., Rodriguez C. L., Rasio F. A., 2017, ApJ, 834, 68Di Carlo U. N., Giacobbo N., Mapelli, M., Pasquato, M., Spera,M., Wang, L., Haardt, F., 2019, MNRAS, 487, 2947Downing J. M. B., Benacquista M. J., Giersz M., Spurzem R.,2010, MNRAS, 477, 1946Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M.,Kalogera V., Holz D. E., 2012, ApJ, 749, 91. doi:10.1088/0004-637X/749/1/91Fujii M. S., Tanikawa A., Makino J., 2017, PASJ, 69, 94Giacobbo N., Mapelli M., 2018, MNRAS, 480, 2011Giersz M., Askar A., Wang L., Hypki A., Leveque A., Spurzem R.,2019, MNRAS, 487, 2412Giersz M., Leigh N., Hypki A., L¨utzgendorf N., Askar A., 2015,MNRAS, 454, 3150Haghi H., Safaei G., Zonoozi A. H., Kroupa P., 2020, ApJ, 904,43. doi:10.3847/1538-4357/abbfb0H´enon M., 1961, AnAp, 24, 369H´enon M., 1975, IAUS, 133, IAUS...69Heggie D. C., 1975, MNRAS, 173, 729Hills J. G., 1975, AJ, 80, 809Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897Hong J., Askar A., Giersz M., Hypki A., Yoon S.-J., 2020, MNARS,498, 4287Hosek M. W., et al., 2019, ApJ, 870, 44Iwasawa M., Tanikawa A., Hosono N., Nitadori K., Muranushi T.,Makino J., 2016, PASJ, 68, 54Iwasawa M., Namekata D., Nitadori K., Nomura K., Wang L.,Tsubouchi M., Makino J., 2020, PASJ, 72, 13Joshi K. J., Rasio F. A., Portegies Zwart S., 2000, ApJ, 540, 969Kinugawa, T., Inayoshi, K., Hotokezaka, K., Nakauchi, D., Naka-mura, T., 2000, ApJ, 540, 969 Kroupa P., 2001, MNRAS, 322, 231Kumamoto J., Fujii M. S, Tanikawa A., 2019, MNRAS, 486, 3942Lu J. R., Do T., Ghez A. M., Morris M. R., Yelda S., MatthewsK., 2013, ApJ, 764, 155Marks M., Kroupa P., Dabringhausen J., Pawlowski M. S., 2012,MNRAS, 422, 2246Mandel, I., de Mink, S. E., 2016, MNRAS, 458, 2634Marchant P., Langer, N., Podsiadlowski, P., Tauris, T. M., Moriya,T. J., 2016, A&A, 588, A50Milone A. P., et al., 2017, MNRAS, 464, 3636Namekata D., Iwasawa M., Nitadori K., Tanikawa A., Mu-ranushi T., Wang L., Hosono N., et al., 2018, PASJ, 70, 70.doi:10.1093/pasj/psy062O’Leary R. M., Kocsis, B., Loeb, A., 2009, MNRAS, 395, 2127Oshino S., Funato Y., Makino J., 2011, PASJ, 63, 881Park D., Kim, C., Lee, H. M., Bae Y.-B., Belczynski K., 2017,MNRAS, 469, 4665Peters P. C., 1964, PhRv, 136, 1224Rizzuto F. P., Naab T., Spurzem R., Giersz M., OstrikerJ. P., Stone N. C., Wang L., et al., 2020, MNRAS.tmp.doi:10.1093/mnras/staa3634Portegies Zwart S. F., McMillan S. L. W., 2000, ApJL, 528, L17Portegies Zwart S. F., Baumgardt H., Hut P, Makino J., McMillanS. L. W.Rodriguez C. L., Chatterjee S., Rasio F. A., 2016, Physical ReviewD, 93, 084029Rodriguez C. L., Haster C.-J., Chatterjee S., Kalogera V., RasioF. A., 2016, ApJL, 824, L8Rodriguez C. L., Morscher M., Wang L., Chatterjee S., Rasio F. A.,Spurzem R., 2016, MNRAS, 463, 2109Rodriguez C. L., Pattabiraman B., Chatterjee S., Choudhary A.,Liao W.-. keng ., Morscher M., Rasio F. A., 2018, ComAC, 5,5. doi:10.1186/s40668-018-0027-3Samsing J., Askar A., Giersz M., ApJ, 855, 124Schneider F. R. N., et al., 2018, Sci, 359, 69Spitzer L., Hart M. H., 1971, ApJ, 164, 399Tanikawa, A., 2013, MNRAS, 435, 1358Tanikawa, A., Susa, H., Yoshida, T., Trani, A. A., Kinugawa, T.,2020, arXiv, arXiv:2008.01890Wang L., Spurzem R., Aarseth S., Nitadori K., Berczik P., Kouwen-hoven M. B. N., Naab T., 2015, MNRAS, 450, 4070Wang L., Kroupa P., Takahashi K., Jerabkova T., 2020, MNRAS,491, 440Wang L., 2020, MNRAS, 491, 2413Wang L., Nitadori K., Makino J., 2020, MNRAS, 493, 3398Wang L., Iwasawa M., Nitadori K., Makino J., 2020, MNRAS, 497,536Weatherford N. C., Fragione G., Kremer K., Chatterjee S.,Ye C. S., Rodriguez C. L., Rasio F. A., 2021, arXiv,arXiv:2101.02217Zhang Z.-Y., Romano D., Ivison R. J., Papadopoulos P. P., Mat-teucci F., 2018, Natur, 558, 260Ziosi B. M., Mapelli M., Branchesi M., Tormen G., 2014, MNRAS,441, 3703MNRAS000