Age Determinations of the Hyades, Praesepe, and Pleiades via MESA Models with Rotation

Of our target clusters, the Hyades is nearest to our solar system, making it a popular object of study among astronomers for many years (e.g., Smart 1939; Wayman 1967; van Altena et al. 1997; Reino et al. 2018). We adopt a distance modulus μ = 3.349 mag (46.75 ± 0.46 pc) from Gaia Collaboration et al. (2017) and an extinction of A_V = 0.0031 (Taylor 2006). Historically, the age of the cluster has been determined to be around 600 Myr, according to nonrotating isochrone fits; for example, Perryman et al. (1998) derive an age of 625 ± 50 Myr from fitting optical CMD data. Martín et al. (2018) measured Li surface abundance in 6 brown dwarf candidates of the Hyades, and using the models of Baraffe et al. (2015) derived an age of 650 ± 70 Myr. The Hyades may have a [Fe/H] = 0.103 ± 0.008 according to Taylor & Joner (2005); recently, Cummings et al. (2017) found 0.146 ± 0.004 from the spectra of 37 Hyads.

Our optical data is comprised of the "high fidelity" stellar members identified by de Bruijne et al. (2001); these are stars with relatively high membership likelihood and evidence supporting their status as single stars. Binary systems are removed from this sample, and so our assumed binary fraction is zero. This catalog is based on the main Hipparcos catalog (Perryman et al. 1997), utilizing derivations of secular parallaxes for cluster candidates to determine membership likelihoods. Using these high fidelity members, we fit stars with V_T < 8 mag, forcing the fitting algorithm to focus on the age-sensitive MSTO.

Due to the close proximity of the Hyades, a significant spread in the apparent magnitudes will exist due to differential distance effects. This effect must be accounted for when fitting observed CMDs of the Hyades to the models. de Bruijne et al. (2001) manage this by considering the absolute magnitude of stars (in their case derived using e.g., the Hipparcos secular parallaxes from the main catalog, and recorded in their data tables). We do the same, but convert back to an apparent magnitude using the average distance modulus, μ, of the Hyades cluster, μ = 3.349 mag. In essence, we place all stars at the mean distance of the cluster.

Our infrared data is sourced by the members of Goldman et al. (2013), who used the members of Röser et al. (2011), where the former researchers made efforts to extend membership down to 0.1 M_⊙. In both instances, membership was determined by the convergent point method. Furthermore, field star contamination has been estimated for these datasets; Goldman et al. (2013) found that contamination is likely, although decreasing to negligible levels (less than 10%) inwards of 18 pc within the catalog. Binaries are present as well, and here we assume a binary fraction of 25%, as suggested by Gunn et al. (1988) (this value was also adopted by Röser et al. 2011).

The Praesepe is the furthest of our target clusters at roughly four times the distance to the Hyades, but also appears similar to the Hyades in both its age and chemical composition. We adopt μ = 6.26 mag (≈179 pc) (Gáspár et al. 2009) and A_V = 0.0837 (Taylor 2006) for modeling Praesepe. The metallicity of Praesepe has been estimated e.g., by Boesgaard et al. (2013), to be [Fe/H] = 0.12 ± 0.04 based on measurements of 11 Praesepe dwarfs; Cummings et al. (2017) found 0.156 ± 0.004 from the spectra of 39 Praesepe members.

Our optical data derives from the 24 members tabulated by Madsen et al. (2002). We have cross-matched these stars to obtain updated Tycho B_T, V_T magnitudes via their Hipparcos ID numbers. This has been checked and cleaned of binary and field stars by the catalog's authors. Here we also impose a V_T < 9 mag cut, for the same reasons that we made this cut in our data of the Hyades (see Section 2.1).

Additionally, we use the near-infrared (NIR) 2MASS J, K_s magnitudes cataloged by Wang et al. (2014). This dataset includes 1040 stellar members in total, ranging from M ≈ 0.11–2.4 M_⊙ acquired via proper motion analysis. A binary fraction of 20%–40% is suggested for these members; correspondingly, we adopt a binary fraction of 30%. Furthermore, as noted by Wang et al. (2014), Praesepe members exhibit distinct proper motions from potentially interfering field stars, virtually eliminating confusion between them. The authors estimate that their member list possesses ∼16% nonmembers (i.e., 862 true members and 168 nonmembers). Here, we also impose a cut to only include stars with V_T < 9; focusing on the MSTO, however, also lessening the probability of contaminating stars, as contamination is estimated to lessen toward the population's bright end.

The Pleiades is closer to us than the Praesepe, at roughly 3 times the distance to the Hyades, possessing an appreciably younger age, as well as a lower metallicity than the previous clusters. Our adopted distance modulus and extinction for this cluster are: 5.64 mag (Gaia Collaboration et al. 2016), or ≈134 pc, and A_V = 0.1054 (Taylor 2008), respectively. Soderblom et al. (2009) found [Fe/H] = 0.03 ± 0.05 from spectroscopic measurements of 20 Pleiads. The age has been estimated by Barrado y Navascués et al. (2004) to be 130 ±20 Myr using the LDB. From CMD analysis, Meynet et al. (1993) determined 100 Myr by fitting nonrotating isochrones to its higher mass stars. Recently, Dahm (2015) found an age of 112 ± 5 Myr via the LDB, using the evolutionary models of Baraffe et al. (1998, 2015).

Our optical data is again found using the Hipparcos IDs of cluster members reported in Madsen et al. (2002), subsequently cross-matched for B_T, V_T magnitudes. Binary and field stars have been removed from this sample by Lindegren et al. (2000) with a maximum likelihood estimation to determine likely (single) cluster members.

We also use the 2MASS magnitudes of members reported by Stauffer et al. (2007). In fitting this data, we only include stars with K_s < 5.0 mag, for the same reasons that we made similar cuts in the Praesepe and Hyades (see Section 2.1). There is no estimate of the binary fraction for this data but stars in this magnitude range are well known, higher mass members of the Pleiades. They are: Alcyone, Electra, Atlas, Maia, Merope, Taygeta, Pleione, and Celeno—many of these are Be stars. Of the stars listed, Atlas is a double line spectroscopic binary system (Zwahlen et al. 2004); evidence from lunar occultation (Richichi et al. 1994) suggests Taygeta may be a binary system; Pleione may be a single line spectroscopic binary (see e.g., Katahira et al. 1996 and Nemravová et al. 2010); Electra's multiplicity is inconclusive thus far, with contradicting evidence (e.g., Abt et al. 1965 and Pearce & Hill 1971); meanwhile Alcyone is a multiple star system, although most of its members are further than 77 arcsec (as listed in the Washington Double Star Catalog). In practice, we assume Atlas and Taygeta are systems whose photometry may be significantly affected by their companions and remove them. We make cuts in our optical data similar to the infrared data, only using stars with V_T < 5.0.

Our models are fundamentally based on the MESA stellar evolution code (Paxton et al. 2011, 2013, 2015, 2018). MESA is a 1D, open source, parallelized code, written with a modular design for flexibility in customizing its incorporated physics. The stellar structure and composition equations are simultaneously solved via a Newton–Raphson solver. Boundary conditions are necessary for solving stellar structure equations; MESA provides a number of options for their choice. In these models, the ATLAS12 code (Kurucz 1970, 1993) model atmosphere tables set these boundary conditions; as in Choi et al. (2016). Furthermore, we use the protosolar abundances of Asplund et al. (2009) where metallicities are scaled to $Z={Z}_{\odot ,\mathrm{protosolar}}=0.0142$. In Section 3.1.1, we give an overview of how rotation manipulates a star and Section 3.1.2 is dedicated to a description of gravity darkening and how it is implemented in our stellar models.

In this work, we have expanded upon the grids presented in Choi et al. (2016) by computing a set of models with fine variation in the initial rotation rates and a denser grid in metallicity. These are not fully evolved, but truncated to the end of core helium burning (CHeB), with a metallicity range more focused on the solar neighborhood and LMC. Our metallicity range has a higher resolution around values appropriate for our clusters of interest (near solar [Fe/H] =0.0). These choices were made to introduce greater variability in the rotation rate as a model parameter while maintaining a reasonable model creation timescale. The mass, [Fe/H], and rotation ranges covered by our models are listed in Table 2.

3.1.1. Stellar Rotation

In this section we provide a brief overview of the key physical processes related to modeling stellar rotation (see Choi et al. (2016) and the MESA papers (Paxton et al. 2011, 2013, 2015, 2018) for details). In MESA, surface magnetic effects are not modeled, thus magnetic braking is neglected presently. To compensate and reproduce the observed slow rotation of solar-like stars, models below 1.2 M_⊙ do not rotate. Between 1.2 and 1.8 M_⊙, rotation is scaled by a factor ranging from 0 to 1 as mass increases, such that models above 1.8 M_⊙ rotate at their full velocity. For reference: MSTO masses range from ∼1.4 to 2.5 M_⊙ at isochrone ages of 650–750 Myr (e.g., the Praesepe and Hyades), and from ∼2.4 to 4.5 M_⊙ at ages of 100–150 Myr (e.g., the Pleiades), in our modeling.

Rotational velocity is commonly characterized by the ratio ${v}_{\mathrm{ZAMS}}/{v}_{{\rm{c}}}={\rm{\Omega }}/{{\rm{\Omega }}}_{{\rm{c}}}$. This is the zero age main sequence (ZAMS) surface equatorial rotation rate, v_ZAMS, compared with the critical rotation velocity, v_c, the velocity that would disrupt the star through centrifugal force. The quantities Ω and Ω_c are the angular counterparts.

The 1D code, MESA, implements rotation through the shellular approximation described by Kippenhahn & Thomas (1970), and employs the diffusive approximation introduced by Endal & Sofia (1978) to model rotationally enhanced chemical mixing. The latter scheme requires that a choice be made for two parameters: f_c, which dictates how closely compositional mixing follows the flow of angular momentum transport, and f_μ, which encodes the efficiency of rotational mixing in the presence of stabilizing molecular weight gradients. Values of f_c = 1/30 (calibrated to reproduce the surface ⁷Li abundance of the Sun by Pinsonneault et al. 1989; Chaboyer & Zahn 1992) and ${f}_{\mu }=0.05$ (found by Heger et al. 2000 to reproduce nitrogen surface enhancement in evolved 10–20 M_⊙ stars from Gies & Lambert 1992; Herrero 1993; Vrancken et al. 2000) are adopted in our models.

Although the shellular approximation is virtually ubiquitous in 1D stellar evolution codes, the diffusive approximation is not. Other codes, such as Geneva (Ekström et al. 2012) take a diffusive-advective approach (as described in Maeder & Meynet 2000a) instead. Consequently, MESA models can exhibit noticeable distinctions from models developed via alternative codes, in some respects (as pointed out by Choi et al. 2016, primarily showcasing the differences of higher mass models). Rotating MIST models tend to be fainter and cooler at the MSTO compared to those same models computed by Geneva, perhaps due to less efficient chemical mixing from rotation effects.

In Figure 1 we reiterate these points, focusing on lower and intermediate masses, as would likely be found in our relatively young target clusters. In these plots we show evolutionary tracks for stellar models of 2, 3, 4, and 7 M_⊙. Our models are displayed as solid lines, while dashed lines correspond to the Geneva models. The rotation rate is encoded by color, wherein red represents Ω/Ω_c = 0.0 (nonrotating) and blue is Ω/Ω_c = 0.5. The two top panels overlay our models and Geneva models with no rotation (top left) and with rotation (top right); these top panels are comparisons across model sets. Our nonrotating models are brighter and hotter than the nonrotating Geneva models; however, the roles are reversed in comparing rotating models. The two bottom panels compare nonrotating to rotating models within each model set to see how rotation affects the models differently according to each respective code. Our models become cooler as the rotation rate increases (due to gravity darkening; Section 3.1.2), whereas Geneva models become hotter and brighter as their rotation rate is increased.

Zoom In Zoom Out Reset image size

Figure 2 shows these effects on a CMD, as they manifest in the MSTO of an isochrone constructed from our models. Figure 2 shows a Ω/Ω_c = 0.0 isochrone (black) compared to a Ω/Ω_c = 0.6 isochrone (red, dashed) at an age near to that which has been classically reported for the Pleiades, ∼100 Myr, (left panel) and an age near classical reports for the Hyades/Praesepe, ∼600 Myr (right panel). The effects are visually quite different between these ages. On the left, a slightly younger isochrone is shown in solid blue to point out how rotation can enhance the MSTO brightness, causing an older rotating population to mimic the morphology of a younger nonrotating one; here our rotating models might predict older ages compared to nonrotating models (albeit, the effect is very slight in our modeling). Whereas on the right, less massive and more slowly rotating stars exist on the MSTO at 600 Myr and the effects of rotation are modest. Thus, the MSTO barely becomes brighter at these older ages, and primarily appears cooler, mimicking the color position of an older population (shown as blue in the right-hand panel); here our models might find a younger age (decreasing age will make the MSTO hotter again) than predicted by nonrotating models, or could cause shifts in derived metallicity. Hence, our models do not behave as straightforwardly as Geneva models, where rotation primarily makes the MSTO brighter and hotter. The behavior exhibited by Geneva models leads to rotating models predicting younger ages than nonrotating models.

Zoom In Zoom Out Reset image size

At 100 Myr, MSTO stellar masses are ∼2.4–4.5 M_⊙, which rotate fully. Meanwhile at 600 Myr, the extant MSTO mass range is roughly ∼1.4–2.5 M_⊙. As discussed at the start of this section, stars below 1.8 M_⊙ do not rotate at full capacity in order to replicate the effects of magnetic braking relevant to lower mass stars with convective envelopes. Thus, the effects of rotation on MSTO morphology do change with age for our models.

That rotating MIST models are cooler and less luminous than their Geneva counterparts (top right panel, Figure 1) may stem from divergent approaches to rotationally induced angular momentum transport (such differences were also pointed out by Choi et al. 2016). Meanwhile, nonrotating MIST models appear more luminous than those of Geneva (top left panel, Figure 1), perhaps due to differing assumptions made for the treatment of convection.

Convective core overshoot (CCO) mixing is a chief aspect of convection that can influence the MSTO's CMD morphology and position. CCO is the idea that convectively driven material should not suddenly stop at the theoretical boundaries (e.g., Schwarzschild or Ledoux) of a convective zone. Rather, it seems feasible that momentum should carry material past these boundaries, allowing it to penetrate into non-convective zones, see e.g., Böhm (1963). The penetration distance, often considered as an extension of the convective zone, is often either a fraction of the pressure scale height (denoted H_P), or is modeled as a diffusive process with an exponentially decaying diffusion coefficient. We adopt the latter formalism, while Geneva adopts the former. Our adoptions (see Section 3.6.1 of Choi et al. 2016) are roughly equivalent to an extended convective core boundary of 0.2H_P (Section 2.1 of Magic et al. 2010), effectively twice the value adopted in Geneva.

Similar to rotationally enhanced mixing, CCO mixing can supply the stellar core with more fuel, enhancing energy production and the MS lifetimes of stars (e.g., Maeder 1975), thereby affecting a population's MSTO on a CMD. Studies presented by Gallart et al. (2003), Woo et al. (2003) and Bertelli et al. (2003) looked at the intermediate age clusters NGC 2173, SL 556, and NGC 2155, with nonrotating Padova (Girardi et al. 2000) and Yonsei–Yale (based on tracks from YREC, Guenther et al. 1992) stellar models. Modifying the efficiency of CCO appeared to provide one avenue for the models to explain the observed CMD features. Although, an alternative perspective on CCO is that material outside the convective core (CC) boundary, once mixed via overshoot or some other process, may drive an increase in opacity or the erasure of composition gradients, subsequently turning sub-adiabatic material super-adiabatic, and expanding the CC (discussed further in Section 2 of Paxton et al. 2018). CCO may partly be the result of a poorly defined CC boundary in our models. The proper treatment of CCO (and convection in general) is another active branch of development in stellar modeling.

Of course, CCO is not a rotational effect, but is discussed here to expound on the behavior shown in Figure 1 and to highlight its important (albeit uncertain) role in stellar evolution. Our models possess a stronger mixing due to CCO (see Choi et al. 2016, Section 3.6.2, and Ekström et al. 2012, Section 2.3). As CCO mixing occurs regardless of rotation, this is a possible cause for our nonrotating models maintaining a higher luminosity than their nonrotating Geneva counterparts.

In addition to the differences shown here on a Hertzsprung–Russell diagram and CMD, there are also differences in the MS lifetime extension afforded by enhanced rotational mixing. These differences were shown in Figure 20 of Choi et al. (2016), where the ratio ${\tau }_{\mathrm{MS},{\rm{R}}}/{\tau }_{\mathrm{MS},\mathrm{NR}}$ (MS lifetime with rotation over the same without rotation) is plotted against the initial stellar mass M_i (in units of M_⊙). Rotating Geneva models with M_i > 2 M_⊙ show a lifetime extension of roughly ∼25%, while corresponding MIST models only garner an increase of ⪅10%.

The physical adoptions made in Geneva do not equate to those made in MIST. The uncertainties present in stellar modeling (convection and rotationally enhanced mixing here; see, e.g., Salaris & Cassisi 2017 for a review) provide enough leeway for inconsistent model behavior. As will be seen, our findings are different from those of Brandt & Huang (2015a, 2015b), likely due to model differences.

3.1.2. Gravity Darkening

In addition to enhanced chemical mixing, rotation also introduces oblateness to the stellar structure. Latitudinally dependent centrifugal forces result in an oblate deformation of the star. The surface gravity, g, of a star is lessened by these forces leading to an effective, latitudinally dependent value, g_eff(θ). Here, θ refers to the polar angle in a spherical coordinate system. The effective gravity of a rotating star is related to its radiative luminosity (e.g., as encapsulated by the von Zeipel theorem; von Zeipel 1924), leading to a relation ${g}_{\mathrm{eff}}\propto {T}_{\mathrm{eff}}^{4}(\theta )$ between the effective gravity and temperature of a star. Hence, the observed colors (temperatures) and magnitudes (luminosities) become dependent on the viewing angle of observation. This effect is commonly termed gravity darkening.

The essential physical arguments of von Zeipel (1924) qualitatively adhere well to observations, but the power law relation presented in the von Zeipel theorem may be too simplistic. As it was derived based on radiative flux relations, the von Zeipel theorem does not directly apply to stars with convective envelopes; Lucy (1967) derived a more general ${g}_{\mathrm{eff}}\propto {T}_{\mathrm{eff}}^{\beta }$, with β = 0.08 for convective stars. More recently, in comparison to interferometric observations of rapidly rotating stars (McAlister et al. 2005; Aufdenberg et al. 2006; van Belle et al. 2006; Zhao et al. 2009), it has been found that ${g}_{\mathrm{eff}}\propto {T}_{\mathrm{eff}}^{4}$ overestimates temperature variation going between the pole and equator. In light of this, Espinosa Lara & Rieutord (2011) were motivated to derive a new formulation that can describe gravity darkening; they do so, deriving a relation that depends only the rotation rate of the star, at a given viewing angle, luminosity, and T_eff.

We use the equations of Espinosa Lara & Rieutord (2011) to translate the model stellar luminosity L, and effective temperature T_eff, across the desired viewing angles. The symbol i denotes the inclination (or viewing) angle in this paper; i = 0° corresponds to an observation directed at the equator, whereas i = 90° is directed at the star's pole.

A demonstration of gravity darkening effects is shown in Figure 3 on our MESA models. Several isochrones are displayed for two scenarios: ∼60 Myr models are shown in the left panel and ∼300 Myr in the right. In both cases, black shows an older isochrone while blue marks a slightly younger isochrone. The solid black line shows i = 0°, dashed is i = 45°, and dotted-dashed is i = 90°. In each case, the MSTO of an older isochrone moves toward mimicking that of a younger isochrone as the inclination angle varies from 0° to 90°, due to the increased luminosity and temperature of the stellar pole. Conceivably, if stars were to host a distribution of various inclination angles, there is the possibility that these effects would create a broadened MSTO.

Zoom In Zoom Out Reset image size

The final models that we fit to the data are constructed in MATCH (Dolphin 1997, 2002), a tool used to study resolved stellar populations (e.g., see de Jong et al. 2008; Gouliermis et al. 2011; Weisz et al. 2011). MATCH uses a given set of stellar models (our MESA models in this case) to create its own library of isochrones. These models have a finite resolution in age and metallicity; we have chosen 0.02 dex for both. Thus, it should be noted that there is an inherent spread to the models, as the isochrone parameters do not pertain to the delta functions in MATCH. Our composite stellar populations are constructed with the effects of gravity darkening via random viewing angles and we have also developed the ability for our synthetic populations to possess a distribution of rotation rates.

CMDs are populated according to an initial mass function (IMF), describing the occurrence of stellar masses, which are subsequently combined into a composite model population. For this purpose, we specify a Kroupa IMF (Kroupa 2001). Additionally, MATCH is able to consider binary systems when drawing its models, given a binary fraction. Binaries are added according to a flat distribution in mass fraction (i.e., from 0 to 1). So long as both stars are alive in the binary, magnitudes are the sum of the fluxes from each star; if the primary has died, the magnitude is of the secondary survivor only. Colors are computed by constructing the magnitude in each filter in this way, and then taking the difference of the magnitudes. MATCH does not model the interactions between binary companions. In the following discussion, we will list our adopted values of binary fraction where appropriate (adopted binary fractions are also listed in Table 1).

Table 1.  Adopted Cluster Parameters

Cluster	μa	AVb	${B}_{T},{V}_{T}$ Binary Fraction	$J,{K}_{s}$ Binary Fractionc
The Hyades	3.349	0.0031	0.0	0.25
	(Gaia Collaboration et al. 2017)	(Taylor 2006)
The Praesepe	6.26	0.0837	0.0	0.30
	(Gáspár et al. 2009)	(Taylor 2006)
The Pleiades	5.64	0.1054	0.0	0.0
	(Gaia Collaboration et al. 2016)	(Taylor 2008)

Notes.

^aDistance modulus. ^bInterstellar reddening. ^cNIR data has not been cleaned of binaries for the Hyades/Praesepe.

Download table as:  ASCIITypeset image

The ability to model gravity darkening was developed and added to MATCH for this project, where the brightness and temperature of our models are altered according to the equations of Espinosa Lara & Rieutord (2011). These alterations are a function of the model's viewing angle (which is drawn randomly as it is added to the composite stellar population under construction) and their rotation rate. In order to demonstrate these effects, we have created artificial datasets using the MATCH program "fake," which generates a set of photometric data for use on a CMD, given some finite period of star formation and metallicity variation. This program is also able to simulate photometric errors, but we take the errors to be zero in all cases, in order to highlight rotation related effects. For this reason, we have also increased model resolution to 0.01 dex in Figures 4 and 5, to lessen the spread due to finite composite population model age and metallicity resolution and create an appearance closer to an SSP.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Model Mass, [Fe/H], and Rotation Rate Range

Mass Range	[Fe/H] Range	$\tfrac{{\rm{\Omega }}}{{{\rm{\Omega }}}_{{\rm{c}}}}$ Range
0.1–8.0 M⊙	−0.75, −0.60, ..., 0.45	0.0, 0.1, ..., 0.6

Download table as:  ASCIITypeset image

The effects from gravity darkening are shown in Figure 4, as incorporated in MATCH. In Figure 4, an artificial dataset that neglects gravity darkening is shown in black, in comparison to an artificial dataset that does include model gravity darkening, which is shown in red. Several scenarios are demonstrated: the left column shows isochrones with Ω/Ω_c = 0.3, while the right shows Ω/Ω_c = 0.6. Meanwhile, the top row compares the effects of gravity darkening at 794 Myr (nearer to the ages of the Hyades and Praesepe), and the bottom row shows 100 Myr (near the age of the Pleiades). The inclusion of gravity darkening has a much stronger effect at higher rotation rates, issuing a greater spread to the MSTO.

Furthermore, stellar populations likely do not have stars rotating at a fixed value, rather they appear to possess a distribution of rates. Studies such as those of Zorec & Royer (2012) and Royer et al. (2007) have found evidence for a bimodal distribution of rotation rates for A and low mass B-type stars, with masses estimated to be around 1.5–3 M_⊙. This mass range is roughly appropriate for the MSTO stars in our target clusters (mentioned at the beginning of Section 3.1.1). However, there is evidence that this distribution may change with stellar type, where e.g., B and O-type stars appear to exhibit a singly peaked asymmetric distribution (e.g., Penny 1996; Huang et al. 2010; Ramírez-Agudelo et al. 2013).

For this work, the ability to model a distribution of rotation rates has been incorporated as a new feature in MATCH. We use a Gaussian distribution as a preliminary choice to model these effects. As the rotation rates of our models range from ${\rm{\Omega }}/{{\rm{\Omega }}}_{{\rm{c}}}\in [0.0,0.6]$, we choose a distribution with mean 0.3 and standard deviation 0.2. This distribution is truncated at the Ω/Ω_c parameter bounds of our model grid; hereafter this distribution is referenced as $P({\rm{\Omega }}/{{\rm{\Omega }}}_{{\rm{c}}}=0.3)$. The rotation rate of a proposed model is drawn from this Gaussian distribution as it is added to the composite population.

Although, as may be seen in Figure 8 of Huang et al. (2010), for the mass range 2.2 ≤ M/M_⊙ ≤ 4.0, stars may possess a distribution of $v/{v}_{{\rm{c}}}$ that is peaked nearer to 0.6 (derived from v sin i measurements acquired from spectra). This mass range roughly corresponds to the MSTO of the Pleiades, for instance, and so our chosen distribution may not fully represent the distribution of rotation rates in this cluster. In order to test the effects of a distribution that includes higher rotation rates, we have also tested a flat distribution of Ω/Ω_c between 0.0 and 0.6. Given that our models are currently limited to Ω/Ω_c = 0.6, we opt for this rather than creating a new distribution centered on Ω/Ω_c = 0.6. Generally, this flat distribution finds ages within 10 Myr and metallicities very similar to those found with our Gaussian distribution. Our usage of a flat distribution is not an extensive test of a distribution including higher rotation rates, for instance, as may be seen in Figure 4 the effects of rotation can vary dramatically toward higher rotation rates, so a lack of very fast rotators may neglect a wide range of model behavior. However, usage of a flat distribution is intended to give some preliminary sense of whether results would change significantly in the presence of a greater number of fast rotators or not. In future work, we plan to include models from Ω/Ω_c = 0.0 to 0.9, allowing us to test a wider range of possible rotation distributions.

In Figure 5, four panels display a progression of the consequences from including gravity darkening (top right), a distribution of rotation rates (bottom left), and finally both together (bottom right), in comparison to what a synthetic population looks like in the absence of these phenomena (black points). Models excluding a distribution of rotation rates are given the fixed value Ω/Ω_c = 0.3. With the relatively modest value of Ω/Ω_c = 0.3, gravity darkening has a weak influence on MSTO morphology, in comparison to that of a distribution of rotation rates (compare the top right and bottom left panels); although bear in mind that the effects of gravity darkening would increase with Ω/Ω_c, as shown in Figure 4.

Here we present our CMD fitting methodology. First, we describe MATCH and give an overview of its operation in fitting stellar models to the data (Section 3.3.1). Next, we discuss mock tests that were performed to demonstrate the accuracy of our results, using simulated observations that have stellar densities similar to our target clusters (Section 3.3.2).

3.3.1. MATCH

3.3.2. Testing Parameter Recovery with Mock Data

Here we perform a number of mock tests in order to demonstrate the level of accuracy that MATCH may achieve in fitting our models to the data. A rotation distribution and the modeling of gravity darkening are new additions to MATCH, and our data is fairly sparse in comparison to what MATCH is typically used for. Thus, we apply our models to an artificial dataset (described below) and test for the ability of MATCH to recover the input age and [Fe/H] values used to construct this artificial data. We vary the number N of artificial data points (i.e., N is the number of stars in our artificial data) to test accuracy as the observations become more sparse.

An artificial dataset with parameters similar to the Hyades, namely: log age = 8.90, [Fe/H] = 0.15 ± 0.01, μ = 3.34, A_V = 0.0031, was fit using our models in conjunction with MATCH. The artificial data is created using the previously mentioned MATCH program "fake." In these tests, our artificial clusters contain no multiplicity, and so our binary fraction is set to zero. Below, we show the results of several trial runs, fitting for the artificial cluster age and metallicity.

Figure 6 shows the results of these mock tests. Here we plot the best-fit age (top row) and metallicity (bottom row) versus N, where each point represents an average of 10 trial runs (to average over stochastic fluctuations). The red points correspond to models assigned a single rotation rate; the fit to artificial data at a single rotation rate. The blue points correspond to models where a Gaussian distribution of rotation rates was used; the fit to artificial data created with a Gaussian distribution of rotation rates. The errors shown are the average errors of the 10 trial runs. These errors were calculated via analysis of resultant posterior probability distribution, marginalized for the corresponding parameter; i.e., from the 16% and 84% percentiles of the posterior. The true values are represented by the light blue horizontal band, bounded by dashed black lines. For these tests, we are using the same models that will be used to fit the observations of our target clusters, which have a resolution 0.02 dex in age and metallicity. We create stars in a single formation episode whose duration in time and spread in metallicity is 0.02 dex (matching the resolution of our model composite stellar populations). Hence, the input age and metallicity bins, represented by the blue bands in Figure 6, span from log age = 8.90 to 8.92, and [Fe/H] = 0.14 to 0.16 dex, respectively.

Zoom In Zoom Out Reset image size

Ranging from artificial clusters comprised of several hundred to several thousand total stellar members, the input cluster parameters are often found, or are off, by at most 0.05 dex or so. This range in stellar numbers is chosen to replicate that of our target clusters. For instance, our sample of the Hyades contains lists of roughly one to several hundred stars in optical and NIR, respectively; the Pleiades lists contain several hundred; our list for the Praesepe consists of only 24 stars at the cluster turn off in the optical (although with ∼1000 total members in infrared for Praesepe). As seen in Figure 6, log age appears to not be recovered well in 2MASS J, K_s for the case N = 100 or so. However our NIR data lists contain nearer to 1000 total members.

Models become more degenerate in age and metallicity with fewer total stars, as both the MSTO and MS become less populous under IMF sampling. Accordingly, it is expected that the recovered parameter uncertainties should become larger as N is decreased, as may be seen in Figure 6. Still, the error bars remain relatively small (of order ≲0.03 dex for log age and [Fe/H] for the worst cases) for results derived with the lowest stellar densities tested here.

We explore the effects that stellar rotation has on derived cluster age and metallicity. With fixed distance and extinction, we fit to the MS and MSTO of the Hyades, Praesepe, and Pleiades. This is performed through Hess diagrams comparing that observed to the synthetic magnitudes and colors. A pair of results is presented sequentially for each cluster below; derived assuming a single population rotation rate, or assuming a Gaussian distribution of rates (referenced as $P({\rm{\Omega }}/{{\rm{\Omega }}}_{{\rm{c}}}=0.3);$ Section 3.3.1).

Figures 7–9 show the best-fitting isochrones overlaid with observed data from the Hyades, Praesepe, and Pleiades. These are representative isochrones from our MESA models; representative in that they exist on model grid points, while the reported best-fit values from MATCH belong to a continuum. The displayed isochrone ages are the closest grid point values to the MATCH reports. We have overlaid the Hess diagram of the best-fit $P({\rm{\Omega }}/{{\rm{\Omega }}}_{{\rm{c}}}=0.3)$ model in these figures to give a sense of what our models truly look like (shown in gray). For instance, refer back to Figure 5 for the general appearance of our models, with the effects of gravity darkening and a distribution of rotation rates broadening the MSTO.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Best-fit isochrones appear to trace the data qualitatively well. Although, in Figures 7 and 8, corresponding to the Hyades and Praesepe, it may be seen that the red clump region in each CMD does not fit very well, perhaps due to complexities in convection (see Section 5.3). Teal isochrones correspond to fits assuming no distribution of rotation rates. In black, the fit utilizing $P({\rm{\Omega }}/{{\rm{\Omega }}}_{{\rm{c}}}=0.3)$ is shown. Solid lines represent isochrones at the best-fit values of age and metallicity. The dashed lines are the best-fit isochrone at its ± uncertainty values in age; these are not displayed with $P({\rm{\Omega }}/{{\rm{\Omega }}}_{{\rm{c}}}=0.3)$, for clarity. The derived errors and best-fit values are collected in Table 3.

For the Hyades, we mostly find ages that are consistent with classical nonrotating CMD analyses, and the LDB age determined by Martín et al. (2018), in both filter sets. Metallicities are within ∼0.05 dex of the spectroscopic values. For instance, in B_T, V_T, the age is found to be ${676}_{-11}^{+67}$ Myr, and $[\mathrm{Fe}/{\rm{H}}]={0.23}_{-0.02}^{+0.02}$ when modeling a distribution of rotation rates. This age is nearer to the results found using nonrotating models, although on the higher end. The ages of ∼680 Myr are found for the Hyades, whether using a distribution of rotation rates or a fixed value. This is true for the Hyades in all cases, except with our age determined via NIR data when modeling a distribution of rotation rates.

For the Hyades in NIR, when using a Gaussian rotation distribution, we see an age of ${741}_{-11}^{+55}$ Myr. This older age may result from the greater color span of MSTO stars. In some sense, this resembles a broadened MSTO; using a distribution of rotation rates also makes the MSTO more broad. Thus, the $P({\rm{\Omega }}/{{\rm{\Omega }}}_{{\rm{c}}}=0.3)$ model can fit both the redder and bluer MSTO stars by sitting in the middle, driving an older age than the fixed Ω/Ω_c model, which is forced to fit either the redward or blueward MSTO stars. On this point, note that our best-fit single Ω/Ω_c model for the Hyades in NIR is high (0.6), compared to the lower value found in the optical data. This is due to the relatively low color span of MSTO stars in optical compared to NIR for the Hyades; the effects of gravity darkening produce a greater color spread in MSTO stars at higher rotation rates, driving the best-fit model to its higher Ω/Ω_c value. However, the MSTO spread in this case may be due to complications in the photometry. Membership and binarity uncertainties could be at play here, i.e., this older age could be spurious; (e.g., Kopytova et al. 2016 find a tighter NIR MSTO morphology, identifying single stars in the Hyades). We excluded several of the redder MSTO stars from this fit, indicated in Figure 7 by the open circles (there are 3 on the upper MSTO) in 2MASS J, K_s. This was done to reduce the influence of redward MSTO stars on the derived age. The Praesepe and Pleiades do not display as large a color span in their MSTOs with NIR data.

It may also be noted that a much higher value of Ω/Ω_c is preferred in fitting the NIR data of the Hyades (0.6 versus 0.3 in optical data). This is likely due to the MSTO of the Hyades having a greater spread in color and more MSTO stars in the NIR data. Adopting a higher Ω/Ω_c grants models a broadened MSTO via gravity darkening (e.g., Figure 4). We speculate that this allows a better fit to the NIR data, but is statistically worse in application to the Hyades optical morphology, which lacks a significant color spread in its MSTO. However, we leave a more thorough explanation of this discrepancy to future work.

Thus, for the Hyades, our ages mostly resemble those of classical nonrotating isochrone determinations, e.g., Perryman et al. (1998). The metallicity results roughly align with the measurement of [Fe/H] ≈ 0.15 from Cummings et al. (2017), although our results from the optical bands tend to favor values nearer to [Fe/H] = 0.20, while in NIR we find them nearer to [Fe/H] = 0.10.

The situation for the Praesepe is similar to that of the Hyades: we see results nearer to the reports from nonrotating models. This is the case across the NIR and optical datasets, whether using fixed Ω/Ω_c or $P({\rm{\Omega }}/{{\rm{\Omega }}}_{{\rm{c}}}=0.3)$. Here, the ages tend toward ∼590 Myr, while the best fits for [Fe/H] are near ∼0.09 in the NIR and in B_T, V_T, we find [Fe/H] ∼ 0.25. Here our optical dataset contains fewer MS stars than our NIR, making metallicity more difficult to determine, driving higher [Fe/H] values.

For the Pleiades, we similarly find ages that agree with classical, nonrotating analyses, as well as LDB results. Our ages fall within the range ∼112–160 Myr, in concordance with values derived from the LDB, and nonrotating isochrones. In the NIR, using a distribution of rotation rates, the ages have large uncertainties. The MSTO of the Pleiades in NIR is relatively more extended, driving this behavior. Values of [Fe/H] range from −0.01 to 0.40 (our [Fe/H] search boundary); these tend to be higher than the established values. Similar to the scenario seen with the optical dataset for the Praesepe, the number of stars here is very low: only 7. Our magnitude cuts of V_T < 5.0 and K_s < 5.0 were made to focus more on the MSTO, at the expense of constraints from the numerous MS Pleiads.

With the Pleiades, we did attempt fits with less severe magnitude cuts, but this in turn led to an increased age. On a Hess diagram, the low stellar density in the MSTO makes its stars appear as outliers in determining a global fit, despite these few stars being perhaps the most significant to extracting an accurate age. With relatively few stars, our uncertainties on age are larger for this cluster than the others. Moreover, many of the derived metallicities are unbounded and not well determined.

Figure 10 shows our best-fit age and [Fe/H] from fitting to optical data in comparison to several literature values; Figure 11 shows the same for our infrared data fits. For the Hyades, Praesepe, and Pleiades, we generally find ages that align with classical results (although see Section 4.1 about the higher age derived for the Hyades in our NIR fit). Our metallicities for the Hyades, Praesepe, in the optical, and the Pleiades in both optical and infrared, appear to be inconsistent with literature values. Metallicity is unconstrained in these cases, leading to a dubious metallicity determination that lies near the end of our search space in metallicity: 0.40 dex. That metallicity appears unconstrained here is probably due to the relatively low stellar number densities present in these datasets, and the issues that this may cause with analysis on a Hess diagram. Additionally, we excluded much of the MS in our fits to the Pleiades, sacrificing its use in constraining metallicity, opting to focus on the sparse, brighter Pleiads for an age determination that is driven more by the MSTO.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The derived ages of the Hyades and Praesepe from CMD analysis using nonrotating isochrones and the LDB congregate near 625–650 Myr (see references in Sections 2.1 and 2.2). For the Pleiades, nonrotating isochrone and LDB analysis have derived ages from about 100 to 130 Myr (see references in Section 2.3). Our ages are consistent with classical nonrotating isochrone CMD fits and LDB results (where available) for the Hyades, Praesepe, and Pleiades.

Our results contrast with the ages derived by Brandt & Huang (2015a) (see e.g., the points labeled BH15 in Figures 10 and 11). The derived ages of the Hyades and Praesepe from CMD analysis using rotating Geneva isochrones, carried out by Brandt & Huang (2015a, 2015c), using a Bayesian approach, are greater than classical reports. An age of 750 ± 100 Myr was derived for Hyades (800 ± 50 if fixing metallicity to [Fe/H] =0.10). The age of the Praesepe was determined to be 790 ±60 Myr. Lastly, an age of ∼95 ± 35 Myr was determined for the Pleiades, consistent with classical nonrotating age determinations and the LDB. Our derived ages agree with these studies for the Pleiades, but we do not find as great of an age increase for the Hyades nor the Praesepe via our models and methods.

The metallicity results for the Hyades and Praesepe are lower by ∼0.05 dex of the literature values in our 2MASS fits; in B_T, V_T, best-fit [Fe/H] for the Hyades and Praesepe are higher by ∼0.05 dex of the literature values. Our metallicity results for the Pleiades are generally inconsistent with the literature values (those being [Fe/H] ≈ 0.0); this is likely due to our magnitude cuts excluding much of the MS in favor of MSTO stars, and thus [Fe/H] becoming unconstrained in this CMD region. For instance, note the similarity in the CMD position of the isochrones in Figure 9 (right panel), although they differ by ∼0.10–0.15 dex in metallicity.

Generally, rotating models are preferred owing to a mildly higher probability over nonrotating models. MATCH uses a Poisson equivalent to χ² as a fit statistic (e.g., see Dolphin 2002, Section 2.3), designated via $-2\mathrm{ln}P$ (a lower value corresponds to higher probability), where P is the cumulative Poisson likelihood ratio, incorporating all Hess diagram bins. For the Hyades, we find $-2\mathrm{ln}P=104$ in optical photometry for the best-fitting model, while the nonrotating model shows $-2\mathrm{ln}P=120$, for instance. The best-fit model with a Gaussian distribution of rotation rates achieved $-2\mathrm{ln}P=101$ in this case. Typically, rotating models are preferred according to the fit statistic by roughly 3%–25% over nonrotating models.

The results derived via models with a Gaussian distribution of Ω/Ω_c or a fixed value are often similar. However, as these clusters contain relatively low stellar densities in their MSTOs, they do not provide a strong distinction between more realistic models that include a distribution of rotation rates, and those that possess a fixed population Ω/Ω_c value. With precise photometry from upcoming datasets like Gaia DR2, and in studying more populous clusters such as those in the LMC and SMC, we will have a better opportunity to assess the role of rotation distribution.

Our models behave differently from the Geneva models, as highlighted in Section 3.1.1, so it may be no surprise that our results differ from work that utilizes Geneva models. Rotation has more modest consequences for stellar evolution under our physical assumptions, so we do not observe a strong affect on derived ages due to stellar rotation. No major differences in derived age appear to manifest from using a Gaussian distribution of rotation rates as opposed to a single value. Indeed, the MSTOs present in these clusters tend to not be especially well populated, and so any spread that may exist due to a distribution of rotation rates may be difficult to observe.

As rotation tends to make the MSTO primarily cooler in our modeling, at least at ages near 600 Myr (see Figure 2), it makes sense that we do not find a significantly older age for the Hyades or Praesepe. Rather, adopting a younger age would increase the brightness and make the MSTO hotter, helping an isochrone compensate for the effects incurred by an increased rotation rate. In the case of the Pleiades at ∼120 Myr it is vice versa, here increasing rotation primarily makes the models brighter, leading models with higher Ω/Ω_c to mimic younger stars at a fixed age.³

This is in contrast to the results found by Brandt & Huang (2015a, 2015c), where the Praesepe and Hyades showed an age increase of roughly 200 Myr as a result of stellar rotation in their modeling. In Geneva models, the changes in isochrone morphology due to rotation are almost completely opposite to the effects seen in our models. Geneva models are hotter, brighter, and live longer on the MS (meaning the population evolves more slowly), as discussed in Section 3.1.1. It seems plausible that this discrepancy may stem from a more modest rotational mixing efficiency in MIST compared to Geneva. Whether this is precisely the case, or perhaps if other model differences come into play more strongly, will require further investigation. It is also important to bear in mind that our model set is limited to Ω/Ω_c = 0.0–0.6 presently, while the Geneva models allow study from Ω/Ω_c = 0.0 to 1. Extending our model grid to include Ω/Ω_c > 0.6 will be fruitful for comparison, aside from perhaps being necessary to study realistic rotation distributions (e.g., those discussed in Huang et al. 2010); we aim to incorporate models with higher rotation in future work.

Furthermore, our models possess a greater mixing due to CCO. This was shown in Figure 1 to make our nonrotating models hotter and brighter than their Geneva counterparts. The effects of CCO mixing are similar to rotationally enhanced mixing, effectively expanding the stellar core, granting it more fuel to burn longer, brighter, and hotter. Our nonrotating age determinations are near to that found by Brandt & Huang (2015a) with their rotating models (see Table 3 and Figures 10 and 11). It seems plausible that this is due to the higher level of CCO mixing in our models, as CCO makes the MSTO hotter and brighter. Proper 1D treatment of convection is yet another complication in stellar evolution theory and CMD-based analyses.

The physical assumptions made in our models (Section 3.1) and in Geneva are both able to simulate observational constraints. Much of our adopted physics follow from the assumptions made in MIST (Choi et al. 2016; see Sections 8 and 9 of that paper for comparisons of the data), while the physics adopted in Geneva (Ekström et al. 2012; see Section 5 of that paper for comparisons of the data) are similar, although there are differences, for example, in the treatment of rotational mixing and the assumed strength of CCO mixing. In our adopted formalism, the rotational mixing parameters, f_c and f_μ, are tuned to match the observed nitrogen enrichment in galactic MS B-type stars (with the observations from e.g., Gies & Lambert 1992; Kilian 1992; Morel et al. 2008; Hunter et al. 2009), following Heger et al. (2000). The Geneva models are capable of reproducing these observed surface abundances without having any parameters calibrated to do so in their formalism, followed from Maeder & Zahn (1998). Thus, although physical assumptions may differ between model sets, neither is ruled out by observational constraints thus far.

Although our fits appear to match the observed MS and MSTO regions relatively well for these clusters, in Figures 7 and 8, it is clear that our best-fit isochrones miss the red clump region. Figure 12 displays the range of convective mixing at various values, to demonstrate how it can affect the red clump region. In particular, we vary the parameter α_MLT, responsible for setting the enhanced range of convective mixing, according to the Mixing Length Theory (MLT) of Böhm-Vitense (1958). This parameter essentially dictates how far a fluid parcel travels, l_MLT, before mixing thermally with its surroundings; this length scale is characterized as some fraction of the pressure scale height, H_P, and is related to the constant α_MLT via ${l}_{\mathrm{MLT}}={\alpha }_{\mathrm{MLT}}{H}_{{\rm{P}}}$.

Zoom In Zoom Out Reset image size

By default α_MLT = 1.82 in our models, based on solar calibration (discussed further in Choi et al. 2016). Figure 12 shows isochrones, at age 649 Myr (similar to our derived ages for the Hyades and Praesepe), [Fe/H] = 0.15, Ω/Ω_c = 0.4, and inclination i = 0°, at various α_MLT = 1.60, 1.82 (default), and 2.00. Altering α_MLT makes the modeled red giant stars hotter or cooler, providing the potential for a better fit to a given dataset. In Figure 12 it is shown that a greater α_MLT may provide a better fit in the case of the Hyades. However, in the Praesepe this may lead to a worse fit, as it would pull the model away from two of the redder red clump stars.

Some 3D modeling efforts have provided evidence that α_MLT may vary with e.g., stellar T_eff, log g, (Trampedach et al. 2015) and metallicity (Magic et al. 2015). Evidence that α_MLT may vary with stellar parameters has also been seen in comparing models to observations, as in e.g., Bonaca et al. (2012), Tayar et al. (2017), Joyce & Chaboyer (2018), Chun et al. (2018), Li et al. (2018), and Viani et al. (2018). However, Choi et al. (2018) found that model T_eff can vary by nearly ±100 K according to the treatment of the surface boundary conditions in 1D codes. Thus, it appears that discrepancies between models and observations may (at least in part) be due to the chosen boundary conditions used in comparing models to the data. Determining precisely how α_MLT may vary, and according to which stellar properties, is an ongoing task. Efforts by Arnett & Moravveji (2017) and Mosumgaard et al. (2017), for example, will aid in having 3D simulations further inform our 1D models. Convection is another uncertain aspect in stellar modeling, alongside the uncertainties of stellar rotation.