# Isochrone fitting to the open cluster M67 in the era of Gaia and improved
model physics

Claudia Reyes1, Dennis Stello2,1, Marc Hon3,4, Regner Trampedach5, Eric
Sandquist6, and Marc H. Pinsonneault7.
1School of Physics, University of New South Wales, NSW 2052, Australia
2Sydney Institute for Astronomy (SIfA), School of Physics, University of
Sydney, NSW 2006, Australia
3 Department of Physics and Kavli Institute for Astrophysics and Space
Research, Massachusetts Institute of Technology, 77 Massachusetts Ave,
Cambridge,
MA 02139, USA.
4Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu,
HI 96822, USA
5Space Science Institute, 4765 Walnut Street, Boulder, CO 80301, USA
6Department of Astronomy, San Diego State University, San Diego, CA 92182, USA
7Department of Astronomy, The Ohio State University, Columbus, OH 43210, USA

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The Gaia mission has provided highly accurate observations that have
significantly reduced the scatter in the colour-magnitude diagrams of open
clusters. As a result of the improved isochrone sequence of the open cluster
M67, we have created new stellar models that avoid commonly used
simplifications in 1D stellar modelling, such as mass-independent core
overshooting and a constant mixing length parameter. This has enabled us to
deliver a precise isochrone specifically designed for M67, available for
download. We follow a commonly-used qualitative approach to adjust the input
physics to match the well-defined colour-magnitude sequence, and we test the
model-predicted masses against a known eclipsing binary system at the main
sequence turnoff of the cluster. Despite using improvements in photometry and
stellar physics we cannot match the masses of both binary components with the
same theoretical isochrone. A $\chi^{2}$-based isochrone fitting approach
using our preferred input physics results in a cluster age of
$3.95^{+0.16}_{-0.15}$ Gyrs.

###### keywords:

Hertzsprung–Russell and colour–magnitude diagrams – open clusters and
associations: individual: M67

††pubyear: 2024††pagerange: Isochrone fitting to the open cluster M67 in the
era of Gaia and improved model physics–Isochrone fitting to the open cluster
M67 in the era of Gaia and improved model physics

##  1 Introduction

The old open cluster M67 (NGC 2682) is known to be an important testbed for
stellar evolution studies because of its low reddening, near-solar
metallicity, and rich stellar population. Its turnoff stars fall right in the
mass range where we can infer how convective core size depend on stellar mass,
making M67 particularly challenging, but also one of the only sites where we
can use the colour-magnitude diagram to map out the convective core effects so
readily. The colour-magnitude diagram of M67, together with stellar models,
has therefore been used to determine the extent of convective core
overshooting, to calibrate colour-temperature relations, and to test solar
abundances. Examples include VandenBerg & Stetson (2004) and VandenBerg et al.
(2007) who used $BV$ photometry from Montgomery et al. (1993) and later Magic
et al. (2010) who used the higher precision $BV$ photometry from Sandquist
(2004).

An unusually large fraction of binaries and blue stragglers have been revealed
through the study of M67’s colour-magnitude diagram (Strom et al., 1971;
Mathieu & Latham, 1986; Mathys, 1991; Gilliland et al., 1993; Mathieu et al.,
2003), and through spectroscopy (Geller et al., 2015; Leiner et al., 2019;
Geller et al., 2021). Various spectroscopic studies have also focused on the
chemical abundances of M67 and have reported the signature of diffusion
processes in the cluster’s stars (Souto et al., 2018; Bertelli Motta et al.,
2018; Souto et al., 2019; Liu et al., 2019; Beeson et al., 2024). In addition,
the cluster’s giants have been the targets of asteroseismic studies with both
ground-based (Gilliland et al., 1993; Stello et al., 2006) and space-based
telescopes (Stello et al., 2016; Li et al., 2024).

Since the first release of data from the Gaia mission (Gaia Collaboration et
al., 2016), the observed colour-magnitude diagrams of clusters like M67 became
significantly less scattered, and as a result new efforts to match theoretical
isochrones to clusters have been made, including to M67 (Gaia Collaboration et
al. (2018b); Choi et al. (2018); Nguyen et al. (2022)). However, none of the
existing Gaia-era isochrones were explicitly created for M67 as the end goal,
and as a result they do not match the Gaia data particularly well (Appendix B,
Figure 12). Furthermore, they use the common approach of adopting a fixed
mixing length parameter, $\alpha_{\mathrm{MLT}}$, usually calibrated to the
Sun, although there is no strong evidence suggesting $\alpha_{\mathrm{MLT}}$
is the same for stars of different evolutionary stages (Song et al., 2020).
This simplification can be expected to influence the temperatures and
luminosities of stars with convective envelopes (Joyce & Tayar, 2023).

In this paper, we aim to develop a new theoretical isochrone that incorporates
the latest research on the physics governing stellar evolution: notably, model
improvements based on hydrodynamical 3D simulations of convection (Trampedach
et al., 2014a, b) while improving the fit to the colour-magnitude diagram.
Such isochrone, validated by Gaia’s M67 data, will provide an isochrone
reference for future studies. Additionally, we address the issue of the
eclipsing binary WOCS11028, whose stellar components have highly accurate
determinations of mass and radius and are found to likely not have interacted
in the past (See Sandquist et al., 2021, for an extensive discussion). The
problem is that no existing isochrone models simultaneously match these
properties for both components.

Our approach is as follows: we begin using a traditional qualitative isochrone
fitting approach (VandenBerg & Stetson, 2004; Pinsonneault et al., 2004;
Barker & Paust, 2018) for empirical calibration of the input physics. Lastly,
inspired by Li et al. (2023), who explored a two-parameter grid to estimate
age and extinction likelihoods in NGC 2516, we estimate age and distance
modulus likelihoods for our models given our preferred input physics using a
$\chi^{2}$ approach.

##  2 The M67 Colour-Magnitude Diagram

###  2.1 Membership sample

We determined cluster membership using Gaia DR3 (Gaia Collaboration, 2022)
from stars within 2 degrees of the centre of the cluster at
$\mathrm{RA}=132.85^{\circ}$, $\mathrm{Dec}=11.81^{\circ}$ (Figure 1a),
selecting the overdensity of stars with radial velocities within $34\pm
13.6\,\mathrm{km/s}$ and parallax $1.15\,\pm 0.12\,\mathrm{mas}$ (Figure 1b),
and proper motions within $0.7\,\mathrm{mas/year}$ from
$\mu_{\alpha}=-11.0\,\mathrm{mas/year}$ and
$\mu_{\delta}=-2.9\,\mathrm{mas/year}$ in the proper-motion space (Figure 1c).
To boost the number of stars along the less populated giant branch, we further
added ten giant stars with secure kinematic membership from ground-based
radial velocities by Geller et al. (2015) that: (1) did not have DR3 radial
velocities available (white circles in Figures 1a and 1c), or (2) narrowly
missed our Gaia-based cuts in parallax or proper motion (orange and/or yellow
circles in Figure 1), but were later confirmed by colour-magnitude diagram
position to belong to the cluster’s standard-evolution single-star sequence.
In total, our selection resulted in 488 cluster members: 22 non-standard
evolution stars such as blue stragglers, yellow giants, and sub-subgiants; 361
main sequence stars, 63 subgiants, 35 red giant branch stars, and 7 red clump
stars. This number includes four new stars that were not identified by Geller
et al. (2015) as M67 members, but met all of our selection criteria for
cluster membership, as shown in Figure 2 with yellow squared symbols. Table 1
lists the stars mentioned above, along with their corresponding IDs from
various catalogues.

Figure 1: In small black symbols are cluster members as per our selection using DR3 values. Large circles are the ten stars we added because they were deemed cluster members by Geller et al. (2015) despite not being in our Gaia-based selection. Of those: (1) white circles show stars with no Gaia DR3 radial velocities available, (2) orange circles are stars that missed our cut in parallax, and (3) yellow circles are those outside our proper motion cut. The star shown in half-orange and half-yellow, WOCS 7004, falls outside of our proper motion and parallax cuts. Figure 2: The three panels show stars considered cluster members according to our Gaia-based selection. Stars in black symbols are also part of the membership sample by Geller et al. (2015), while orange squares show the four red giant stars that were not considered cluster members by Geller et al. (2015). Sanders111Catalogue | WOCS222Catalogue Yadav et al. (2008) | EPIC333Huber et al. (2017) | Gaia DR3444Gaia Collaboration (2022) | In DR3 selection | Geller’15 member | Comment
---|---|---|---|---|---|---
S 1016 | 2004 | 211410817 | 604917594995372544 | No | Yes | No RV available
S 978 | 1008 | 211407537 | 604911375882674560 | No | Yes | No RV available
S 1319 | 5021 | 211420451 | 604923470510324224 | No | Yes | No RV available
S 1557 | 1065 | 211380313 | 604696730596495744 | No | Yes | No RV available
S 1264 | 3004 | 211411629 | 604917835513458688 | No | Yes | Missed cuts in Parallax
S 806 | 5016 | 211417812 | 604969061587592064 | No | Yes | Missed cuts in Parallax
S 1221 | 2014 | 211406541 | 604904503934969856 | No | Yes | Missed cuts in Proper Motion
S 961 | 2017 | 211403555 | 604910001493155584 | No | Yes | Missed cuts in Proper Motion
S 1463 | 4017 | 211413064 | 604918522708195584 | No | Yes | Missed cuts in Proper Motion
S 1000 | 7004 | 211409644 | 604917320117398528 | No | Yes | Missed cuts in PLX + PM
S 258 | 1054 | 211414351 | 604965664269158656 | Yes | No | Red Giant, single star
S 1135 | 2059 | 211443624 | 605015309794935552 | Yes | No | Red Giant, single star
S 247 | 5059 | 211406144 | 598959032246073216 | Yes | No | Red Giant, binary star
S 2000 | - | 211384259 | 604689308892956416 | Yes | No | Red Giant, No WOCS ID
1Sanders (1977), 2Yadav et al. (2008), 3Huber et al. (2017), 4Gaia
Collaboration (2022)
Table 1: The 10 stars from Figure 1 and the 4 stars from Figure 2, along with
their commonly used names.

###  2.2 Observed Colour-Magnitude Diagram

####  2.2.1 Photometry

We compared M67 photometry from Gaia DR2 (Gaia Collaboration et al., 2018a)
and DR3 (Gaia Collaboration, 2022) and found a similar scatter in the colour-
magnitude sequence of the two systems. Therefore, either system could be used
for isochrone fitting without precision loss. We decided to use Gaia DR2
because the deconvolution photometry of the eclipsing binary WOCS 11028
(Sandquist et al., 2021) at the cluster’s turnoff is not directly available in
Gaia DR3 magnitudes, and we want to use the binary system during the isochrone
fitting (section 3).

####  2.2.2 Isochrone fitting sample

We exclude from our isochrone-fitting sample those stars identified as binary,
likely binary, stragglers, or possibly contaminated by a nearby source
according to Geller et al. (2015), and stars from the Gaia DR3 non-single star
catalogue Gaia Collaboration et al. (2022). Following Khan et al. (2019), we
also remove stars with a renormalised unit weight error (RUWE) that exceeds
1.2 in DR2, all of which are likely to fall far from the tight single-star
cluster sequence in the colour-magnitude diagram. This selection results in a
sample of 369 stars (small grey and black dots in Figure 3). The confidence
level that this sample is free of multiple star systems mainly depends on the
accuracy of the radial velocities from Geller et al. (2015), which degrade
toward stars fainter than $12^{\mathrm{th}}$ magnitude in $V$ (see their
Figure 4). To take care of the remaining stars that clearly belong to the
binary main sequence by their position in the colour-magnitude diagram, we
made a further manual cut that brings the isochrone fitting sample to the 309
stars shown in black in Figure 3.

####  2.2.3 Reddening

For differential de-reddening of the photometry, we sought the latest version
of the Bayestar 3D dust maps (Green et al., 2019). However, we found that
these maps lacked resolution and did not behave smoothly in the direction of
M67 (See Appendix C). Thus, we performed de-reddening on a star-by-star basis
using 2D dust maps by Planck Collaboration et al. (2016), and calculated
individual extinction coefficients using the formulae from Danielski et al.
(2018) with polynomial coefficients as calculated by Gaia Collaboration et al.
(2018b) assuming $R_{V}=3.1$. The averages of the full membership sample of
488 stars are $E(BP-RP)=0.059\pm 0.007$ and $A_{G}=0.104\pm 0.012$ magnitudes.
Assuming $E(B-V)=0.041$ as in Taylor (2007), our values are in agreement with
Stassun et al. (2019) who found the relations $E(BP-RP)=1.31\cdot E(B-V)$ and
$A_{G}=2.72\cdot E(B-V)$. Even in a low reddening cluster like M67, applying
individual corrections for dereddening and extinction reduced the scatter of
the colour-magnitude sequence. This was most noticeable around the hook and
the lower red giant branch, but also along the main sequence (compare Figures
3a and b).

Figure 3:  Black symbols represent stars in the isochrone fitting sample, grey
‘x’ symbols show presumed binaries or presumed products on non-standard
evolution, and grey points show stars we later removed from the isochrone-
fitting sample (see Section 2.2.2) (a) M67 colour-magnitude diagram (absolute
magnitudes) after uniform reddening (2.72 x 0.041 G-magnitudes) and extinction
(0.06 BP-RP magnitudes) corrections. (b) Applying differential corrections
based on dustmaps and stellar coordinates.

.

####  2.2.4 Distance

We derived the distance modulus from our sample using Gaia DR3 parallaxes,
which are shown to have increased precision by 30% over DR2 parallaxes (Gaia
Collaboration et al., 2021), and we applied zero point offsets calculated for
each star following Lindegren et al. (2021), who found that zero-points depend
on stellar magnitude, colour, and ecliptic latitude. The mean zero-point
corrected parallax derived this way is $\bar{\omega}=1.195\pm 0.027$ mas
corresponding to a distance modulus $=9.614\pm 0.049$ magnitudes. This is a
smaller distance modulus than the values of 9.726 and 9.730, derived by Gaia
Collaboration et al. (2018b) and Choi et al. (2018), respectively, and also
smaller than 9.690 adopted by Nguyen et al. (2022), but close to 9.630 by
Sandquist et al. (2021). These four studies base their distance modulus
calculation on Gaia DR2 parallaxes, but the latter two apply a parallax zero-
point offset. Nguyen et al. (2022) correct their value applying directly the
offset of $-30\mu\mathrm{as}$ found by Lindegren et al. (2018), while
Sandquist et al. (2021) estimate an offset of $-54\mu\mathrm{as}$ considering
the results by Stassun & Torres (2018), Zinn et al. (2019), and Schönrich et
al. (2019) in addition to Lindegren et al. (2018).

####  2.2.5 Metallicity

Figure 4: Metallicities and $\log(g)$ by APOGEE DR17 and GALAH DR3 for stars
overlapping with our 488 cluster members. The orange lines show fits to the
metallicities of the stars with $\log(g)<3.0$ for both APOGEE (solid line) and
GALAH (dashed line). In the latter, two extreme outliers at low $\log(g)$ and
low [Fe/H] were excluded from the fit. All but one of the stars in the GALAH
sample are also in the APOGEE sample

.

M67 has been observed extensively by spectroscopic surveys and it has even
been used as a calibration cluster for APOGEE (Spoo et al., 2022). Naturally,
we turned to this survey for updated M67 metallicity values. Figure 4 shows
[Fe/H] from APOGEE DR17 for our cluster members (black dots). We noticed a
decreasing metallicity trend with decreasing $\log(g)$, affecting the sample’s
most evolved giants. This was also observed in other open clusters from DR17
(See Appendix D). For comparison, we looked at GALAH DR3 results of M67 (grey
dots) (Buder et al., 2021). Even though the GALAH data has higher dispersion
and uncertainties, it also showed a descending trend with $\log(g)$, albeit
more pronounced. This common trend is likely systematics from not using a full
3D-NLTE approach in the spectroscopic analysis, affecting stars with
$\log(g)\approx 3$-2.5 or lower. No physical mechanism intrinsic to the stars
can explain such metallicity decline in stars undergoing convective mixing in
the atmosphere. This issue could lead to an error in [Fe/H] of up to 1 dex in
some cases (for example, NGC 2632 in Figure 15). Therefore, we caution against
relying on giant star metallicities from non-3D-NLTE spectroscopic results,
and we turn to recent M67 literature to find an appropriate [Fe/H] value for
our models.

The average metallicity across all SDSS DR12 (Alam et al., 2015) M67
observations was reported as $\mathrm{[Fe/H]}=0.08\pm 0.03$ by Stello et al.
(2016), with solar reference from Asplund et al. (2005). Two recent studies
have used high-resolution spectroscopic observations to examine the chemical
abundance trends of stars along the cluster isochrone to investigate the
effects of atomic diffusion. Liu et al. (2019) used Keck/HIRES spectra to
obtain high-precision atmospheric parameters and elemental abundances, and
Souto et al. (2019) based their work on APOGEE DR14 data (Holtzman et al.,
2018). They both found that stellar models needed to have a [Fe/H] between 0.0
and 0.1 to match the detailed chemical abundance patterns, and both used MIST
tracks with diffusion (Choi et al., 2016), and the solar mixture by Asplund et
al. (2009) (the AGSS09 scale). Souto et al. (2019) found no significant
abundance difference in M67 stars of the same evolutionary class, but found
abundance differences between stars in different evolutionary stages of up to
$\sim 0.5\mathrm{dex}$. They found that the observed abundance trends can be
explained by diffusion processes, and estimated a reduction in the efficiency
of gravitational settling of 15%.

##  3 Isochrone Models

We produce custom isochrones tailored to M67 using MESA version 23.05.1
(Paxton et al., 2011, 2013, 2015, 2018, 2019; Jermyn et al., 2023). To
transform the theoretical isochrones into magnitudes that allow for direct
comparisons with observations, we use bolometric correction tables from the
MIST project555https://waps.cfa.harvard.edu/MIST/model_grids.html, which are
based on the grid of stellar atmospheres and synthetic spectra described by
Choi et al. (2016). Our approach for selecting model parameters to match M67,
resulting in our final preferred isochrone, named Isochrone A666Isochrone A
and Inlists available on Zenodo https://zenodo.org/records/12616441, is
described in the following. Table 2 summarises Isochrone A’s parameters.

Parameter | Isochrone A
---|---
Age | 3.95 Gyrs
Distance modulus | 9.614
Nuclear reaction network | h1, he3, he4, c12, n14, o16, ne20, mg24
| (MESA’s basic.net)
Solar mixture | AGSS09 (Asplund et al., 2009)
Z | 0.016
Y | 0.267
$\Delta$Z/$\Delta$Y | 1.2
EoS | HELM (Timmes & Swesty, 2000)
| \+ Skye (Jermyn et al., 2021)
| \+ FreeEOS (Irwin, 2004)
| \+ OPAL (Rogers & Nayfonov, 2002)
| \+ SCVH (Saumon et al., 1995)
Interior opacities | OPAL (Iglesias & Rogers, 1993, 1996)
Opacities | AESOPUS (Marigo & Aringer, 2009)
Atmospheres | T($\tau$), varying, Trampedach solar
| (Trampedach et al., 2014a; Ball, 2021)
Overshooting | Exponential, increasing with mass.
$\alpha_{\mathrm{MLT}}$ | Varying with $T_{\mathrm{eff}}$ and $\log(g)$,
| follows 3D grid (Trampedach et al., 2014b)
Scaling factor $k_{\alpha_{\mathrm{MLT}}}$ | 1.11
Mass loss | No
MLT_option | Henyey
Convection | Schwarzschild criterion
Diffusion | Yes

Table 2: Summary of parameters in MESA models for the isochrone shown in
Figure 7. The isochrone is presented in 5

.

###  3.1 Initial chemical composition

Guided by the results detailed in 2.2.5 we decided to adopt the AGSS09 scale,
$\mathrm{[Fe/H]}=0.05$, and $Z=0.016$, implying a helium enrichment law of
$\Delta Y/\Delta Z=1.2$ when primordial helium is taken as
$Y_{\mathrm{p}}=0.248$. While helium is not solar-scaled, similar values of
$\Delta Y/\Delta Z$ have been used to match models to data in other open
clusters: 1.4 was found for NGC 6791 (Brogaard et al., 2012), 1.2 for the
Hyades (Brogaard et al., 2021), and an average of 1.25 for NGC 6811 (Sandquist
et al., 2016). Souto et al. (2019) predicted [Fe/H] $\sim 0.05$ (our chosen
primordial metallicity) for M67 solar analogues using their estimated
diffusion factor, which is consistent with our models if we consider our lower
than solar initial helium, and their solar scaled models.

Our choice of $Z$ achieves good colour agreement with the main sequence for
ages within the range quoted in literature: 3.5 - 4.3 Gyrs (VandenBerg &
Stetson, 2004; Sarajedini et al., 2009; Barnes et al., 2016), while the chosen
helium content, $Y=0.267$, results in slightly more massive stars for a given
luminosity compared to solar-scaled models. This becomes relevant when taking
into account the eclipsing binary WOCS 11028 (Section 3.6).

###  3.2 Model atmospheres

One-dimensional stellar evolution codes like MESA include simplifications in
their representation of the outer layers of a star, such as using a single and
fixed value of the convective mixing length $\alpha_{\mathrm{MLT}}$, which is
commonly anchored to the Sun. In our models, we seek to incorporate more
realistic values of $\alpha_{\mathrm{MLT}}$. This, we take from interpolations
to the grid of hydrodynamical 3D simulations of stellar convection – the 3D
grid – by Trampedach et al. (2014b) evaluated at every step according to
$\log(g)$ and $T_{\mathrm{eff}}$ with a scaling factor
$k_{\alpha_{\mathrm{MLT}}}$, following Mosumgaard et al. (2018), except that
we calibrate $k_{\alpha_{\mathrm{MLT}}}$ to M67 rather than to the Sun
(Section 3.5.1).

To ensure the correct photospheric transition to the outer part of the models,
we follow the formulation of generalised Hop-functions for atmospheres by
Trampedach et al. (2014a) from 3D simulations of convection in deep stellar
atmospheres: ($T(\tau)$=Trampedach_solar) as implemented by Ball (2021), and
we modify the radiative gradient so that regions of low optical depth have a
temperature that follows this $T(\tau)$ relation. The results of applying the
improvements described above to a base form of Isochrone A can be observed in
Appendix B, by comparing Figures 13a and b.

While studying the effects of varying model parameters, we saw that the
selection of low-temperature opacities had a significant impact on the
temperatures and luminosities of the models. As a result, we have opted to use
low-temperature opacities from the AESOPUS (Accurate Equation of State and
Opacity Utility Software) tables because the AESOPUS code (Marigo & Aringer,
2009) accounts for the processing of CNO in giant stars while the MESA default
tables FA05 (Ferguson et al., 2005) do not. The effect of additionally
changing from FA05 to AESOPUS opacity tables can be observed in Appendix B, by
comparing Figures 13b and c.

###  3.3 Overshooting

We implement exponential diffusive overshoot (Herwig, 2000) at the Hydrogen
core with a strength set to increase gradually with mass, following the
expression proposed by VandenBerg et al. (2006):

$f_{\mathrm{ov}}=f\left[1-\mathrm{cos}\left(\pi\frac{M_{\star}-M_{\mathrm{full\\_off}}}{M_{\mathrm{full\\_on}}-M_{\mathrm{full\\_off}}}\right)\right]$
(1)

where $M_{*}$ is the mass along the isochrone,
$M_{\mathrm{full\\_off}}=1.10\,\mathrm{M}_{\odot}$,
$M_{\mathrm{full\\_on}}=1.80\,\mathrm{M}_{\odot}$, which provides a smooth
transition between 0 and $f_{\mathrm{ov}}$. We set the parameter $f_{0}$ as
$f_{\mathrm{ov}}/2$, where $f_{0}$ is a distance from the edge of the
convective zone where the diffusion mixing coefficient is measured. This
coefficient is then multiplied by $e^{-2z/f_{\mathrm{ov}}H\\!p}$, where
$H\\!p$ is the pressure scale height, so that the mixing efficiency declines
exponentially into the radiative zone with distance $z$ from the convective
boundary.

While it is common to define overshoot in terms of $H\\!p$, this poses a
problem when $f_{\mathrm{ov}}$ is kept fixed for all masses because
$H\\!p\\!\rightarrow\\!\infty$ when the convective core radius
$\rightarrow\\!0$. By default, MESA limits $f_{\mathrm{ov}}H\\!p$ to the size
of the convective zone, however, this could still result in unphysically large
convective cores in low mass stars, as, according to Roxburgh (1992)
overshooting should not expand the convective core more than 18% of its ‘no-
overshoot’ size. A mass-dependent overshoot like the one from Equation (1)
helps overcome this problem. Another common criterion that avoids the issue is
to express the overshoot extent as a fraction of the core radius. The effect
of the change from a fixed overshoot to our adopted mass-dependent values can
be seen in Appendix B, by comparing the black and green curves in Figure 13c.

###  3.4 Mass Loss

We do not have a complete understanding of the mass-loss processes operating
during the red giant phase (Karakas, 2017). Choi et al. (2016) selected a
Reimers parameter $\eta=0.1$ for use in their MIST models based on the
initial-final mass relation in the Magellanic Clouds. Miglio et al. (2012)
found from asteroseismology that the mass difference between red giant branch
and red clump stars of the metal-rich open cluster NGC 6791 can be described
with a Reimers parameter in the range $0.1<\eta<0.3$, lower than predicted by
the commonly used prescription by Reimers (1975). They also found that the
near-solar-metallicity open cluster NGC 6819 is compatible with no mass loss.
Stello et al. (2016) measured the masses of M67 giants using asteroseismology
and found no significant mass difference between the red clump and the red
giant branch. Considering all these results, and especially the evidence from
M67 asteroseismology, we decided to run models without mass loss. Regardless,
mass loss is unlikely to significantly affect our M67 isochrone fit because
nearly all mass loss occurs at the very end of the red giant branch, and red
clump stars will be excluded from the final fit (see Section 4).

###  3.5 Empirical Calibrations

Bringing together our model prescriptions as described above and the
extinction-corrected M67 colour-magnitude diagram from Section 2.2, we can now
fine-tune the mixing length and overshooting.

Figure 5: Isochrone A and M67 presumed single stars in absolute magnitude and
colour-coded by (a) $\alpha_{\mathrm{MLT}}$, (b) $\log(g)$, (c) overshooting
parameter $f_{\mathrm{ov}}$, and (d) stellar mass. The dashed line in b)
crosses the isochrone where the $\alpha_{\mathrm{MLT}}$ 3D grid stops at
$\log(g)$ 2.4. After that point in the evolution up the RGB, we set our models
to maintain the last measured $\alpha_{\mathrm{MLT}}$ value and after helium
ignition $\alpha_{\mathrm{MLT}}$ is set to 1.93.

####  3.5.1 Mixing Length Parameter

Figure 5a shows Isochrone A colour-coded by our adopted mixing length,
resulting from applying the scaling factor $k_{\alpha_{\mathrm{MLT}}}$= 1.11
to the 3D grid. Because the 3D grid stops below $\log(g)\sim$2.4 –segmented
line in Figure 5b–, we keep the last interpolated grid value for the rest of
the evolution up the giant branch, but doing so puts the isochrone’s red clump
at lower temperatures than observed from the cluster’s clump stars. To resolve
this issue, we decided to set $\alpha_{\mathrm{MLT}}$ to 1.93 after the helium
flash, which makes Isochrone A agree with the observations of the red clump.
However, this does not necessarily mean that 1.93 is the correct mixing length
value for stars with $T_{\mathrm{eff}}$ and $\log(g)$ similar to the M67 clump
stars. It merely expresses our ignorance of the behaviour of
$\alpha_{\mathrm{MLT}}$ after the 3D grid stops and also compensates for other
shortcomings that our helium core burning models might have that would affect
$T_{\mathrm{eff}}$.

####  3.5.2 Overshooting

The effect of the overshooting parameter is seen primarily near the
isochrone’s turnoff due to the turnoff’s sensitivity to the size of the
convective core and how the size of the convective core depends on stellar
mass. Figures 5c and d show Isochrone A colour-coded by $f_{\mathrm{ov}}$ and
stellar mass, respectively. Figure 6 shows our mass-dependent overshooting
parameter $f_{\mathrm{ov}}$, compared to the reported $f_{\mathrm{ov}}$ by
Claret & Torres (2017, 2018) in orange symbols, where models using exponential
overshooting with diffusion are fitted to eclipsing binary masses.

There have been several recent studies on overshooting, but they often use
different methods that cannot be directly compared to our calibration. For
example, Lindsay et al. (2024) applies ‘step’ overshoot instead of
exponential, and maintains a fixed value of $f_{0}$ while they vary
$f_{\mathrm{ov}}$ and take care of the issue of large cores at low masses
mentioned in Section 3.3 by switching from units of $H\\!p$ to units of core
radius if $H\\!p$ is larger than the core radius. For an in-depth review of
their results in the context of results from various overshooting studies,
including Claret & Torres (2017, 2018), we refer the reader to Lindsay et al.
(2024).

Figure 6: Comparing our mass-dependent overshoot used in Isochrone A with
independent modelling of eclipsing binaries (Claret & Torres, 2017, 2018). The
dark section of the curve indicates the mass range along Isochrone A.

###  3.6 Eclipsing Binary WOCS 11028

The discovery of the eclipsing binary WOCS 11028 (Sandquist et al., 2021),
with one component near the cluster turnoff, provides an opportunity to test
how well stellar models can describe the cluster. Sandquist et al. (2021)
found discrepancies of at least 0.17 mag ($\sim 10\sigma$) between the
absolute magnitude of the primary component of the binary system and
predictions by various stellar evolution codes, or equivalently a lower than
predicted mass by $\delta m=0.05M_{\odot}$, even for isochrones as young as
3.2 Gyr. Nguyen et al. (2022) with their updated PARSEC isochrones report
similar discrepancies.

Figure 7: The de-reddened colour-magnitude diagram of our M67 sample and the 3.95 Gyr old best fitting isochrone from this work, Isochrone A. In small black symbols: stars with RUWE>1.2, stars in multiple star systems, and stellar products of non-standard evolution, as well as suspected binary members manually removed from the sample (end of Section 2.2). Several blue straggler members are bluer than 0.5 and are not shown in the figure. In grey symbols: presumed single members. Inset: the eclipsing binary system WOCS 11028 in large off-white circles resulting from the deconvolution of the binary’s photometry, see Table 3. Dashed lines show the isochrone’s predicted magnitudes for the masses of the two binary components. $\Delta G_{\mathrm{mag}}$ is $-0.05\pm 0.09$, or $-0.67\sigma$ for the secondary component and $+0.21\pm$ 0.06, or $+12\sigma$ for the primary component, where $\sigma$ values correspond to those shown in Table 3. Binary Member | BP-RP∗ | $G^{*}$ | 2$\sigma$(G) | Mass M⊙ | Radius R⊙
---|---|---|---|---|---
WOCS 11028a | 0.686 | 13.336 | 0.035 | 1.222 $\pm$ 0.006 | 1.430 $\pm$ 0.030
WOCS 11028b | 0.916 | 14.871 | 0.150 | 0.909 $\pm$ 0.004 | 0.904 $\pm$ 0.015
Table 3: The eclipsing binary results from Sandquist et al. (2021). ∗BP-RP and
G are from their main sequence fit in Gaia DR2 magnitudes and have been
dereddened and corrected for extinction using $A_{V}$ from dust maps.

Our Isochrone A-predicted G-magnitudes for the eclipsing primary and secondary
(Figure 7, inset) are respectively, $+0.21\pm 0.06$ and $-0.05\pm 0.09$
magnitudes relative to the eclipsing binary estimates in Table 3. This
includes uncertainties from our distance modulus and from the magnitudes and
masses of the eclipsing binary. These offsets are comparable to the ones
obtained by Sandquist et al. (2021) and Nguyen et al. (2022). Like us, Nguyen
et al. (2022) found isochrones that agreed with the secondary within
uncertainties but not with the primary. This is concerning because Sandquist
et al. (2021) reported lower confidence in estimating the secondary star’s
magnitude compared to the primary.

####  3.6.1 Alternative Isochrone EB

To address the eclipsing binary issue, we attempted to anchor an alternative
isochrone on the mass and magnitude of the primary (WOCS 11028a) within
$1\sigma$, as shown in the inset of Figure 8.

Figure 8: Alternative Isochrone EB. Dashed lines in the inset show predicted
magnitudes for the masses of the two binary components. $\Delta
G_{\mathrm{mag}}$ is $0.01\pm 0.06$ ($0.3\sigma$) for the primary component
and $-0.27\pm 0.09$, or $-4\sigma$ for the secondary component, where $\sigma$
values correspond to those shown in Table 3. Symbols follow Figure7.

To match the mass of the primary, we needed an isochrone with a main-sequence
absolute brightness that is 0.2 magnitudes fainter compared to Isochrone A for
a given mass. This can be achieved by making specific parameter adjustments,
such as lowering the helium content, changing the reference element mixture,
reducing isochrone age, increasing the distance modulus, or adding more
reddening. Although all modifications will also move the secondary, applying
one or more of these modifications could reconcile discrepancies with the
primary, but it may also result in inconsistencies in other parts of the
isochrone. For instance, a lower isochrone age leads to a subgiant branch
morphology that does not match the M67 observed data, but that, in turn, could
be resolved with increased shell overshoot values. Considering the above, we
developed a few alternative isochrones that could be anchored on WOCS 11028a.
None of them, however, could match both binary components simultaneously, and
all of them violated at least one other constraint. We consider that some
constraints are more reliably established than others. The cluster’s reddening
has been consistently measured at $E(B-V)=0.04\pm 0.04$ magnitudes, or its
equivalent in other bands (Nissen et al., 1987; Schlegel et al., 1998; Taylor,
2007), which agrees with our calculations from Planck Collaboration et al.
(2016) dust maps. Additionally, the helium content should not be much lower
than the Sun’s, given the near-solar age and metallicity of the cluster, and
it should certainly be higher than the primordial helium abundance. In
contrast, historical discrepancies in the parallax measurements for M67 exist
(see Section 2.2.4), and convective overshoot is known to be a significant
source of uncertainty in stellar models (Kippenhahn et al., 2013; Salaris &
Cassisi, 2017; Viani & Basu, 2020). Therefore, we present here our best-fit
alternative model, Isochrone EB, which varies these less-established
constraints. We achieve a good match between Isochrone EB and the data by
using: (a) the solar mixture from Grevesse & Sauval (1998) with initial
$Z=0.0196$, $Y=0.2710$ ($\Delta Y/\Delta Z=1.18$), (b) a distance modulus of
9.73, $2.4\sigma$ larger than our derived value from zero-point corrected Gaia
DR3 parallaxes, (c) and a sharp and significant increase in overshoot at and
beyond $1.32\,\mathrm{M}_{\odot}$, with manually adjusted values of
$\alpha_{\mathrm{MLT}}$. The model parameters of Isochrone EB are summarised
in Table 4, and a comparison of its overshoot and mixing length to those of
Isocrone A are shown in Figures 9a and b, respectively.

Parameter | Isochrone EB
---|---
Age | 3.55 Gyrs
Distance modulus | 9.73
Nuclear reaction network | h1, he3, he4, c12, n14, o16, ne20, mg24
| (MESA’s basic.net)
Solar mixture | GS98 (Grevesse & Sauval, 1998)
Z | 0.0196
Y | 0.2710
$\Delta$Z/$\Delta$Y | 1.18
EoS | HELM (Timmes & Swesty, 2000)
| \+ Skye (Jermyn et al., 2021)
| \+ FreeEOS (Irwin, 2004)
| \+ OPAL (Rogers & Nayfonov, 2002)
| \+ SCVH (Saumon et al., 1995)
Interior opacities | OPAL (Iglesias & Rogers, 1993, 1996)
Low temperature opacities | FA05 (Ferguson et al., 2005)
Atmospheres |  $T(\tau)$, varying, Eddington
Overshooting | Exponential, increasing with mass,
| sharp increase at M=1.32M⊙
$\alpha_{\mathrm{MLT}}$ | Variable, adjusted manually
MLT_option | Henyey
Mass loss | No mass loss
Convection | Schwarzschild criterion
Diffusion | Yes
Table 4: Summary of parameters in MESA models for Isochrone EB shown in
Figure8. Figure 9: (a) The overshooting parameter $f_{\mathrm{ov}}$ and the
mixing length parameter $\alpha_{\mathrm{MLT}}$ used in Isochrone EB (dashed)
shown along the values used in Isochrone A (solid grey).

We note that in Isochrone EB (Figure 9a), the values of overshooting in stars
more massive than 1.35M⊙ make the models exceed the upper limit given by the
Roxburgh criterion for small convective cores, where overshooting is limited
to a level that would prevent the core from growing to more than 18% of its
‘no-overshoot’ size (Roxburgh, 1992). Furthermore, the study of convective
overshooting remains a very active field, with various stellar signatures
being used to calibrate models in different studies. Still, all recent studies
find either a gradual increase with mass (Lindsay et al., 2024) or no increase
at all (Constantino & Baraffe, 2018).

From Figure 9b, we see some correspondence between the shapes described by
both $\alpha_{\mathrm{MLT}}$ curves if we consider that the masses of
Isochrone EB are close to $0.5\,$M⊙ larger than Isochrone A’s at any given
point in the colour-magnitude or $\log(g)$-$T_{\mathrm{eff}}$ diagrams.
However, the $\alpha_{\mathrm{MLT}}$ values we applied to stellar models of
different masses are not based on any known physics but are solely adapted by
eye to fit the CMD diagram. From all the above, we conclude that even though
Isochrone EB matches the mass and magnitude of WOCS11028a, the isochrone is
unlikely to be a good model for the cluster. Also, preliminary results
indicate that the luminosities of the giants predicted by the alternative
isochrone EB are in tension with asteroseismic observations (Reyes et al., in
preparation).

####  3.6.2 Possible causes of the model-mass conflict

In addition to potential remaining shortcomings in stellar modelling, possible
explanations for the conflict with the eclipsing binary masses include non-
standard evolution of at least one component of WOCS11028, or systematics in
the estimations of the eclipsing binary parameters. However, the presence of
lithium in the spectra of the binary makes it unlikely that these stars have
interacted in the past or resulted from a merger event, as, in such scenarios,
lithium would have been consumed by the original lower mass stars with deeper
surface convection zones, before merging (Sandquist et al., 2021). Possible
grounds for an altered eclipsing binary mass or luminosity would be the non-
detected presence of a third star in the WOCS 11028 system. Such has been the
case in the past in the young Hyades open cluster and the multiple system
HD27130 where a previous mass overestimate due to an unseen third component
led Lebreton et al. (2001) to estimate a helium content of $Y=0.255$, while
more recent studies indicate a higher than solar helium abundance for the
Hyades (see Brandner et al., 2023, and references therein). For a discussion
on the effects of HD27130 on the derived properties of the Hyades see Brogaard
et al. (2021). However, even if there were a faint third star hiding in the
photometry of WOCS11028 due to a particularly pathological setup of the
multiple system, it would mean a similar and small mass shift to both known
components, with no significant change in their relative positions in the
colour-magnitude diagram. This opens up the possibility that models do not
accurately represent the speed at which stars evolve during the main sequence,
at least for stars around $1.2\mathrm{M}_{\odot}$, which could make us
reconsider our understanding of stellar evolution. In summary, the case of
WOCS11028 is still very much open and requires further investigation.

##  4 Age uncertainty through chi-square analysis

Although we followed the traditional qualitative approach of matching the
input physics of Isochrone A to the M67 colour-magnitude diagram, we will in
the following look to estimate the age uncertainty using a $\chi^{2}$
approach. We take our isochrone-fitting sample from section 2.2 and further
remove the stars that have evolved beyond the point where the red giant branch
leaves the $\log(g)$-$T_{\mathrm{eff}}$ area covered by the 3D grid, marked by
the segmented line in Figure 5b.

We vary the model age between 3.5 and 4.3 Gyrs, covering the age range quoted
in literature (VandenBerg & Stetson, 2004; Sarajedini et al., 2009; Barnes et
al., 2016), and the distance modulus within 2$\sigma$ of our Gaia DR3
parallax-derived value of $9.614\pm 0.049$. Figure 10 shows the extremes of
varying age and distance modulus like this.

Figure 10: Isochrones of 3.50 Gyrs (solid) and 4.3 Gyrs (dashed) using the
same physics as Isochrone A. (a) M67 isochrone fitting sample with a distance
modulus is 9.714. (b) same as panel a, but for a distance modulus of 9.514.

For each combination of age and distance modulus, we used the nearest matching
point on the isochrone as the expected value ($\vec{x}_{\mathrm{exp},i}$) for
each star ($\vec{x}_{\mathrm{obs},i}$) in a $\chi^{2}$-like goodness of fit
analysis. Following Sandquist et al. (2021), we define a $\chi^{2}$-like
parameter involving fractional colour and magnitude differences, of the form:

$\displaystyle\frac{1}{N}\sum_{i=0}^{N}\frac{(\vec{x}_{\mathrm{obs},i}-\vec{x}_{\mathrm{exp},i})^{2}}{\vec{x}_{\mathrm{exp},i}^{2}},\;\;\vec{x}=(\mathrm{col}_{i},\mathrm{mag}_{i})$
(2)

that we later rescale so that the lowest $\chi^{2}$ value equals 1. We do not
include a $\sigma$ in this expression because Gaia does not provide magnitude
uncertainties. Figure 11 shows the age versus distance modulus space colour-
coded by this (rescaled) parameter. Figure 11a presents the combined
statistics of both main-sequence and giant stars, which we will use in the
following analysis. We find that the lowest $\chi^{2}$ (black triangle) is at
a distance modulus about 1$\sigma$ below the parallax-based value (black
circle), and with a corresponding age of $4.04\,\mathrm{Gyrs}$. If this were
to be caused by any systematics in Gaia DR3 parallaxes, it would correspond to
about $25\,\mu\mathrm{as}$. However, there is no indication that such
systematics exist beyond the Lindegren et al. (2021) corrections (Zinn, 2021),
which we have already included. Furthermore, a smaller distance modulus would
increase discrepancies between models and WOCS 11028. We, therefore, prefer to
quote the age as $3.95\,\mathrm{Gyrs}$ –the black circle in Figure 11a – and
our confidence interval as $[3.80,4.11]\,\mathrm{Gyrs}$, the point where the
$\chi^{2}$ has risen by 1 from its minimum value along the adopted parallax-
based distance modulus (black error bars). For completeness, we also present
the $\chi^{2}$ results restricted to main sequence stars in Figure 11b, and
giants in Figure 11c, where the split into those groups is defined by core
hydrogen exhaustion (the hottest point along the isochrone). For reference, we
copy in the latter panels with grey symbols the location of the lowest
$\chi^{2}$ points from panel a.

Figure 11: $\chi^{2}$-like statistics from M67 main sequence and giant stars
for isochrone ages between 3.5 and 4.3 Gyrs and distance modulus within
2$\sigma$ of our Gaia parallax-derived value of 9.614 $\pm$ 0.049. The
triangle shows the overall lowest $\chi^{2}$ for panel a, and the circle,
including errorbars, shows our quoted age of $3.95^{+0.16}_{-0.15}$ Gyrs,
corresponding to the lowest $\chi^{2}$ at distance modulus = 9.614. The
position of the circle and the triangle are copied in the other panels for
completeness.

##  5 Conclusions and Future Work

We have presented a detailed isochrone, Isochrone A, specifically designed for
the low-reddening cluster M67 using state-of-the-art stellar input physics.
Our models match the morphology of the tight cluster sequence better than
previously reported in the literature in Gaia colours and magnitudes. Our
results indicate that adopting a mixing length $\alpha_{\mathrm{MLT}}$ that is
dependent on $T_{\mathrm{eff}}$ and $\log(g)$, such as the one obtained from
the 3D-grid, can improve the fit of models to observations. Therefore, this
approach should be considered in future standard stellar models. Our best age
estimate is $3.95^{+0.16}_{-0.15}\,\mathrm{Gyrs}$, in agreement with the
estimates from VandenBerg & Stetson (2004) ($4.0\pm 0.4\,\mathrm{Gyrs}$ from
$BV$ photometry), Sarajedini et al. (2009) ($3.5-4.0\,\mathrm{Gyrs}$ from
2MASS photometry), (Barnes et al., 2016) ($4.2\pm 0.2\,\mathrm{Gyrs}$ from
gyrochronology), Bressan et al. (2012) ($3.7\,\mathrm{Gyrs}$ from $BV$
photometry using PARSEC isochrones v1.2), and Nguyen et al. (2022)
($3.9\,\mathrm{Gyrs}$ from Gaia DR2 photometry using PARSEC isochrones v2.0).

Despite the improvements in the match to the cluster colour-magnitude diagram,
our isochrone shows some disagreement with the measured mass for the primary
component of the eclipsing binary WOCS11028 near the cluster turnoff. The only
alternative isochrone we could anchor on the primary eclipsing binary is based
on improbable input physics and fundamental cluster properties. Therefore, we
conclude that the literature values of the masses and luminosities for the
eclipsing binary system are incompatible with the rest of the cluster in a
scenario where both binary components have evolved independently. Of course,
this is under the perhaps unrealistic assumption that there are no
shortcomings in our current model input physics.

To help resolve the issues surrounding WOCS 11028, one could potentially use
detailed asteroseismic modelling to confirm masses and radii of the cluster
stars. In addition, high-resolution imaging of the binary might be able to
detect any undiscovered component that could affect the mass and radius
estimates.

## Software

Extinction correction performed using the dustmaps python tool (Green, 2018).
https://dustmaps.readthedocs.io/en/latest/

Interpolation to bolometric correction tables was performed using the python
package isochrones v2.1 (Morton, 2015).
https://isochrones.readthedocs.io/en/latest/

Isochrones are constructed using the fortran code iso described in Dotter
(2016a, b). https://cmasher.readthedocs.io/user/introduction.html

## Acknowledgements

D.S. is supported by the Australian Research Council (DP190100666). This work
has made use of data from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/Gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/Gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

This research includes computations using the computational cluster Katana
supported by Research Technology Services at UNSW Sydney,
https://doi.org/10.26190/669x-a286.

We would like to thank the participants of the MIAPbP Program "Stellar
Astrophysics" in Garching August 2023 for the discussion on artificial
metallicity trends with $\log(g)$ when not using 3D-NLTE spectroscopic
analyses.

## Data Availability

The Isochrone and the inlists for the models presented in this work can be
found on the Zenodo online software platform,
https://zenodo.org/records/12616441.

##  Appendix A Isochrone A

Our final 3.95 Gyr isochrone, Isochrone A, is partially presented in Table 5,
and the full isochrone from main sequence to core helium burning is available
for download through the online version of the paper. The column names
correspond to standard MESA parameters, except mixing_length_alpha and
f_overshoot, which correspond to $\alpha_{\mathrm{MLT}}$ and
$f_{\mathrm{ov}}$, respectively, as denoted in the text.

isochrone_age_yr | star_mass | log_L | log_Teff | log_R | log_g | … | mixing_length_alpha | f_overshoot
---|---|---|---|---|---|---|---|---
3.9500e9 | 0.8506 | -0.4003 | 3.7133 | -0.1038 | 4.5755 | … | 1.9849 | 0.0000
3.9500e9 | 0.8591 | -0.3778 | 3.7165 | -0.0989 | 4.5700 | … | 1.9774 | 0.0000
3.9500e9 | 0.8676 | -0.3557 | 3.7195 | -0.0940 | 4.5644 | … | 1.9691 | 0.0000
3.9500e9 | 0.8759 | -0.3339 | 3.7226 | -0.0891 | 4.5587 | … | 1.9637 | 0.0000
3.9500e9 | 0.8841 | -0.3127 | 3.7255 | -0.0843 | 4.5532 | … | 1.9602 | 0.0000
Table 5: First few rows of selected columns of Isochrone A. The full
isochrone, including key internal structure parameters and photometric bands,
is available for download from the online version of the paper, or on Zenodo
(See Data Availability.)

##  Appendix B Other Isochrones

Figure 12 shows isochrones from two recent works that have been fitted to the
M67 colour-magnitude diagram in Gaia colours. These are models designed for
general use and, therefore, not specifically for M67. Nguyen et al. (2022) use
The PAdova and tRieste Stellar Evolutionary Code v2.0, and similarly to
Isochrone A from this work, incorporates a mass-dependent overshoot and the
AESOPUS opacity tables. Their best fit to M67 is shown with our stellar
fitting sample in Figure 12a. Choi et al. (2018) -the MIST models- use the
MESA evolutionary code, and a fixed core overshoot $f_{\mathrm{ov}}=0.016$.
Figure 12b shows their fit to M67. In each panel, the fitting sample is
corrected to absolute magnitudes using the distance modulus adopted by the
corresponding authors.

Figure 12: Best fit isochrones to the M67 Gaia colour magnitude diagram from
a) Nguyen et al. (2022), and b) Choi et al. (2018).

Figure 13 shows the qualitative effects of our atmosphere and core-overshoot
choices from sections 3.2 and 3.5.2 in our models. In all panels, Isochrone A
is shown in green and: Panel (a) shows a base-model from where we start:
similar to the MIST models, but with our chosen initial composition (2.2.5),
age, and the carefully chosen (fixed) values of overshoot $f_{\mathrm{ov}}$
and $\alpha_{\mathrm{MLT}}$ to best fit M67; Column (a) in the table just
below the figure indicates five key parameter values. We see that the
simplified model can still be made to fit the red giant branch well with
$\alpha_{\mathrm{MLT}}$=1.91, but the model is too hot compared with the main-
sequence. In (b) we show the base model from (a), but now applying the
3D-grid-based varying $\alpha_{\mathrm{MLT}}$ and associated atmospheric
improvements from Trampedach et al. (2014a, b) (see column (b) in table).
These changes move the isochrone towards colder temperatures, making the model
fit the main sequence and subgiants better, however the red giant branch is
now "too cold", at least with the age and [Fe/H] of Isochrone A. In (c) we now
also add low opacity tables that account for CNO processes in red giants (see
column (c)). This brings the model back into agreement with the red giants.
The only remaining simplified parameter in (c) is the fixed overshoot
$f_{\mathrm{ov}}$, which has a subtle, but noticeable, effect on the
morphology of the main sequence turnoff, where the even smaller values of
$f_{\mathrm{ov}}$ used in Isochrone A at these masses (Figure 5c, d) appear to
fit the stellar sequence marginally better.

Parameter | (a) | (b) | (c) | Isochrone A
---|---|---|---|---
$\alpha_{\mathrm{MLT}}$ | 1.91 | 3D grid | 3D grid | 3D grid
T$(\tau)$ | Eddington | Trampedach solar | Trampedach solar | Trampedach solar
photospheric transition | default | improved | improved | improved
low T opacities | FA05 | FA05 | AESOPUS | AESOPUS
core overshoot | fixed $f_{\mathrm{ov}=0.004}$ | fixed $f_{\mathrm{ov}=0.004}$ | fixed $f_{\mathrm{ov}=0.004}$ | mass-dependent

Figure 13: In b) and c), the isochrones in black show the effects of
cumulatively applying parameter ‘improvements’ to the simplified form of
Isochrone A shown in (a). The table below the figure shows the model
parameters we vary and their values, along with the values used in Isochrone
A.

##  Appendix C Extinction Maps

We present extinction maps of the region around the M67 cluster, centred at
$\mathrm{RA}=215.69^{\circ}$, $\mathrm{Dec}=+31.92^{\circ}$ (galactic
coordinates), that we considered for differential de-reddening of M67
photometry, aiming to reduce the scatter of the cluster sequence in the
colour-magnitude diagram. We compared the 2D extinction map by Planck
Collaboration et al. (2016) with the Bayestar 3D reddening map (Green et al.,
2019) and used the coefficients from Table 6 from Schlafly & Finkbeiner (2011)
to convert from Bayestar reddening to $A_{V}$ extinctions. We found that the
2D map results in a smoother map of extinctions that works well to reduce the
scatter in the M67 colour-magnitude diagram (Figure 3), while the 3D Bayestar
map does not behave well, at least in this area of the sky, and at the small
scales we require. The 3D extinction map looks pixelated and shows sharp
variations, as seen in Figure 14. This could be related to the systematic
trends in reddening identified by Green et al. (2019) at low reddenings. This
is a note of caution for future cluster studies.

Figure 14: Extinction maps in galactic coordinates of the M67 cluster region
(left) and of our M67 stellar sample (right). The colour expresses the
extinction $A_{V}$ at each point of the sky.

##  Appendix D Metallicity trend in Red Giants from APOGEE DR17

The trend of decreasing metallicity with decreasing $\log(g)$ observed in M67
was observed in almost all of the clusters included in the APOGREE DR17 Open
Cluster Chemical Abundances and Mapping Catalogue
777https://www.sdss4.org/dr17/data_access/value-added-catalogs/?vac_id=open-
cluster-chemical-abundances-and-mapping-catalog. In Figure 15 we show the
clusters in the catalogue with more than 80 members, and fit the metallicity
trend for $\log(g)<3.0$ in each cluster.

Figure 15: Clusters from the APOGEE DR17 Open Cluster Chemical Abundances and
Mapping Catalogue with more than 80 cluster members (alphabetically ordered).
The orange lines show a fit to the metallicities of all stars with $\log(g)$<=
3.0 within each cluster.

## References

  * Alam et al. (2015) Alam S., et al., 2015, ApJS, 219, 12
  * Asplund et al. (2005) Asplund M., Grevesse N., Guedel M., Sauval A. J., 2005, arXiv e-prints, pp astro–ph/0510377
  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
  * Ball (2021) Ball W. H., 2021, Research Notes of the American Astronomical Society, 5, 7
  * Barker & Paust (2018) Barker H., Paust N. E. Q., 2018, PASP, 130, 034204
  * Barnes et al. (2016) Barnes S. A., Weingrill J., Fritzewski D., Strassmeier K. G., Platais I., 2016, ApJ, 823, 16
  * Beeson et al. (2024) Beeson K. L., et al., 2024, MNRAS, 529, 2483
  * Bertelli Motta et al. (2018) Bertelli Motta C., et al., 2018, MNRAS, 478, 425
  * Brandner et al. (2023) Brandner W., Calissendorff P., Kopytova T., 2023, MNRAS, 518, 662
  * Bressan et al. (2012) Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127
  * Brogaard et al. (2012) Brogaard K., et al., 2012, A&A, 543, A106
  * Brogaard et al. (2021) Brogaard K., et al., 2021, A&A, 645, A25
  * Buder et al. (2021) Buder S., et al., 2021, MNRAS, 506, 150
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Choi et al. (2018) Choi J., Conroy C., Ting Y.-S., Cargile P. A., Dotter A., Johnson B. D., 2018, ApJ, 863, 65
  * Claret & Torres (2017) Claret A., Torres G., 2017, ApJ, 849, 18
  * Claret & Torres (2018) Claret A., Torres G., 2018, ApJ, 859, 100
  * Constantino & Baraffe (2018) Constantino T., Baraffe I., 2018, A&A, 618, A177
  * Danielski et al. (2018) Danielski C., Babusiaux C., Ruiz-Dern L., Sartoretti P., Arenou F., 2018, A&A, 614, A19
  * Dotter (2016a) Dotter A., 2016a, ISO: Isochrone construction, Astrophysics Source Code Library, record ascl:1601.021 (ascl:1601.021)
  * Dotter (2016b) Dotter A., 2016b, ApJS, 222, 8
  * Ferguson et al. (2005) Ferguson J. W., Alexander D. R., Allard F., Barman T., Bodnarik J. G., Hauschildt P. H., Heffner-Wong A., Tamanai A., 2005, ApJ, 623, 585
  * Gaia Collaboration (2022) Gaia Collaboration 2022, VizieR Online Data Catalog, p. I/355
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2018a) Gaia Collaboration et al., 2018a, A&A, 616, A1
  * Gaia Collaboration et al. (2018b) Gaia Collaboration et al., 2018b, A&A, 616, A10
  * Gaia Collaboration et al. (2021) Gaia Collaboration et al., 2021, A&A, 649, A1
  * Gaia Collaboration et al. (2022) Gaia Collaboration et al., 2022, arXiv e-prints, p. arXiv:2206.05595
  * Geller et al. (2015) Geller A. M., Latham D. W., Mathieu R. D., 2015, AJ, 150, 97
  * Geller et al. (2021) Geller A. M., Mathieu R. D., Latham D. W., Pollack M., Torres G., Leiner E. M., 2021, AJ, 161, 190
  * Gilliland et al. (1993) Gilliland R. L., et al., 1993, AJ, 106, 2441
  * Green (2018) Green G., 2018, The Journal of Open Source Software, 3, 695
  * Green et al. (2019) Green G. M., Schlafly E., Zucker C., Speagle J. S., Finkbeiner D., 2019, ApJ, 887, 93
  * Grevesse & Sauval (1998) Grevesse N., Sauval A. J., 1998, Space Sci. Rev., 85, 161
  * Herwig (2000) Herwig F., 2000, A&A, 360, 952
  * Holtzman et al. (2018) Holtzman J. A., et al., 2018, AJ, 156, 125
  * Huber et al. (2017) Huber D., Bryson S. T., et al. 2017, VizieR Online Data Catalog, p. IV/34
  * Iglesias & Rogers (1993) Iglesias C. A., Rogers F. J., 1993, ApJ, 412, 752
  * Iglesias & Rogers (1996) Iglesias C. A., Rogers F. J., 1996, ApJ, 464, 943
  * Irwin (2004) Irwin A. W., 2004, The FreeEOS Code for Calculating the Equation of State for Stellar Interiors, http://freeeos.sourceforge.net/
  * Jermyn et al. (2021) Jermyn A. S., Schwab J., Bauer E., Timmes F. X., Potekhin A. Y., 2021, ApJ, 913, 72
  * Jermyn et al. (2023) Jermyn A. S., et al., 2023, ApJS, 265, 15
  * Joyce & Tayar (2023) Joyce M., Tayar J., 2023, Galaxies, 11, 75
  * Karakas (2017) Karakas A. I., 2017, Low- and Intermediate-Mass Stars. Springer International Publishing, Cham, pp 1–21, doi:10.1007/978-3-319-20794-0_117-1, https://doi.org/10.1007/978-3-319-20794-0_117-1
  * Khan et al. (2019) Khan S., et al., 2019, A&A, 628, A35
  * Kippenhahn et al. (2013) Kippenhahn R., Weigert A., Weiss A., 2013, Stellar Structure and Evolution, doi:10.1007/978-3-642-30304-3.
  * Lebreton et al. (2001) Lebreton Y., Fernandes J., Lejeune T., 2001, A&A, 374, 540
  * Leiner et al. (2019) Leiner E., Mathieu R. D., Vanderburg A., Gosnell N. M., Smith J. C., 2019, ApJ, 881, 47
  * Li et al. (2023) Li G., et al., 2023, arXiv e-prints, p. arXiv:2311.16991
  * Li et al. (2024) Li T., Bi S., Davies G. R., Bedding T. R., Li Y., Stello D., Reyes C., 2024, MNRAS, 530, 2810
  * Lindegren et al. (2018) Lindegren L., et al., 2018, A&A, 616, A2
  * Lindegren et al. (2021) Lindegren L., et al., 2021, A&A, 649, A4
  * Lindsay et al. (2024) Lindsay C. J., Ong J. M. J., Basu S., 2024, arXiv e-prints, p. arXiv:2402.12461
  * Liu et al. (2019) Liu F., Asplund M., Yong D., Feltzing S., Dotter A., Meléndez J., Ramírez I., 2019, A&A, 627, A117
  * Magic et al. (2010) Magic Z., Serenelli A., Weiss A., Chaboyer B., 2010, ApJ, 718, 1378
  * Marigo & Aringer (2009) Marigo P., Aringer B., 2009, A&A, 508, 1539
  * Mathieu & Latham (1986) Mathieu R. D., Latham D. W., 1986, AJ, 92, 1364
  * Mathieu et al. (2003) Mathieu R. D., van den Berg M., Torres G., Latham D., Verbunt F., Stassun K., 2003, AJ, 125, 246
  * Mathys (1991) Mathys G., 1991, A&A, 245, 467
  * Miglio et al. (2012) Miglio A., et al., 2012, MNRAS, 419, 2077
  * Montgomery et al. (1993) Montgomery K. A., Marschall L. A., Janes K. A., 1993, AJ, 106, 181
  * Morton (2015) Morton T. D., 2015, isochrones: Stellar model grid package, Astrophysics Source Code Library, record ascl:1503.010 (ascl:1503.010)
  * Mosumgaard et al. (2018) Mosumgaard J. R., Ball W. H., Silva Aguirre V., Weiss A., Christensen-Dalsgaard J., 2018, MNRAS, 478, 5650
  * Nguyen et al. (2022) Nguyen C. T., et al., 2022, A&A, 665, A126
  * Nissen et al. (1987) Nissen P. E., Twarog B. A., Crawford D. L., 1987, AJ, 93, 634
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton B., et al., 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton B., et al., 2019, ApJS, 243, 10
  * Pinsonneault et al. (2004) Pinsonneault M. H., Terndrup D. M., Hanson R. B., Stauffer J. R., 2004, ApJ, 600, 946
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Reimers (1975) Reimers D., 1975, in , Problems in stellar atmospheres and envelopes.. Springer-Verlag, pp 229–256
  * Rogers & Nayfonov (2002) Rogers F. J., Nayfonov A., 2002, ApJ, 576, 1064
  * Roxburgh (1992) Roxburgh I. W., 1992, A&A, 266, 291
  * Salaris & Cassisi (2017) Salaris M., Cassisi S., 2017, Royal Society Open Science, 4, 170192
  * Sanders (1977) Sanders W. L., 1977, A&AS, 27, 89
  * Sandquist (2004) Sandquist E. L., 2004, MNRAS, 347, 101
  * Sandquist et al. (2016) Sandquist E. L., et al., 2016, ApJ, 831, 11
  * Sandquist et al. (2021) Sandquist E. L., et al., 2021, AJ, 161, 59
  * Sarajedini et al. (2009) Sarajedini A., Dotter A., Kirkpatrick A., 2009, ApJ, 698, 1872
  * Saumon et al. (1995) Saumon D., Chabrier G., van Horn H. M., 1995, ApJS, 99, 713
  * Schlafly & Finkbeiner (2011) Schlafly E. F., Finkbeiner D. P., 2011, The Astrophysical Journal, 737, 103
  * Schlegel et al. (1998) Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525
  * Schönrich et al. (2019) Schönrich R., McMillan P., Eyer L., 2019, MNRAS, 487, 3568
  * Song et al. (2020) Song N., Alexeeva S., Sitnova T., Wang L., Grupp F., Zhao G., 2020, A&A, 635, A176
  * Souto et al. (2018) Souto D., et al., 2018, ApJ, 857, 14
  * Souto et al. (2019) Souto D., et al., 2019, ApJ, 874, 97
  * Spoo et al. (2022) Spoo T., et al., 2022, AJ, 163, 229
  * Stassun & Torres (2018) Stassun K. G., Torres G., 2018, ApJ, 862, 61
  * Stassun et al. (2019) Stassun K. G., et al., 2019, AJ, 158, 138
  * Stello et al. (2006) Stello D., et al., 2006, MNRAS, 373, 1141
  * Stello et al. (2016) Stello D., et al., 2016, ApJ, 832, 133
  * Strom et al. (1971) Strom S. E., Strom K. M., Bregman J. N., 1971, PASP, 83, 768
  * Taylor (2007) Taylor B. J., 2007, AJ, 133, 370
  * Timmes & Swesty (2000) Timmes F. X., Swesty F. D., 2000, ApJS, 126, 501
  * Trampedach et al. (2014a) Trampedach R., Stein R. F., Christensen-Dalsgaard J., Nordlund Å., Asplund M., 2014a, MNRAS, 442, 805
  * Trampedach et al. (2014b) Trampedach R., Stein R. F., Christensen-Dalsgaard J., Nordlund Å., Asplund M., 2014b, MNRAS, 445, 4366
  * VandenBerg & Stetson (2004) VandenBerg D. A., Stetson P. B., 2004, PASP, 116, 997
  * VandenBerg et al. (2006) VandenBerg D. A., Bergbusch P. A., Dowler P. D., 2006, ApJS, 162, 375
  * VandenBerg et al. (2007) VandenBerg D. A., Gustafsson B., Edvardsson B., Eriksson K., Ferguson J., 2007, ApJ, 666, L105
  * Viani & Basu (2020) Viani L. S., Basu S., 2020, ApJ, 904, 22
  * Yadav et al. (2008) Yadav R. K. S., et al., 2008, A&A, 484, 609
  * Zinn (2021) Zinn J. C., 2021, AJ, 161, 214
  * Zinn et al. (2019) Zinn J. C., Pinsonneault M. H., Huber D., Stello D., 2019, ApJ, 878, 136