Combined Effects of Rotation and Age Spreads on Extended Main-Sequence Turn Offs

The intermediate age cluster NGC 2203, located in the LMC is the oldest of our target clusters according to the literature, at about 1.55 Gyr (Goudfrooij et al. 2014; Rosenfield et al. 2017). Photometry for this cluster is the same as that used in Rosenfield et al. (2017), where the data reduction process is also described. In brief, this data is part of a larger set that was re-reduced from two HST programs: GO-9891 (PI: Gilmore) and GO-12257 (PI: Girardi). So, this photometry comes from archival Advanced Camera for Surveys (ACS) and WFC3 data, re-reduced with the University of Washington data reduction pipeline, designed to reduce the HST treasury programs ANGST (Dalcanton et al. 2009) and PHAT (Dalcanton et al. 2012); see Williams et al. (2014) for further details. ASteCA (Perren et al. 2015) was used to determine the cluster center, found via the maximum spatial density using a 2D Gaussian kernel density estimator. The cluster radius corresponds to where the radial density profile becomes indistinguishable from the background stellar density. ASteCA was also used to determine contamination, utilizing a nonparametric Bayesian decontamination algorithm based on the method of Cabrera-Cano & Alfaro (1990). The final membership was limited to stars within the cluster radius with >70% membership probability.

We adopt the cluster parameters cited in Goudfrooij et al. (2014), i.e., their age, [Fe/H], A_V, and distance estimated via best-fit isochrones from the Padova group (Marigo et al. 2008). Respectively, these are values of 1.55 Gyr, −0.30, 0.16 dex, and 18.37. We also adopt the binary fraction of 0.18 used by Goudfrooij et al. (2014). Taking into account that slight differences in fits to CMD features can arise due to different input physics between our MIST and those Padova models. Although, see that in Figure 1 that these parameters still provide a reasonable fit with our models.

On the boundary of the intermediate age regime, NGC 2249 is a star cluster also located in the LMC. Membership for this cluster was determined in an identical way to what was described for NGC 2203. Isochrone analysis in Correnti et al. (2014) has estimated the cluster to be about 1 Gyr old. Data reduction for this cluster is the same as described for NGC 2203, with further details in Rosenfield et al. (2017). In this work, we adopt the same mean age, [Fe/H], A_V, and distance modulus for NGC 2249 as cited by Correnti et al. (2014); namely, 1 Gyr, −0.46, 0.07 dex, and 18.2 mag, respectively. We adopt the binary fraction of 0.30 cited by Correnti et al. (2014).

NGC 2818 enters the regime of young cluster ages for our target clusters at about 700 Myr and it is located in the Milky Way. Our NGC 2818 data is taken from the publicly available data gathered by Bastian et al. (2018), originating from Gaia DR2 (Gaia Collaboration et al. 2016, 2018). Details of the member selection for this cluster are given in Bastian et al. (2018); briefly, members were selected via proper motion and parallax cuts. NGC 2818 is one example of a number of young Galactic clusters with eMSTOs that have been revealed with Gaia DR2 data in the last year or so, thanks to Gaia's enhanced photometric precision. The adopted cluster parameters are based on values used in Bastian et al. (2018): with an age of 700 Myr, a solar [Fe/H] of 0.0, A_V = 0.90, distance modulus of 12.76. We adopt a binary fraction of 0.29 from Cordoni et al. (2018).

NGC 1831 is also an approximately 700 Myr old cluster as well, located in the LMC. For this cluster, we carried out PSF photometry on the flat-field corrected, and bias-subtracted HST "flc" images (Program ID: GO-14688) using the WFC3 module of DOLPHOT, a modified version of HSTphot (Dolphin 2000) and following the procedure described in Balbinot et al. (2009). For our analysis, and for the CMD shown in Figure 1, we chose stars within a half-light radius (from McLaughlin & van der Marel 2005), r_h = 3385 of the cluster's center. We have adopted the parameters based on those used by Goudfrooij et al. (2018) for NGC 1831, determined from PARSEC isochrones. Thus, an age of 700 Myr, [Fe/H] of −0.25, A_V of 0.14, distance modulus of 18.35, and a binary fraction of 0.20.

NGC 1866 is the youngest cluster in this study, at about 200 Myr old, and it is located in the LMC. This photometry was obtained in the same manner as described for NGC 1831 (except from HST Program ID: GO-14204). Here we have taken stars within the half-light radius r_h = 429 of the cluster's center, based on McLaughlin & van der Marel (2005). Our adopted cluster parameters for NGC 1866 are adopted from Milone et al. (2017) for the metallicity and binary fraction. These parameters are an age of 200 Myr, [Fe/H] = −0.36, A_V = 0.34, a distance modulus of 18.31, and a binary fraction of 0.25.

The starting point for our models is the MESA stellar evolution code (Paxton et al. 2011, 2013, 2015, 2018, 2019), version r7503, which is a modular and open source 1D stellar evolution code. We closely followed the physics used for the MIST database (Choi et al. 2016), adopting the protosolar abundances of Asplund et al. (2009) and using boundary conditions from ATLAS12, while SYNTHE is used for bolometric corrections (Kurucz 1970, 1993). Our models are evolved to the end of core helium burning. Hereafter, we will refer to our models as "MIST models" or "MIST based."

The MIST models are set rotating at the zero-age MS (ZAMS) with a given velocity denoted by the ratio of equatorial angular velocity ${{\rm{\Omega }}}_{\mathrm{ZAMS}}$ at ZAMS, over the critical Ω_c. The critical velocity Ω_c is a property intrinsic to the star that depends on its mass (see, e.g., Maeder & Stahler 2009); it represents the limit where centrifugal force overcomes the star's gravity. The ratio ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}$ is equivalent to the linear velocity form $v/{v}_{{\rm{c}},\mathrm{ZAMS}}$ in the MESA formalism. Previously, our models were limited to ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}=0.6$ in Gossage et al. (2018); we now include models ranging up to ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}=0.9$, in steps of 0.1. The initial MIST models released by Choi et al. (2016) only included nonrotating and ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}=0.4$ models.

MIST models rotation under the shellular approximation developed by Kippenhahn & Thomas (1970), with chemical and angular momentum transport described by the equations of diffusion equations of Endal & Sofia (1978). This diffusive formalism is also adopted in the stellar evolution codes STERN (Brott et al. 2011) and the recent version of PARSEC (Costa et al. 2019b). The shellular approximation is standard in 1D stellar evolution codes. The treatment of angular momentum and chemical transport varies between codes. For instance, GENEC uses a diffusive-advective formalism, described in Zahn (1992), Maeder & Zahn (1998), and Maeder & Meynet (2000). The differences in these two formalisms have significant effects on the models, e.g., leading to different MS lifetime extensions (by up to 20% or so) and color–magnitude variations (see Choi et al. 2016; Gossage et al. 2018 for examples). Our models possess a stronger convective mixing with weaker rotation mixing, whereas GENEC features the opposite. Consequently, our rotating models are primarily affected by the structural changes of gravity darkening when they rotate, they do not see a dramatic MS lifetime extension or luminosity enhancement from rotational mixing, as is seen in the GENEC models.

Gravity darkening is handled by the equations of Espinosa Lara & Rieutord (2011) (recently adopted by Paxton et al. (2019) in MESA as well) in determining the (surface averaged) luminosity and temperature of a given stellar model at viewing angle i. The viewing angle corresponds to $i=90^\circ$ when viewing is equator-on, versus 0° when viewing pole-on. Gravity darkening is the effect of centrifugal force reducing the surface gravity of a rotating star. This effect is stronger at the equator than at the poles, due to the greater centrifugal force at the equator. Thus, gravity darkening causes the equator of a rotating star to become cooler and dimmer than the poles, introducing a viewing-angle dependence on the apparent magnitude and color of a rotating star. The effects can be substantial; examples for our models exist in Gossage et al. (2018). The formalism that we adopt from Espinosa Lara & Rieutord (2011) is similar to the gravity-darkening formalism used for the GENEC-based SYCLIST stellar population models.

Another important aspect of rotation, at least for masses $\lesssim 1.8\ {M}_{\odot }$ is magnetic braking. Modeling this process is an active area of research (e.g., Garraffo et al. 2016, 2018; Sadeghi Ardestani et al. 2017; Garraffo et al. 2018; Fuller et al. 2019). Niederhofer et al. (2015) predicted the Δ(Age)–age trend should stop after magnetic braking becomes effective. This limit is expected to be reached by TO stars in clusters with ages older than about 1.5 Gyr, depending on the metallicity. In their recent work, Georgy et al. (2019) used models developed with the STAREVOL code (Amard et al. 2016), for masses between 1 and 2 ${M}_{\odot }$, including a prescription for magnetic braking according to Matt et al. (2015). This mass range is roughly where surface convection zones develop, leading to surface magnetic fields that can act on extended stellar material, braking the star. We take a crude approach to simulate this, in absence of a proper model of the effects of magnetic braking. Below $M=1.3\ {M}_{\odot }$, models are forced to be nonrotating; from $M=1.3\,\mathrm{to}\,1.8\ {M}_{\odot }$ models have ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}$ scaled up to the full value. The SYCLIST models do not model magnetic braking, but also exclude stellar masses below $M=1.7\ {M}_{\odot }$ (Georgy et al. 2014). The MSTOs of NGC 2203 and 2249 have TO masses that are low enough for magnetic braking to become important, so our results for these clusters in particular will be affected by uncertainties due to magnetic braking.

From our MIST-based stellar models, we compute synthetic stellar populations, as in Gossage et al. (2018), using the code MATCH (Dolphin 2002). Specifically, we use MATCH to compute Hess diagrams of CMDs, including unresolved binaries, at fixed values of Z, age, and ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}$. Photometric errors are simulated with MATCH via artificial star tests. The populations are created at distinct ages, covering log(age) = 8.0 to 9.5 (in 0.02 dex steps), each of which is also created at ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}=0.0,0.1,\ldots$, up to 0.9. Our synthetic populations include the effect of gravity darkening via randomly drawn viewing angles for constituent stars. Stellar models are drawn according to a Kroupa (2001) initial mass function (IMF). We combine these synthetic populations (weighting them as described in Section 3.2) to form a composite stellar population that may possess stars from a range of ages and rotation rates.

Colored points in the top row of Figure 2 shows several examples of these models. In this figure, models have fixed age pertaining to the representative cluster, from left to right: NGC 2818, 2249, and 2203. The colors map the full range of rotation rates (i.e., ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}=0.0$–0.9), according to a flat distribution of ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}$, to clearly show the full effects of stellar rotation on the MSTO in these clusters. Red lines show nonrotating MIST-based isochrones, dashed at the mean age, and solid at ±Δ(age). This Δ(Age) was chosen so that the isochrones roughly covered the full extent of the eMSTO. The bottom row of panels shows the same nonrotating isochrones overlaid on the data of NGC 2818, 2249, and 2203. Broadly, both stellar rotation and a range of ages can cover the extent of the eMSTO, but in different ways, hence the contention between the theories.

Zoom In Zoom Out Reset image size

Figure 2 shows how increasing age reddens and decreases the luminosity of TO stars in a similar manner to the effects of gravity darkening. The effects of stellar rotation manifest in a confined region of CMD space at a given age, in comparison to an age spread that can modify the luminosity and temperatures of stars over as wide a space as is useful for covering the eMSTO. In this sense, the effects of an age spread have relatively more freedom than the effects of a rotation-rate distribution. MATCH was used to generate the models displayed in the top row, but the code currently lacks the ability to perform a straightforward fit for rotation-rate distributions. We opted for a separate fitting method outside of what is provided by MATCH (outlined below in Section 3.2) to handle this.

We build composite stellar populations according to three scenarios (below "weights" measure the total number of stars in a Hess diagram whose stars may be observations or simulations):

• 1.  
  Gaussian log(age) spread (σ_t model): This has a Gaussian SFH. More complex SFH are imaginable, but would expand parameter space further, and so are omitted for now. Here ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}$ is restricted to 0.0 (nonrotating). This model has three parameters describing the overall weight, or number of stars, in the composite ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}=0.0$ population (i.e., amplitude), plus the Gaussian log(age) distribution's mean and standard deviation.
• 2.  
  Nonparametric rotation-rate distribution (Ω model): This model considers 10 free weights ranging from 0 to the total weight of observed stars. Each free weight corresponds to one of the 10 possible ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}$ populations. All stars are assumed to have a single age, in this model, which is also fit as an 11th free parameter.
• 3.  
  Age spread with rotation (σ_tΩ model): This model combines the σ_t and Ω models. A Gaussian age distribution is allowed (mean age and standard deviation are free parameters), as are the 10 free weights for a nonparametric rotation distribution, giving a total of 12 parameters. Like the σ_t model, the age distribution here is in terms of log(age).

We measure the probability of models matching the data using Hess diagrams and a Poisson likelihood as described in Dolphin (2002). Our fitting considers up to 10 independent weights for the density of stars at our 10 values of ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}$, plus up to two more parameters describing the Gaussian age distribution's mean and standard deviation. We take this standard deviation of the Gaussian age distribution to represent the "age spread." So, given a derived mean log(age) μ_τ, we take the age spread to be ${10}^{({\mu }_{\tau }+{\sigma }_{\tau })}-{10}^{({\mu }_{\tau }-{\sigma }_{\tau })}$, with σ_τ being the standard deviation. We use Markov chain Monte Carlo (MCMC) to sample the probability distributions, determining the most likely parameter values for rotation-rate weights, mean ages, and age spreads. We employ emcee's affine-invariant ensemble sampling algorithm (Foreman-Mackey et al. 2013). To initialize our ensemble of walkers, we chose randomized positions from a uniform distribution within ±0.2 dex of the chosen mean log(age) for a cluster, and within ±0.05 dex of the arbitrarily chosen initial age spread of Δlog(age) = 0.05, if applicable. In initializing the walker positions for the various ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}$ weights was done using a Dirichlet distribution. The reasoning behind this choice being that we desired these random initial positions to lie within the solution plane, such that all rotation-rate weights add up to the total combined weight of all bins for the data. In other words, we set these positions in a way that preserves a cluster's total number of stars, rather than initializing in an invalid portion of parameter space where the total counts is far off from that of the data. Thus, our model for rotation distributions is a nonparametric model consisting of 10 free parameters; our eSF model is a Gaussian model described by its mean and standard deviation.

Here we present the results of mock tests examining accuracy in parameter recovery. These tests were carried out for our three scenarios of eMSTO presence under consideration: population age spread (σ_t model), stellar rotation distribution (Ω model), or both (σ_tΩ model). We generated mock data according to each of these scenarios using MATCH and applied our models to check that the input age and relative weights of populations at various ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}$ were recovered. The mock data is generated with a mean log(age) 9.0 in all cases, metallicity of [Fe/H] = −0.40, A_V = 0.07, akin to NGC 2249, as determined in Correnti et al. (2014). We generate mock data according to each scenario and the log(age) (measured in yr) spread is 0.05, while the input rotation distribution is a Gaussian centered at ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}=0.4$ with a standard deviation of 0.2 dex in applicable cases.

The recovered weights shown in Figure 3 find values near the inputs (black, hatched), and the truth is contained within errors in most cases, except the two shown in the bottom-left and -right panels (d) and (f). These two cases correspond to fitting the Ω model to mock data created with either the σ_t model in panel (d), or σ_tΩ model in panel (f). In both of these cases, the model forms roughly a bimodal distribution of rotation rates; both of these cases are the Ω model fit to mock data containing an age spread. Thus, we expect that our Ω model sees an age spread as a bimodal distribution of rotation rates. The reason for this is that fast rotators can enhance the eMSTO spread with more dramatic gravity-darkening effects, as shown in, e.g., Bastian & de Mink (2009), Brandt & Huang (2015), and Gossage et al. (2018). Additionally, our MIST-based rotating models mostly become redder as ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}$ rises, populating a reddened MSTO, allowing a greater TO spread but leaving behind a depopulated blue MSTO. Low rotation rates refill the bluer side of the MS when added to the full ensemble of stars. The bimodal nature of rotation seems to arise from fast rotators being used in an effort to fit the eMSTO spread, leaving a depopulated blue MSTO, alongside slow rotators counteracting this offset and populating the red MSTO. So, the properties of fast and slow rotators appear to drive our models to favor the presence of both in explaining an age spread, or otherwise broad MSTO. Unless the underlying rotation distribution in the data is something specific, like a Gaussian, we expect to see a bimodal distribution arise in order to model the width of the eMSTO. This is shown in panel (e) for Ω model fit to itself when the input rotation distribution is a Gaussian. The full model (i.e., the σ_tΩ model) is capable of recovering the input rotation-rate distribution in all cases, as may be seen in panels (a)–(c). The additional degrees of freedom allowed by the Gaussian age spread in the σ_tΩ model removes the necessity for a bimodal rotation-rate distribution.

Zoom In Zoom Out Reset image size

Figure 4 shows residuals for our mock tests, again omitting 1:1 comparisons and only showing cases where the comparisons produced interesting residuals. In panel (b) of Figure 4, one may see that the Ω model does not match the smooth variation of stellar densities created by a Gaussian age spread (as in the σ_t model). It does better matching the σ_tΩ mock, seen in panel (a) of Figure 4, (i.e., it acquires a higher likelihood), but similarly misses the continuous morphology of a Gaussian age spread; the placement of rotating models on the CMD imposes a combination of relatively specific morphologies, discretized according to a corresponding rotation rate. Thus, there is a low chance that each of these morphologies matches the more ambiguous morphology of a Gaussian age spread, creating distinct features in the residuals and relatively poor fits when the two try to match each other. Finally, the bottom row of Figure 4 shows the σ_t model fit to the σ_tΩ model in panel (c) and to the Ω model in panel (d). The σ_t model is generally capable of achieving higher likelihoods than the Ω model shown in the top row, but it still has difficulty reproducing the densities of stars at non-zero rotation rates, leading to the features shown in the residuals. This indicates that the σ_t model possesses enough ambiguity to smooth out inconsistencies and produce a higher fit statistic with disregard to the presence of a distribution of rotation rates. In comparing panels (a) and (c), it is noticeable that the sub-giant branch (SGB) of the σ_tΩ mock data is better matched by the Ω model.

Zoom In Zoom Out Reset image size

Derived Gaussian age distributions for the σ_t (blue-shaded region) and σ_tΩ models are shown in Figure 5, with the best-fit age of the Ω model as a solid cyan line. Black dashed lines show the literature ages described in Section 2. For NGC 2203 and 2249, black solid curves indicate the "pseudo-age" distributions determined by Goudfrooij et al. (2014) and Goudfrooij et al. (2017), respectively. It is also mentioned in those works that the pseudo-age distributions are broader than what photometric errors allow, and so are not spurious in that manner.

Zoom In Zoom Out Reset image size

Briefly, the pseudo-age distribution is one method of determining the age spread on the eMSTO. It is created with a parallelogram that encloses the width of the eMSTO. The colors and magnitudes of stars within this parallelogram are translated into an age distribution. This is done by taking the ages predicted by stellar models of these stars. Like our Gaussian age distributions shown in Figure 5, the pseudo-age is not reflective of the SF history of the cluster on its own (Goudfrooij et al. 2014), but rather emulates the distribution of ages that may be present in the eMSTO at the time of observation.

In comparison to our derived age distributions, the pseudo-age distributions have multiple peaks, with the strongest amplitude at younger ages than we find. This may be due to the different models used (Padova in Goudfrooij et al. 2014 and SYCLIST in Goudfrooij et al. 2017, while ours are MIST based). Additionally, it could be due to our age-spread model only allowing a single peak, causing it to compromise for a peak in between the multiple peaks found by Goudfrooij et al. (2014, 2017). In spite of these offsets in mean age, the widths of the pseudo-age and Gaussian age distributions are comparable, suggesting that both methods predict similar age spreads.

However, a common trend seen with all clusters is that the ${\sigma }_{t}{\rm{\Omega }}$ model predicts a smaller age spread. This is expected, as stellar rotation also contributes to the eMSTO morphology with this model. Gravity darkening can introduce substantial color variations and eMSTO width, as demonstrated in Figure 2. Furthermore, it may be seen that the inclusion of stellar rotation tends to reduce the predicted mean age.

Due to the reddening effect of gravity darkening, it is also expected that stellar rotation would reduce the mean age. Age is primarily determined by the CMD position of the MSTO, which is fixed by the data. As stellar rotation tends to redden stellar models, selecting a younger age counteracts this effect by making the rotating stars blue again. Hence, stellar rotation is seen to derive a younger age than nonrotating models (the σ_t model) in Figure 5. The ages found by the Ω and σ_tΩ models are either similar or coincide.

Although the σ_tΩ and Ω model agree on cluster ages, the derived ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}$ distributions shown in Figure 6 highlight where these scenarios disagree. The top row shows the distributions found by the σ_tΩ model, while results for the Ω model are on the bottom row. In all cases, the Ω model finds a more distinct population of fast (e.g., ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}\geqslant 0.5$) and slow rotators (${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}\lt 0.5$). In contrast, the ${\sigma }_{t}{\rm{\Omega }}$ model generally finds a smaller presence, or lack of slow rotators.

Zoom In Zoom Out Reset image size

The Ω model's rotation-rate distributions agree with observations more than the σ_tΩ model in this context. They qualitatively reproduce the observed fast and slow rotators found by, e.g., Dupree et al. (2017) and Bastian et al. (2018). In the observations, slow rotators reside blueward, while fast rotators lie redward. The stellar models capture this behavior as well, as seen in Figure 2; gravity darkening tends to redden fast rotators, causing the two populations to occupy distinct color spaces in the CMD.

The two populations are necessary in the stellar rotation scenario so that the full blue to redward extent of the eMSTO is reproduced.

The two populations are not required in the σ_tΩ model. The Gaussian age spread can compensate for a lack of slow rotators. This does lead to a clear lack of slow rotators in NGC 2203 and 2249 (the oldest clusters). For these clusters, the predicted lack of slow rotators is not in line with recent findings for younger clusters.

Figure 7 shows the residuals of the data compared to the best-fit models for the σ_tΩ (left column), Ω (middle), and σ_t model (right). Dashed black lines show nonrotating MIST isochrones placed to trace the eMSTO width (similar to the isochrones in Figure 2). Each row corresponds to a cluster, and the global likelihoods (-2lnP) are kept in the upper right corner of each panel. While the σ_tΩ model may not reproduce the rotation-rate distributions one might expect, it is the best-fit model overall, considering the global likelihoods.

Zoom In Zoom Out Reset image size

In NGC 2249, 2818, and 1831, the likelihoods of all three scenarios are comparable, and the best-matched clusters on the basis of the residuals. Meanwhile, NGC 2203 and 1866 (the oldest and youngest clusters) are the worst matched. In the latter two cases, performance between the three scenarios is more disparate. The σ_tΩ model formally achieves the best fit with all clusters (perhaps unsurprisingly given more free parameters), but clearly so with these two clusters in particular. None of the residuals are clean, showing that all scenarios produce imperfect fits, although in different ways.

Inclusion of an age spread does tend to deliver a higher likelihood. At the same time, the residuals show that the ${\sigma }_{t}$ and σ_tΩ models often create eMSTOs that are broader than the data. Pixels overfit by these models tend to lie outside of the bounding isochrones that roughly trace the eMSTO. The Ω model tends to overpredict within the isochrones, without the extended behavior of the Gaussian age-spread models.

The Ω model can find comparable likelihoods to the age-spread models but also finds mismatches. In fact, it is formally the worst fit in all cases. Yet, it does appear to match the eMSTO extent well qualitatively, suggesting that it fails more in getting the precise stellar densities correct. In other words, the Ω model appears to reproduce the data morphologically, though it does not show a strong statistical advantage over the other two models.

The Δ(Age)–age trend, such as the one derived here and shown in Figure 8, has been a central point in eMSTO studies such as Niederhofer et al. (2015). The points in Figure 8 mark age spreads and ages found by the σ_t (blue) and σ_tΩ models (the latter in red). For a series of ages, Niederhofer et al. (2015) determined the effective Δ(Age) that a nonrotating GENEC-based, SYCLIST isochrone would need in order to match a ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}=0.5$ isochrone. They selected distinct points on each isochrone, and determined how much the nonrotating isochrone's age needed to shift from its original value in order to match the CMD position of the rotating isohrone. Remarkably, they found a Δ(Age)–age trend that agreed with measured eMSTO widths from a range of studies. The lines in Figure 8 show such model predicted trends. These Δ(Age)–age trends imply that the effect of stellar rotation can mimic an age spread.

Zoom In Zoom Out Reset image size

The magenta line in Figure 8 is the trend formed by their second isochrone point, "M_V at MSTO." We performed the same analysis to see what the MIST models predict and this trend is shown as the solid black line. The MIST models produce a shallower trend. In Gossage et al. (2018), we found that MIST models predict a smaller apparent age spread as a result of rotation when compared to SYCLIST, due to our weaker rotational mixing. This causes the luminosities of our rotating and nonrotating isochrones to separate less than SYCLIST models, leading to a smaller effective "age spread." We also found that gravity darkening appeared to be the strongest effect of stellar rotation in the MIST models.

The black and magenta lines in Figure 8 do not show the inclination angle dependence of rotationally induced eMSTO width. These lines were measured with the luminosity and temperature of the stellar models averaged over all inclination angles. Thus, this trend does vary slightly depending on the chosen inclination angle, but it remains that on average the MIST models predict a shallower trend than SYCLIST when measuring the eMSTO spread with this isochrone-based method. In a synthetic population (as opposed to single isochrones), the combined effects of a distribution of rotation rates and inclination angles allow the rotating MIST models to achieve comparable eMSTO widths to the observed Δ(Age)–age trends.

In the derived ages and age spreads of the σ_t and σ_tΩ models, we see that the inclusion of stellar rotation does not reduce the age spread to zero. The age spreads derived by the σ_tΩ model are smaller than those found by the σ_t model. This is expected, as stellar rotation and age spread compete to explain the eMSTO morphology in the σ_tΩ model. The fact that the σ_tΩ model's age spreads still correlate with cluster age suggests that they trace a stellar evolution effect, rather than SF. The fact that the σ_tΩ model finds a reduced, but non-zero age spread may be indicative of missing physics in the rotating models, compensated for by an age spread.

NGC 1866 has a split MS, in addition to its eMSTO, which warrants further analysis. This additional complexity may contribute to the relatively poor fit achieved for NGC 1866 (Figure 7). The split MS has been argued to imply that star clusters host both eSF and a distribution of rotation rates. The fact that the σΩ model achieves the best fit here, and the Ω model the worst, serves to demonstrate how eSF may be compensating for an incomplete modeling of stellar evolutionary effects (e.g., binary mergers and decretion disks). In spite of the lower statistical likelihood of the Ω model (Figure 7), we find that a coeval distribution of rotation rates is capable of reproducing the split MS and the eMSTO simultaneously, albeit with caveats.

Schematically highlighted in Figure 9 is the debate between whether the split MS is due to a rotation-rate distribution or whether eSF may be present. Figure 9 displays content similar to Figure 10 from Milone et al. (2017). A range of ages clearly does not reproduce the split MS (right panel). As also noted by Milone et al. (2017, 2018), modeling the split MS appears to require fast and slow rotators in the cluster. This is shown for the simple case of a nonrotating and ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}=0.5$ isochrone (${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}=0.7$ is included as well to show how they aid in extending the MSTO). The split MS may be driven by bimodal rotation-rate distributions, but $V\sin i$ confirmations still need to be obtained. As plotted here, and as seen in Figure 6, the fast rotators are mainly stars rotating at ${\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{crit},\mathrm{ZAMS}}=0.5$, in contrast to the predictions of Geneva shown in Milone et al. (2017, 2018) where the red MS appears to be comprised of stars rotating near critical velocity.

Zoom In Zoom Out Reset image size

The right panel of Figure 9 shows that a range of stellar ages can aid in spanning the bluemost and redmost regions of the eMSTO. This was also shown in the corresponding Figures 10 and 11 of Milone et al. (2017), where isochrones at several rotation rates and multiple ages are plotted together. This is one example of how eSF in models (e.g., σΩ) may optimize a fit to the data, but it is uncertain whether the invoked age spread is physical. As mentioned, there is sparse empirical evidence for eSF, and alternatives are known (although, see recent work by Costa et al. (2019a), finding possible evidence for eSF in NGC 1866 via modeling of its Cepheid stars).

Our best-fit Ω model (which obtains the worst likelihood) for NGC 1866 is shown in Figure 10, with an inset focusing on the split MS. The recovered rotation-rate distribution for NGC 1866 is indeed bimodal (for both the σΩ and Ω model), and the model qualitatively matches the split MS. The eMSTO is matched in some respects, e.g., that it predicts a mixed population of slow and fast rotators, as observed spectroscopically by Dupree et al. (2017). In other areas, a rotation-rate distribution is inconsistent with the data as modeled here. Particularly, it misses stars blueward (i.e., about F336W–F814W < −0.65 mag) and redward (about F336W–F814W > −0.35) in the eMSTO, but these areas of the CMD are known to be affected by complex stellar evolution effects, e.g., interacting binaries and Be stars, respectively.

Zoom In Zoom Out Reset image size

For instance, our models do not include properties such as decretion disks nor binary merger products. The population of fast rotators extending blueward may not reproduce the findings of Dupree et al. (2017), and could be a spurious result. Increasing the binary fraction to 50% and removing the possible "blue stragglers" does not influence the presence of these near-critical rotators in the best fit; their presence appears to aid in matching stellar densities along the MS, rather than the eMSTO, with these near-critical rotators extending blueward of the TO as a possible side effect. On the matter of reddened, fast rotators, Bastian et al. (2017) found evidence of a high number of Hα emitters (suspected Be stars) throughout the eMSTOs of these young clusters, a number of which lie in this redward region. Though, Correnti et al. (2017) examined whether Hα emitters account for the redward extension of the MSTO in the young cluster NGC 1850; in comparing their data with SYCLIST stellar models, they found that Hα emitters did not in NGC 1850. See the modeling done by Granada et al. (2018) for more on Be stars and further examples in young clusters. As for the blueward extent, Li et al. (2019) show some examples of the blue straggler populations that might exist in Magellanic Cloud clusters, but Dalessandro et al. (2019) warn that such stars may also be unaccounted for field contaminants. Presently for our NGC 1866 data, these stars do not reside in the field, and so appear to be cluster members.

Rather than eSF, it is possible that the redmost and bluemost populations of eMSTO stars (where the Ω model does not reproduce the data well, thus calling to question if this is evidence of eSF) are affected by decretion disks and binary interaction. These features are also missed by the Ω model in NGC 1831, 2818, and 2249. The σΩ model may compensate for these missing factors by invoking an age spread. We will need to improve our stellar modeling, or obtain further data on what sort of stars these are before we can say for sure.