© The Authors 2025

1 Introduction

The role of stellar rotation in determining the observed features in the color-magnitude diagram (CMD) of young and intermediate-age clusters is now well established. Clusters younger than 500–600 Myr show split main sequences (in near-UV and/or blue bands; for a recent review see Li et al. 2024) that are due to a bimodal stellar rotational distribution (e.g., Kamann et al. 2023). Stars on the blue main sequence (bMS) are slow rotators, while the red main sequence (rMS) is populated by rapidly rotating stars. Additionally, clusters with ages up to ∼2 Gyr show main sequence turnoffs (MSTOs) that are significantly broader compared to expectations of simple stellar populations; these broad MSTOs are also caused by the spread in stellar rotation rates (e.g., Dupree et al. 2017; Kamann et al. 2020). Beyond ∼2 Gyr, the extended MSTOs disappear abruptly as the turnoff mass drops below the threshold where stars become magnetically braked and, hence, all stars become slow rotators (Martocchia et al. 2018; Georgy et al. 2019).

The origin of the bimodal rotational distribution, however, is still uncertain. To date, three hypotheses have been put forward in the literature. The first is that tidal interactions within binary systems cause the stars to become tidally locked and subsequently their rotation slows (D’Antona et al. 2017). In this scenario, the slowly rotating stars (i.e., the bMS), would be made up of binary stars that have been braked through binary interactions. In this scenario, we would expect the bMS to have a much higher binary fraction than the rMS, in conflict with observations that have shown that they are similar (Kamann et al. 2021); see Appendix A.

The second scenario is also based on binary stars. Wang et al. (2022) suggest that stellar mergers can result in a merger product that is a slowly rotating star. They suggest that the mergers happen mostly during the pre-main sequence (PMS) or shortly thereafter; however, the process can continue at a lower rate during the subsequent few tens or even hundreds of megayears. Here, the bMS is formed due to two effects. The first is that the stars are slowly rotating (and hence hotter on average for a given luminosity), and the second is that the stars are rejuvenated, i.e., the merger product resets at the zero-age main sequence. It should be noted, however, that a key assumption of this model is that all or most merger products result in a slowly rotating star. This is not a given, as theoretically, depending on the orientation of the orbit or collision, the merger product can result in a rapid rotator (e.g., de Mink et al. 2013). Simulations of merging stars have suggested that the merger product can lose a large amount of angular momentum during a “bloated phase” directly after the merger (Schneider et al. 2019).

The final scenario, put forward by Bastian et al. (2020), is that the rotation rate of a star is dictated by whether a star has been able to retain its circumstellar disk throughout its PMS phase (slow rotators), or whether the disk gets disrupted due to X-ray/UV photoionization or dynamical encounters (fast rotators). This scenario is based on the idea that circumstellar disks transport angular momentum away from the collapsing PMS star, so stars that retain their disks become slow rotators and stars that lose their disks become fast rotators. This scenario is best tested in very young clusters whose stars may still retain their disks. The prediction is that there should be a strong anticorrelation between the presence of a disk and the rotation rate of the stars. Bu et al. (2024a) report this trend in the young Galactic star-forming region NGC 2264. Further studies with larger samples of stars in more regions are required to confirm this initial finding.

In the binary interaction model (the first scenario discussed above), the expectation is that the fraction of binaries amongst the bMS stars should be near unity due to the fact that in order for a star to become a slow rotator, a close and tidally interacting companion must be present. However, no evidence of an elevated binary fraction amongst bMS stars has been found in spectroscopic observations of bMS and rMS stars (Kamann et al. 2021; He et al. 2023; Bu et al. 2024b).

In the binary merger model (the second scenario), which is the focus of the current study, the blue straggler stars located blueward of the extended MSTO are of particular relevance. Wang et al. (2022) considered them as a natural extension of the bMS, hence making the implicit assumption that they are merger products. We note that while detailed studies of high-mass blue stragglers are still scarce, studies of low-mass stars in old star clusters have shown that many blue stragglers form via stable mass transfer rather than mergers (e.g., Geller & Mathieu 2011; Gosnell et al. 2019)¹. In this work, we investigate the validity of the assumption that blue stragglers form an extension of the bMS, i.e., a blue upper main sequence (bUMS), rather than being a separate phenomenon. For clarity, we highlight specific regions discussed in the current work in Fig. 1.

The binary merger model can be tested directly by comparing the rotation rates (V sin i) of stars on the bUMS with those observed on the bMS. In the Wang et al. (2022) scenario, the rotation rates should be statistically the same, given that the same process is responsible for both populations. If, however, the blue stragglers and the bMS are formed by different processes, then the rotation distribution may be very different.

We used a combination of Hubble Space Telescope (HST) photometry and Very Large Telescope (VLT) Multi Unit Spectroscopic Explorer (MUSE) spectroscopy of large samples of stars in three young clusters in the Large Magellanic Cloud (LMC), namely NGC 1850 (90 Myr, 2184 stars), NGC 1866 (200 Myr, 2531 stars), and NGC 1856 (300 Myr, 4406 stars), to locate stars in the CMD and measure their rotation rates. The data are presented in Sect. 2, and we discuss the results and their implications in Sect. 3.

2 Data

We used two independent datasets, both of which have been presented in detail in other works. Firstly, we used the published catalogs of HST photometry for all three clusters from Niederhofer et al. (2024), which include the F336W, F438W, and F814W bands, along with proper motion (PM) measurements. Proper motions were used to separate stars belonging to the clusters from the surrounding field stars. The NGC 1856 data have been corrected for differential extinction, while no such corrections were needed for, or applied to, NGC 1850 and NGC 1866. We note that the extinction correction has a negligible impact on the results of this study.

NGC 1850 has a projected nearby neighbor, NGC 1850b, whose young age (<10 Myr; Gilmozzi et al. 1994) results in its main sequence stars overlapping with the bUMS and blue straggler stars of NGC 1850. Fortunately, the two clusters have PMs that can be used to separate the two populations. We applied a PM cut of μ_α cos δ>1.9 to select NGC 1850 stars along with a spatial cut of RA < 77.172. No such complications were seen in the NGC 1866 or NGC 1856 data, as these clusters are clearly separated from their surroundings.

Secondly, we used VLT/MUSE integral field spectroscopy to measure the rotation rates (V sin i) of stars in all three clusters. Measurements for stars in NGC 1850 were taken from Kamann et al. (2023), while those for NGC 1866 and NGC 1856 were taken from Kamann et al. (2025). We refer the interested reader to those works for details of the observations, spectral extraction, and data analysis. Stellar rotation rates were measured for each star following the procedures presented in Kamann et al. (2021) and Kamann et al. (2023).

To make a quantitative comparison between the rotation rate distribution of stars on the bUMS and on the bMS, we selected stars in the regions shown in Fig. 2. The bMS selection region was chosen based on the split main sequence. The bUMS region was selected to find classical blue straggler stars (i.e., stars that are brighter and bluer than the MSTO region of each cluster).

Fig. 1 Observed CMD of NGC 1866 with specific regions and features highlighted.

3 Results and discussion

Wang et al. (2022) highlighted that a significant difference in the predictions between the binary merger model and the disk evolution model can be found in stars on the bUMS (see Fig. 1). In the binary merger model, these stars are effectively the more massive counterparts to the bMS stars (see Fig. 3 and Supplementary Fig. 14 from Wang et al. 2022). In the model, all stars on the bMS and the bUMS are merger products. As such, they should display the same rotational velocity distributions. Specifically, both populations are expected to be made up of slow rotators.

Alternatively, in the disk evolution model, the bMS is made up of stars that were able to retain their circumstellar disks during their early evolution and have been braked to low V sin i values. The bUMS is an unrelated phenomenon, effectively just made up of blue straggler stars, which indeed are likely to be merger products or mass transfer systems. In this case, we might expect that the two populations (bMS and bUMS) would have different V sin i distributions.

We used the datasets for the three clusters to search for such differences and similarities between the bUMS and the bMS by comparing the median, mean, and standard deviation of their V sin i distributions. The results are shown in Table 1. As can be seen in Fig. 2, the rotational distribution in the bUMS and the bMS are very different, with the bUMS having a relatively high fraction of rapid rotators, which are unseen on the bMS.

We also used a Kolmogorov-Smirnov (KS) test to determine the probability that each of the populations was drawn from the same parent distribution. These numbers are also shown in Table 1. In all three clusters we find a low probability of the two distributions (bUMS and bMS) being drawn from the same parent populations (p<0.025 in all clusters). This is in direct conflict with the predictions and assumptions of the Wang et al. (2022) merger scenario

To further quantify the differences, we combined populations on the bUMS and the bMS for each cluster in order to improve our number statistics. The cumulative distributions are shown in Fig. 3. A KS test finds a low (P=0.014) probability that they were drawn from the same parent population.

From these tests we find that, contrary to the predictions of the Wang et al. (2022) binary merger scenario, the bUMS and the bMS have different V sin i distributions. This suggests that the origin of the two populations is different. We conclude that the bMS is not made up of binary merger products, and that the blue straggler population within young clusters is not related to the observed bMS.

A corollary of these results is that not all blue straggler stars in young clusters are slow rotators. Hence, whatever the cause of blue straggler stars, be it stellar mergers, interactions, or mass transfer, it does not necessarily lead to a slowly rotating star.

Our results qualitatively agree with investigations into the spin distributions of blue stragglers in old globular clusters, where the formation of blue stragglers is widely believed to result from a mixture of mass transfer and stellar collisions. The V sin i distribution of blue stragglers found in Galactic globular clusters is typically characterized by a dominant population of slow rotators and a tail that extends to larger values (e.g., Simunovic & Puzia 2014; Mucciarelli et al. 2014; Billi et al. 2023), with the contribution of the latter decreasing with increasing cluster density (Ferraro et al. 2023). The dashed line in Fig. 3 resembles such a distribution, although for a higher V sin i on average. In globular clusters, the fast rotating blue stragglers can be linked to recent interactions because of the fast spin-down times of lowmass stars (e.g., Leiner et al. 2018). However, such a link cannot be readily adopted for our targets, because of the young cluster ages and the lack of braking experienced by early-type stars. Instead, one might speculate that the fast rotators in our sample are the result of mass transfer, while the slow rotators predominantly originate from stellar collisions (e.g., Bailyn 1992). Future studies of the binary properties of both types of blue stragglers may provide more definite conclusions on this point.

We only tested the specific model of Wang et al. (2022) and not a general binary merger model. Our results suggest that, in analogy to older clusters, bUMS/blue stragglers in younger clusters are formed through multiple physical processes, i.e., mass transfer systems and/or stellar mergers. These processes could be included in future models. However, we note that such gainer stars would be expected to also be found on the main sequence², as being a gainer is independent of a star’s evolutionary status. Hence, in this case we can still expect similar V sin i distributions between the bMS and the bUMS. Additionally, such models would need to predict the relative fractions of bUMS and bMS stars. Finally, we note that the estimated age distributions of the binary mergers in the Wang et al. (2022) model are highly sensitive to the exact shape of the isochrones adopted, which are known to vary between different stellar evolutionary models (see Appendix B).

Fig. 2 Observed CMDs of NGC 1850 (left), NGC 1866 (middle), and NGC 1856 (right). Only stars that were considered cluster members after PM and radial velocity cuts are shown. Boxes highlight the regions in each cluster that were selected for bUMS stars (upper rectangles) and bMS stars (lower parallelograms).

Fig. 3 Cumulative V sin i distribution of the combined bMS (solid blue line) and bUMS (dashed blue line) from all clusters. A KS test reveals that the two populations differ significantly (p=0.014). We also show the V sin i distribution of the rMS (dot-dashed red line) of the three clusters combined, selected from stars with similar magnitudes as the bMS.

Acknowledgements

C.W., S.d.M. and N.L. are thanked for useful discussions on the binary merger model. S.K. acknowledges funding from UKRI in the form of a Future Leaders Fellowship (grant no. MR/T022868/1). F.N. acknowledges funding from DLR grant 50 OR 2216. S.S. acknowledges funding from the European Union under the grant ERC-2022-AdG, “StarDance: the non-canonical evolution of stars in clusters”, Grant Agreement 101093572, PI: E. Pancino. Based on observations made with ESO Telescopes at La Silla Paranal Observatory. Based on observations of the NASA/ESA Hubble Space Telescope, obtained from the data archive at the Space Telescope Science Institute. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS5-26555.

Appendix A A reanalysis of the binary fractions of the blue and red main sequences

In addition to the stellar merger hypothesis for the origin of the bimodal stellar rotation distribution, D’Antona et al. (2017) have suggested that tidal forces between binary components could slow down rapid rotators. The authors suggest that if all or most stars in clusters are born as fast rotators, tidal locking may slow down a significant fraction of these stars, leading to the bMS. A prediction of this model is that all or most of the bMS stars should be in binary systems, specifically in close binary systems where tidal forces can lock the stars into a synchronous orbit, effectively slowing each star down. Kamann et al. (2021) tested this scenario by spectroscopically measuring the binary fraction through radial velocity variations of both the bMS and rMS in the young massive cluster, NGC 1850. They found that both sequences had a low (and statistically identical) binary fraction (∼5%), in conflict with predictions of the tidal forces scenario.

Recently, Muratore et al. (2024) have revisited this problem by estimating the binary fractions in the bMS and rMS in three young massive clusters in the LMC. These authors used multicolor HST photometry to identify binaries in certain filter combinations and then associate them with the different populations, using other combinations. The method is described in detail in Milone et al. (2020).

Muratore et al. (2024) report a much higher binary fraction amongst bMS stars (23 ± 4%) relative to rMS stars (5 ± 4%) in NGC 1850, and a similar trend in the other two clusters. This result broadly supports the tidal locking scenario; however, it is in conflict with the spectroscopic binary estimates (Kamann et al. 2021). The authors offer a possible explanation, that tidal locking is more efficient than previously thought; essentially, binary stars whose orbital separation is large enough to avoid spectroscopic detection might be able to still tidally lock (i.e., soft binaries might be able to tidally lock). While this is a possibility, it remains unclear why the two populations (rMS and bMS) would have similar hard binary populations.

In this section we reevaluate the Muratore et al. (2024) method and subsequently their results using synthetic clusters. Effectively, the method adopted is to look at a CMD (specifically m_F275W - m_F336W vs. m_F336W) where the bMS and rMS overlap, in order to select high mass ratio binaries. By then looking at the position of the binary systems in a different CMD projection (namely, m_F336W - m_F18W vs. m_F336W), where the two populations are separated, the relative binary fractions can be inferred by looking at the resulting distribution compared to equal mass binary sequences. This is done by looking at the binary systems positions in color space (after verticalization of the main sequence) in combination with simulated populations. Effectively, the cumulative color distribution of the observations is compared to that of synthetic cluster populations with different input parameters (binary fractions of the red and blue populations, fraction of rMS/bMS stars, mass distributions, etc) in order to select the best fitting model.

A potential cause for concern in the analysis of Muratore et al. (2024) is the relatively low number of N≲50 binaries in their photometric selection boxes. When stochastically sampling N data points from two distributions, with probabilities of F_bin,red and F_bin,blue = 1−F_bin,red, the number of actually drawn rMS and bMS binaries will be subject to Poisson noise. This limits the precision with which F_bin,red can be recovered from the numbers drawn. If we adopt F_bin,red and N=38, based on Fig. 3 in Muratore et al. (2024), we find that F_bin,red can be recovered to an accuracy of 8%, which is already higher than the uncertainties provided for NGC 1850 in Muratore et al. (2024). Contrary to this simple estimate, for the real data, we cannot tell for each individual binary whether it belongs to the bMS or the rMS sample, because the two sequences overlap. This will likely substantially increase the uncertainties compared to 8% stated above. We note that Muratore et al. (2024) used mock samples containing about 10× the number of stars as their real samples to estimate their uncertainties. This implies that stochastic effects play a much smaller role in the mock samples than they do in the real data, potentially leading to unrealistically low uncertainties.

To directly test the method, we created a number of synthetic cluster datasets with the nonrotating MIST isochrones at LMC metallicity (Choi et al. 2016) and a Kroupa initial mass function (Kroupa 2001). We adopted an age of 100 Myr (e.g., Bastian et al. 2016). Each cluster contains two populations, a bMS (containing one third of the stars) and a rMS (two thirds of the stars), offset by ∼ 0.1 mag in color in m_F336W - m_F814W (shifted by hand to match observations) but overlapping in m_F275W - m_F336W. We assumed that the binary fraction of the red population is constant at 0.1(F_bin,red). For the blue population, we adopted different binary fractions (F_bin,blue) of 0.1, 0.2, or 0.4. For all binaries, we adopted a flat mass-ratio distribution. For each combination of parameters, we created ten synthetic clusters that are stochastically sampled.

In a first experiment, we created stochastic clusters with ∼ 10 times more stars than observed (i.e., for NGC 1850 there are 38 binary systems in the selection box used by Muratore et al. 2024), meaning that the synthetic clusters have 350 to 450 binary systems in the selection box. This is meant to broadly reproduce the simulations used in Muratore et al. (2024). The resulting cumulative distributions are shown in the left panel of Fig. A.1.

From this figure we can see that the overall method is sensitive to the binary fraction, and that the cumulative distributions are significantly different for different F_bin,blue values. Of course, all other parameters in the model are fixed, leading to an idealized simulation, but the results are clearly very sensitive to F_bin,blue.

In a second set of simulations, we attempted to reproduce the observations more accurately, specifically the number of binary systems that fall within the selection box. Muratore et al. (2024) found 38 binary systems (that they used for the analysis) within their selection box. For this second set of simulations, the synthetic clusters have the same set of parameters as previously adopted, but now they have only 35-45 binary systems that fall into the selection box.

The results are shown in the right panel of Fig. A.1. As can be seen, in the low-N limit (and we note that NGC 1850 had the highest number of binary systems in the selection box of the three clusters studied by Muratore et al. 2024) the diagnostic power of this method becomes significantly worse. In this case, even synthetic clusters with radically different F_bin,blue values overlap, meaning that no significant constraints can be put on the binary fractions of the two populations. Again, we note that in these idealized simulations, all other parameters are known and kept fixed, limiting any additional uncertainties.

We conclude that the method adopted by Muratore et al. (2024) is potentially very powerful in the high-N regime. However, in the space occupied by young massive clusters, the method is not able to provide useful constraints. Hence, we conclude that the uncertainties cited by Muratore et al. (2024) are underestimated due to the lack of inclusion of accurate number statistics. This is likely the reason for the discrepancy between the estimated binary fractions between that work and Kamann et al. (2021).

Fig. A.1 Left panel: Cumulative distribution of the synthetic clusters, separated by color for the three different Fbin,blue values. Each synthetic cluster has 350 to 450 binary systems that fall into the selection box, representing the high-number regime. Right panel: Same but now for synthetic clusters with only 35 to 45 binary systems in the selection box, in order to be more consistent with the observed clusters.

Appendix B Magnitude-dependent merger histories

Wang et al. (2022) derive the merger histories of binaries in their sample of clusters by fitting the position of stars along the bMS with rejuvenated stellar models. Under the assumption that each star is a merger product (in this case from an equal mass merger) the position of the star between the current isochrone that fits the bMS (i.e., the age of the cluster) and the zero-age main sequence can be used to determine when the merger happened. This is based on the idea that a merger would reset the main sequence evolution of a star.

Wang et al. (2022) note that toward the MSTO more stars are to the blue of the best fitting bMS isochrone compared to fainter magnitudes. The authors use this to conclude that the binary merger history is magnitude- (i.e., mass-) dependent.

An alternative explanation is that uncertainties in the models lead to small offsets between the observations and the isochrones. This is to be expected to some extent, given the uncertainties and parameterizations within stellar evolutionary models (especially when rotation is considered; see, e.g., Ekström et al. 2018).

As an example, we can look at NGC 330 as studied in Fig. 3 in Wang et al. (2022). At a magnitude of m_F814W ∼ 18.5 we see a cloud of stars that lie m_F336W − m_F814W ∼ 0.05−0.1 mag to the blue of the adopted bMS isochrone (translating into the box in Fig. B.1).

To get a feel on the level of model uncertainties, we compared the adopted Wang et al. (2022) isochrones to those of other groups. Specifically, we adopted the MIST (Choi et al. 2016) and Parsec models (Bressan et al. 2012). We removed the extinction and distance modulus applied to the isochrones used in Wang et al. (2022), adopted the same metallicity for all models (namely Z_SMC = 0.00218), and used the nonrotating models. A comparison of the isochrones are shown in Fig. B.1.

While all three models show a very similar shape and at fainter magnitudes the agreement between the models is excellent. However, we do see differences at brighter magnitudes. Specifically, in the box highlighted in Fig. B.1 (which is where the main offset between the models and data in Wang et al. (2022) is found) the Parsec and MIST models are not as steep as the Wang et al. (2022) models, leading to bluer models. Within the region of interest, the offset is m_F336W − m_F814W ∼ 0.05 mag. This is similar to the reported offset between the Wang et al. (2022) isochrones and the observed stars, highlighting that caution is necessary when comparing small offsets between data and isochrones.

Fig. B.1 Comparison of the isochrones from three different stellar evolutionary codes.

While it is impossible to make any statements about which model is best, this comparison demonstrates that caution is needed when comparing the predictions of stellar isochrones with actual data, as offsets may not have a physical origin.

References

All Tables

All Figures

	Fig. 1 Observed CMD of NGC 1866 with specific regions and features highlighted.
In the text

	Fig. 2 Observed CMDs of NGC 1850 (left), NGC 1866 (middle), and NGC 1856 (right). Only stars that were considered cluster members after PM and radial velocity cuts are shown. Boxes highlight the regions in each cluster that were selected for bUMS stars (upper rectangles) and bMS stars (lower parallelograms).
In the text

	Fig. 3 Cumulative V sin i distribution of the combined bMS (solid blue line) and bUMS (dashed blue line) from all clusters. A KS test reveals that the two populations differ significantly (p=0.014). We also show the V sin i distribution of the rMS (dot-dashed red line) of the three clusters combined, selected from stars with similar magnitudes as the bMS.
In the text

Fig. A.1 Left panel: Cumulative distribution of the synthetic clusters, separated by color for the three different Fbin,blue values. Each synthetic cluster has 350 to 450 binary systems that fall into the selection box, representing the high-number regime. Right panel: Same but now for synthetic clusters with only 35 to 45 binary systems in the selection box, in order to be more consistent with the observed clusters.

In the text

	Fig. B.1 Comparison of the isochrones from three different stellar evolutionary codes.
In the text