Binary-driven stellar rotation evolution at the main-sequence turn-off in star clusters

Extended main-sequence turn-offs (eMSTOs; e.g., Mackey & Broby Nielsen, 2007; Goudfrooij et al., 2011) and split main sequences (MSs; e.g., D’Antona et al., 2015; Li et al., 2017a) pose a fundamental challenge to our traditional understanding of star clusters as ‘simple stellar populations’. Both features are thought to be driven by differences in stellar rotation rates (e.g., Niederhofer et al., 2015), supported by an increasing body of spectroscopic evidence in both Magellanic Cloud clusters (Dupree et al., 2017; Kamann et al., 2020) and Galactic open clusters (OCs; Bastian et al., 2018; Sun et al., 2019a, b). Fast rotators appear redder than their slowly rotating counterparts because of the compound effects of gravity darkening and rotational mixing (Yang et al., 2013).

As the effects of stellar rotation in defining the morphology of cluster colour–magnitude diagrams (CMDs) are now well-understood, the rotation distributions of MSTO stars and their origin in star clusters represent the next key open question. A cluster mostly composed of initially rapidly rotating stars may reproduce the observed ‘converging’ subgiant branch in NGC 419 (Wu et al., 2016). D’Antona et al. (2015) argued that split MSs are composed of two coeval populations with different rotation rates: a slowly rotating blue-MS (bMS) and a rapidly rotating red-MS (rMS) population. This bimodal rotational velocity distribution was confirmed for NGC 2287 (Sun et al., 2019b), where the well-separated double MS in this young OC is tightly correlated with a dichotomous distribution of stellar rotation rates. Kamann et al. (2020) also discovered a bimodal distribution in the rotation rates of MSTO stars in the 1.5 Gyr-old cluster NGC 1846. However, the peak rotational velocities in this cluster are slower than those in the younger cluster NGC 2287, thus offering a hint of stellar rotation evolution.

This naturally raises the question as to the origin of such a bimodal rotation distribution. D’Antona et al. (2015, 2017) suggested that bMS stars in young clusters might have resulted from tidal braking of initially fast rotators on time-scales of a few $\times 10^{7}$ yr. Alternatively, the bimodal rotation distribution could have been established within the first few Myr, regulated by either the disc-locking time-scale (Bastian et al., 2020) or the accretion abundance of circumstellar discs (Hoppe et al., 2020). In this paper, we discover a tight correlation between the stellar rotation rates and binary fractions in five intermediate-age Galactic OCs with similar dynamical ages, favouring a binary-driven model. This is further supported by the greater central concentration of the bMS stars relative to the rMS population in at least two of our clusters. Intriguingly, our observational results are inconsistent with the prevailing theoretical models, which thus invites more detailed future investigations.

This paper is organised as follows. We present our photometric and spectroscopic data reduction in Section 2. Our analysis to unravel the OCs’ rotation distributions and its validation are described in Section 3. A discussion about the formation mechanism determining the stellar rotation rates in clusters and our conclusions are summarised in Section 4.

Several recent attempts have been made to unravel cluster rotation distributions (Gossage et al., 2019; de Juan Ovelar et al., 2020). However, most of these estimates were based solely on their CMD morphologies (see also Kamann et al., 2020), i.e., based on the probability of models matching the observations using Hess diagrams. Whereas this may be suitable for MC clusters, where the numbers of member stars are sufficiently large and where it is difficult or impossible to obtain direct velocity information, such analyses of OCs inevitably suffer from stochastic sampling effects (see Section 3.1.2). Spectroscopic observations can ameliorate these effects.

We used the SYCLIST stellar rotation models (Georgy et al., 2013, 2014) to generate synthetic clusters, assuming solar metallicity ($Z=0.014$), an age of $\log(t\mbox{ yr}^{-1})=8.95$ and a 50 per cent binary fraction. The model also considers limb darkening (Claret, 2000). We adopted the gravity-darkening law of Espinosa Lara & Rieutord (2011). The inclination angles follow a random distribution. The model suite is limited to a minimum stellar mass of $1.7$ M_⊙ and it does not account for the evolution of interacting binary systems.

Since the SYCLIST models are limited to high masses ($\geq 1.7$ M_⊙; suitable for covering the MSTO) and it is difficult to constrain a cluster’s binary fraction based on the morphology of its MSTO, we exploited the less massive MS stars (see Fig. 3) to derive total binary fractions, $f_{\mathrm{b}}$, at masses where the effects of stellar rotation are negligible. Binary models were generated from the MIST isochrones by adding unresolved binaries with different mass ratios, assuming a flat mass-ratio distribution. Next, we compared the pseudo-colour distributions of stars in our synthetic clusters with those in the observed clusters, using the same MS selection criteria. The goodness of the comparison is given by the Anderson–Darling $p$ value. We added a second parameter, $\sigma$, to characterise the observational scatter in the pseudo-colour, combining the effects of a possible internal age spread, photometric uncertainties and differential extinction. A minor colour shift among the clusters was also taken into consideration. We employed the Markov-chain Monte Carlo method (emcee; Foreman-Mackey et al., 2013) to determine the best model and estimate the uncertainties in the resulting parameters.

We adopted this binary fraction as our input parameter and applied a similar method to the MSTO stars to derive their rotation distribution. In the SYCLIST models, on the zero-age MS (ZAMS) this distribution is controlled by the ratio of the equatorial angular velocity, $\Omega$, to the critical velocity $\Omega_{\mathrm{crit}}$. The latter is the rotational angular velocity where the surface gravity can no longer maintain equilibrium with the centrifugal motion. We binned the resulting rotation rates into 10 bins, from $\Omega/\Omega_{\mathrm{crit,ZAMS}}=0.0$ to 1.0 in steps of 0.1, and assigned 10 weights, $F_{\Omega_{i}}$. Each weight corresponded to the fraction of stars with a given rotation rate.

We first generated a model cluster of 200,000 stars with a flat rotation distribution. For each set of input parameters ($f_{\mathrm{b}}$, $F_{\Omega_{i}}$), we randomly selected the $N$ most appropriate stars from the total pool to generate a new synthetic cluster which satisfied the rotation distribution required. Subsequently, we added scatter ($\sigma$) to the synthetic cluster and selected the remaining MSTO stars for comparison with the observed cluster. In practice, when selecting MSTO stars and calculating their pseudo-colours, we shifted the ridge line towards bluer colours by 0.3 mag to include MSTO stars which have scattered away from the MSTO region.

As for the MS region, a comparison of the pseudo-colour distributions of MSTO stars in the synthetic and observed clusters yields the corresponding $p_{\mathrm{phot}}$ values. We also used the rotational velocities to address the degeneracy between the pseudo-colour and rotation-velocity distributions (see Section 3.1.2). For each MSTO star with spectroscopic measurements, we selected the closest 100 stars in the synthetic cluster’s CMD. We then computed the probability of deriving the same $v\sin i$ distribution as observed. We first resampled the measured projected rotational velocities according to the error distribution and estimated the probabilities for this alternative realisation. We iterated each run 100 times and adopted the median value to minimise stochastic sampling effects. Finally, we combined the corresponding $p_{\mathrm{spec}}$ values from the spectroscopic data with $p_{\mathrm{phot}}$ and used emcee to find the rotation model which best reproduces both the pseudo-colour and the $v\sin i$ distributions. The uncertainties were derived from the samples’ 16th and 84th percentiles in the marginalised distributions.

The best-fitting rotation rates’ weights are shown in Fig. 6 (small panels), whereas the rotation rate probability distributions of the MSTO stars for all clusters are shown in the top left-hand panel. All five clusters contain a significant fraction of fast rotators. However, the locations and fractions of fast rotators vary among the clusters. Whereas some clusters, like NGC 6134 and NGC 2818, appear to host three rotational populations (including two slowly rotating populations), it is difficult to reliably separate between populations of slow rotators because of their minor colour differences. Therefore, we exercise caution and suggest that these clusters may only contain two distinct populations (i.e., a bimodal distribution) in stellar rotation space.

The number fraction of slow rotators, $N_{\mathrm{slow}}/N_{\mathrm{tot}}$, among synthetic MSTO stars shows a positive correlation with the clusters’ total binary fraction $f_{\mathrm{b}}$: see Fig. 7. We verified that our approach to selecting slowly and rapidly rotating subsamples, within reasonable ranges, minimally affects this correlation by changing the velocity adopted for our subsample selection to fixed values. In all tests, the correlation between the number fraction of slow rotators and the binary fractions of the clusters remains almost the same, modulo minor shifts in the absolute value of $N_{\mathrm{slow}}/N_{\mathrm{tot}}$ In fact, one slow rotator ($v\sin i=123.1\,\mathrm{km\,s^{-1}}$) in NGC 6134 (Gaia DR2 ID 5941409684301615360) exhibited variation in its radial velocity during two observation epochs, spanning three days. The observed variation is around $80\,\mathrm{km\,s^{-1}}$, which is four times larger than the uncertainty associated with the velocity measurements. The estimated rotational velocities remained unchanged, with differences of less than the $v\sin i$ error ($7\,\mathrm{km\,s^{-1}}$). If this radial-velocity variation were induced by binaries, the star is likely a member of a close binary system with a separation of a few tens of solar radii. To arrive at this estimate, we assumed that its luminosity is not severely affected by the companion and we adopted the observed variation as the amplitude of radial velocity curve. Follow-up time-series observations are required to confirm this star’s nature.

The approximately constant number ratios of bMS stars found in massive MC clusters (Milone et al., 2018) may result from the similar binary fractions ($\sim 0.3$) which are prevalent in MC clusters (Milone et al., 2009) if the correlation of Fig. 7 also pertains to more massive clusters. For instance, Milone et al. (2017) estimated that the number ratio of bMS stars in NGC 1866 is 35 per cent; the cluster’s overall binary fraction is 0.28. These properties will place it comfortably on the apparent correlation. However, some younger clusters, e.g., NGC 1818 ($\sim 40\,\mathrm{Myr}$) and NGC 330 ($\sim 40\,\mathrm{Myr}$), do not follow the correlation if we adopt their binary fractions from Li et al. (2017a) and their number ratios from Milone et al. (2018). This may suggest that these young clusters are still in the early stages of their stellar rotation evolution, which may last for several tens of Myr. In NGC 1846 ($\sim 1.5\,\mathrm{Gyr}$), Kamann et al. (2020) found that $\sim 45$ per cent of MSTO stars are slow rotators. They derived a binary fraction of $\sim$7 per cent based on radial velocity variations from four-epoch spectra. However, because their binary fraction measurement is limited by their temporal coverage and most sensitive to the detection of tight binaries, a large fraction of the cluster’s binaries are likely missing. Therefore, we did not include this cluster in Fig. 7.

To date, none of the prevailing theoretical models (e.g., D’Antona et al., 2015; Bastian et al., 2020) can naturally explain the observed correlation. Our observational result is in apparent conflict with the model of Bastian et al. (2020), which suggested that binaries may destroy stellar discs, thus resulting in faster rotational velocities on the MS. Moreover, this interpretation cannot explain the cluster-to-cluster variations in the MSTO pseudo-colour distributions in Fig. 3 (top right), which we verified using a two-sample Kolmogorov–Smirnov test. This is so, because our result is not an artefact of requiring a higher binary fraction to compensate for the potential narrowing of the MSTO owing to a lower fraction of fast rotators.

Neither does our correlation follow the scenario of D’Antona et al. (2015), where all or most stars form as fast rotators with a fraction undergoing close-binary tidal braking. Although this model predicts a larger fraction of slow rotators for a higher binary fraction, everything else being equal, the fraction of slow rotators can never exceed the binary fraction. Since only close binaries may become tidally locked, only a small subset of binaries could contribute. However, as shown in Fig. 7, the slow rotators’ number ratios are comparable to or even higher than the binary fractions. Moreover, the slope between $N_{\mathrm{slow}}/N_{\mathrm{tot}}$ and $f_{\mathrm{b}}$ is close to unity regardless of the subsample selection, thus suggesting an alternative binary-driven mechanism. However, we cannot rule out a possible origin associated with the initial phase. The slow rotators we observe at the present time originate from both the initial population of slow rotators and may also have evolved from the fast rotators. Huang et al. (2010, their figure 7) selected a young subpopulation of B-type stars that just evolved from the zero-age main sequence and showed that a large fraction of them were formed as slow rotators. So far, it is unclear which physical property of (close and wide) binaries determines the rotation distributions. This may explain why Kamann et al. (2020) did not find any significant differences among the binary fractions of slow and fast rotators. If we consider the possibility of missing binaries from their sample, combined with the disruption of wide binaries, there could be distinct differences between slow and fast rotating populations.

If binaries are the dominant drivers of slow rotators, one would expect slow rotators to be more centrally concentrated in a cluster because they are generally more massive than their rapidly rotating counterparts, and they would thus sink more easily to the cluster centre owing to two-body relaxation (dynamical mass segregation; Binney & Tremaine, 1987). However, this scenario appears at odds with recent observations in young MC clusters. Dupree et al. (2017) and Milone et al. (2017) found that in NGC 1866 ($\sim 400\,\mathrm{Myr}$), bMS stars (slow rotators) are less centrally concentrated than rMS stars (fast rotators). Meanwhile, in NGC 1856 ($\sim 300\,\mathrm{Myr}$-old), the number ratios of bMS and rMS stars remain unchanged at different radii, thus suggesting that they are spatially homogeneously distributed (Li et al., 2017b).

We used the best-fitting synthetic cluster to infer the rotational velocities for all member stars (Sun et al., 2019a) and employed the spatial locations of the observed stars to estimate their degree of mass segregation. Because of their low number densities, it is hard to robustly determine the centres of OCs. Therefore, we adopted minimum spanning trees (MSTs) to quantify a cluster’s degree of mass segregation (Allison et al., 2009). The mass segregation ratio, $\Lambda_{\mathrm{MSR}}$, of a given population is defined as the ratio of the average random path length, $l_{\mathrm{random}}$, to that of the entire population, $l_{\mathrm{pop}}$,

where $l$ represents the length of the shortest path connecting all data points and $\langle l_{\mathrm{random}}\rangle\pm\sigma_{\mathrm{random}}$ quantifies the length distribution. We determined the OCs’ MSTs using MiSTree (Naidoo, 2019), in celestial coordinates. Given our OCs’ close proximity and the low stellar densities, the impact of sampling incompleteness of the low-mass stars in estimating $\Lambda_{\mathrm{MSR}}$ is negligible.

To link our results with the observations of blue and red MSs, we adopted bMS stars as those stars with $v\sin i<150\,\mathrm{km\,s^{-1}}$, whereas rMS stars have $v\sin i>150\,\mathrm{km\,s^{-1}}$. This velocity threshold was selected based on the velocity dip at $v\sin i\approx 100\,\mathrm{km\,s^{-1}}$ observed in NGC 1846 (Kamann et al., 2020). Given the mass differences of MSTO stars in clusters with different ages, this corresponds to $v\sin i\approx 150\,\mathrm{km\,s^{-1}}$ for our sample. We assumed that the rotation rate ($\Omega/\Omega_{\mathrm{crit}}$) of the dip does not change for clusters of either 1 Gyr or 1.5 Gyr (NGC 1846). The critical rotation velocity only marginally decreases as a star evolves, whereas the rotation velocity strongly depends on stellar mass (Bastian et al., 2020, their figure 1). The typical MSTO stellar mass is around $1.5-1.6\,\mathrm{\mathrm{M}_{\odot}}$ and $1.8\,\mathrm{\mathrm{M}_{\odot}}$ for NGC 1846 and our cluster sample, respectively. Based on Georgy et al. (2014), the velocity dip for bMS and rMS should be $\sim$1.5 times larger in a 1 Gyr-old cluster than in NGC 1846.

The corresponding $\Lambda_{\mathrm{MSR}}$ values are listed in Table 1. In all clusters, except for NGC 2818, the bMS stars exhibit significant spatial segregation, $\Lambda_{\mathrm{MSR}}>1$. Moreover, $\Lambda_{\mathrm{MSR,b}}>\Lambda_{\mathrm{MSR,r}}$ in NGC 3860 and NGC 5822, with differences $>2\sigma$. This suggests that the bMS stars in these clusters are more centrally concentrated than their rMS counterparts. This, hence, offers promising supporting evidence of a more massive origin of the bMS stars in these OCs. However, this is at odds with MC cluster results (e.g., Milone et al., 2017; Kamann et al., 2020). The difference could be attributed to binary disruption in the cluster centre. Milone et al. (2017) found the binary fraction to increase towards the cluster outskirts, following the same radial distribution as the bMS stars in NGC 1866, whereas this effect is not obvious in OCs.

In summary, we have discovered a strong observational correlation between the number ratios of the slow rotators and the binary fractions in five Galactic OCs and shown support for a marked concentration of bMS stars. Future work will extend the survey to younger clusters, cover a larger parameter space (e.g., in $f_{\mathrm{b}}$, stellar and cluster mass, metallicity, etc.) and study the detailed history of stellar rotation in cluster environments to gain additional important insights into the star cluster formation processes.

W.S. thanks Xiao-Wei Duan for insightful discussions. He is grateful for financial support from the China Scholarship Council. L.D. acknowledges research support from the National Natural Science Foundation of China through grants 11633005, 11473037 and U1631102. This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/web/gaia/dpac/consortium).

The data analysed in this paper will be shared upon request.