The Long-term Evolution of Main-sequence Binaries in DRAGON Simulations

Among four models of DRAGON simulations, model D4-R3-IMF01 was initialized with a more compact density profile. This model was evolved for only 1 Gyr, which is too short for a comparison with observed globular clusters. Cluster model D3-R7-ROT included global rotation. However, it did not show any significantly different results in the evolution of the binary population, which is the focus of our study. The other two models, D1-R7-IMF93 and D2-R7-IMF01, have many similar initial conditions and a few differences. The small differences between both models result in a different evolution that deserves a comparison. Therefore, in this work we select D1-R7-IMF93 and D2-R7-IMF01 as our target models to study the evolution of the binary population in detail. Both models aim to simulate the dynamical evolution of globular cluster NGC 4372. In both models, the globular cluster is placed in an external tidal field corresponding to that of the Milky Way, at NGC 4372's Galactocentric distance of 7.1 kpc. The star cluster is placed in a circular orbit around the Milky Way's center. The Galactic center is represented with a point mass, corresponding to the enclosed mass of the Milky Way of 8 × 10¹⁰ M_⊙ and exerting the tidal force on the star cluster. In both models, the clusters are initially tidally underfilled.

As direct N-body simulations, the DRAGON simulations represent each particle as an individual star in the cluster. These globular cluster models are initialized with 950,000 single stars and 50,000 primordial binary systems. We list the main properties of the two models in Table 1. Both globular cluster models have a King (1966) initial density profile with a King dimensionless parameter W₀ = 6, and are assigned initial half-mass radii, R_h,0, of 7.5 pc for D1-R7-IMF93 and 7.6 pc D2-R7-IMF01. Both globular clusters are initialized in virial equilibrium, Q = 0.5, where the virial ratio Q = ∣T/U∣ is the ratio between the total kinetic energy (T) and the total potential energy (U) of the star cluster.

One major difference between the initial conditions of D1-R7-IMF93 and D2-R7-IMF01 is the initial mass function (IMF). Model D1-R7-IMF93 adopts the IMF of Kroupa et al. (1993), while D2-R7-IMF01 uses that of Kroupa (2001). The IMF is initialized following f(m) ∝ m^−α in the mass range 0.08–100 M_⊙. Both simulations adopt α = 1.3 for the mass range 0.08 < m/M_⊙ ≤ 0.5. The value of α differs in 0.5 < m/M_⊙ ≤ 1, with α = 2.2 for D1-R7-IMF93, and α = 2.3 for D2-R7-IMF01. Finally, for the m > 1 M_⊙ range, D1-R7-IMF93 uses α = 2.7 and D2-R7-IMF01 uses α = 2.3. Due to different slopes in the IMFs, D2-R7-IMF01 has a higher total mass than D1-R7-IMF93. Therefore it also has a slightly larger tidal radius (R_t,0), as shown in Table 1. The initial metallicity of both models, Z = 0.00016, was chosen to be identical to that of the globular cluster NGC 4372 (Geisler et al. 1995; Kacharov et al. 2014; San Roman et al. 2015).

Each model is initialized with a binary fraction of 5%. The start of the simulations (t = 0) represents the time after which the star clusters have expelled all their interstellar gas, have smoothed out their initial substructure, and have obtained virial equilibrium. Moreover, soft binaries are assumed to have been destroyed; the binary population at t = 0 Myr is thus modeled to represent the hard binary population at that epoch. In the analysis below, the first snapshot (0.1 Myr) after the start of the simulations is defined as t = 0 Myr. The initial semimajor axes of the primordial binary system follows a log-normal distribution in the range 0.005–50 au. The lower limit is a consequence of the physical size of the stars, while the upper limit is near the hard-soft boundary for stars in these globular clusters. The initial eccentricity distribution is thermal: f(e) ∝ e (see, e.g., Heggie 1975; Goodman & Hut 1993).

The initial mass-ratio distribution of the binary stars, f(q), is also different for the two simulations. Here, the mass ratio is defined as q = m₂/m₁, where m₁ is the mass of the primary star, and m₂ is the mass of the companion star (i.e., the least massive of the two stars in the binary system). In D1-R7-IMF93, the mass of primordial binaries were generated by random pairing from the IMF. In D2-R7-IMF01, the mass of the primary star (m₁) is randomly chosen from the IMF, and the secondary mass m₂ = qm₁ is subsequently assigned after drawing the mass ratio from the probability distribution f(q) ∝ q^−0.4 (Kouwenhoven et al. 2007) in the mass range 0.08 M_⊙ ≤ m₂ ≤ 100 M_⊙. The latter lower limit ensures that both primary stars and companion stars are in the mass range 0.08–100 M_⊙. This major difference will lead to variation in binary evolution (see Section 3.1).

Stellar evolution is modeled following the prescriptions of Hurley et al. (2000) for single stars, and those of Hurley et al. (2002) for interacting binary systems. The gravitational kicks exerted on neutron stars and black holes follow a Maxwellian velocity distribution. The variance of the Maxwellian distribution in model D1-R7-IMF93 differs from that in model D2-R7-IMF01 (see Table 1). The kick velocity determines the number of remaining bound neutron stars and the retained black hole subsystems in the cluster, which affects the global dynamics of the cluster.

Primordial binaries are present immediately after the star formation process has ended (e.g., Kouwenhoven et al. 2003), and have therefore, in an idealized scenario, not yet been affected by dynamical evolution. Over time, existing binary systems are disrupted, and new binary systems are formed. Dynamical binary systems are formed as the simulation proceeds, through dynamical close encounters between single stars and/or binary systems. The two stars that form the new binary system may or may not have been part of a primordial (and/or dynamical binary) at earlier times. In this paper, we consider the binaries present in the initial conditions of D1-R7-IMF93 and D2-R7-IMF01 as primordial binaries, which can be identified with the unique IDs. Any binaries formed at later times are considered as dynamical binaries, whose IDs are different from the primordial ones.

We identify all binary systems in each snapshot of both models, for a period of 12 Gyr, and identify the evolution of each binary system. The evolution of total number of MS stars, N, and the properties of MS binaries for both simulations are listed in Tables 2 and 3, for snapshots at times 0 (more precisely, 0.1), 100, 300, 600, 1000, 3000, 6000, and 12,000 Myr. The detailed binary samples for each snapshot are provided in the Appendix.

Table 2. Numerical Evolution of the Number of MS Binaries in D1-R7-IMF93, at Different Times

Property		0	100	300	600	1000	3000	6000	12,000
All particles	N	1,050,000	1,046,536	1,045,610	1,044,067	1,041,297	1,020,529	986,958	923,860
Binary fraction (%)	${ \mathcal B }$	4.73	4.63	4.50	4.36	4.21	3.84	3.55	3.31
All binaries	B	49,628	48,458	47,036	45,489	43,885	39,169	35,007	30,568
Primordial binaries	BP	49,628	48,446	47,018	45,454	43,833	39,059	34,860	30,359
Dynamical binaries	BD	0	12	18	35	52	110	147	209

Note. Binary fraction is defined as the total number of binaries divided by the total number of binaries and single stars. Time, in megayears, is indicated in the first row.

Download table as:  ASCIITypeset image

Table 3. Numerical Evolution of MS Binaries in D2-R7-IMF01, at Different Times

Property		0	100	300	600	1000	3000	6000	12,000
All particles	N	1,050,000	1,041,914	1,039,032	1,032,165	1,021,608	966,365	887,621	729,738
Binary fraction (%)	${ \mathcal B }$	4.76	4.72	4.61	4.48	4.39	4.13	3.94	3.77
All binaries	B	50,003	49,151	47,880	46,252	44,798	39,894	34,935	27,530
Primordial binaries	BP	49,997	49,136	47,857	46,223	44,756	39,831	34,840	27,438
Dynamical binaries	BD	6	15	23	29	42	63	95	92

Note. The parameters are the same as those in Table 2.

Download table as:  ASCIITypeset image

Figure 1 shows the secular evolution of MS binary mass ratio in model D1-R7-IMF93. Throughout this work, we define the mass ratio, q = m₂/m₁, as the ratio of the masses of the secondary star, m₂, and the primary star, m₁, where m₁ ≥ m₂. Binary systems in model D1-R7-IMF93 are generated through random pairing of the two components, which are independently drawn from the IMF. The mass-ratio distribution for the combined sample of primordial MS binaries generated through random pairing is highest in the range 0.2 < q < 0.4, and remains so during the entire simulation. The relative frequency of binary systems with q ≤ 0.2 decreases over time. At early times, primary stars of binaries with q ≲ 0.2 are typically more massive than those of binaries with q ≳ 0.2. The former evolve faster and experience higher mass loss, resulting in a tendency of the average mass ratio of the MS binary population to increase over time.

Zoom In Zoom Out Reset image size

The relative frequency of MS binaries q ≤ 0.2 in D2-R7-IMF01 also decreases over time (see Figure 2), similar to in model D1-R7-IMF93. However, the decrease is less prominent due to the different initial mass-ratio distribution of Kouwenhoven et al. (2007) that is assigned to the population. For MS binaries in model D2-R7-IMF01, the initial fraction of MS stars with q ≤ 0.2 in the combined sample is lower than that for model D1-R7-IMF93. Model D1-R7-IMF93 has a higher fraction of MS binaries with lower mass ratios, and therefore more massive MS binaries.

Zoom In Zoom Out Reset image size

Throughout the entire simulation, the number of dynamically formed MS binaries remains more than two orders of magnitude smaller than the number of primordial MS binaries, for both model D1-R7-IMF93 and model D2-R7-IMF01. The bottom panels of Figures 1 and 2 show the mass-ratio distributions for the dynamical binary systems in both models. Note that while the top and middle panels in these figures we show mass-ratio distributions, f(q), the bottom panels show cumulative mass-ratio distributions, F(q), which are more suitable for representing the much smaller number of dynamical binary systems. Note that the mass-ratio distribution of a population of binary systems with a restricted primary star mass range may look substantially different from that of the entire population of binary systems (see, e.g., Kouwenhoven et al. 2009, and references therein).

Figure 3 shows the secular evolution of the median primary masses of the MS binaries in models D1-R7-IMF93 and D2-R7-IMF01. The median primary mass resulting from random paring (D1-R7-IMF93) is generally larger than that resulting from the Kouwenhoven et al. (2007) recipe (model D2-R7-IMF01), fully consistent with Figures 1 and 2. Over time, the median primary mass in both models decreases as a consequence of stellar evolution.

Zoom In Zoom Out Reset image size

Numerous observational studies have made efforts to obtain parameter distributions of binary systems in star clusters (e.g., Sollima et al. 2010; Milone et al. 2012, 2016). In addition to the binary fraction, the most sought-after parameter distributions are that of the mass ratio and that of the semimajor axis (or orbital period). The mass-ratio distribution of massive MS binaries (B-type), in the OB-association Sco OB2, is well described by the power law f(q) ∝ q^−0.5 (Shatsky & Tokovinin 2002), similar to Kouwenhoven et al. (2007). The mass-ratio distribution in the regime q ≤ 0.1 of these massive binaries in Sco OB2, is similar to our binary samples at the time before 1 Gyr. On the other hand, a flat distribution of the mass ratio (for q > 0.5) has been observed in 59 Galactic globular clusters (see Milone et al. 2012).

Raghavan et al. (2010) carried out a systematic study for binary companions among 454 solar-type stars in the Galactic neighborhood. The mass-ratio distribution of these field binaries peaks at q ≈ 1 (Raghavan et al. 2010), with companion star masses in the range 0.05–0.98 M_⊙. D2-R7-IMF01 appears to be consistent with Raghavan et al. (2010). However, if the high-order binaries (i.e., triples or quadruples) are removed in Raghavan et al. (2010), the mass-ratio distribution flattens.

The semimajor axis, a, is another important parameter in binary evolution. We show the secular evolution of the semimajor axis distribution, f(a), in Figures 4 (for D1-R7-IMF93) and 5 (for D2-R7-IMF01). The initial semimajor axis has a log-normal distribution, ranging from a = 0.005 au to a = 50 au for both simulations. Binary systems with large semimajor axes are easily disrupted, since they are vulnerable to close encounters with neighboring stars. Therefore, values of the semimajor axis are limited to a = 50 au at the start of the simulation. Generally, the evolution pattern of a in primordial binaries is similar for both simulations, in that the relative frequency of small a increases with time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The case is different in dynamical binaries. In D1-R7-IMF93 before 1 Gyr, a high fraction of dynamical binaries (Figure 4) have a ≲ 10 au. The fraction of dynamical binaries with a ≲ 10 au drops from 0.4 at t = 100 Myr to 0.1 at t = 12 Gyr. On the contrary, in D2-R7-IMF01, the fraction of dynamical binaries with semimajor axes a ≲ 10 au increases with time. As compared to the MS field binaries in Raghavan et al. (2010), the semimajor axis ranges from 10⁻²–10⁵ au and peaks at 10–100 au; D1-R7-IMF93 and D2-R7-IMF01 produced a much larger number of close binaries with a ≲ 10 au, as wide binaries are less likely to survive in a star-cluster environment.

Although the evolution of the semimajor axis is closely linked to that of the eccentricity, e, the secular evolution of eccentricity is not as apparent as that of the semimajor axis and the mass ratio. We show the secular evolution of the eccentricity distribution, f(e), in Figures 6 (for D1-R7-IMF93) and 7 (for D2-R7-IMF01). Both simulations D1-R7-IMF93 and D2-R7-IMF01 initially adopted a thermal eccentricity distribution, f(e) ∝ e. To avoid immediate collisions between the two components of binary systems at the start of the simulations, binary systems with very small periastron distances are not included in the DRAGON simulations. This is the reason why there is a depression at e ≈ 1. Through dynamical encounters, dynamical binaries reach an energy equipartition state, which leads to a thermal eccentricity distribution for dynamical binaries, similar to the earlier discovery by Jeans (1919).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our models predict that the mean eccentricity at t = 12 Gyr is 〈e〉 = 0.623 and 〈e〉 = 0.627 for models D1-R7-IMF93 and D2-R7-IMF01, respectively. Milone et al. (2016) find that the mean eccentricity of binaries within the half-mass–radius (excluding the core region) of Galactic globular clusters is e = 0.35 ± 0.18 (Milone et al. 2016), assuming a binary fraction of 33%. The discrepancy between the two may be explained by the choice for the initial conditions of D1-R7-IMF93 and D2-R7-IMF01, notably the initial thermal eccentricity distribution, and the initial binary fraction of 5%. However, it should be noted that in Milone et al. (2016) only one of the binary systems is a MS–MS binary, and therefore a direct comparison has to be carried out with caution.

In this section, we search for the correlations between the parameters of MS binaries at 12 Gyr. We analyze the evolution of the mass-ratio distributions (q), semimajor axis (a), eccentricity (e), period (p), primary mass, secondary mass, and cluster-centric distances (r) of binaries in the star clusters under consideration.

We carry out this analysis for both models D1-R7-IMF93 and D2-R7-IMF01, and we consider both primordial binary systems and dynamical binaries. One of the major differences between D1-R7-IMF93 and D2-R7-IMF01 is the initial binary mass ratio q, which generates a statistically different distributions, as shown in Figures 1 and 2. We display q and primary/secondary mass in Figure 8. Red dots are primordial binaries, while blue dots are dynamical binaries. There is an asymptotic boundary for the distribution of q and the primary mass, which is due to the lower limit of secondary mass 0.08 M_⊙, equal to the IMF lower boundary, as q is the ratio between the secondary and primary mass. The comparison between q and the secondary mass shows a nearly straight curve with a slope of 1.25 as the lower boundary, which is also the maximum possible mass of the primary star. At the end of the simulation, at 12 Gyr, all massive stars evolve over the MS. The highest mass of MS stars reaches 0.8 M_⊙.

Zoom In Zoom Out Reset image size

The number of dynamical binaries (blue dots) is statistically too small. The trend of total binaries (black dots) is dominated by primordial binaries (red dots). The mean mass ratio declines with increasing primary mass (black dotted curves), and as the secondary mass decrease.

Comparing black curves in D2-R7-IMF01 and D1-R7-IMF93, the D2-R7-IMF01 has larger average q value because of the different initial mass-ratio distribution (see Section 2.1). The latter generated relatively more high-q binaries, and thus a larger average q value.

Figures 9 and 10 show another set of correlations between parameters in D1-R7-IMF93 and D2-R7-IMF01 at 12 Gyr. Despite the different mass-ratio distribution of binaries on the initial conditions, the MS binary parameter correlations in these two simulations show very similar behaviors. We find no statistical difference of these six correlations in two simulations. The evolution of these parameters are not influenced by mass ratio. Although the considered binaries are only MS stars, they form the vast majority of the binary population at all times. Therefore, the behavior of the general binary population is well described with the MS binaries.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Generally, binaries located in the outskirts of the star clusters typically have a larger semimajor axis (upper left panels in Figures 9 and 10), and longer orbital periods (bottom left panels) than their counterparts in the inner part. The linear relation between the semimajor axis and the period is simply a consequence of Kepler's third law. The dynamical binaries (blue dots) preferentially have a larger semimajor axis, longer period, and more eccentric orbits than the average primordial binaries. The encounter probability in the inner region of the star cluster is far larger than that in the outer region, therefore the dynamical binaries form more easily in the inner region. The blue dots in the left panels shows that the dynamical binaries are statistically more often located in the inner regions than primordial binaries. The eccentricity of the system (middle left panels) does not depend on the position of the binary system in the star cluster. The values of eccentricity are related to the strength of encounters, which do not depend on the position either (e.g., Spurzem et al. 2009; Flammini Dotti et al. 2019).

Binary stars with a small semimajor axis or short orbital periods do not have large eccentricities (middle right panels), because two very close stars are more likely to have mass transfer through the Roche-lobe overflow (Eggleton 1983), and circularize and/or merge shortly after. Due to initial thermal eccentricity distribution, the mean eccentricity is e ≈ 0.618, and the distribution in e² is flat. The uniform number density in the distribution of e² versus the semimajor axis indicates a tendency toward energy equipartition among MS binaries.

As a consequence of the dependence of dynamical processes on the local stellar density, the binary fraction in the star-cluster core tends to evolve differently from that in the outskirts. We show the radial binary fraction, normalized to the value within the core binary fraction, in the upper panels of Figures 11 and 12. For both models D1-R7-IMF93 and D2-R7-IMF01, before 1 Gyr, the binary fraction is larger within the core radius, and slightly drops outside the core. Only in D1-R7-IMF93 does a significant drop of binary fraction outside the core occur after 1 Gyr. The radial decreasing trend becomes more and more steeper with time. In model D2-R7-IMF01, the binary fraction drops outside the core radius, more profoundly so at 1 Gyr. Binary evolution is often closely linked to the long-term dynamical evolution of star clusters. The radial binary distribution can therefore be used as a probe of notable dynamical events that have occurred in the history of the star cluster.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to determine the related dynamical process at 1 Gyr, we show the evolution of the Lagrangian radii in the bottom panels of Figures 11 and 12. The significant decreases in the 0.1% and 1% Lagrangian radii are the evidence of core collapse, which terminated at t ≈ 1 Gyr in both simulations. The other Lagrangian radii curves show a smooth expansion, since the core collapse affects only the innermost regions of the star cluster. The classical core-collapse phase is caused by two-body relaxation processes, which transfer energy from the inner region to the outskirts of a star cluster. In the DRAGON simulations, the center of the cluster is occupied by a black hole subsystem after 1 Gyr, which is quite different from the classical paradigm. As time passes, the radial gradient of the binary fraction increases, with the normalized radial binary fraction dropping at larger radii as the cluster evolves approaching the time of ∼1 Gyr when the core-collapse phase for both simulations is terminated. A steep radial binary fraction distribution may thus indicate a post-core-collapse phase in globular clusters, and vice versa.

The decreasing radial binary fraction at 12 Gyr appears to be a common trend, and is similar to the observed globular clusters (Milone et al. 2016). We fit the radial binary fraction distribution at 12 Gyr of D1-R7-IMF93 with the function proposed by Milone et al. (2016):

$$\begin{eqnarray}&&f({R}_{* })=\displaystyle \frac{{a}_{1}}{{\left(1+{R}_{* }/{a}_{2}\right)}^{2}}+{a}_{3},\end{eqnarray} \tag{ 1 }$$

where R_* is the cluster-centric distance in units of the core radius, the same as in Figures 11 and 12. This function generally fits the D1-R7-IMF93 binary sample (Figure 13), with fitted parameters of a₁ = 0.72 ± 0.10, a₂ = 2.75 ± 0.71, and a₃ = 0.62 ± 0.02. These values are different from the parameters of Galactic globular clusters, i.e., a₁ = 1.05, a₂ = 3.5, and a₃ = 0.2 (Milone et al. 2016). The smaller value of a₁ in D1-R7-IMF93 might be due to the sample only including MS binaries. A larger value of a₂ indicates a steeper decrease of the radial binary fraction in our sample. When all types of binaries in the simulation are included in the fit, we obtain a₁ = 0.82 ± 0.10, a₂ = 3.04 ± 0.71, and a₃ = 0.56 ± 0.02, which show no significant differences from the values obtained from MS binaries when fitting errors are taken into account.

Zoom In Zoom Out Reset image size

The major dynamical process in a star cluster is two-body relaxation, which segregates massive stars to the cluster center, and low-mass stars to the outskirts (Khalisi et al. 2007). Even a cluster as young as 1 Myr can exhibit mass segregation (Pang et al. 2013), since the segregation timescale is faster for massive stars. The 12 Gyr simulation time in D1-R7-IMF93 and D2-R7-IMF01 is much longer than the two-body relaxation timescale, and long enough to segregate low-mass members. We found significant mass segregation within the core in both simulations (see Figure 14). The mean stellar mass in the core region (within the 0.1%–1% Lagrangian radii; see Figures 11 and 12) is more than an order of magnitude larger than in the outskirts of each cluster. Gaburov & Gieles (2008) investigated the integrated color of simulated star clusters, and found that the process of mass segregation will eventually change them. There was a V − I color difference of roughly 0.1–0.2 mag between the center and outer part of the mass-segregated clusters.

Zoom In Zoom Out Reset image size

To investigate the radial color difference in the segregated DRAGON clusters, we derived the integrated color F330W − F814W (Hubble Space Telescope (HST) filters) of two modeled clusters via GalevNB (Pang et al. 2016) and show the color value in each annulus in Figure 15. A color gradient is found in all star and binary samples of D1-R7-IMF93, and only in all star samples of D2-R7-IMF01. The color is bluer in the center and redder in the outskirts, with a difference of ∼0.1 mag. The result is in agreement with the findings of Gaburov & Gieles (2008). The color gradient of binary stars may reflect the steep radial profile of the binary fraction in D1-R7-IMF93, which is not found in D2-R7-IMF01. Thus no color gradient is found in binaries of D2-R7-IMF01 for the same reason. Therefore, mass segregation and the binary radial distribution are related. Both have an imprint in the star clusters' photometry.

Zoom In Zoom Out Reset image size

The color difference of V − I or F330W − F814W = 0.1 mag requires a photometric uncertainty of less than 0.1 mag. This corresponds to V (F555W) ∼ 23–24 mag in HST depending on the exposure time. HST has been the major instrument for extragalactic star-cluster observations, which usually cannot be resolved into single stars.

The Chinese Space Station Telescope (CSST) will be an essential piece of equipment of observing extragalactic objects after HST, with a spatial resolution of ∼015 (Cao et al. 2018; Gong et al. 2019). CSST will observe down to g = 26 mag. With the ultra-deep-field observation (down to 30 mag), much better photometry is expected. In the right panels in Figure 15, we show radial color NUV − y distribution for both models. The color gradient becomes more significant in NUV − y color, with a color difference of ∼0.2 mag, which is larger than HST filters. Therefore, CSST will be an excellent instrument to use to search for the dynamical signature through radial color difference in extragalactic star clusters.