EXPLANATION OF A SPECIAL COLOR–MAGNITUDE DIAGRAM OF STAR CLUSTER NGC 1651 FROM DIFFERENT MODELS^*

Many intermediate-age star clusters, e.g., NGC 1651, show complicated color–magnitude diagrams (CMDs), which include blue stragglers, a broad main sequence, an extended main sequence turn-off (eMSTO), and an extended red clump (eRC). Such clusters are found mainly in the Large Magellanic Cloud (LMC; Mackey & Broby Nielsen 2007; Mackey et al. 2008; Milone et al. 2009; Piatti 2013; Kerber et al. 2007; Rubele et al. 2011), and Small Magellanic Cloud (SMC; Girardi et al. 2009, 2013). There has been a long history of explanation for such special CMDs. Factors considered as possible causes for CMDs with eMSTO and eRC features include the spread of chemical abundance (Piotto et al. 2005, 2007; Mackey et al. 2008; Goudfrooij et al. 2009), the capture of field stars (Mackey et al. 2008; Goudfrooij et al. 2009), the merger of existing star clusters (Mackey & Broby Nielsen 2007), the formation of a second-generation of stars from the ejecta of first-generation asymptotic giant branch stars (D'Ercole et al. 2008; Goudfrooij et al. 2009), binary stars (e.g., Milone et al. 2009), observational selection and uncertainty effects (Keller et al. 2011), the mixture of stars with and without overshooting (Girardi et al. 2011), differential reddening (Platais et al. 2012), age spread (Girardi et al. 2011; Richer et al. 2013), stellar rotation (Bastian & de Mink 2009), and the combination of binaries and rotation (Li et al. 2012b). Meanwhile, some other studies, including Mucciarelli et al. (2008), Goudfrooij et al. (2011a, 2011b), Glatt et al. (2008), Rubele et al. (2010), and Girardi et al. (2011, 2013), pointed out some challenges to these assumptions. Finally, a spreading in age (e.g., Girardi et al. 2011) and stellar rotation (e.g., Bastian & de Mink 2009; Li et al. 2012b) are thought to be the most probable causes for special CMDs with eMSTO and eRC, because only such models can reproduce the eMSTO part sufficiently well.

Some special CMDs of clusters have been compared with theoretical models, but no work compares the results from multiple different stellar population models simultaneously. Such an approach can, however, provide new insight into the relative importance of various factors involved in the observed CMD. This work aims to supply such an attempt via a typical cluster, NGC 1651. This is a red globular cluster in the LMC, which is situated 3 kpc southwest of the Bar in a region that appears fairly free of recent star formation. The CMD of NGC 1651 has been studied before, e.g., by Mould et al. (1986, 1997), Brocato et al. (2001), Sarajedini et al. (2002), and Grocholski et al. (2007). Although not very clear, special CMD structures with eMSTO and eRC, which have first been described by Mackey & Broby Nielsen (2007) and Girardi et al. (2009), can be seen in the results of Brocato et al. (2001). The special CMD of NGC 1651 was recently confirmed by the present authors using data from the Hubble Space Telescope (HST) archive, in which blue stragglers, broad MS, eMSTO, and eRC structures seem very clear. We intend to study this CMD via various model assumptions. The study will help to clarify the roles of various factors in explaining the observed CMDs and in determining the cluster parameters. In the most complex model, the effects of binaries, age spread, rotating stars, and star formation history are taken into account, in addition to those of distance, extinction, and metallicity. The structure of this paper is as follows. Sections 2 and 3 introduce the observed and theoretical CMDs. Section 4 introduces the technique of CMD fitting. Section 5 presents the best-fit results from different models. Finally, Section 6 contains a summary and discussion.

2.1. Data and Photometry

2.2. Photometric Errors and Completeness

A kind of photometric uncertainty has been reported by HSTphot, but it is obviously less than the real error because simulated CMDs are much narrower than the observed one when the error reported by HSTphot is taken into account. Furthermore, the sample seems incomplete because the stars of NGC 1651 are in a crowded field. This will surely lead to some extra errors. We therefore perform artificial star tests (ASTs) to characterize the photometric errors and completeness of stars. Via this technique, most uncertainties caused by crowding and the photometry process can be estimated. The ASTs are performed as follows. First, a large number of images that include more than 10$^{5}$ artificial stars in the same color and magnitude ranges as the observed ones are generated. For convenience, the real magnitudes of artificial stars are recorded as input magnitudes. Then the generated images are processed by HSTphot to find out stars and measure their magnitudes. The measured magnitudes are recorded as recovered magnitudes. Finally, the input and recovered magnitudes of artificial stars are compared, to characterize the photometric errors, and the input and recovered star numbers are compared, to provide star completeness. Figure 1 shows the photometric errors as a function of V and I input magnitudes, and Figure 2 shows the completeness of different CMD parts of NGC 1651. We see that photometric errors mainly depend on recovered magnitudes. The fainter the magnitudes, the larger the spread of photometric errors. In addition, we find that the completeness of various CMD parts are different. Fainter areas have lower completeness on average. Thus we applied position-dependent completeness to the observed CMD. Note that the results for where observed stars are located are highlighted in Figure 2, in order to help readers understand the observed CMD better. Figure 3 shows some compensating stars that are randomly added into the observed CMD according to star completeness. Because the added stars are far fewer than the observed stars, they do not affect the final results too much. The final observed CMD of NGC 1651, which has included compensating stars, is shown by Figure 4. This CMD consists of 4803 stars and is much clearer than what is used by previous works (Mould et al. 1997; Sarajedini et al. 2002). We can see broad MS, eMSTO, and eRC (see Section 5 for details) structures clearly from this CMD. It is therefore a good example for our detailed study.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We construct synthetic CMDs following previously cited works (e.g., Zhang et al. 2004; Li & Han 2008a, 2008b; Li et al. 2012a). Simple and composite stellar populations (SSPs and CSPs) are built by taking the initial mass function (IMF) of Salpeter (1955) with lower and upper mass limits of 0.1 and 100 ${{M}_{\odot }}$, respectively. We build both single star stellar populations (ssSPs) and binary star stellar populations (bsSPs). The difference between the two kinds of populations is that bsSPs contain some binaries but ssSPs do not. In each binary, the mass of the primary component is generated following the selected IMF, and the mass of the secondary component is calculated by taking a random secondary-to-primary mass ratio (q), which obeys a uniform distribution within 0–1. The eccentricity (e) of each binary is given randomly within 0–1. The separations (a) of two components are given by

$$\begin{eqnarray}an(a)=\left\{ \begin{array}{ccccccccccccccc} {{\alpha }_{{\rm sep}}}{{(a/{{a}_{0}})}^{m}} & a\leqslant {{a}_{0}}; \\ {{\alpha }_{{\rm sep}}}, & {{a}_{0}}\lt a\lt {{a}_{1}}, \\ \end{array} \right.\end{eqnarray} \tag{ 1 }$$

where ${{\alpha }_{{\rm sep}}}\approx 0.070$, ${{a}_{0}}=10{{R}_{\odot }}$, ${{a}_{1}}=5.75\times {{10}^{6}}{{R}_{\odot }}=0.13{\rm pc},$ and $m\approx 1.2$ (Han et al. 1995). This leads to about 50% of the stars in binaries having orbital periods of less than 100 yr. Although the typical binary fraction of LMC clusters is about 35% (Elson et al. 1998; Bastian & de Mink 2009), the values for various clusters may be different. Thus some population models with different binary fractions are built in this work. We remove some random binaries to make the binary fraction equal to what we need.

After the generation of the star sample, we calculate the evolution of stars using the rapid stellar evolution code of Hurley & Tout (1998) and Hurley et al. (2002) (hereafter the Hurley code), which can evolve both single stars and binaries. In binary evolution, most binary interactions such as mass transfer, mass accretion, common-envelope evolution, collisions, supernova kicks, angular momentum loss mechanism, and tidal interactions are taken into account. The typical uncertainty of this code is about 5% in stellar luminosity, radius, and core mass, which affects the results slightly.

Because the Hurley code does not take stellar rotation into account, we add the effect of rotation on effective temperature, luminosity, and MS lifetime to massive (>1.7 ${{M}_{\odot }}$) stars, using a recent result of evolution of rotating stars (Georgy et al. 2013). The database of Georgy et al. 2013 is particularly useful for constructing synthetic populations of stars, accounting for mass, rotation, and metallicity distributions. It includes an accurate computation of the angular-momentum and stellar-wind anisotropy. The changes in surface temperature and luminosity, which are caused by stellar rotation, are calculated by comparing the evolutionary tracks of rotating stars to those of non-rotating ones. The changes can be described as a function of stellar metallicity, mass, and initial rotation rate. Besides changing the evolutionary track, rotation significantly affects the MS lifetime of stars. Because of centrifugal forces, rotating stars behave like non-rotating stars of lower mass (see also Meynet & Maeder 2000; Ekström et al. 2008). We therefore calculate the MS lifetime change of rotating stars using the result of Georgy et al. (2013). Stars with rotation rate (the ratio of the actual angular velocity to the critical one, ω; see Georgy et al. 2013) larger than 0.3 usually have a lifetime change greater than 15%, and the most rapid rotators may reach values of 30%–35%. Because of the lack of evolutionary tracks of stars with Z = 0.008, stellar metallicity is interpolated to the typical value of LMC clusters to fit the needs of this work. In addition, some random values are assigned to the stars of a population following the observed results of Royer et al. (2007; similarly, Zorec & Royer 2012) because the rotation rate of stars has obvious distributions. Finally, the stellar evolutionary parameters ([Fe/H], T_eff, ${\rm log} g$, ${\rm log} L$) are transformed into colors and magnitudes using the atmosphere library of Lejeune et al. (1998). The photometric errors derived from ASTs are applied to the simulated stars when comparing the synthetic CMDs with the observed one.

In this work we want to determine the distance modulus, color excess, binary fraction, rotational star fraction, and star formation mode of clusters. In order to determine these effectively, we apply a new CMD fitting method. The new technique searches for best-fit parameters by comparing the star fraction in every part (or grid) of a CMD. An observed CMD is divided into many grids via selected color and magnitude intervals. The goodness of fit is judged by the average difference of grids that contain observed or theoretical stars (hereafter effective grids), when taking the weight of each grid into account. For convenience, we defined a parameter, weighted average difference (WAD), to denote the goodness of fit. It can be calculated by

$$\begin{eqnarray}&&{\rm WAD}=\frac{{\Sigma }[{{\omega }_{i}}.|{{f}_{{\rm ob}}}-{{f}_{{\rm th}}}|]}{\sum {{\omega }_{i}}}\;,\end{eqnarray} \tag{ 2 }$$

where ${{\omega }_{i}}$ is the weight of the ith grid, and f_ob and f_th are the star fractions of observed and theoretical CMDs in the same grid. ${{\omega }_{i}}$ is calculated as

$$\begin{eqnarray}&&{{\omega }_{i}}=\frac{1}{|1-{{C}_{i}}|}=\frac{1}{{{\sigma }_{i}}}\;,\end{eqnarray} \tag{ 3 }$$

where C_i is the completeness of the ith grid that is given by AST, and ${{\sigma }_{i}}$ means star fraction uncertainty.

Although various indicators are used for the goodness of CMD fitting by different works, there are obvious advantages for using WAD in this work, which counts stars but is different from the approach used by Bertelli et al. (2003). First, WAD can give us a rough estimation of the difference between observed and synthetic CMDs. In other words, we believe that WAD has more physical meaning than other approaches, such as the widely used ${{\chi }^{2}}$, which is suitable for Gaussian observational error distributions, and Poisson equivalent ${{\chi }^{2}}$ (hereafter $\chi _{e}^{2}$), which is suitable for Poisson error distributions (Dolphin 2002). second, WAD has considered most information about the observed CMD, including distributions of magnitude errors, completeness of stars, and star fraction difference of observed and theoretical CMDs. Third, unlike the $\chi _{e}^{2}$ method of Dolphin (2002), in which the calculation of goodness parameters requires observed stars in every grid (in the calculation of $\chi _{e}^{2}$, i.e., Equation (5) in Dolphin (2002), ${\rm log} ({{f}_{{\rm th}}}/{{f}_{{\rm ob}}})$ is meaningless when f_ob = 0), our method, allows grids to have an absence of observed stars when calculating the goodness indicator. In order to help readers to understand the relation between WAD and other indicators, Figure 5 compares WAD with two other indicators (${{\chi }^{2}}$ and $\chi _{e}^{2}$). Note that the data are taken from the CMD fitting of NGC 1651 using a few different models. We see that WAD increases with ${{\chi }^{2}}$ as a whole, which suggests that they may give similar results. However, WAD and $\chi _{e}^{2}$ will possibly give different results because there is no correlation between them. Note that ${{\chi }^{2}}$ is calculated by

$$\begin{eqnarray}&&{{\chi }^{2}}={\Sigma }\frac{{{({{f}_{{\rm ob}}}-{{f}_{{\rm th}}})}^{2}}}{{{(1-C)}^{2}}},\end{eqnarray} \tag{ 4 }$$

and $\chi _{e}^{2}$ by

$$\begin{eqnarray}&&\chi _{e}^{2}=2{\Sigma }[({{f}_{{\rm ob}}}-{{f}_{{\rm th}}})+{{f}_{{\rm th}}}.{\rm log} ({{f}_{{\rm th}}}/{{f}_{{\rm ob}}})].\end{eqnarray} \tag{ 5 }$$

Zoom In Zoom Out Reset image size

When we test the above three kinds of fitting methods using a few artificial clusters (Z = 0.008; age = 0.9, 1.2, 1.5, 1.8, 2.1 Gyrs; ${{f}_{b}}$ = 0.0, 0.3, 0.6), it is shown that the WAD and ${{\chi }^{2}}$ methods usually give correct results, but the $\chi _{e}^{2}$ method does not report correct results for some cases. In particular, WAD fitting can determine the stellar ages of artificial clusters accurately, although there is an uncertainty of about 0.1 in ${{f}_{b}}$. When we apply WAD fitting to two simple clusters (NGC 1261 and NGC 2257), which do not have an eMSTO structure, the best-fit models agree well with observed CMDs as a whole (Figure 8). We see that a few stars, in particular those on the horizontal branch (HB), are not fitted well. This may result from the still large uncertainties in the late phases of stellar evolution. Indeed, HBs are very tricky to model because of mass loss. The test therefore suggests that the WAD method is reliable for CMD studies. Note that the data of two clusters are also taken from HST, and one can check their magnitude versus error relations via Figures 6 and 7. In fact, magnitude errors affect the shapes of CMDs visibly.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this section we use eight kinds of stellar population models to study the observed CMD of NGC 1651 and compare the results. These models are different from each other by binary fraction, fraction of rotational stars, star formation types (single or multiple bursts), and star formation histories. The eight types of models and their abbreviations are listed in Table 1. They are used to derive the best-fit model based on WAD fitting. Note that the observed CMD possibly contains a few field stars, e.g., those below the red giant branch, but they did not affect the best-fit result as they enlarge the WAD values of all models simultaneously.

In CMD fitting, the magnitude uncertainties of a simulated star are taken from those of AST stars with similar (difference less than 0.05 mag) magnitudes. The distance modulus, color excess, binary fraction, rotational star fraction, and star formation mode of a cluster are determined simultaneously via the WAD fitting method because it was shown to be reliable. A new CMD fitting code, i.e., "Binary Star to Fit for CMD" (BS2fit for CMD), which was developed by Dr. Zhongmu Li, is used for the study. The code is part of research project BS2fit, and it will be available to the public in the future. We use five star formation modes to describe the SFH of clusters in a simple way. These are the uniform, linearly increasing, linearly decreasing, Gaussian, and decreasing–increasing ("V" mode) formation modes. In each mode, the fraction of stars of each single burst is given as a function of stellar age.

The ranges and steps of parameters for the final run of the NGC 1651 CMD fitting are shown in Table 2. The values are estimated by some runs that took larger ranges and steps for these parameters in advance.

Table 2.  Parameter Ranges and Steps for the Last Run of the NGC 1651 CMD Fitting 

Parameter	Range	Step	Unit
Z	0.008	...	...
${{(m-M)}_{0}}$	19.10–19.20	0.01	mag
$E(V-I)$	0.00–0.25	0.01	mag
${{N}_{{\rm sf}}}$	1–5	1	...
${{M}_{{\rm sf}}}$	1–5	1	...
Age	1.0–2.5	0.1	Gyr
${{f}_{{\rm bin}}}$	0.0–0.8	0.1	...
${{f}_{{\rm rot}}}$	0.0 –100.0	0.1	...

Note. ${{N}_{{\rm sf}}}$, ${{M}_{{\rm sf}}}$, ${{f}_{{\rm bin}}}$, and ${{f}_{{\rm rot}}}$ are the number of star bursts, mode of star formation, binary fraction, and rotational star fraction, respectively. ${{M}_{{\rm sf}}}$ from 1 to 5 means uniform, linearly increasing, linearly decreasing, gaussian, and decreasing–increasing modes, respectively.

Download table as:  ASCIITypeset image

Table 3 shows the best-fit results from various models. As we see, when a relationship of ${{A}_{v}}=2.253*E(V-I)$ is adopted, the best-fit true distance modulus of cluster NGC 1651 is possibly between 18.68 and 18.70, which is consistent with previous results (18.5–18.7, see Mould et al. 1986; Sarajedini et al. 2002). The SSP-fit age (1.5–1.6 Gyr) of this cluster is different from Sarajedini et al. 2002 (∼1.8 Gyr), but is similar to that of Mould et al. 1997 and Mould et al. 1986 (∼1.6 Gyr). The total color excess $E(V-I)$ is about 0.19 mag. In addition, this cluster possibly contains about 50%–60% binaries. Note that "binaries" here mean the ones with orbital period less than 100 yr today. Because interacting binaries are only a part of such binaries, the binary fraction (${{f}_{{\rm bin}}}$) shown by this work is different from some other works (e.g., Davis et al. 2008), in which the binary fraction usually means the fraction of interacting binaries.

Table 3.  Best-fit Parameters of NGC 1651 from Various Models

Parameter	Model 1	Model 2	Model 3	Model 4	Model 5	Model 6	Model 7	Model 8
Z	0.008	0.008	0.008	0.008	0.008	0.008	0.008	0.008
${{(m-M)}_{0}}$	18.72	18.70	18.70	18.68	18.71	18.70	18.69	18.70
$E(V-I)$	0.19	0.19	0.19	0.18	0.19	0.18	0.18	0.19
${{N}_{{\rm sf}}}$	1	1	6	6	1	1	6	5
${{M}_{{\rm sf}}}$	⋯	⋯	5	4	⋯	⋯	5	3
Age (Gyr)	1.5	1.6	1.4–1.9	1.4–1.9	1.5	1.5	1.4–1.9	1.4–1.8
${{f}_{{\rm bin}}}$	0	0.5	0	0.5	0	0.6	0	0.5
${{f}_{{\rm rot}}}$	0	0	0	0	1.0	1.0	0.3	0.7
WADx1000	3.7381	2.4689	3.0112	2.1692	3.3273	2.2674	3.2933	2.0486

Note. The meanings of the models correspond to Table 1 via model numbers, and the units of parameters are the same as in Table 2. The value of ${{(m-M)}_{0}}$ has been corrected for galactic extinction. WAD is the difference of the star fraction in an effective grid of 0.047 mag for color and 0.020 mag for magnitude (the total effective grid number is about 2230). The less the WAD, the better the CMD fit.

Download table as:  ASCIITypeset image

Panel $(a)$ of Figure 9 shows the best-fit CMD from a traditional kind of model, i.e., non-rot single star simple stellar population (ssSSP). Complementary information is also presented in Figure 10. As we see, there are obvious differences between the best-fit (red) and observed (black) CMDs. The theoretical model does not reproduce well the extended MS, turn-off, and red clump. These three parts are labeled as parts 1, 2, and 3, respectively, to improve the readability of the figures. Note that the RC of a non-rot ssSSP is like a simple curve, as shown by red points.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Then, panel (b) shows the result from another kind of simple stellar population model, in which the effects of binaries have been taken into account. We find that although broad MS and blue stragglers are reproduced, there is significant difference between the observed and best-fit CMDs, in particular for the turn-off ("1") and red clump ("3") parts. This indicates that binaries are not the main cause of eMSTO (see also the review of Girardi et al. 2013) and eRC, although this contributes a few (about six) stars to eMSTO. It suggests that the stellar population of NGC 1651 is not an SSP without rotational stars. Our result is therefore different from that of Yang et al. 2011, which takes binaries as a main reason for eMSTO because an overly large star number (50,000 binaries) and a fixed binary fraction was used. In their work, the number of binaries are considerably overestimated.

After that, panel (c) shows the result from non-rot single star composite stellar population (ssCSP) fits. In the best-fit model, single stars formed by a few bursts, in fact all stars, have the same metallicity. It is shown that the best-fit model reproduces eMSTO and eRC, but it does not reproduce well the broad MS ("2") and blue stragglers (upper left). This implies that multiple star bursts may be the reason for eMSTO and eRC.

Next, panel (d) shows the result of a similar model, but taking binaries into account. We see that the best-fit population (non-rot bsCSP) fits the observed CMD much better than previously tested models. The main observational features, including eMSTO, broad MS, eRC, and blue stragglers, are reproduced here. The WAD of this model is much less than that of non-rot ssSSP, non-rot bsSSP, and non-rot ssCSP models (Table 3). Thus, among all non-rotational models, a bsCSP (parameters can be found in Table 3) is the best explanation for the observed CMD of NGC 1651.

If we use another kind of simple stellar population model, which takes the effects of stellar rotation into account, some new results are shown. Panels (e)–(g) give the comparison of best-fit and observed CMDs.

Panel (e) shows the result from std-rot ssSSP fits. In the best-fit model, all stars are found to be rotators. Such models can reproduce the eMSTO part ("1"), but cannot reproduce broad MS ("2"), eRC, and blue stragglers, as we can see. According to this figure, there is an obvious difference between fitted and observed CMDs. The stellar population of NGC 1651 is not likely to be an std-rot ssSSP.

When we use rotational binary population models (std-rot bsSSP) for CMD fitting, the result seems obviously better. Although the best-fit population is a simple one that formed its stars via a starburst, all observed features are reproduced by it, and its WAD is small, at only 0.00227. We see that stellar rotation can lead to eMSTO and eRC because it changes the luminosity, surface temperature, and evolutionary speed of stars. One can read panel (f) for a detailed comparison of best-fit and observed CMDs. Comparing the results of non-rot bsCSP and std-rot bsSSP fits (panels (d) and (f)) shows that both simple and CSP models can explain the CMD of NGC 1651. Therefore, it is not clear whether the eMSTO and eRC result from an age spread (i.e., multiple bursts), as the combination of stellar rotation and binarity can also lead to the eMSTO and eRC.

Because it is possible that a cluster contains a fraction of rotational stars with different ages, we also test such models. Panel (g) compares the best-fit std-rot ssCSP model and the observation. We find that broad MS ("2") and blue stragglers are not well fitted, while we find that std-rot bsCSP models fit best to the observed CMD (panel (h)) among eight kinds of models. This suggests that the stellar population of NGC 1651 may contain 50% binaries and 70% rotators. The ages of stars can vary from 1.4 to 1.8 Gyr.

In order to make the comparison between fitted and observed CMDs clearer, Figure 11 shows the difference between the observed CMD and best-fit models (those shown in Figure 9). By comparing two panels at the same line (e.g., panels (a) and (b)), we can conclude from the figure that stellar population models with binaries fit to the observed CMD better than those without binaries. Meanwhile, the comparisons of panels (b) and (f) and panels (d) and (h) indicate that models with rotational stars reproduce the observed CMD features better than those without rotators. This suggests that both binaries and rotational stars possibly contribute to the special CMD of NGC 1651. Because the fit of the eMSTO part from the std-rot bsSSP (panel (f) is better than from non-rot bsCSP (panel (d)), while it is not as good as that from the std-rot bsCSP (panel (h)), the eMSTO of NGC 1651 could be caused by not only age spread, but also stellar rotation and binarity. It is also possible that the observed CMD results from all three reasons. There is a degeneracy between their effects on the CMD. In addition, we find from panels (f) and (h) of Figure 10 that there are only a few grids with differences larger than 0.4 $f{{0}_{{\rm diff}}}$. This implies that most parts of the observed CMD can be reproduced well by both the std-rot bsSSP and std-rot bsCSP models.

Zoom In Zoom Out Reset image size

Finally, we check the color and magnitude distributions of the observed and fitted CMDs in Figures 11 and 12. It can be seen from Figure 11 that for the part near eMSTO ($V-I$ = 0.4 ∼ 0.6 mag), CSP models (non-rot bsCSP and std-rot bsCSP) have color distributions closer to the observed result, compared to simple stellar population models (non-rot bsSSP and std-rot bsSSP). This is in agreement with the result of Figure 9. At the same time, the distribution of magnitude has a similar trend for all models, as shown in Figure 12 (V = 20.5 ∼ 22.0 mag). We can thereby conclude that composite population models can better reproduce the eMSTO structure.

Zoom In Zoom Out Reset image size

From the comparison in Figures 9–12, we find that all models have slight differences from the observed CMD of NGC 1651 because all WAD values are greater than zero and there is an obvious difference in the color and magnitude distributions between the observed and fitted CMDs. A main part of this difference may result from the uncertainties in the IMF and the distribution of the initial binary separation of the stellar populations. Meanwhile, the WAD of std-rot bsSSP is close to that of non-rot bsCSP. This implies that the stellar population of NGC 1651 is equally possible to be a simple stellar population of rotating stars and a composite population without rotators. When comparing the WAD values of best-fit std-rot bsSSP and std-rot bsCSP models, we find a small difference of 0.0002188. This is similar to the uncertainty in CMD synthesis, which can be caused by the uncertainties in rotation rate distribution and in the treatment of the effect of stellar rotation on the magnitudes of stars. If we change the distribution of the stellar rotation rate, the WAD of std-rot bsSSP can possibly be smaller. Therefore, it is difficult to reject the model of bsSSP with rotating stars from our results. In this case, to study whether clusters like NGC 1651 contain a large number of rotational stars is crucial for understanding the stellar populations of such clusters.

In this work, we used eight kinds of stellar population models to fit the special CMD of LMC cluster NGC 1651. Our results suggest that stellar binarity, rotation, and age spread can demonstrably improve the goodness of CMD fit, but there is a degeneracy among their effects. It is therefore difficult to judge the stellar population feature (e.g., age and star formation history) of this cluster using traditional stellar population (TSP) models such as non-rotational ssSSP and ssCSP models. For future studies we strongly suggest using those advanced stellar population (ASP) models that take into account stellar binarity, rotation, and star formation history.

The result of this study shows that cluster NGC 1651 has a true distance modulus between 18.68 and 18.70 mag, $E(V-I)$ around 0.19 mag, and binary fraction around 50%. We can conclude that both bsSSP and bsCSP models can reproduce the main shape of the CMD of NGC 1651, if the effect of stellar rotation is taken into account. A bsCSP model with 50% binaries and 70% rotational stars can fit the observed CMD best among all the test models. The star formation history types (one or a few star formations), stellar ages, and star formation histories (or modes) from different kinds of models are significantly different. Although our result prefers that NGC 1651 be a bsCSP with rotational stars, it cannot exclude the bsCSP model without rotators or bsSSP model that consists of rotational stars. Further study (e.g., spectral study) is therefore needed to confirm the results. For instance, we can check how many and how fast the massive stars of a cluster are really rotating from the stellar spectra to judge the role of stellar rotation. Note that TSP models such as ssSP models seem significantly worse than bsSPs. Thus NGC 1651 certainly contains some binaries, and we need only to check how many rotators are included in a cluster like NGC 1651.

In this work, we use the calculation of Georgy et al. (2013) to estimate the effects of stellar rotation. Although the value is not accurate enough, the treatment is reasonable. First, in the work of Georgy et al. (2013), the variation in the angular momentum content is precisely tracked as it changes under the influence of stellar winds and mechanical mass loss. The effects of initial rotation on the Hertzsprung–Russell diagram, the evolution of surface rotation and abundances, and main sequence lifetime are computed simultaneously in their paper. This makes their results more reliable than some other works (e.g., Bastian & de Mink 2009). Second, the observed CMD of NGC 1651 is well reproduced when we add the result of Georgy et al. (2013) to our bsSPs models. Another important reason is that different works show that the trend of CMD change is correct when taking stellar rotation into account (Bastian & de Mink 2009; Li et al. 2012b; Yang et al. 2013), although observational work by Platais et al. (2012) and theoretical work by Girardi et al. (2011) showed that rotation works in the opposite sense in a star cluster. When we use the fitting formulae of Bastian & de Mink (2009) to calculate rotational effects, we also get the same conclusion. Third, as pointed out above, the effects of the fraction of rotational stars, distribution of rotation rate, and fraction of binaries are degenerated in some cases. Finally, although only the effects of massive rotational stars have been taken into account by rotational population models because of the lack of evolutionary data of low-mass rotators, the treatment does not affect the final result too much, as low-mass stars usually rotate much more slowly than massive ones. Therefore, the treatment for the rotational effects of stars is reasonable at this moment, and the main conclusion is accordingly reliable.

Furthermore, although Girardi et al. (2013) argue that a spread in the rotation rate cannot lead to eMSTO, because rotation affects the lifetime of stars, our work implies that the conclusion may strongly depend on stellar evolutionary models. It shows that when taking the new result of Georgy et al. (2013) for rotational effects, the MS lifetime change not only reproduces our previous result from fitting formulae (Li et al. 2012b), but also makes the fit of std-rot bsCSP to eRC better. Compared to our previous work, the treatment in this work can reproduce the observed CMD better because most rotational stars leave the MS later than non-rotating ones. This changes the positions of the components of binaries in a CMD when they leave the MS and reach the RC. Because the evolutionary track and MS lifetime of a rotating binary depend on the masses, separation, and rotation rates of its component stars, rotation cannot lead to the same CMD position (see Girardi et al. 2011 for comparison) for a given age. The special CMD structure of a rotating-star population results from both the changes in the evolutionary tracks and MS lifetime. Therefore, the MS lifetime change also contributes to the final CMD, and rotation plays different roles in populations with and without binaries. This is different from Yang et al. (2013), who claim that the increase in MS lifetime does not compensate for the changes in the evolutionary tracks.

Finally, although the AST technique was used for estimating the uncertainties in magnitudes of NGC 1651, more uncertainties exist and may affect the results. In fact, certain other factors, including uncertainties in metallicity, IMF, binary assumptions, calculation of stellar evolution, atmosphere library, and contamination of field stars, can also affect the CMDs of stellar populations. Because the study of these uncertainties is not the purpose of this work, we did not test their effects. However, further detailed studies on these uncertainties will be very helpful. One can read many examples of this kind of detailed study, such as Dolphin (2012, 2013).

The authors thank the referee for constructive comments, Prof. Yan Li for suggesting the use of evolutionary tracks of rotational stars, Mr. Chengyuan Li for help with the use of HSTphot, and Mr. Thomas Dodd for English proofreading. This work has been supported by the Chinese National Science Foundation (Grant No. 11203005), the Open Project of Key Laboratory for the Structure and Evolution of Celestial Objects, the Chinese Academy of Sciences (OP201304), and the Foundation of Yunnan Education Department (No. 2010Z004).