Discovery of Extended Main-sequence Turnoffs in Four Young Massive Clusters in the Magellanic Clouds

The origin of the multiple stellar populations that have been discovered in almost all Galactic globular clusters (GCs) remains an open question (Gratton et al. 2004; Piotto et al. 2015). Since it is not possible to observe the young GCs in the early universe, studying young GC candidates in the Milky Way and its satellites therefore plays a fundamental role in addressing this issue. However, almost all young Galactic clusters are located in the Galactic disk. They are affected by severe foreground extinction, which may mask the multiple stellar populations in their observed color–magnitude diagrams (CMDs).

An alternative way is to study the young clusters in the Large and Small Magellanic Clouds (LMC and SMC). Current observations of star clusters in the LMC and SMC mainly reveal two features: (a) almost all intermediate-age (1–2 Gyr old) star clusters exhibit extended main-sequence turnoffs (eMSTOs; e.g., Bertelli et al. 2003; Mackey & Broby Nielsen 2007; Milone et al. 2009; Girardi et al. 2013; Li et al. 2014a), while some even exhibit extended or dual red clumps (Girardi et al. 2009; Li et al. 2016a); (b) some young massive clusters (YMCs; with ages $\leqslant 400$ Myr) in the LMC exhibit bimodal or extended main sequences (MSs) (Milone et al. 2015, 2016, 2017; Bastian et al. 2016).

The discovery of eMSTO regions in intermediate-age clusters has triggered a surge of interest in explaining this feature, starting from the age-spread scenario, which attributes the observed eMSTOs to primordial age spreads of up to 700 Myr (Goudfrooij et al. 2011, 2014; Girardi et al. 2013; Girardi 2016). Next came the scenario of rapid stellar rotation (Bastian & de Mink 2009; Yang et al. 2013; Li et al. 2014b, 2016a; Brandt & Huang 2015; Bastian & Niederhofer 2015; Wu et al. 2016). Recently, a new scenario has been suggested, i.e., that stellar variability may play a potential role in shaping the eMSTO regions as well (Salinas et al. 2016; de Grijs 2017). The debate is indeed heating up.

Further excitement was added by the discovery of bimodal or extended MSTOs in YMCs. These features are usually explained by a coeval stellar population with a dispersion in stellar rotation rates (Li et al. 2014a, 2014b; Bastian et al. 2016) or a coeval, rapidly rotating population combined with a slowly/non-rotating population exhibiting an age spread (Milone et al. 2017). It has been suggested that the blue MS stars may hide a binary component, where binary synchronization is responsible for their small rotational rates (D'Antona et al. 2015). A recent study based on the YMCs NGC 1866 and NGC 1850 showed that a coeval stellar population featuring a distribution of stellar rotation rates can only partially explain the broadening of the MSTO (Correnti et al. 2017; Milone et al. 2017).

In this paper, we study four YMCs in the LMC and SMC, i.e., NGC 330 (SMC), NGC 1805, NGC 1818, and NGC 2164. Their isochronal ages do not exceed 40 Myr, 40 Myr, 35 Myr, and 100 Myr, respectively. Most are younger than the recently studied YMCs (e.g., Milone et al. 2015, 2016, 2017; Li et al. 2017). We analyzed their CMDs using high-resolution Hubble Space Telescope (HST) photometry, obtained using ultraviolet and visual filters. We find that they all exhibit extended MSTOs, although at different levels of significance. We test whether the observations can be reproduced assuming an age spread, a dispersion in stellar rotation rates, or a combination of these scenarios.

This article is arranged as follows. Section 2 presents our data reduction, which is followed by the main results of our analysis in Section 3. In Section 4 we present a discussion. Section 5 summarizes and concludes the paper.

The data used in this work come from the HST Ultraviolet and Visual Channel of the Wide Field Camera 3 (WFC3/UVIS) and Wide Field and Planetary Camera 2 (WFPC2) images. The WFC3/UVIS images were observed through the F225W and F336W passbands and the HST/WFPC2 images were observed in the F555W and F814W bands. Relevant information pertaining to the data is summarized in Table 1.

As regards the WFC3/UVIS and WFPC2 images, we performed point-spread function (PSF) photometry on the "_flt" or "_c0m" frames using the WFC3 and WFPC2 module of the dolphot 2.0 photometry package (Dolphin 2011a, 2011b, 2013). To obtain a stellar catalog with high-quality photometry, we adopted a filter employing the sharpness and "crowding" parameters calculated by dolphot. The sharpness is a measure comparing the profile of an object relative to the PSF. A perfect star would have a sharpness of zero. A negative sharpness that is too small usually indicates a cosmic ray, while a very large positive sharpness means that the detected object is extended, probably a background galaxy. The crowding parameter quantifies how much brighter a star would have been measured had nearby stars not been fitted simultaneously (in units of magnitudes). For an isolated star, the crowding is zero. High crowding usually means that the star is poorly measured. We confirmed that most high-crowding objects are faint, which would not dramatically affect our analysis. We only selected objects with −0.2 ≤ sharpness ≤ 0.2 and crowding ≤ 0.5 in both frames, which left us with ∼70% of the objects with high-accuracy photometry. The numbers of objects before the selection (and after, in brackets) for the clusters NGC 330, NGC 1805, NGC 1818, and NGC 2164 were 9996 (7208), 3808 (3102), 4478 (3346), and 4284 (3322), respectively. dolphot can automatically flag "good stars" and centrally saturated objects. We kept only objects that were indicated as good in our final sample for further analysis.

For each cluster, we combined stellar catalogs of different exposure times into a deep catalog. If one star appeared in two stellar catalogs, we selected the longest exposure time as the best representation. Next, we obtained two stellar catalogs, derived from the WFC3/UVIS and WFPC2 images. We transferred each star's CCD coordinates to equatorial coordinates ($X,Y\to {\alpha }_{{\rm{J}}2000}$, ${\delta }_{{\rm{J}}2000}$). We selected only stars located in the areas covered by both CCD fields for each cluster. This process sacrifices a large field but provides us with deep, multiband stellar catalogs.

We used the method introduced by Milone et al. (2012, their Section 3.1) to calculate differential reddening maps. We corrected all stars for differential reddening by assuming ${A}_{\mathrm{UVIS},{\rm{F}}336{\rm{W}}}=1.658{A}_{{\rm{V}}}$ and ${A}_{\mathrm{WFPC}2,{\rm{F}}555{\rm{W}}}=1.017{A}_{{\rm{V}}}$. We found that differential reddening does not dramatically change the morphology of our CMDs.

We generated figures showing number-density contours to determine the center coordinates of our clusters. We simply assigned the coordinate where the number density reaches the highest value as the cluster center. We then used the center of the number density to define a circular region within twice the half-light radius (2${r}_{\mathrm{hl}}$) as the cluster region. The ${r}_{\mathrm{hl}}$ values were taken from McLaughlin & van der Marel (2005) and were based on fits to a Wilson model (Wilson 1975).⁶ In Figure 1 we show the stellar spatial distributions, as well as their corresponding cluster centers and number-density contours. We confirmed that for all clusters, the stellar number density decreased to the field level at $2{r}_{\mathrm{hl}}$. The clusters' structural parameters are included in Table 2.

Zoom In Zoom Out Reset image size

Table 2.  Structural Parameters

Cluster	${\alpha }_{{\rm{J}}2000}$	${\delta }_{{\rm{J}}2000}$	${h}_{{\rm{c}}}$ (pc)	${h}_{\mathrm{hl}}$ (pc)
			(1)	(1)
NGC 330	${00}^{{\rm{h}}}{56}^{{\rm{m}}}31\buildrel{\rm{s}}\over{.} 00$	$-72^\circ 27^{\prime} 55\buildrel{\prime\prime}\over{.} 80$	${2.47}_{-0.19}^{+0.14}$	${6.95}_{-0.43}^{+0.96}$
NGC 1805	${05}^{{\rm{h}}}{02}^{{\rm{m}}}22\buildrel{\rm{s}}\over{.} 14$	$-66^\circ 06^{\prime} 41\buildrel{\prime\prime}\over{.} 40$	${1.34}_{-0.10}^{+0.09}$	${2.86}_{-0.08}^{+0.13}$
NGC 1818	${05}^{{\rm{h}}}{04}^{{\rm{m}}}13\buildrel{\rm{s}}\over{.} 44$	$-66^\circ 26^{\prime} 02\buildrel{\prime\prime}\over{.} 40$	${2.36}_{-0.09}^{+0.08}$	5.05±0.08
NGC 2164	${05}^{{\rm{h}}}{58}^{{\rm{m}}}55\buildrel{\rm{s}}\over{.} 80$	$-68^\circ 30^{\prime} 59\buildrel{\prime\prime}\over{.} 40$	${1.69}_{-0.05}^{+0.09}$	${4.37}_{-0.09}^{+0.05}$

Note. (1) McLaughlin & van der Marel (2005).

Download table as:  ASCIITypeset image

The final step is field-star decontamination. For NGC 330, NGC 1805, NGC 1818, and NGC 2164, we adopt stars that are located beyond 3.75${r}_{\mathrm{hl}}$, $6.00{r}_{\mathrm{hl}}$, $4.00{r}_{\mathrm{hl}}$, and $4.00{r}_{\mathrm{hl}}$, respectively, as our reference field stars. This choice results in numbers of field stars close to ∼500 for all clusters. We decontaminated the cluster CMDs using a similar method to that employed by Li et al. (2013a). We divided both the cluster and field CMDs into a carefully considered number of cells and counted the number of stars in each. The adopted numbers of grid cells were used to determine a reasonable average size for the grid cells. They should contain enough field stars to statistically correct the cluster CMDs, but they should not be too large, since that would prevent us from distinguishing detailed features along the MSs. Finally we randomly removed stars corresponding to those in the area-corrected field-star CMDs from the cluster CMDs. We found that our method tends to oversubtract field stars in the faint tail of the MS. This is because the field regions usually have lower detection limits than the clusters' central regions because of crowding and background effects (which results in higher completeness levels for the field samples at the bottom of the MSs; see Figure 7 of Li et al. 2013a). In this paper, we focus only on the MS range covering F555W ≤ 21 mag, which is ∼3 mag brighter than the magnitude of the detection limits. Our field-star contamination is thus reliable for the bright MSs.

The decontaminated YMC CMDs are presented in Figure 2. Here we present only the (F555W versus F336W – F555W) CMDs, because the broadening of the MS is most obvious in this parameter space. Nevertheless, this feature can be detected in other ultraviolet–visual CMDs as well, e.g., in the (F814W versus F336W – F814W) CMDs. As shown in Figure 2, all clusters exhibit wide MSs. We show the distribution of reference field stars in the same CMD. Contamination by field stars is minimal because the clusters' broadened MS regions are bright.

Zoom In Zoom Out Reset image size

We next used different isochrones to fit the observations. Specifically, we used two different stellar evolution models to compare with the observations, i.e., the Geneva SYCLIST code (Ekström et al. 2012; Georgy et al. 2013, 2014)⁷ and the MESA Isochrone & Stellar Tracks (MIST) models (Paxton et al. 2011, 2013, 2015; Choi et al. 2016; Dotter 2016).⁸ We first used the Geneva SYCLIST database to test whether a coeval population of stars characterized by different rotation rates could reproduce the observations. We overplotted nine Geneva isochrones with rotation rates of ω = 0.0–0.95⁹ to the observed CMD, where ω represent the ratio of the stellar rotation rate to the critical, breakup rate. All isochrones have the same age. The data of the isochrones were transformed into observational planes using the model atmospheres by Castelli & Kurucz (2004), convolved with the HST filter transmission. The effects of extinction are not included in the Geneva models; we simply shifted the isochrones to match the MS ridge-lines. The input metallicities in the Geneva models are limited to Z = 0.002, 0.006, and 0.014 (Z_⊙). The physical parameters adopted for our fits are included in Table 3 (rows 1–4). Our fits were based on initial estimates taken from previous publications. The typical age of NGC 330 is 20 Myr or $\mathrm{log}(t\,{\mathrm{yr}}^{-1})=7.30$ (McLaughlin & van der Marel 2005). For NGC 1805 and NGC 1818, the relevant numbers are 25 Myr to 40 Myr, i.e., $\mathrm{log}(t\,{\mathrm{yr}}^{-1})=7.40$–7.60 (Johnson et al. 2001), and for NGC 2164 it is 80 Myr, $\mathrm{log}(t\,{\mathrm{yr}}^{-1})=7.91$ (Mucciarelli et al. 2006). The best-fitting metallicities are [Fe/H] = −0.82 dex (0.15Z_⊙; McLaughlin & van der Marel 2005), −0.3 dex ($0.5{Z}_{\odot };$ Li et al. 2013a), 0.0 dex (Z_⊙; Li et al. 2013a), and −0.4–0.0 dex (0.4Z_⊙ to Z_⊙; Sagar & Richtler 1991) for NGC 330, NGC 1805, NGC 1818, and NGC 2164, respectively. The extinction, $E(B-V)$, for these clusters varies from 0.08 mag to 0.12 mag (Bessell 1991; Vallenari et al. 1991; Li et al. 2013a). The canonical distance moduli to the LMC and SMC are ${(m-M)}_{0}=18.49$ mag and 18.96 mag, respectively (de Grijs et al. 2014; de Grijs & Bono 2015).

Table 3.  Physical Parameters for the Isochrones or Synthetic Clusters in Figures 3 and 4 (Rows 1–4), Figure 5 (Rows 5–8), and Figures 6–9 (Rows 9–12)

Cluster	$\mathrm{log}(t\,{{\rm{yr}}}^{-1})$	Z	$E(B-V)$ (mag)	${(m-M)}_{0}$ (mag)	ω
NGC 330	7.30	0.002	N/A	19.00	0.00–0.95
NGC 1805	7.50	0.006	N/A	18.50	0.00–0.95
NGC 1818	7.55	0.014	N/A	18.55	0.00–0.95
NGC 2164	8.00	0.014	N/A	18.55	0.00–0.95
NGC 330	6.00, 7.40, 7.60	0.003	0.10	18.90	0.40
NGC 1805	6.00, 7.30, 7.60	0.006	0.12	18.40	0.40
NGC 1818	6.00, 7.25, 7.55	0.011	0.08	18.60	0.40
NGC 2164	7.70, 7.85, 8.00	0.011	0.10	18.55	0.40
NGC 330	6.00–7.60	0.002	N/A	19.00	0.00–0.95
NGC 1805	6.00–7.60	0.006	N/A	18.50	0.00–0.95
NGC 1818	6.00–7.55	0.014	N/A	18.55	0.00–0.95
NGC 2164	7.70–8.00	0.014	N/A	18.55	0.00–0.95

Download table as:  ASCIITypeset image

We show our fits in Figure 3. We find that for three young star clusters (NGC 330, NGC 1805, and NGC 1818), the non- and fast-rotating isochrones converge in the MSTO region. Although isochrones with different rotation rates create a broadened region near the MSTO region of NGC 2164, the broadening still seems too narrow to fully explain its extended MSTO region. For all our clusters, the observed MS widths gradually increase toward the MSTO region. Our fits seem to show minor disagreements between the models and the observations.

Zoom In Zoom Out Reset image size

Isochrone fitting alone cannot precisely describe the similarities between the models and the observations. In order to better demonstrate the reliability of our fits, we constructed a synthetic CMD for each cluster for comparison. Our synthetic CMDs are based on the adopted isochrones, to which we added the same photometric uncertainties as pertaining to the real data. We used the method of Milone et al. (2012). We determined MS–MS binary fractions of (53.8 ± 4.9)%, (46.1 ± 3.5)%, (39.9 ± 2.9)%, and (53.3 ± 4.5)% for NGC 330, NGC 1805, NGC 1818, and NGC 2164, respectively.¹⁰ The average binary fraction of these four clusters is ∼48.3%. We thus include 50% of unresolved MS–MS binaries in each simulated CMD. If the mass of the secondary component of a simulated binary system is below 1.7 M_⊙ (corresponding to the low-mass limit of the grid of B-type stars with different rotation rates in the Geneva models), we interpolate the magnitudes of the secondary star using the isochrones calculated based on the large grid (M = 0.8–500 M_⊙, only two rotation rates; Ekström et al. 2012; Georgy et al. 2013). Where necessary, we extrapolate outside the grid's boundaries down to 0.08 M_⊙ (i.e., the stellar hydrogen-burning limit).

For each CMD, we generated more than 3 × 10⁶ artificial stars. The Geneva models can only generate rapidly rotating isochrones for stellar masses down to 1.7 M_⊙. The numbers of the corresponding stars in the CMDs of NGC 330, NGC 1805, NGC 1818, and NGC 2164 are 1827, 493, 1206, and 1234, respectively. The number of artificial stars is thus at least 1600 times larger than the numbers of stars in the observations. For each observed star, we selected the 10 nearest artificial stars from the synthetic CMD. Finally, we randomly selected one of these 10 fake stars as representative of the observed star. Using this procedure, if a simple stellar population (SSP) characterized by different stellar rotation rates and a realistic fraction of unresolved binaries could fully cover and match the observed MS region, then the synthetic CMD should be almost identical to the observed CMD. We present our results in Figure 4.

Zoom In Zoom Out Reset image size

Figure 4 shows that the use of only a coeval stellar population characterized by different stellar rotation rates cannot fully reproduce the observed wide MSs. There is an excess of blue MS stars that cannot be reproduced by our simulations. Our method shows that many simulated stars are non-rotating (i.e., ω = 0.0). This is not surprising, because most fast-rotating stars would appear redder than their non-rotating counterparts. Since a fraction of the observed MS stars is too blue to be reproduced by our SSP models, our method is forced to select numerous non-rotating stars to try to cover those blue MS stars in the CMD. This confirms our speculation based on Figure 3: all sample clusters exhibit populations of very blue stars that cannot be explained by a coeval stellar population, not even when adopting a range in stellar rotation rates. The synthetic CMDs have the same distributions of the photometric uncertainties as the observations, and the resulting synthetic CMDs are based on large samples, each containing more than 3 × 10⁶ artificial stars, thus these additional blue MS stars cannot be explained by large scatter in the measurements. As shown in Figure 2, only a few background stars could contaminate the bright part of the MS; residual background contamination thus cannot explain these blue MS stars either. The additional blue stars may suggest the presence of young stellar populations.

We subsequently explored the age-spread scenario. We used the MIST models to describe the observations. The MIST models cover a large grid of single-star stellar evolution models, extending across all evolutionary phases for different stellar masses and metallicities (Paxton et al. 2011, 2013, 2015; Choi et al. 2016; Dotter 2016). The MIST 1.0 version includes a default rotational rate of $V/{V}_{\mathrm{crit}}=0.4\,\,(\omega \sim 0.6)$ in the output isochrones. Huang et al. (2010) studied the rotational velocities of ∼530 B-type stars. They found that the highest probability density occurs around $V/{V}_{\mathrm{crit}}=0.49$, which is close to the MIST default value. Using the Geneva models, we confirmed that the color difference between isochrones of $V/{V}_{\mathrm{crit}}=0.4$ and 0.49 is small (less than 0.02 mag in U − V color). Given that the broadened sections of our clusters' MSs are expected to be mainly populated by B-type stars, an average rotational rate of $V/{V}_{\mathrm{crit}}=0.4\,\,(\omega \sim 0.6)$ is thus a good approximation.

In Figure 5 we show the outcome of our fits including age spreads. For each cluster, we have overplotted a young, a median, and an old isochrone onto the observations. The best-fitting young and old isochrones were determined by visual inspection, by comparing their loci to the blue and red edges of the MSs. The median isochrones were determined arbitrarily; they will be used to indicate the trends of the turnoff (TO) regions for different ages. The metallicities, extinction values, and distance moduli for these differently aged isochrones are the same: see Table 3 (rows 5–8).

Zoom In Zoom Out Reset image size

As shown in Figure 5, as regards our fits to the colors, the age-spread scenario works better than the stellar rotation scenario. We found that in order to fit the additional blue MS stars, the age range of the stellar population needs to extend to almost zero for NGC 330, NGC 1805, and NGC 1818. We simply adopted an age of 1 Myr to represent the zero-age population in our fits (but see below). For the relatively old cluster NGC 2164, an isochrone with age down to 50 Myr is required to fit the blue MS boundary. Our fits indicate that the age spreads in NGC 330, NGC 1805, NGC 1818, and NGC 2164 are 40 Myr, 40 Myr, 35 Myr, and 50 Myr, respectively.

However, the age-spread model is not perfect either. The MIST models suggest that if the age spread were to reflect an extended star formation history, a 1 Myr old stellar population should contain a fraction of massive O-type and pre-MS stars. However, there is no obvious evidence of the presence of such objects. In NGC 330, a small number of stars are detected around the pre-MS locus, but we confirmed that they are background residuals, based on inspection of their spatial distribution. The absence of massive O-type and pre-MS stars indicates that the origin of such young stellar populations must be reconsidered. We will return to this issue in Section 4.

Based on their exploration of the TO region in the cluster NGC 1850, Correnti et al. (2017) suggested that the best solution may be a combination of stellar rotation and an age spread. We similarly explored the promise of the combination of an age spread and stellar rotation. The input parameters are presented in Table 3 (rows 9–12). The method we used is similar to that used for Figure 4. This time, we generated artificial stars that are different in both age and rotation properties. For each observed star, we randomly selected an artificial star from among the 10 nearest candidates as its best representation. Our simulation was based on the Geneva database, because the Geneva models provide isochrones for different stellar rotation rates. Each time when we generated an SSP for a different stellar rotation rate, we compared its CMD with the observational counterpart, similarly to our analysis of Figure 4. Again, we confirmed that the reproduction is not satisfactory even if we vary the isochronal age. Specifically, once we adopted a young isochronal age for the synthetic CMD, a fraction of red MS stars could not be reproduced. This is a similar conclusion to that reached by Milone et al. (2016, their Figure 10). Synthetic CMDs of stellar populations with different ages and rotation rates are presented in the right panels of Figures 6–9. For a more direct comparison, we also included the observed CMDs (gray dots). In the synthetic CMDs, colors represent the ages of the simulated stars.

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 first examined NGC 330. We found that if we use a combined model including an age spread and a distribution of stellar rotation rates, only a small number of blue MS stars would be poorly reproduced. We find that for the low-mass end of the MS, the contributions from young and old populations are not very different. However, for the MSTO region, these contributions are indeed very different. Almost all young stars reproduce the observed stars located in the blue part of the MSTO region, while the old stars mainly contribute to the red part of the MSTO. This result further supports the notion that a fraction of young stars is required in order to reproduce the bluest MS stars.

The same result for NGC 1805 is shown in Figure 7. The synthetic CMD of multiply aged stellar populations is almost identical to the observations. Only three MS stars around F555W ∼ 18 mag are not well reproduced. This is not surprising, because the synthetic CMD fully covers the region of the observed MS; for each observed star we should be able to find a corresponding artificial star with similar color and magnitude. This result stands in sharp contrast to the synthetic CMD of a coeval stellar population (see Figure 4), where numerous additional blue MS stars in the range 16 ≤ F555W ≤ 19 mag are not well reproduced.

This result also holds for NGC 1818, as shown in the right panel of Figure 8: almost all blue MS stars that appear in Figure 4 are reproduced. Finally we show the same result for NGC 2164 in Figure 9. Again, the simulated MS is significantly broadened compared with that of a single-aged stellar population. In summary, the cluster CMDs can be reproduced well by a combination of SSPs that cover an age range of 35–50 Myr and a wide variety of rotation rates.

At face value, the successful reproduction of the observations through a combination of an age range and a dispersion in rotation rates seems to indicate the need for an extended star formation history for the clusters. However, the adoption of using a "nearest" star to represent an observed star is complicated. This treatment is similar to the application of a mass truncation to the young stellar population, which thus avoids the appearance of massive O-type stars. This is indeed a contrived approach.

As we already emphasized, the origin of any ongoing star formation in these clusters is suspicious. In Figures 6–9, we also included the young isochrones adopted in Figure 5, as well as the equal-mass binary locus of the corresponding old isochrone. Although these two lines seem to adequately describe the boundaries of the MSs, the young isochrone clearly predicts the presence of bright O-type and pre-MS stars. An absence of O-type stars would indicate an unphysical mass function, while the lack of pre-MS stars reflects the notion that star formation has long been terminated.

4.1. Stellar Rotation?

4.2. Extended Star Formation?

As we showed in Figures 6–9, an age spread of 35 Myr (NGC 1818) to 50 Myr (NGC 2164) is required to reproduce the observed broad MSTOs. An age spread of ∼10 Myr (${\rm{\Delta }}\mathrm{log}(t\,{{\rm{yr}}}^{-1})\,\sim$ 0.1–0.2 for the typical ages of NGC 330, NGC 1805, and NGC 1818, roughly equal to the timescale of initial gas expulsion; Krause et al. 2016) is insufficient. As Figures 6–9 clearly illustrate, only a stellar population with an age as young as 1 Myr could reach the positions of the bluest MS stars. Our result is in agreement with the conclusions of Milone et al. (2017) and Correnti et al. (2017).

Does this imply that all of our sample clusters have extended star formation histories? One should be cautious as regards such speculations. We emphasize once again that the absence of O-type and pre-MS stars contradicts the ongoing star formation hypothesis. It seems that the most massive stars in the blue stellar population are not more massive than those of the bulk stellar population. As we showed in Section 3, to reproduce the observed MSs, we have to assume a mass truncation for the young stellar population. This seems a contrived approach to reproduce the observations, which cannot be naturally explained by a invoking a scenario of extended star formation.

A possible explanation is that a few million-year-old stars may still be embedded in their natal dust cocoons, which would prevent us from observing them in the UV band. However, this cannot explain the presence of the very blue MS stars. Why would the very bright O-type stars be obscured while the B-type stars (and thus those blue MS stars) are discernible in the UV band? In addition, Bastian & Strader (2014) studied the gas and dust contents in NGC 330, NGC 1818, and NGC 2164; no significant gas residuals nor any dust were detected in these clusters. All these arguments thus challenge the continuous star formation hypothesis.

We can estimate the minimum masses for our clusters to retain their gas using the equation from Georgiev et al. (2009),

$$\begin{eqnarray}&&{M}_{\mathrm{cl}}\approx 100{v}_{\mathrm{esc}}^{2}{r}_{{\rm{h}}},\end{eqnarray} \tag{ 1 }$$

where ${M}_{\mathrm{cl}}$ is the total mass of the cluster expressed in units of M_⊙, ${v}_{\mathrm{esc}}$ is the escape velocity of the initial gas in km s⁻¹, and ${r}_{{\rm{h}}}$ is the cluster's half-mass radius in parsecs. We assumed that the clusters' half-mass radii are equal to their half-light radii, ${r}_{\mathrm{hl}}$.¹¹ If star formation in these clusters can last for several tens of millions of years, numerous Type II supernova explosions should have taken place. The escaping gas can be accelerated to several hundred km s⁻¹. The corresponding minimum masses required to retain the initial gas in our sample clusters are log$({M}_{\mathrm{cl}}/{M}_{\odot })=6.84$, 6.46, 6.70, and 6.64 (assuming a minimum escape velocity of 100 km s⁻¹). However, the current masses of NGC 330, NGC 1805, NGC 1818, and NGC 2164 are only log$({M}_{\mathrm{cl}}/{M}_{\odot })=4.61$, 3.70, 4.41, and 4.18, respectively (McLaughlin & van der Marel 2005). Clearly, it is difficult for these clusters to sustain long star formation episodes by accreting their initial gas.

Another explanation to explain the occurrence of an extended star formation episode that seems viable is that after the initial gas expulsion phase (∼10 Myr), these clusters may still have been sufficiently massive to accrete the subsequent stellar ejecta of the stars on the asymptotic giant branch (AGB). The minimum velocity for the AGB stellar ejecta is about ${v}_{\mathrm{esc}}=10$ km s⁻¹ (Renzini 2008; Li et al. 2016b). The minimum masses to retain the gas of the AGB ejecta for these clusters are $\mathrm{log}(M/{M}_{\odot })=4.84$, 4.46, 4.70, and 4.64, which is still 1.7 (NGC 330) to 5.8 (NGC 1805) times the current cluster masses. Although the clusters would have lost their stellar mass through dynamical evaporation, this would not dramatically affect our clusters, because the typical timescale for such mass loss is expressed in units of billions of years (McLaughlin & Fall 2008; Li et al. 2016b).

In summary, although the observed broad MSs in our clusters indicate the possible presence of an age spread, it seems that such age spreads are unlikely to originate from extended star formation histories.

4.3. Blue Straggler Stars?

An alternative explanation is that the population of puzzling blue MS stars are blue straggler stars (BSSs). Because BSSs are produced through binary mass transfer, mergers, or stellar collisions, their maximum mass does not exceed twice the mass of TO stars. This may explain the absence of O-type and pre-MS stars.

D'Antona et al. (2015) studied the split MS in the cluster NGC 1856. They speculated that all observed blue MS-component stars may hide a binary component. They suggested that the periods pertaining to these binary systems range from 4 to 500 days. Based on Kepler's Third Law, the relationship between the binary period, P, and its semimajor axis, a, is

$$\begin{eqnarray}&&P=2\pi \sqrt{\displaystyle \frac{{a}^{3}}{G(M+m)}}=2\pi \sqrt{\displaystyle \frac{{a}^{3}}{{GM}(1+q)}},\end{eqnarray} \tag{ 2 }$$

where G is the usual gravitational constant, M and m are the masses of the primary and secondary stars, and q = m/M is their ratio. Rearranging Equation (2), we can calculate a in terms of P, M, and q:

$$\begin{eqnarray}&&a=\sqrt[3]{\displaystyle \frac{{P}^{2}{GM}(1+q)}{4{\pi }^{2}}}.\end{eqnarray} \tag{ 3 }$$

Assuming that q ranges from 0 to 1, P from 4 to 500 days, and that the typical mass of the young-population stars in our clusters is $2\,{M}_{\odot }\leqslant M\leqslant 9\,{M}_{\odot }$, the resulting distribution of the semimajor axes ranges from a = 0.06 to 3.23 au. Using a Monte Carlo approach, we derived that the average length of the semimajor axis of the binary population is $a=1.63\,\mathrm{au}$, with 23% of the binaries having semimajor axes $a\leqslant 1\,\mathrm{au}$ (we adopted flat distributions for q, P, and M). This indicates the presence of a significant fraction of compact binaries, which may form a potential population for mass-transfer BSS candidates. This was also suggested by Yang et al. (2011). Specifically, the latter authors suggested that the eMSTO regions in intermediate-age star clusters include significant contributions from binary interactions and mergers.

However, the fraction of BSSs with respect to the bulk stellar population is usually very small (e.g., for the old LMC GCs, see Mackey et al. 2006; Li et al. 2013b). It is unclear how a cluster environment could produce such a large number of BSSs. To our knowledge, there are no direct observations of BSSs in YMCs. Xin et al. (2007) studied a large sample of ($\gt 1$ Gyr old) Galactic open clusters. They found that the specific frequency of BSSs (i.e., the ratio of the number of BSSs to that of stars spanning a 2 mag range below the MSTO) ranges from 1% to 20%. Numerical simulations also show that the specific frequency of BSSs should be small during the early stages of a cluster's evolution (e.g., at an age of 50 Myr, the number of BSSs is fewer than 50 for a cluster that initially contains 100,000 stars; Lu et al. 2011; Hypki & Giersz 2013).

Despite this challenge, if the observed young stars are indeed BSSs, most should have been formed through binary mass transfer. Future studies should focus on their carbon and oxygen abundances to test this suggestion (Ferraro et al. 2006).

In this paper, we have studied the CMDs of the clusters NGC 330, NGC 1805, NGC 1818, and NGC 2164. They all exhibit broad MSs, which cannot be explained by an SSP with unresolved binaries and photometric uncertainties. We suggest that it is likely that most YMCs may exhibit wide MSs in their ultraviolet–visual CMDs.

We found that the gravity darkening caused by stellar rotation plays a very limited role in hot, massive, early-type MS stars. In the meantime, rotational mixing would cause an isochrone to have a TO position that is almost indistinguishable from that of a non-rotating isochrone. Therefore, a dispersion of stellar rotation velocities in coeval ensembles of stars cannot reproduce the observed wide MSs.

The failure of the stellar rotation scenario to fully reproduce the observations indicates that the age-spread scenario may still be viable. An age spread of 35–50 Myr is required to explain the observations. We confirm that a combination of an age spread and stellar rotation can adequately reproduce the observations for all sample clusters. Similar conclusions were derived previously for the ∼100–200 Myr old clusters NGC 1850 and NGC 1866 (Correnti et al. 2017; Milone et al. 2017). However, we argue that the apparent age spread is unlikely, owing to continuous star formation. Indeed, the clusters' masses seem too small to sustain extended star formation episodes. Moreover, the absence of O-type and pre-MS stars also contradicts the hypothesis of extended star formation.

We suggest that the young stars may be BSSs, which reduces the need for very massive O-type and pre-MS stars. However, if this were correct, it is not clear why the number of BSSs in the clusters should be comparable to the number in the bulk stellar population. To understand the origin of these young stars, details about their chemical composition and rotation rates are required.

We thank Choi Jieun (Harvard–Smithsonian Center for Astrophysics) for providing us with the relevant model database. C.L. is supported by the Macquarie Research Fellowship Scheme. R.d.G. and L.D. acknowledge research support from the National Natural Science Foundation of China through grants U1631102, 11373010, and 11633005. A.P.M. acknowledges support from the Australian Research Council through a Discovery Early Career Researcher Award, number DE150101816.