The impact of primordial binary on the dynamical evolution of intermediate massive star clusters

In this work, we use the $N$-body code petar (Wang, Nitadori & Makino, 2020a) to perform the numerical simulations of the star clusters. petar is a high-performance $N$-body code that combines the particle-tree particle-particle method (Oshino, Funato, & Makino, 2011) and the slow-down algorithmic regularization method (SDAR; Wang, 2020). The latter is the algorithm for accurately evolving the dynamical evolution of multiplicity. By using a hybrid parallel method based on the fdps framework (Iwasawa et al., 2016, 2020; Namekata et al., 2018). The code can handle the simulation of massive star clusters with a binary fraction up to 100 percent. This unique feature allows us to carry out this study.

We also use the single and binary stellar evolution packages, sse and bse, in the simulations (Hurley, Pols, & Tout, 2000; Hurley, Tout, & Pols, 2002; Banerjee et al., 2020b). These population synthesis codes can follow the stellar wind mass loss, the stellar type changes and the mass transfer between interacting binaries. The semi-empirical stellar wind prescriptions from Belczynski et al. (2010) is assumed. We adopt the “rapid” supernova model for the remnant formation and material fallback from Fryer et al. (2012), along with the pulsation pair-instability supernova (PPSN; Belczynski et al., 2016). Because supernova explosion is asymmetric, the compact objects forming after supernovae gain natal kick velocities. Based on the measurement of the proper motions of field neutron stars (NSs) (Hobbs et al., 2005), we adopt a Maxwellian distribution of the kick velocity with a velocity dispersion of 265 km s^-1. The randomly sampling of this distribution result in high velocities that most exceed the escaping velocity of star clusters. The fallback scenario assumes that a part of BHs do not have less or zero kick velocity due to mass fallback or failed supernovae, thus these BHs can be retained in star clusters. Meanwhile, electron-capture supernovae that result in NSs also have a low kick velocity (assuming a velocity dispersion of 3 km s^-1), thus a part of NSs also can be retained. The gravitational wave driven mergers of compact object binaries are also included in the bse following the approximation from Peters (1964).

It is complicated to set up realistic initial conditions of star clusters including the gas component and its large-scale influence. In this study, we focus on how the different binary fractions affect the long-term dynamical evolution of star clusters and GW events. Thus, we ignore the complexity of the gas embedded phase and start the simulation in a gas-free and virialised initial condition which is commonly adopted in previous studies (e.g., Wang et al., 2016).

For all models in this work, the positions and velocities of stars and center-of-the-masses of binaries are generated by randomly sampling from the Plummer density profile (Plummer, 1911) with an initial half-mass radius of 2 pc. The masses of stars are assigned by randomly sampling from the Kroupa (2001) initial mass function (IMF) with the minimum and the maximum masses of $0.08$ and $150~{}M_{\odot}$, respectively. To keep the shape of the IMF, we firstly sample the masses of all stars, and then, pair the two components in binaries based on the mass ratio distribution. No tidal field is assumed. We adopt a common metallicity of $Z=0.001$. The modified mcluster code¹¹1https://github.com/lwang-astro/mcluster is used to generate the initial conditions.

We have investigated 7 different initial fractions of primordial binaries ($f_{\mathrm{b,0}}$), as listed in Table 1. The B0 set has no primordial binary; other sets except the B1.0 adopt the flat distribution of semi-major axes ($a$) and the thermal distribution of eccentricities ($e$). The minimum and the maximum semi-major axes are $3~{}R_{\odot}$ and $100$ AU, respectively. This range covers $a$ of the mass-transfer binaries and all tight binaries that can survive for a long term in the star clusters. Due to the Heggie-Hills law (Heggie, 1975; Hills, 1975), wide (soft) binaries tend to be disrupted after a few strong encounters with intruders while tight (hard) binaries tend to be tighter after encounters. The maximum $a$ of hard binaries can be estimated as

where $G$, $m_{1}$, $m_{2}$, $\langle m\rangle$ and $\sigma$ are gravitational constant, masses of two binary components, the local average mass and the local velocity dispersion, respectively. In our model, the average $a_{\mathrm{h}}\approx 32$ AU. Thus, a part of primordial binaries are in the wide region.

The random sampling from the distributions of $a$ and $e$ can form binaries where two components touch each other at the peri-center position. Such binaries should already have mass transfer or merged during the star formation and should not be included as the primordial binaries. We apply the eigenevolution method suggested by Kroupa (1995b) to adjust the orbital parameters of these binaries to avoid the forbidden region of $a$-$e$.

The pairing of masses of the two components in binaries is different for massive OB ($>5M_{\odot}$) and low-mass stars. We adopt the uniform mass ratio distribution between 0.1 and 1.0 for massive binaries following Sana et al. (2012) and the random pairing of masses for the low-mass binaries. In addition, the observation shows a high multiplicity frequency of OB-stars (Moe & Di Stefano, 2017), thus we adopt 100 percent fraction of OB binaries which is independent of $f_{\mathrm{b,0}}$.

The B1.0 set follows the properties of primordial binaries of low-mass stars ($<5M_{\odot}$) derived by Kroupa (1995a, b) and Belloni et al. (2017) and of high-mass stars based on Sana et al. (2012). Belloni et al. (2017) updated the model in order to be consistent with the observed color distribution of GCs.

Figure 1 compares the $a$-$e$ distribution for the B0.9 set (also representing the sets with $0<f_{\mathrm{b,0}}<1$) and the B1.0 set. For the B0.9 set, most of binaries have $a>10$ AU with the $e$ peaked at 0.8. The B1.0 set has a much larger region of $a$ with the maximum value being approximately $1\times 10^{4}$ AU. The distribution of $e$ is flatter between 0.6 and 0.9 and has a clump near zero.

We create a group of models for each set listed in Table 1 to investigate the long-term dynamical evolution. The group is named as "L-group". The names of models combines the set names in Table 1 and the suffix "L", such as B0L. Simulations of the B0L, B0.1L, B0.3L and B0.5L models have reached 12 Gyr, while others have reached 7 Gyr. After 7 Gyr, most of BHs have escaped from the systems, thus it is sufficient for the purpose of this research.

The random sampling of the initial masses, positions and velocities of stars introduces a stochastic effect. To investigate how this affects the long-term evolution, we create the second "S" group. In this group, we create 5 models with different random seeds for each $f_{\mathrm{b,0}}$. The suffixes "S0, "S1", …, "S5" are used in the names to distinguish the models. These models have been evolved to 500 Myr.

In the L- and S-groups, the binary fractions and the mass pairing of massive and low-mass binaries are different. For a more general study, we also create a group of models without the special treatment of massive binaries, i.e., the property of massive binaries follow that of the low-mass binaries. This group is named as (randomly pairing) "R-group". Similar to the S-group, the R-group contains 5 models with different random seeds for each $f_{\mathrm{b,0}}$, and each model has beend evolved to 500 Myr. Table 2 summarizes the difference of the three groups.

Hereafter, when we refer to all the models with the same $f_{\mathrm{b,0}}$, we simply call them a model name without a suffix, such as B0, or B0.1.

The principle of Hénon (1971, 1975) showed that in a gravitationally bound stellar system, the energy flux demanded by the global system should be balanced by the energy generated by the central engine (binary heating). This principle links the dynamics of central binaries to the long-term evolution of the global structure. Breen & Heggie (2013) showed that when the BH subsystem exists, BBHs are the major source of the binary heating. The balance of the energy in the virial equilibrium can be reflected by the half-mass radii ($R_{\mathrm{h}}$) of the global system and of the BH subsystem (see Eq. 3 in Breen & Heggie, 2013):

where the suffix "BH" reflects the properties of the BH subsystems. When the BH subsystem exists, the core radius $R_{\mathrm{c}}$ is linked to $R_{\mathrm{h,BH}}$. Thus, in Figure 2, we compare the evolution of $R_{\mathrm{h}}$ and $R_{\mathrm{c}}$ for all models in the L-group and investigate that how $f_{\mathrm{b,0}}$ affects the structural evolution.

A scatter of $R_{\mathrm{h}}$ among the models appears after around 100 Myr. However, there is no monotonic relation between $R_{\mathrm{h}}$ and $f_{\mathrm{b,0}}$. The B0L and the B0.7L models have the largest and the smallest $R_{\mathrm{h}}$, respectively, while the B0.3L, B0.5L and B0.9L models have a similar $R_{\mathrm{h}}$. Thus, such a difference is not likely to be caused by different $f_{\mathrm{b,0}}$.

We can understand the reason why the $R_{\mathrm{h}}$ evolution is independent of $f_{\mathrm{b,0}}$ as follows. $R_{\mathrm{h}}$ expands because of binary heating at the cluster center. The binary heating is controlled by the global energy requirement, according to Henon’s principle. The global energy requirement is similar among the models in Figure 2. This is because the initial global structure is the same for all the models, although only the number of objects differs. Note that a (tight) binary star is counted as one object, since it is not resolved in dynamics unless its orbit is significantly perturbed by a close encounter. The global energy requirement is met by only the massive primordial binaries (here, BBHs and their progenitors). The massive primordial binaries have identical properties among the L-group models except for the B0 regardless of $f_{\mathrm{b,0}}$. Thus, their $R_{\mathrm{h}}$ evolution is independent of $f_{\mathrm{b,0}}$.

Binaries are relatively heavier than single stars. Thus, binaries become more centrally concentrated in the cluster due to the mass segregation. The bottom panel of Figure 2 shows the ratio of the half-mass radii between binaries and all stars ($f_{\mathrm{rh}}$). The value of 1.0 indicates that binaries have the same density profile as singles, while the value below 1.0 indicates the mass segregation of binaries. The B1.0L model has $f_{\mathrm{b,0}}>1.0$ initially due to the low-number statistic of single stars. There is a trend that the models with higher $f_{\mathrm{b,0}}$ show larger $f_{\mathrm{rh}}$ (smaller degree of mass segregation), because the mass segregation firstly appears in the sub-group of massive binaries, and the larger fraction of low-mass binaries in these models smooths out the feature of mass segregation for all binaries.

There is an rise of $f_{\mathrm{rh}}$ around 100 Myr for all models, while models with lower $f_{\mathrm{b,0}}$ show a more pronounced feature. There are two mechanisms contributing to this rise. (1) the natal kicks of compact objects formed after supernovae of massive stars invert the mass segregation; these kicks mostly occurred before 100 Myr in massive binaries; (2) the binary heating from BBHs after core collapse ejects BHs from the center.

The evolution of $R_{\mathrm{c}}$ is roughly identical for all models except the B0L. Comparing the B0L and the other models, they show two major differences:

(1) The B0L model shows a pronounced feature of core collapse at approximately 50 Myr (see the middle right panel of Figure 2) while the others do not. The binary heating is the major mechanism to suppress core collapse. In the B0L model, such the process starts when the first binary forms after the core collapse. In the other models, the primordial binaries can provide the heating. The mass segregation starts to appear around 10 Myr as shown in the evolution of $f_{\mathrm{rh}}$ (see the bottom right panel of Figure 2). It is expected that the binary heating of massive binaries can already start at this time, and thus, no clear sign of core collapse appears around 50 Myr.

(2) After about 7 Gyr, $R_{\mathrm{c}}$ of the B0L model continues increasing while those of the other models (B0.1L, B0.3L and B0.5L) show the signal of core collapse. To explain the reason, we show the properties of the BH subsystems in Figure 3. The upper panels compare the evolution of the half-mass radius of BHs ($R_{\mathrm{h,BH}}$). For the B0L model, $R_{\mathrm{h,BH}}$ continues increasing till 12 Gyr, while for all the other models, $R_{\mathrm{h,BH}}$ decrease with sharp peaks appearing between 4 to 8 Gyr. The evolution of the number of BHs inside $R_{\mathrm{c}}$ ($N_{\mathrm{c,BH}}$; the middle panels of Figure 3) explains the reason for this behaviour. After about 100 Myr, the B0L model contains much more BHs than the other models. At the end of simulation, the B0L model still keeps about 10 BHs inside $R_{\mathrm{c}}$. But BHs are already depleted between 4 to 8 Gyr in the other models. Thus, the sharp peaks appear in the evolution of $R_{\mathrm{h,BH}}$ during this period because of the ejection of BHs and the low-number statistics. After BHs are depleted, the leak of binary heating drives the core collapse of light stars as shown in Figure 2. This is also consistent with the theory of Breen & Heggie (2013). Therefore, the initial difference of $N_{\mathrm{c,BH}}$ significantly affects the long-term structural evolution of star clusters.

The bottom panel of Figure 3 shows the evolution for the mass ratio of BHs and all objects inside $R_{\mathrm{c}}$ ($f_{\mathrm{c,BH}}$). After the BH subsystem forms, BHs contribute to about 40 and 60 percent of the mass inside $R_{\mathrm{c}}$. $f_{\mathrm{c,BH}}$ decreases as BHs escape via few-body interactions. At the end, BHs in the B0L model still occupy about 10 percent of the core mass.

To further investigate the reason for the initially different number of BHs, we collect all BH progenitors with the initial stellar masses $>20.3~{}M_{\odot}$ and check how many of them have escaped from the star clusters before 100 Myr and how many have merged in binaries. Such events reduce the remaining number of BHs in the cluster. The result is shown in Table 3. The numbers of formed BHs ($N_{\mathrm{BH}}$) differ among the models due to the random effect of initial condition and the mergers. The major mechanism to remove BHs are from single escapers (E-S), where most of them are driven by the high natal kick velocities after asymmetric supernovae (E-SK). The numbers of E-SK are similar for all L-group models, while the numbers of E-S and binary escapers (E-B) from the B0L model are about 10 less than those of the other models, respectively. Meanwhile, there is no natal kick driven binary escapers (E-BK) for all models. There are also a few less mergers in the B0L model. Therefore, the escaping of singles and binaries driven by dynamical interactions and mergers are the major mechanism that causes the less number of BHs in the models with primordial binaries.

To validate that the massive binaries (BBHs) dominate the binary heating, we investigate the evolution of semi-major axes ($a$) of all compact binaries inside $R_{\mathrm{c}}$ for the L-group models (see Figure 4). Firstly, the tracks of $a$ for individual BBHs show a pronounced decrease of $a$ down to 1 – 10 AU. This is the sign of binary heating process as described by the the Heggie-Hills law, where tight binaries become tighter after encounters with others and transfer the binding energy outwards. When $a$ becomes small enough, the BBH either merge or escape from the core. Generally, only a few BBHs can provide heating simultaneously. Once tight binaries in the core are ejected, new tight binaries start to evolve.

The B0L model continues producing new BBHs up to 7 Gyr, while the formation rate of tight BBHs significantly decreases in other models after 4 Gyr due to the lack of BHs. The total number of dynamically formed BBHs inside $R_{\mathrm{c}}$ before 7 Gyr is about 570, where only 62 BBHs are tight with the minimum semi-major axis below 200 AU. The BBHs formed via the exchange of members after three- or four-body encounters are also counted as new ones.

Meanwhile, there are also a large group of low-mass binaries with white dwarfs (WDs) and NSs in all models except the B0L. They form from the low-mass primordial binaries. Unlike the BBHs, these low-mass binaries do not show the decrease of $a$, and thus, they do not contribute to the binary heating. A few exceptions with small $a$ are probably due to the binary stellar evolution, which leads to the mergers of WD-WD binaries at the end. This explains why the evolution of $R_{\mathrm{h}}$ and $R_{\mathrm{c}}$ is independent of $f_{\mathrm{b,0}}$.

Figure 4 only shows the result of the first 7 Gyr, when BBHs dominate the binary heating. After most BHs escape, the core collapse of light objects including NSs and WDs occur as shown in Figure 2. Then, low-mass binaries start to dominate the binary heating.

As shown in Figure 2, the structural evolution in the L-group does not monotonically depends on $f_{\mathrm{b,0}}$. However, it can be possible that the random sampling of masses, positions and velocities of individual stars and binaries cause the stochastic scatter, which may override the general $f_{\mathrm{b,0}}$-dependent trend. We investigate this by comparing the structural evolution of the S-group models in the left two panels of Figure 5. By showing the 5 models with different random seed for each $f_{\mathrm{b,0}}$ together (same color), we can identify that the evolution of $R_{\mathrm{h}}$ for the models with primordial binaries indeed has a wide stochastic scatter due to the random sampling (see the first and the second rows in Figure 5). Including such a scatter, the models with primordial binaries show an identical general trend, while the B0S models show a pronounced difference. This again suggests that only the massive binaries (BBHs) affect the long-term structural evolution.

To further confirm this, we also compare the structural evolution of the R-group models, where the primordial binary fraction follows $f_{\mathrm{b,0}}$, in the right panels of Figure 5. Unlike those in the S-group, the models excluding the B1.0R show a closer evolution of $R_{\mathrm{h}}$ and $R_{\mathrm{c}}$ to the B0R models, The low-$f_{\mathrm{b,0}}$ models also show a similar timescale of core collapse to that of the B0R model.

Figure 3 suggests that $N_{\mathrm{c,BH}}$ is the key parameter to determine the structural evolution. We also compare this for the S- and R-group models in the bottom panels of Figure 5. Similar to the L-group, the S-group shows an initially large difference of $N_{\mathrm{c,BH}}$ between the B0S models and the others. But the R-group has a different behavior: except the B1.0R models, the others show a similar $N_{\mathrm{c,BH}}$ that weakly depends on $f_{\mathrm{b,0}}$. This is consistent with the behaviour of $R_{\mathrm{h}}$ and $R_{\mathrm{c}}$. It is interesting to see that the large-$f_{\mathrm{b,0}}$ models in the R-group do not show a pronounced difference of $N_{\mathrm{c,BH}}$ compared to that of the B0R models. This suggests that the fraction of massive primordial binaries is not the major parameter that determines $N_{\mathrm{c,BH}}$, the orbital properties, including the mass ratio, the period distribution and the eccentricity distribution, have a stronger impact. In summary, we can conclude that only massive binaries are important on the dynamics evolution of the system.

The binary heating process is also the mechanism to produce GW sources, especially for the dynamically formed BBH mergers. Thus, we also investigate how the GW mergers depend on $f_{\mathrm{b,0}}$ based on the L-group models. The numbers of mergers for different combination of WD, NS and BH binaries are summarized in Table 4. For a self-consistent comparison, we only select mergers before 7 Gyr for all models. For each model, we also distinguish the origins of mergers by two ways, which are represented as the suffixes and prefixes in the type names, respectively.

Suffix: we detect the mergers from primordial binaries (PB) and from dynamically formed binaries (DY). The former can be triggered either by the binary stellar evolution or by the stellar dynamics. Thus, we also use the standalone bse code to directly evolve all primordial binaries from the initial conditions of each model and detect the GW mergers (ISO). Some primordial binaries merged in both ISO and $N$-body models, we name it as common mergers (C). In Table 4, we show the numbers of common mergers from $N$-body models. But the component stellar types can be different for ISO and $N$-body models. For example, the E-C in B0.3L has one BBH merger, but in the ISO group, it is a double-WD binary (BWD) merger.

Prefix: meanwhile, mergers in the $N$-body models can be inside star clusters (the prefix "I") or already escape from the clusters (the prefix "E").

Table 4 shows that there is no BBH merger from primordial binaries in both ISO and $N$-body models (see the numbers of I-C and E-C in the BBH column). We have investigated and found that this is due to a bug in the interface between petar and bse. This bug appears when the main function (evolv2 in the source file) of bse is called multiple times. bse is originally designed as a standalone code so that evolv2 is called only once to evolve a binary to a given time. When the code is plugged in a $N$-body code, the evolv2 function needs to be called many times to couple the binary stellar evolution and the dynamical perturbation to the binary. We found that when evolv2 is called to evolve a binary in the Roche-overflow phase, a parameter is not properly initialized so that an artificial expansion of the orbit happens. Such an expansion changes the following evolution and finally reduce the merger rate. Unfortunately it is extremely computationally expensive to redo the $N$-body models (a few months), and thus, we can only redo the binary stellar evolution by using the fixed standalone bse code, similar to the ISO groups in Table 4. The comparison between the old (ISO) and the bug-fixed versions are shown in Table 5. The bug-fixed version of bse results in much more BBH mergers, especially for the B1.0L model. We expect that if the correct version is used in $N$-body models, the BBH mergers from primordial binaries would also increase. Since we focus on the dynamically formed BBHs, although there is a systematical bias of BBH merger numbers, it does not change our major conclusion. In the following analysis, we take into account the impact from this bug.

Table 4 shows that most of the GW mergers are from BWDs, and no BWD merger is from dynamically formed binaries. This is expected since BWD is not strongly affected by the stellar dynamics when BBHs dominate the binary heating as shown in Figure 4. We expect that after the core collapse of light objects when most BHs have escaped, the dynamically formed BWDs will start to appear. Because BWDs formed from low-mass binaries, a larger $f_{\mathrm{b,0}}$ results in a larger number of BWD mergers.

The numbers of BWD mergers are comparable in the ISO and in the $N$-body models and most of them occurred in both cases. There are a few ISO BWD mergers that do not appear in the $N$-body models and vice versa. We check the history of the binary stellar evolution and dynamical encounters of these mergers, and find three reasons that cause the diverged evolution:

(1) The binary stellar evolution of the same binary is not fully reproducible when the calling frequency of the evolution function (evolv2) is different. The different calling frequency can cause the numerical variation. For example, the mass loss rate of the two components may be slightly different case by case. Then, the difference grows up and the final result diverges. Particularly, if the evolution track of one component in a binary is close to the boundary that can either evolve to a WD or a NS, the binary may become a BWD merger, a BWN merger or not a merger. The bug we found in the bse code is also related to the calling frequency and also contribute to the divergency.

(2) The asymmetric stellar wind and supernovae result in a natal kick when a WD, NS or BH forms. In bse, the kick velocity is randomly sampled based on an assumption of the velocity distribution function. Such a random effect can significant change the final state of a binary, especially for the mergers including NSs or BHs.

(3) In star clusters, the dynamical interaction can also perturb the binary orbit and affect the final result. But this effect is weak for most BWD mergers. It is mostly important for the BBH mergers.

The number of DYN mergers (mainly from BBHs) are much less than that of the PB mergers and show a weak dependence on $f_{\mathrm{b,0}}$. As shown in Figure 3, the evolution of $a$ for BBHs is similar for all models with primordial binaries. Thus, although it may be affected by the low number statistics, it is reasonable to see a similar number of DYN mergers.

There are a few escaped mergers, where most are also PB mergers (see E-PB in Table 4). In $N$-body models, these mergers are driven by a mix effect of binary stellar evolution and stellar dynamics.

In Figure 6, we compare the initial (zero-age) orbital parameters of the GW mergers. The dynamical mergers show the initial orbital parameter when the binaries form. The ISO and PB mergers have the $a$-eccentricity ($e$) distribution close to the boundary where the peri-center separation is close to the stellar radii of two components. The escaped mergers from primordial binaries (E-PB) have very different $a$ and $e$ from the boundary, thus these mergers are probably driven by the dynamical interactions.

Generally, we find that the low-mass binary (BWD) mergers do not show pronounced difference in the ISO and in the $N$-body models. When massive BBHs exist, they dominate the binary heating process by which the dynamical impact on the low-mass binaries is strongly suppressed, as shown in Figure 4. We expect that the dynamical channel to form low-mass GW mergers can only be important when most BHs are kicked out from the host star clusters.