ON USING THE COLOR–MAGNITUDE DIAGRAM MORPHOLOGY OF M67 TO TEST SOLAR ABUNDANCES

Most of the stellar model calculations presented in this paper were done with the GARSTEC (Weiss & Schlattl 2008) program. Here, we briefly present aspects or modifications of the code, which are relevant for this work. The standard assumptions and settings can be found in the quoted reference.

Convective overshooting is treated as a diffusive process according to the approach described by Freytag et al. (1996). The exact implementation can be found in Herwig et al. (1997). Based on hydrodynamical simulations, the diffusive constant therein (D_os) is assumed to decay exponentially beyond the formal Schwarzschild border as

$$\begin{equation} D_{\rm os}\left(z\right) = D_0 \ \exp \ \frac{-2z}{f H_P},\ D_0 = \frac{1}{3} v_0 \cdot H_P. \end{equation} \tag{ 1 }$$

This is equivalent to exponentially declining velocities of convective elements. In Equation (1), z is the radial distance from the Schwarzschild border and H_P is the pressure scale height taken there. The constant D₀ sets the scale of diffusive speed and depends on the convective velocity v₀ inside of the convective border. f is a free parameter defining the scale of the overshooting. It is known that fitting the CMD of young open clusters usually leads to an overshooting region extending for about 0.2 H_P in the classical local prescription. We have computed some test models for stars between 2 and 6  M_☉ and found that this overshooting value is well reproduced by the diffusive approach, during the hydrogen core burning phase, with a value of about f = 0.018, a confirmation of previous results by Herwig et al. (1997).

For small convective cores, the amount of overshooting has to be limited. In the local description this is done, for example, by a multiplicative factor, which increases linearly with mass in the range of 1–2 M_☉ (see Pietrinferni et al. 2004, for an example). In our code, this is achieved by a geometrical cutoff, where H_P in Equation (1) is replaced by

$$\begin{equation} \widetilde{H_P}{:}\,{=}\, H_P \cdot \min \left[1,\left(\frac{\Delta R_{\rm cz}}{H_P}\right)^2 \right], \end{equation} \tag{ 2 }$$

where ΔR_cz is the thickness of the convective zone. In this way, the geometric cutoff has the advantage that the overshooting region is always limited to a fraction of the extension of the convective region. It is our standard limiting procedure for overshooting. A similar cut is used in Ventura et al. (2005, 1998).

The alternative approach, namely the use of a ramp function for the f parameter, is also applied for some of our calculations. For masses lower or equal than 1.1 M_☉ we use f = 0 (i.e., no overshooting), whereas for masses equal to or higher than 2.0 M_☉ we set f = 0.018. In the intermediate mass range, the efficiency of overshooting varies linearly with mass

$$\begin{equation} f = (0.13\,M/M_{\odot } - 0.098)/9. \end{equation} \tag{ 3 }$$

The small discontinuity in the overshooting efficiency is of no consequence for our study. 1.1 M_☉ models with f = 0.005, as would result from application of the above equation to this mass, lead, given the low overshooting efficiency, to the same evolutionary tracks as models without overshooting. This prescription for the overshooting efficiency leads to evolutionary tracks that reproduce well those by Pietrinferni et al. (2004).

Atomic diffusion is treated within the same diffusive numerical scheme. In this paper, hydrogen, helium, and heavier elements (including iron) are diffusing. Radiative levitation is not taken into account.

To agree with the physical assumptions in VG07, the default set of nuclear reaction rates is those of the NACRE collaboration (Angulo et al. 1999). To investigate their influence for individual reactions alternative rates were used. For the ¹⁴N(p, γ)¹⁵O reaction, we also used the rate from Marta et al. (2008), the newest result from the LUNA collaboration, which is lower by about a factor of 2 with respect to the NACRE rate at the relevant temperatures. This has a considerable effect on the TO morphology. Another reaction that turned out to be of surprisingly strong influence is the ¹⁷O(p, α)¹⁴N reaction. We tested this by employing either the recent measurements by Moazen et al. (2007), which for T > 4 × 10⁸ K is very similar to the NACRE recommendation, or the one given by Caughlan & Fowler (1988). A detailed discussion on the reaction rates can be found in Section 4.2.

For both solar composition choices (AGS05, GS98), consistent Rosseland mean opacity tables were prepared following the procedure described in Weiss & Schlattl (2008).

The final step in comparing isochrones with observed CMDs is the transformation to colors. For this we used that by VandenBerg & Clem (2003, VC03) which, together with a choice for the distance modulus and reddening of M67, results in satisfying CMD fits on the main sequence and subgiant branch. Alternatively, we used the transformations by Cassisi et al. (2004) for testing purposes.

The initial stellar parameters (Y_in, Z_in, mixing length parameter) are obtained from solar model calibrations. This will be discussed in more detail in Section 3.2. With these, we computed the evolution from the zero-age main sequence to the tip of the red giant branch (RGB) for mass values from 0.6 to 1.5 M_☉ in steps of 0.1 M_☉. We recall here that, as in VG07, we are not interested in considering the RGB of M67, where deficiencies in our treatment of the outer layers of the star could have an impact on the location of the models in the Hertzsprung-Russell diagram (HRD). We confirmed by tests starting from the pre-main sequence that for stars with mass below the critical mass for the occurrence of a convective core, M_ccc, the transient convective core at the end of the pre-main sequence phase also appears here. This convective core is the result of the short phase of CN conversion, but may be sustained if convective overshooting is included. Finally, for constructing isochrones, the tracks are normalized to the so-called equivalent points (Bergbusch & Vandenberg 1992; Pietrinferni et al. 2004) and the interpolation to an isochrone is done between the normalized tracks. After an isochrone age was fixed, we calculated an additional model with the TO mass and recalculated the isochrone to make sure that the TO morphology does not depend on the isochrone interpolation scheme.

To compare our models with M67 we used the photometric data by Sandquist (2004). These are more accurate than the older data from Montgomery et al. (1993), which were used by VG07, and contain bona fide single stars only. Thus, this CMD has a narrower main-sequence band. However, we also tested some of our isochrones with the data by Montgomery et al. (1993). We used both (B − V) and (V − I) colors.

The dereddened distance modulus to M67 is (m −  M)_V = 9.70 according to VG07, and 9.72 ± 0.05 (Sandquist 2004) based on subdwarf fitting to the lower main sequence following Percival et al. (2003). The reddening is E(B − V) = 0.038 (VG07) in good agreement with Sarajedini et al. (1999), who gave 0.04 ± 0.02. Similar values have been recently obtained by Twarog et al. (2009). For the metallicity we assumed a value of [Fe/H] = 0.0, which is well within spectroscopically determined errors, for example by Gratton (2000), who gave [Fe/H] = 0.02 ± 0.06.

To be consistent with VG07, we used GS98 and AGS05 for the old and new solar abundances, although the former are simply an update of GN93 taking into account additional literature values, and AGS09 would be the most recent and complete re-analysis. Since AGS09 abundances are slightly higher than AGS05, our choice is to test the more extreme case. The initial abundances of the stellar models for the M67 isochrones are taken to be the same as those resulting from solar model calibrations.

The idea of VG07 to use M67 for testing the effect of the new solar abundances rests on the fact that the TO mass of this open cluster is very close to the critical mass for the onset of core convection, M_ccc, or "transition mass" (VandenBerg et al. 2007). This transition mass depends on the average exponent of the energy generation rate in the core (Kippenhahn & Weigert 1990, chap. 22), which increases with the contribution from the CNO cycle, for which _CNO ∼ T¹⁵, in contrast to _pp ∼ T⁵ for the pp-chains at the relevant temperatures. The importance of the CNO cycle depends on the nuclear reaction rates (see Section 4.2) as well as on the amount of "catalysts," i.e., the sum of the CNO-abundances. These are, incidentally, the elements with the largest reduction in their abundance according to AGS05. Therefore, the morphological change in the CMD, i.e., displaying the characteristic hook at the TO, or a gap in star density, can be used to determine whether the TO mass is below or above M_ccc.

However, the sensitivity of the M67 TO morphology implies that other effects may also influence it, apart from the metallicity. The nuclear reaction rates were already mentioned. Other possible aspects could be the amount of overshooting, atomic diffusion, and the pre-main sequence history. Last but not least, technical details of the stellar evolution codes may play a role. It is therefore necessary to show that our results do not depend on the particular numerical code. As a first and crucial step, we have successfully attempted to reproduce the key result by VG07. This will also show where changes to their procedure are indicated.

The first step concerns the solar model calibration. Here, as for the M67 models, VG07 used NACRE nuclear reaction rates, the OPAL (Iglesias & Rogers 1996) and Ferguson et al. (2005) opacities (as is done in our code), but ignored diffusion. This point is crucial, as we will show later on. Some amount of convective overshooting (see Section 4.1) was included for the M67 models, but is not relevant for the solar calibration. For the atmospheres they used MARCS model (Gustafsson et al. 2008). In this respect we differ since we use standard Eddington gray atmospheres. However, as shown by VandenBerg et al. (2008), this has in our present case no significant influence on the tracks on the main-sequence and subgiant branch.

Table 1 contains in the first two rows the resulting solar model parameters by VG07 and in rows 3 and 4 the equivalent ones obtained with our code. Additionally, we also show the results when using the Dartmouth and LPCODE codes (Section 4.3). Note that the mixing length parameters can never be compared in their absolute values due to the different formulations of Mixing Length Theory (MLT). Overall, the agreement with VG07 for the initial abundances is within 1%–2%, with our codes returning systematically lower initial metallicities (by 3%) for the AGS05 mixture. Although a small effect, it additionally disfavors the appearance of a convective core in the M67 TO.

Using these initial composition values we calculated stellar tracks and isochrones for M67, which we show in Figure 1. The numerical values for distance and reddening are identical to VG07. The best-fitting isochrones have somewhat higher ages (by 0.3 Gyr). The TO mass is 1.229 and 1.200  M_☉, respectively, for the old and new composition. One recognizes the same basic result as in VG07: for the AGS05 mixture, the isochrone does not show the characteristic hook. However, it displays a slight inclination indicating the presence of a very small convective core at the TO mass. Consistent with this morphological sign, the TO mass is marginally higher than M_ccc, which is 1.175  M_☉, while in the case of the GS98 mixture it is clearly above M_ccc = 1.139  M_☉. These values appear to be very similar to those found by VG07. For these models, as well for other important models presented later in this work, the M_ccc values are listed in Table 2. This exercise demonstrates that we are able to reproduce correctly VG07 and that their result is independent of the stellar evolution code used.

Zoom In Zoom Out Reset image size

The effects of diffusion on the value of M_ccc has been studied by Michaud et al. (2004) for solar compositions with high metallicity, and applied to fitting M67 and NGC 188 CMDs. They found that atomic diffusion (effectively sedimentation) reduces M_ccc by increasing the metallicity of the core over time. As emphasized by VG07, atomic diffusion could then modify the results presented in the previous section. As explained in the previous subsection, an increase of the core metallicity, particularly of CNO elements, will favor the occurrence of a convective core and reduce M_ccc. The solar model calibration, however, should also include diffusion, as it was shown (e.g., Bahcall et al. 1995; Christensen-Dalsgaard et al. 1996) that only with this physical effect the best agreement with the seismic Sun can be achieved. Only in this case the solar parameters, in particular the initial composition, are determined as accurately as possible.

The corresponding results are listed in Table 1 as well. Obviously, both helium and in particular metallicity are now much higher as compared to the previous case ignoring diffusion in the calibration; the latter values are increased by 10%–13% (similar values for the comparison codes). We are using the same solar initial abundances for the M67 isochrones. This is not completely consistent as M67 currently displays a solar metallicity, but due to its slightly younger age diffusion it should have started with a somewhat lower metallicity than the Sun. However, this effect amounts to a few parts in 10⁻⁴ in Z only, such that we did not iterate further the initial composition.

The final fits to M67 are shown in Figure 2. We still used the same distance and reddening values, but the isochrone ages are now lower by 0.3 Gyr compared to Figure 1. The agreement with the ages by VG07 is purely incidental. Instead, the age for the GS98 mixture of 3.9 Gyr should be compared to that by Michaud et al. (2004) of 3.7 Gyr. Now, for both compositions the CMD morphology is reproduced, and AGS05 is no longer disfavored. The fit itself is of the same quality as before.

Zoom In Zoom Out Reset image size

We first repeated the models of Section 3.1, i.e., diffusion was completely ignored, but included overshooting according to our standard prescription outlined in Section 2.1, Equations (1) and (2). The parameter f can be varied, but we keep it below the standard value of 0.02 (Herwig et al. 1997). In this work, we take f = 0.018 because with this value, in our alternative approach (Equation (3)), we are able to reproduce well the evolutionary tracks from Pietrinferni et al. (2004), as mentioned in Section 2.1. In Figure 3, the resulting CMD for f = 0.018 and the standard cluster parameters are displayed. We recall that the amount of overshooting is reduced due to our geometrical restriction (Equation (2)); this reduction can be quite restrictive for stars with lower masses, around 1.2 M_☉, and be crucial in this study. However, even if the fit itself is not as good as before, now both mixture cases have TO masses with a convective core being present.

Zoom In Zoom Out Reset image size

The situation seems to be typical for cases without diffusion in both the solar calibration and the cluster models. In VandenBerg & Stetson (2004), the authors concluded that a small amount of overshooting is needed to best reproduce the CMD of M67. Overshooting was treated following the integral criterion by Roxburgh (1989), and the adaption of a varying amount of overshooting for very small convective cores close to M_ccc determined by VandenBerg et al. (2006). We stress the fact that this procedure was calibrated using open clusters with solar abundances on the old GS98 scale. For complete consistency, this procedure should be repeated for the new AGS05 composition, and this might help explaining why the same small amount of overshooting used by VG07 did not lead to a persistent convective core in case of the new abundances.

The case of Section 3.2, i.e., the one with full consideration of diffusion, was also extended by including overshooting in the tracks for M67. We refrain from showing the result here, as for both mixtures the convective core, as expected, persisted, since the fit is similar to Figure 3, with overshooting increasing the TO luminosity, and leading to a slightly degraded fit. This coincides with the isochrone fit shown by Pietrinferni et al. (2004), where the CMD could be better reproduced if overshooting was completely ignored. It agrees also with Michaud et al. (2004), who also found no necessity for overshooting, and who also included diffusion in their models.

As mentioned above, the nuclear reaction rates of the CNO cycle are equally important for the occurrence of a convective core as is the abundance of CNO-nuclei. So far we have presented models employing the NACRE-library reaction rates for these reactions. The bottleneck reaction of the CNO cycle, ¹⁴N(p, γ)¹⁵O, that determines the overall cycle rate has been measured at stellar energies in the laboratory by the LUNA collaboration (Formicola et al. 2004; Marta et al. 2008) and was found to be lower by about 50% with respect to the NACRE rate. The consequences for globular cluster age determinations and for some sensitive phases of low- and intermediate-mass star evolution have been investigated by Imbriani et al. (2004) and Weiss et al. (2005).

The rate being lower, we expect that the transition mass to harboring a convective core increases. Therefore, isochrones using the lower CNO-abundances of the AGS05 mixtures will be even less likely to show the TO hook. We repeated the case of Section 3.1, (i.e., ignoring diffusion completely) with the updated and most likely more accurate rate by Marta et al. (2008). The resulting CMD for M67, again using our standard distance modulus and reddening, is shown in Figure 4. In agreement with Imbriani et al. (2004), the isochrone age had to be increased (by 0.3 Gyr for both mixtures) and M_ccc increased by ≈0.08 M_☉ in both cases. As a consequence the GS98 case also now lacks the characteristic hook. We add briefly that the inclusion of overshooting does not alter this result because the geometric cutoff is restrictive enough that the convective region formed at the end of the pre-main sequence phase cannot be maintained during the main sequence evolution. The values for M_ccc are 1.215 (GS98) and 1.258 (AGS05), as compared with TO masses of 1.202 and 1.196 M_☉, respectively.

Zoom In Zoom Out Reset image size

This illustrates the fundamental problem with such tests: imagine that we would assume that the GS98 solar composition is the correct one, but want to test which reaction rate is to be preferred. From Figures 1 and 4, we would clearly conclude that the older one is to be preferred.

We now repeat the computation of Section 3.2, i.e., the case with diffusion included, but with the newer and lower ¹⁴N(p, γ)¹⁵O reaction rate. In this case, we recognize that with the older solar composition a small convective core is present (TO mass/M_ccc = 1.214/1.172 M_☉), which is completely absent in the AGS05 case (1.201/1.241 M_☉). While the CMD fit is not very good even with the GS98 isochrone, it can be improved by including overshooting. The results, including overshooting with the geometric cut, are shown in the left panel of Figure 5. As mentioned in the previous section, the geometric cutoff limit for overshooting is very restrictive particularly for lower masses. Therefore, we recompute the tracks with the alternative limiting algorithm, the ramp function described by Equation (3). Changes in TO morphology for the GS98 composition are only minor. On the other hand, for the AGS05 composition, we now get a similarly good CMD fit as with the old mixture (right panel of Figure 5), as the convective core that appears in the pre-main sequence can be maintained during main sequence evolution.

Zoom In Zoom Out Reset image size

We found a second reaction that is also influencing the efficiency of the CNO-double cycle, though indirectly. This is the ¹⁷O(p, α)¹⁴N reaction that closes the CNO-II cycle and which is together with the competing ¹⁷O(p, γ)¹⁸F the slowest one of that subcycle. For this reaction new cross section measurements by Moazen et al. (2007) and Chafa et al. (2007) are available. The resulting two new reaction rates are quite similar to each other and vary in the interesting temperature range of 15–22 MK between ∼50% and ∼135% of the NACRE rate. Additionally, we have also considered for this reaction the rate given by Caughlan & Fowler (1988, CF88). The comparison of all rates is shown in Figure 6.

Zoom In Zoom Out Reset image size

In Figure 7, we show the evolutionary track of a 1.2 M_☉ stellar model with GS98 composition for different choices of the two key nuclear reactions mentioned above. The track corresponding to our standard choice of nuclear rates, NACRE, is shown by the dotted line. Now, if the CF88 rate for the ¹⁷O(p, α)¹⁴N reaction is used, which is less than 10% of the NACRE rate throughout the relevant temperature regime, the track for this crucial mass is strongly modified and no longer shows a sign of a convective core, as depicted by the long-dashed line. If the Moazen et al. (2007) rate is used, the differences with respect to the track using NACRE rates are very small; for this reason, we do not show this track. However, it is interesting to note that the effect that changing the ¹⁷O(p, α)¹⁴N rate has on the tracks depends on the rate used for the ¹⁴N(p, γ)¹⁵O rate. Indeed, as stated before, using the LUNA rate for this reaction (and NACRE rates for all others) leads to the absence of a convective core in the evolution of a 1.2 M_☉ stellar model. This is shown in the short-dashed line in Figure 7. However, when in addition the Moazen et al. (2007) rate is used, a small convective core is apparent in the track, as shown by the solid line in the same figure.

Zoom In Zoom Out Reset image size

We have identified the reason for this to be the following: the low ¹⁷O(p, α)¹⁴N rate from CF88 slows down the CNO-II subcycle, reducing thereby the flux of ¹⁴N back into the CNO-I cycle, which produces the overwhelming majority of the energy. While the loss of energy from CNO-II is rather unimportant, the drainage of available ¹⁴N by up to 90% from the CNO-I cycle is significant for its efficiency. One may consider this as a storage of nitrogen in the form of useless ¹⁷O. This makes the appearance of a convective core more difficult, and thus raises M_ccc. Note that for the higher, older, ¹⁴N(p, γ)¹⁵O rate from NACRE this effect is not important; it appears that in this case the branching is favoring the CNO-I cycle in any case, and the CNO-II cycle is always negligible.

Since the importance of the ¹⁷O(p, α)¹⁴N was rather surprising, and to make sure that our H-burning network included all necessary reactions, we did a test calculation with a variant of our code (M. Alves-Cruz 2009, private communication) that includes an extensive p-capture network up to silicon. In particular, the ¹⁷O(p, γ)¹⁸F reaction, which is the branching reaction to the next higher cycle, is included. We saw no difference in the evolution of stars in this mass range up to post main-sequence phase. Our treatment of the CNO cycles thus appears to be complete and correct.

We have seen that the behavior of the convective core for masses around the TO mass in M67 is very sensitive to the constitutive physics included in the stellar models. Given this sensitivity, it is also important to determine if calculations done with different stellar evolution codes lead to the same conclusions. Here, we have used, in addition to GARSTEC, results from two additional codes: LPCODE and the Dartmouth code. We summarize our findings below.

LPCODE was originally developed at La Plata Observatory. Extensive descriptions of the code have been presented in Althaus et al. (2002, 2003). Here, we only describe updates and changes done in the code that were specifically implemented for this work. The equation of state has been updated to the latest release by OPAL.⁵ Radiative opacities from OPAL were calculated for both GS98 and AGS05 compositions and complemented at low temperatures with those from Ferguson et al. (2005). Conductive opacities, although not relevant for this work, are now those from Potekhin and collaborators, presented in Cassisi et al. (2007). For the nuclear reaction rates, NACRE is the standard choice but the LUNA rate for the ¹⁴N(p, γ)¹⁵O from Marta et al. (2008) has been adopted when necessary. Overshooting is treated in a diffusive approach, but the geometric cutoff (Equation (2)) has been implemented to allow better comparison with GARSTEC results.

We have computed a large set of models covering the ranges of mass, composition, and input physics relevant to our work. In all cases, agreement between GARSTEC and LPCODE has been very good. Here, in Figure 8 we show some evolutionary tracks that include the most complex input physics, as used for the CMD in the left panel of Figure 5: NACRE rates updated by the ¹⁴N(p, γ)¹⁵O reaction rate from Marta et al. (2008), diffusion in both the solar calibration and in the tracks, and overshooting limited by our geometrical restriction (Equation (2)). Left and right panels show results for AGS05 and GS98 compositions, respectively. In both cases, tracks for 1.1, 1.2, and 1.3 M_☉ are shown that bracket the TO mass of M67 and, consequently, determine the morphology of the isochrones around the TO. By inspection of Figure 8, it can readily be seen that models with both codes show a very similar evolution along the HRD; differences are hardly noticeable. Results are of a similar quality for other choices of constitutive physics and for both composition options. In relation to systematic uncertainties affecting the conclusions of our work, this result is particularly encouraging because GARSTEC and LPCODE have been developed completely independently, and do not share any numerical algorithm for solving the equations of stellar evolution.

Zoom In Zoom Out Reset image size

The Dartmouth stellar evolution code is derived from the Yale code (Guenther et al. 1992), with modifications and updates described in Chaboyer et al. (2001), Bjork & Chaboyer (2006) and Dotter et al. (2008). For this project, the NACRE nuclear reaction rates were implemented, with the exception that the LUNA rate (Marta et al. 2008) for the ¹⁴N(p, γ)¹⁵O reaction was used. Convective core overshoot in the Dartmouth code is parameterized as a multiple of H_P. Normally, the prescription of Demarque et al. (2004) is used in which the amount of core overshoot is small (0.05 H_P) for stars with small convective core masses, and gradually increased to 0.20 H_P for stars with masses 0.2 M_☉ above the critical mass for the turn-on of convection in the core. This prescription for convective core overshoot was found to yield good agreement with open cluster CMDs, when the Grevesse et al. (1996) solar mixture was used. This mixture has (Z/X)_☉ = 0.0244. For this project, a small constant amount of convective core overshooting regardless of the stellar mass was used. The amount of convective core overshoot used in specific models is identified in the figure captions, and was typically 0.07 H_P.

In the case of the comparison between results of GARSTEC and the Dartmouth code we choose to present CMD fits to M67 computed with both AGS05 and GS98 compositions. Models for GARSTEC use overshooting constrained by Equation (3). Results for the best-fit isochrones for both codes are shown in Figure 9, where the left (right) panel shows results for the AGS05 (GS98) composition. As throughout this paper, our goal is to determine whether isochrones based on stellar models with one or the other composition can reproduce the occurrence of the observed hook in M67. As can be seen in Figure 9, isochrones from both codes reproduce well the TO morphology for both compositions. In the case of AGS05, the age for M67 obtained from both codes is 4.5 Gyr and for GS98 is 4.2 Gyr. Although the TO morphologies are not identical for both codes, some differences are likely to be present since we have not attempted, even within the uncertainties given by the input physics and observational data, to obtain the best possible agreement between the codes and with the data.

Zoom In Zoom Out Reset image size

As a last comparison, we state that the results determined in VG07 (see Section 3.1) could also be reproduced by both codes.

From comparing results from GARSTEC with those from LPCODE and Dartmouth code, we conclude that, even if numerical schemes for solving stellar structure and evolution equations and implementation of physics are different in the different codes, our conclusions are robust; they do not depend on the stellar code used. It implies that systematics between the codes are not an important source of uncertainty in our conclusions.