The K2 M67 Study: Precise Mass for a Turnoff Star in the Old Open Cluster M67

M67 was observed during Campaigns 5, 16, and 18 of the K2 mission. WOCS 11028 was observed during all of the campaigns in a custom aperture with long cadence (30 min) exposures. Single eclipses were observed during campaigns 5 and 16, having about a 2.7% decrease in flux. (No eclipse was observed in campaign 18 data, and the ephemeris indicates that one was not expected.) The eclipse times were measured to be at BJD $2457178.9531\pm 0.00088$, and $2458117.8789\pm 0.00075$ using the bisector method of Kwee & van Woerden (1956). The radial velocity curve shows that these are secondary eclipses (of the cooler, less massive star). In the K2 light curves, there is no sign of primary eclipses at the four epochs predicted by the orbital ephemeris (see Figure 1).

Because the K2 camera has relatively large pixels ($3\farcs 98$ on a side), we briefly discuss the possibility of contaminating light from unrelated stars. Although WOCS 11028 is a member of M67, it is in the outskirts of the cluster. As shown in Figure 2, the nearest bright star is WOCS 6030 / Sanders 619 at a distance of $54\arcsec$, and about 0.3 mag fainter in Gaia $G$ band. There is a Gaia source (indicated in Figure 2) $8\arcsec$ away, but nearly 8 $G$ magnitudes fainter. There are three other Gaia sources between 27 and $34\arcsec$ away, but all are at least 5.3 mag fainter in $G$. In the light curves described below, the photometric apertures did not extend more than 5 pixels away from the photocenter, and so avoided the brighter stars. As a result, we assume contamination of the K2 apertures used for WOCS 11028 is negligible.

Because of the incomplete gyroscopic stabilization during the K2 mission, systematic effects on the light curve are known to be substantial. We used two different light curves for WOCS 11028 as part of our investigation of systematic error sources. The first light curve we experimented with utilized the K2SFF pipeline (Vanderburg & Johnson, 2014; Vanderburg et al., 2016) that involved stationary aperture photometry along with correction for correlations between the telescope pointing and the measured flux. The light curve that resulted still retained small trends over the long term, so these were removed by fitting the out-of-eclipse points with a low-order polynomial and dividing the fit. The second light curve came from version 2.0 of the EVEREST pipeline (Luger et al., 2018), which uses pixel-level decorrelation to remove instrumental effects. We used the detrended and co-trending basis vector-corrected light curve. However, the eclipse depth is about 4% deeper in the EVEREST light curve than in the K2SFF version, so we need to gauge the effects on measurements of the binary characteristics. Because the K2 light curves have high signal-to-noise, they have a substantial influence on best-fit binary star parameter values (see section 3.4). The scatter among the best-fit parameter values when comparing runs with different K2 light curve reductions is sometimes larger than the statistical uncertainties. Because there is no clearly superior method for processing the K2 light curves, we will take the scatter in the binary model parameter values as a partial indicator of systematic error resulting from the processing.

Because the eclipses of the WOCS 11028 binary provide minimal information on the sizes of the individual stars, we need other well-measured characteristics of the binary’s stars in order to extract age information for the cluster. Fortunately, M67 has been heavily observed over a wide wavelength range, making it possible to construct well-sampled spectral energy distributions for the binary and similarly bright stars in the cluster. In this section, we describe the spectroscopy and photometry we have assembled, and the efforts to put the observations on a consistent flux scale.

Ultraviolet: We obtained photometry and a spectrum from the Galaxy Evolution Explorer (GALEX; Martin et al. 2005) archive for the NUV passband ($1771-2831$ Å). WOCS 11028 was imaged twice, for 1691.05 s (GI1 proposal 94, P.I. W. Landsman) and 5555.2 s (GI1 proposal 55. P.I. K. Honeycutt). As the archived magnitudes are based on count rates with minimal background contributions, we computed a final magnitude and flux based on the average count rate for both observations. Morrissey et al. (2007) describes the characteristics of the GALEX photometry and its calibration to flux. GALEX magnitudes are on an AB system (Oke & Gunn, 1983), and we used the zero point magnitude ($m_{NUV}=20.08$) and reference flux ($2.06\times 10^{-16}$ erg s^-1 cm² Å^-1) to convert to flux.

A NUV grism spectrum was taken as part of GI1 proposal 94 (PI: W. Landsman) on 2005 Feb. 18 with a total exposure time of 27653 s. The NUV spectrum covers a wider wavelength range than the NUV photometry filter, but there was effectively no signal detected for WOCS 11028 at wavelengths less than about 2000 Å. The spectrum was obtained from MAST, and the flux calibration from the pipeline reduction was used.

Observations were taken of a smaller portion of M67 (including WOCS 11028) using the UVOT telescope on the Swift satellite. We collected UV fluxes in the $uvw1$, $uvm2$, and $uvw2$ bands from the Swift UVOT Serendipitous Source Catalogue (version 1.1; Page et al. 2014). The three bands mostly cover the same wavelength range as the GALEX NUV filter, but with somewhat finer resolution. Flux correction factors from Poole et al. (2008) were applied to go from the gamma-ray burst spectrum calibration given in the archive to one utilizing Pickles (1998) library stars.

Near-ultraviolet, Optical and Near-infrared: M67 has been frequently observed from the ground for the purposes of photometric calibration. For the purposes of an SED, narrow-band filters are particularly useful, and we discuss these first. Balaguer-Núñez et al. (2007) presented Strömgren $uvby$ photometry for the cluster. Because M67 stars are commonly used as standards in Strömgren photometry (Nissen et al., 1987), we can be fairly assured that the observations are tied to the standard system. We employed reference fluxes from Gray (1998) to convert the magnitudes to fluxes.

Fan et al. (1996) conducted wide-field observations of M67 in a series of narrow-band filters (“BATC”) covering from 3890 to 9745 Å  using bandpasses avoiding most important sky lines. The wavelength coverage of the filters is better in the near-infrared, making them a good complement to Strömgren photometry. The BATC survey goal was spectrophotometry at the 1% level, and the study largely achieved that, judging from the low scatter in their color-magnitude diagrams. Their reported magnitudes are on the Oke & Gunn AB system. We found that fluxes from these magnitudes were systematically higher than similar observations in other systems. We collected magnitudes in a larger set of filters from BATC Data Release 1¹¹1http://vizier.u-strasbg.fr/viz-bin/VizieR-3?-source=II/262/batc, and found that these were more consistent, probably due to improved calibration (Zhou et al., 2003).

Narrow-band photometry is also available in the Vilnius filter system for many stars in M67 (although not WOCS 11028) in the study of Laugalys et al. (2004). We converted the reported magnitudes to fluxes using the zeropoints from Mann & von Braun (2015). The filters in the Vilnius system are somewhat denser in the blue portions of the optical, and these again complement the spectrophotometry in the BATC system.

We have Johnson-Cousins photometry in $BVI_{C}$ from Sandquist (2004) and Yadav et al. (2008), in $UBVR_{C}I_{C}$ from Montgomery et al. (1993), in $BVRI$ from Nardiello et al. (2016b), and in $BV$ from Data Release 9 of the AAVSO Photometric All-Sky Survey (APASS; Henden et al. 2015). These measurements are on a Vega magnitude system, and so have been converted to fluxes using reference magnitudes from Table A2 of Bessell et al. (1998), and accounting for the known reversal of the zero point correction rows for $f_{\lambda}$ and $f_{\nu}$. The absolute calibrations of each of these studies are unavoidably different for the same filter bands, and this will contribute to the noise in the SEDs. However, this does not affect their use in determining the relative contributions of the two stars in the WOCS 11028 binary (see §2.2.1). (The $R$ and $I$ filter observations given in Nardiello et al. (2016b) were actually taken in SDSS $r$ and $i$ filters, and have been recalibrated to the Sloan DR12 system (D. Nardiello, private communication). These data were not used in fits to SEDs.)

There are several additional surveys that provide calibrated broad-band photometric observations. We used PSF magnitudes from Data Release 14 of the Sloan Digital Sky Survey (SDSS; York et al. 2000), calibrated according to Finkbeiner et al. (2016). The SDSS is nearly on the AB system, with small offsets in $u_{SDSS}$ and $z_{SDSS}$ that we have also corrected for here. The Pan-STARRS1 survey (Kaiser et al., 2010) contains photometry in 5 filters, and we use their mean PSF magnitudes here. Zero points for its AB magnitude system are given in Schlafly et al. (2012). The APASS survey also observed the cluster in Sloan $g^{\prime}r^{\prime}i^{\prime}$ filters, and their photometry was flux calibrated from the AB system.

Finally, Gaia has already produced high-precision photometry extending far down the main sequence of M67 as part of Data Release 2. We obtained the fluxes in the $G$, $G_{BP}$, and $G_{RP}$ bands from the Gaia Archive.

Infrared: We obtained Two-Micron All-Sky Survey (2MASS; Skrutskie et al. 2006) photometry in $JHK_{s}$ from the All-Sky Point Source Catalog, and have converted these to fluxes using reference fluxes for zero magnitude from Cohen et al. (2003). Many stars also have photometry in three bands from the Wide Field Infrared Explorer (WISE; Wright et al. 2010), which were also converted to fluxes using tabulated reference fluxes at zero magnitude.

We have measured radial velocities for the WOCS 11028 binary using three spectroscopic datasets. The first and largest set we employ here comes from the CfA Digital Speedometers (Latham, 1985, 1992) on the 1.5 m Tilinghast reflector at Fred Lawrence Whipple Observatory and the MMT. These observations were recordings of a single echelle order covering 516.7 to 521 nm around the Mg I b triplet using an intensified photon-counting Reticon detector. The spectra were taken as part of a larger monitoring campaign of M67 stars.

We made two observations of the binary using the HARPS-N spectrograph (Cosentino et al., 2012) on the 3.6 m Telescopio Nazionale Galileo (TNG). HARPS-N is a fiber-fed echelle that has a spectral resolving power $R=115000$, covering wavelengths from 383 to 693 nm. The spectra were processed, extracted, and calibrated using the Data Reduction Software (version 3.7) provided with the instrument.

We also obtained five archival spectra from the APOGEE database (Holtzman et al., 2015) that were taken between January and April 2014. These are $H$-band infrared ($1.51-1.70\mu$m) spectra with a spectral resolution $R\sim 22500$ taken on the 2.5 m Sloan Foundation Telescope. We used APOGEE flags to mask out portions of the spectrum that were strongly affected by sky features, and continuum normalized the spectra using a median filter.

The radial velocities were measured using the spectral separation algorithm described in González & Levato (2006). In the first iteration step, a master spectrum for each component is isolated by aligning the observed spectra using trial radial velocities for that component and then averaging. This immediately de-emphasizes the lines of the non-aligned stellar component, and after the first determination of the average spectrum for each star, the contribution of the non-aligned star is subtracted before the averaging in order to better clean each spectrum. The radial velocities can also be remeasured from spectra with one component subtracted. This procedure is repeated for both components and continued until a convergence criterion is met. We measure the broadening functions (Rucinski, 2002) to determine the radial velocities using narrow-lined synthetic spectra as templates. Radial velocity measurements from broadening functions improve accuracy in cases when the lines from the two stars are moderately blended (Rucinski, 2002).

The CfA spectra covered a wide range of orbital phases that made it possible to use a large number of spectra in calculating average spectra for the two stars. In most cases, spectra were left out due to low signal-to-noise ratio. We used 17 out of 28 spectra in determining the average spectrum for the primary, but only used the 13 spectra with the most clearly detected secondary component to determine its average spectrum. With good average spectra, it was possible to subtract each component out of observed spectra even when they were taken close to crossing points (phases $\phi\approx 0.05$ and 0.6). After some experimentation, we found that synthetic templates (from the grid of Coelho et al. (2005)) of 5250 K optimized the detection of broadening function peaks for the secondary star. A template with 6250 K was used for the primary star. For the CfA velocities only, we used run-by-run corrections derived from observations of velocity standards. Run-by-run offsets measured for the CfA spectra ranged from $-1.52$ to $+2.21$ km s^-1. We initially assigned velocity uncertainties derived from the spectral separation analysis, but scaled these to ensure that the scatter around a best-fit model for this dataset was consistent with the measurement errors. (In other words, we forced the reduced $\chi^{2}$ value to 1.) The average uncertainty for the primary star velocities was 1.25 km s${}^{-1}\;$, and for the secondary star it was 4.41 km s^-1.

The APOGEE spectra were analyzed separately, and detection of the secondary star features were considerably more secure. The synthetic spectral templates were taken from ATLAS9/ASS$\epsilon$T models (Zamora et al., 2015) calculated for the APOGEE project. The uncertainty estimate for each APOGEE velocity was generated from the rms of the velocities derived from the three wavelength bands in the APOGEE spectra. This was typically around 0.15 km s${}^{-1}\;$for the primary star and 0.65 km s${}^{-1}\;$for the secondary star. The velocity was corrected to the barycentric system using values calculated by the APOGEE pipeline. For all five spectra, the broadening function peaks for the two stars were resolvable thanks to the high resolution of the spectra and fortuitous phases of observation.

Because we only had two observations using the HARPS-N spectrograph and one of them was taken at a phase very near a crossing point, we determined velocities using broadening functions alone. To get an internal measure of the velocity uncertainties, we measured the velocities in four 30 nm subsections of the spectra and computed the error in the mean.

Because the synthetic spectra used in our broadening function measurements do not account for gravitational redshifts of the stars, our velocities should have an offset relative to velocities derived using an observed solar template. For the CfA spectra, the run-by-run corrections were computed using dawn and dusk sky exposures, and this defines the native CfA velocity system used in Geller et al. (2015), incorporating the gravitational redshift of the Sun. For consistency, we subtracted the combined gravitational redshift of the Sun and gravitational blueshift due to the Earth (0.62 km s^-1) from the APOGEE and HARPS-N velocities we tabulate. Based on the derived masses and radii for the stars in the binary, they should have slightly larger gravitational redshifts than the Sun (0.66 and 0.69 km s^-1, versus 0.64 km s${}^{-1}\;$for the Sun) but we did not correct for these smaller differences. In our later model fits, we did, however, allow the system velocities for the two stars to differ to account for effects like this and convective blueshifting that could produce small systematic shifts. These could affect the measured radial velocity amplitudes $K_{1}$ and $K_{2}$ (and the stellar masses) if not modeled.

The measured radial velocities are presented in Table 1, and the phased radial velocity measurements are plotted in Fig. 9. The model fits are described in Section 3.4, but the parameters that were used to fit the velocity data were the velocity semi-amplitude of the primary star $K_{1}$, mass ratio $q=M_{2}/M_{1}=K_{1}/K_{2}$, eccentricity ($e$), argument of periastron ($\omega$), and systematic radial velocities $\gamma_{1}$ and $\gamma_{2}$. In the end, we do not see noticeable differences between the systematic velocities.

It is always prudent to check for systematic differences in velocity when utilizing datasets from more than one observational set-up. Among the measured primary star velocities, mean residuals (observed minus computed) were $-0.14$ km s${}^{-1}\;$for APOGEE, $+0.02$ km s${}^{-1}\;$for CfA, and $0.02$ km s${}^{-1}\;$for HARPS-N. For the secondary star, we found mean residuals of $+0.78$ km s${}^{-1}\;$for CfA, $+0.10$ km s${}^{-1}\;$for APOGEE, and $+0.09$ km s${}^{-1}\;$for HARPS-N (one measurement). Of these, only the average residual for the CfA secondary velocities was more than one standard devation (0.20 km s^-1) away from zero. Because the APOGEE and HARPS-N spectra have higher signal-to-noise than the CfA spectra, they are the main contributors to measurement of the radial velocity amplitude of the secondary star, which critically affects the calculated mass of the primary star.

WOCS 11028 is a fairly large distance (about $13\farcm 7$) from the center of M67, but Carrera et al. (2019) find that half of the cluster turnoff stars fall within $12\farcm 5$ of the cluster center. As a result, its position is not strong evidence of being a field star.

All ground-based proper motion studies (Sanders, 1977; Girard et al., 1989; Zhao et al., 1993; Yadav et al., 2008; Nardiello et al., 2016a) indicate a high probability ($\geq 90$%) of cluster membership for WOCS 11028, with the exception of Krone-Martins et al. (2010), who give 10%). Gaia observations from Data Release 2 (Gaia Collaboration et al., 2018a) have produced much more precise proper motions recently, and the binary still has a proper motion vector ($\mu_{\alpha}=-11.08\pm 0.04$ mas yr^-1, $\mu_{\delta}=-3.14\pm 0.04$ mas yr^-1) safely residing among other likely members (centered at $\mu_{\alpha}=-10.97$ mas yr^-1, $\mu_{\delta}=-2.94$ mas yr^-1; Gaia Collaboration et al. 2018b). Similarly, the Gaia parallax ($\omega=1.111\pm 0.024$ mas) is within $1\sigma$ of the cluster mean ($\bar{\omega}=1.132$ mas).

Geller et al. (2015) published membership probabilities based on their radial velocity survey by comparing the velocity distributions of cluster members and field stars. They classified WOCS 11028 as a binary member (98% probability) based on their system velocity ($32.91\pm 0.17$ km s^-1) measured from the CfA observations and carefully placed on the zeropoint of other stars in the field. The best-fit system velocity from our measurements ($\gamma=32.37$ km s^-1) is slightly off from the mean cluster velocity of Geller et al. (33.64 km s${}^{-1}\;$with a radial velocity dispersion of $0.59^{+0.07}_{-0.06}$ km s^-1). However, even with our lower $\gamma$ velocity, the radial velocity membership probability is still 95%.

Based on the three-dimensional kinematic evidence, the binary is a high probability cluster member.

With the release of Gaia DR2 parallaxes, the M67 distance modulus should be revisited, as it will play a large role in comparisons with models later in the article. The mean parallax derived for cluster members from Gaia DR2 measurements ($\bar{\omega}=1.1325\pm 0.0011$ mas; Gaia Collaboration et al. 2018b) is very precise, but possibly affected by a zeropoint uncertainty. Lindegren et al. (2018) finds an offset of about $-30\>\mu$as (with tabulated parallaxes being smaller than reality), while systematic variations with position, magnitude, and color were below $10\>\mu$as. At the other extreme, Stassun & Torres (2018) found an offset of $-82\pm 33\>\mu$as using a sample of eclipsing binaries. Zinn et al. (2019) find offsets of $-52.8$ and $-50.2\>\mu$as using comparisons with asteroseismic giants and red clump stars in the Kepler field. Schönrich et al. (2019) derived an offset of $-54\pm 6\>\mu$as from an analysis of all stars (more than 7 million) in the radial velocity sample of Gaia DR2.

These indications are consistent with the situation for M67. Without correction, the Gaia cluster mean distance modulus would be $(m-M)_{0}=9.730\pm 0.002$, but this is comparable to what has been derived previously for the distance modulus with extinction. For example, Pasquini et al. (2008) derive $(m-M)_{V}=9.76\pm 0.06\pm 0.05$ (statistical and systematic uncertainties) from ten solar analogs, Sandquist (2004) derived $9.72\pm 0.05$ using metal-rich field dwarfs with Hipparcos parallaxes, and Stello et al. (2016) found $9.70\pm 0.04$ from K2 asteroseismology of more than 30 cluster red giants. With a $V$-band extinction $A_{V}\approx 0.12$ expected, this is a significant discrepancy. For subsequent calculations in this article, we will assume a $-54\>\mu$as offset in the Gaia parallaxes and include a systematic uncertainty of $30\>\mu$as to conservatively account for disagreement in the parallax offsets in the literature. The resulting distance modulus is $(m-M)_{0}=9.63\pm 0.06$.

Even though the WOCS 11028 binary has a relatively large orbital separation, our analysis of another wide binary in M67 (WOCS 12009; Sandquist et al., 2018) indicated that the brighter star was likely the product of an earlier merger. A key piece of evidence for that binary was the lack of detectable Li. Single main sequence stars of similar brightness in the cluster inhabit the Li plateau, a broad maximum in the Li abundance. The lack of detectable Li was an indication of nuclear processing, consistent with surface material having originally been in a lower-mass main sequence star with a deeper surface convection zone.

To have greater assurance that the stars of WOCS 11028 have evolved as isolated single stars and can be used to constrain the cluster age, we can also use the Li abundance here. If the brighter star in WOCS 11028 was unmodified since the cluster’s formation, its mass ($\sim 1.2M_{\odot}$) would place it toward the bright end of the cluster’s Li plateau with abundances $A(\mbox{Li})\approx 2.5$ (Pace et al., 2012). On the other hand, the secondary star (mass near $0.9M_{\odot}$) should not have detectable Li. Stars in the plateau have outer convective zones that do not transport Li nuclei down to temperatures where nuclear reactions can burn them, and additional mixing processes must have little, if any effect. Even with the diluting effects of the secondary star’s light on the Li lines, Li should be detectable if it is present with the plateau abundance.

We are able to detect Li in our HARPS-N spectra as shown in Fig. 10. If the star had a history of interaction (mass transfer, or a merger involving lower mass stars), it is unlikely that the Li would be detectable. If the brighter star formed in a merger, the more massive merging star could not have been more than about $0.2M_{\odot}$ less massive than the current star or else a surface convection zone would have had a chance to consume its Li during its main sequence evolution. In addition, the less massive merging star would take up residence in the core of the merger remnant because of its lower entropy, rebooting its nuclear evolution with gas having higher hydrogen abundance. After some adjustment on a thermal timescale, it would have been bluer than other cluster stars of the same mass. A precise determination of the Li abundance is complicated by the secondary star’s light, and is beyond the scope of this work. However, the strong Li detection coupled with strong indications that WOCS 11028 A has a color consistent with its fellow cluster members effectively rules out anything except a very early merger and/or a merger with a very low mass object. Effectively, this would make today’s star the same as an unmodified star of the same mass.

We do not see evidence of rotational spot modulation in the K2 light curves due to spots, but our HARPS-N spectra allow us to constrain the rotational speeds of the stars in the binary via the widths of the broadening function peaks. We are unable to detect rotation with an upper limit $v_{\rm rot}\sin i<5$ km s^-1. If the primary was pseudo-synchronized, its rotation velocity would be around this limit (rotation period near 14 d). However, the timescale calculated for this (Hut, 1981) is far longer than the cluster age, thanks to the substantial separation even at periastron. Considered as a single star, the primary has probably spun down somewhat, although not to the same degree as lower mass stars with deeper surface convection zones. Angus et al. (2015) calibrated their gyrochronological models against the rotation of asteroseismic targets in the Kepler field, and their sample probably gives the best indication of typical rotation rates for old stars more massive than the Sun. Their isochrones predict a rotation period near 16 day, and this would imply a rotation velocity below our ability to detect it spectroscopically. As an old cluster, M67 is expected to have slowly rotating stars, but this does put a limit on how long ago mass transfer or a merger could have happened if either star was somehow modified. A rapidly rotating star with the mass of our primary star would require approximately 3 Gyr to reach our detection limit according to the rotational isochrones of Angus et al. (2015), and this is the majority of the cluster’s age.

We used the ELC code (Orosz & Hauschildt, 2000) to simultaneously model the ground-based radial velocities and Kepler K2 photometry for the WOCS 11028 binary. The code is able to use a variety of optimizers to search the complex multi-dimensional space of parameters used in the binary star model. In our particular case, we used a set of 15 parameters. Two of the parameters were the orbital period $P$, and the reference time of periastron $t_{P}$. Six additional parameters mostly characterize the orbit: the velocity semi-amplitude of the primary star $K_{1}$, mass ratio $q=M_{2}/M_{1}=K_{1}/K_{2}$, systematic radial velocities⁴⁴4We allow for the possibility of differences for the two stars that could result from differences in convective blueshifts or gravitational redshifts. $\gamma_{1}$ and $\gamma_{2}$, eccentricity ($e$), and argument of periastron ($\omega$). For a typical double-lined spectroscopic and eclipsing binary, both the radial velocities and eclipse light curves constrain $e$ and $\omega$, with the phase spacing of the eclipse playing a large role. But because this binary only has one eclipse per cycle, these parameters are constrained by the radial velocities almost exclusively.

The light curve has a large role in constraining the inclination parameter $i$, but there are potentially correlations between inclination and choices for radius and temperature parameters. For radii, we used the sum of the stellar radii $R_{1}+R_{2}$ and their ratio $R_{1}/R_{2}$ as parameters. For this binary, the sum of radii is constrained by the measurements of the single grazing eclipse in each cycle along with the well-determined spectroscopic parameters of the orbit. The radius ratio can be constrained by luminosity ratios we have derived from our SED analysis (§2.2) in concert with temperature information. For a typical eclipsing binary, the temperature ratio would be constrained by the relative depths of the eclipses, but that is not measurable here. We have separate constraints on the temperatures of the stars from SED fits to the proxy stars, and so we use these as fit parameters: temperatures of the primary $T_{1}$ and secondary $T_{2}$.

Finally, there is some dependence of the fits on the limb darkening, although relatively little because only the near-limb of the secondary star is probed by the grazing eclipses we see. So we fit for two quadratic limb darkening law coefficients ($q_{1B},q_{2B}$) for the secondary star (the only star that is eclipsed) in the Kepler bandpass. Our runs of the ELC code employed the Kipping (2013) algorithm, which involves a search over a triangular area of parameter space of physically realistic values. This forces the model star to a) darken toward the limb, and b) have a concave-down darkening curve. In the final analysis, we find that the limb darkening coefficients are virtually unconstrained, but by allowing the code to systematically explore potential values, the uncertainties will be incorporated in our uncertainties of the other stellar parameters. Owing to the substantial orbital separation, we assumed that the stars are spherical. Given this, we used the algorithm given in Short et al. (2018) as implemented in ELC to rapidly compute the eclipses. Finally, because there is little or no out-of-eclipse light modulation, we did not see a need to model spots.

For the K2 light curve modeling, we only included observations that were within about $\pm 0.01$ in phase of the eclipse, and a similar section around where the other eclipse would have been if the binary was closer to edge-on. In this way, our models “see” that there is only one eclipse per orbital cycle. In doing this, we leave out the majority of observations, which may be affected by systematics due to spacecraft pointing corrections and other instrumental signals. Because the K2 light curves use long cadence data with 30 minute exposures, they also effectively measure the average flux during that time. To account for this, the ELC code integrates the computed light curves in each observed exposure window.

The quality of the model fit was quantified by an overall $\chi^{2}$ derived from comparing the radial velocity and light curve data to the models, as well as from how well a priori constraints were matched. For this binary, we used stellar temperatures derived from SED fits and luminosity ratios in $BVI_{C}JHK_{s}$ derived from the modeling of the binary SED with single cluster stars. To illustrate the effect of these constraints, we conducted a run without the luminosity ratios, as shown in Table 2.

To try to ensure that different datasets were given appropriate weights in the models, we empirically scaled uncertainties on different datasets to produce a reduced $\chi^{2}$ of 1 relative to a best-fit model for that particular dataset. For spectroscopic velocity measurements, this meant scaling differently according to the source spectrograph and the star (primary and secondary). The K2 photometric light curve uncertainties were scaled independently.

We used a differential evolution Markov Chain optimizer (Ter Braak, 2006) for seeking the overall best-fit model and exploring parameter space around that model to generate a posterior probability sampling. Approximate $1\sigma$ parameter uncertainties were derived from the parts of the posterior distributions containing 68.2% of the remaining models from the Markov chains. Gaussians were good approximations to the posterior distributions for all parameters except limb darkening coefficients in the run where luminosity constraints were applied. The results of the parameter fits are provided in Table 2, and the posterior distributions are shown in Figure 11. A comparison of the K2 light curve with the best-fit model is shown in Figure 12, and a comparison of the radial velocity measurements with models are shown in Figure 9.

As can be expected, when we do not enforce luminosity ratio constraints from SED fitting, there is no information on the ratio of radii in the data, and as a result, the radii of the stars are very poorly constrained. The mass determinations are essentially unaffected, however, because that information is contained in the radial velocity curves and the orbital inclination (derived from the light curve).

As is usually the case, the chemical composition is an important component of interpreting the measurements and understanding the implications for the cluster age and for the stellar physics in the models we use. We discussed high-precision cluster [Fe/H] measurements relative to the Sun in Sandquist et al. (2018) in order to minimize systematic errors due to composition. Since that paper though, it has become clearer that diffusion is playing a role in the measured surface compositions of M67 stars. Significant trends in abundance are seen as a function of the evolutionary state of the stars, with subgiants having higher measured abundances than turnoff stars (Önehag et al., 2014; Bertelli Motta et al., 2018; Liu et al., 2019; Souto et al., 2019). Because subgiants should have deep enough convection zones to homogenize surface layers that were affected by heavy element settling, the subgiant abundances should be more representative of the initial bulk heavy element composition of the stars. The spectroscopic studies generally come to the conclusion that initial composition of the cluster members was [Fe/H]$=+0.05$ to $+0.10$. In the model comparisons below, we will consider initial metal compositions between solar and [Fe/H]$=+0.10$.

The verticality of the cluster CMD at the position of WOCS 11028 A constrains the present direction of its evolution track, and we illustrate the discussion with MESA model tracks (version r11701; Paxton et al. 2011, 2013, 2015, 2018, 2019) in the largest panel of Figure 14. The models predict that for the first 2 Gyr of the star’s life, it was evolving mostly upward in luminosity, with slowly increasing temperature before hitting a maximum and starting to decrease. The star must have started evolving more rapidly toward cooler temperature in order to change the local slope of the main sequence isochrone from increasing temperature to constant temperature with increasing mass. Comparing the markers showing common age for three stars of slightly different mass, the triplet becomes vertical near the midpoints of the tracks during this coolward movement. For greater age, the triplets have decreasing temperature with higher mass, and in M67 this corresponds to stars between the lower turnoff and the faint end of the gap.

The precise point where the vertical slope is produced will depend on details of the fuel consumption in this phase, but we can be assured that the star must still be evolving faster toward lower temperature than slightly lower mass stars (in other words, accelerating) or else it would not be able to match the temperature of those stars that should be evolving at a very similar speed. This in turn assures us that convection must be present in the core of the star — low-mass stars that do not establish convective cores on the main sequence only accelerate toward lower temperature when they are past central hydrogen exhaustion while on the subgant branch. So the mass measurement for WOCS 11028 A gives us an observational lower limit for stars that establish convective cores on the main sequence.

In the evolution tracks, the turn toward lower temperature occurs when the extent of the convective core increases in the latter half of core hydrogen burning. The CNO cycle overtakes the $pp$ chain as the dominant energy generating reaction network while central temperature increases and hydrogen abundance declines. So not only do we have good measurements of the WOCS 11028 A star (described below), but we have an indication of physical conditions occuring within the star. A complete model of the star would match all of the measurements and evolutionary constraints with a unique set of parameters, while mismatches would indicate systematic errors.

With the new data on the massive star in the WOCS 11028 binary, there is an opportunity to re-examine the cluster turnoff and the constraints it places on the age and stellar physics. We can make a more precise estimate of star masses at the turnoff than has been possible before by employing models that are capable of reproducing the characteristics of WOCS 11028 A, reducing the importance of systematic errors by making a relative comparison. The observed properties of the CMD gap and the upper turnoff are strongly influenced by the extent of the convective core, which in turn affects the amount of accessible fuel and the main sequence lifetimes.

A major uncertainty about convective core extent is overshooting of the convective boundary. The stars at M67’s turnoff are uniquely interesting because they are close to the minimum mass for having convective cores, and the amount of overshooting could play an outsized role in stars with small cores. Some studies of eclipsing binaries in the $1.2-2.0M_{\odot}$ mass range (Claret & Torres, 2017, 2018) indicate that convective overshooting is consistent with zero at the mass of WOCS 11028 A, but rapidly increasing to a plateau for masses greater than about $2M_{\odot}$. On the other hand, the models of Higl & Weiss (2017) also show little or no need for overshooting in binaries with main sequence stars having masses near $1.2M_{\odot}$, but some weak evidence for overshooting in BG Ind, having a somewhat more massive ($1.42M_{\odot}$) and more evolved primary star. Constantino & Baraffe (2018) found little evidence for a mass-dependent overshooting for masses between 1.2 and $2.0M_{\odot}$. A major problem with the binary star samples in most of these studies is that overshooting reveals itself most clearly very close to core hydrogen exhaustion, and unless stars are in a very specific evolutionary phase, their conventional stellar properties (like $R$, $L$, and $T_{\rm eff}$) will not be sensitive to the overshooting. This disagreement in the literature reinforces the value of using the information contained in M67, with its larger sample of stars and mass inferences from the WOCS 11028 binary.

More recently, a number of asteroseismic studies have also attempted to constrain convective core overshooting with low-degree p modes. Studies like Viani & Basu (2020) find positive trends in overshooting amount with mass, while others like Angelou et al. (2020) find overshooting amounts that are consistent with what is found from binary and cluster calibrations but without evidence of a clear mass dependence. These issues highlight that there are systematics and degeneracies in the fitting of parameters in asteroseismic studies, and that they have not fully clarified the overshoot issue. Asteroseismic masses have significantly lower precisions than binary system masses as well. All of these factors call out for additional precise constraints on the overshooting in this mass range.

Overshooting did not play a significant role in our interpretation of the WOCS 11028 binary because the primary star is well before core hydrogen exhaustion, but it does more strongly affect the CMD positions of stars around the turnoff gap and beginning of the subgiant branch. Models indicate that a star like WOCS 11028 A reaches a maximum in surface temperature when the convective core starts expanding as the core temperature increases and the CNO cycle becomes more dominant. After the star passes its surface temperature maximum, models show that the convective core is maintained and the star evolves mostly in $T_{\rm eff}$ until shortly before core hydrogen exhaustion. At that time the convective core rapidly shrinks and disappears, and the star executes a quick movement to higher temperatures and luminosities in the HR Diagram. There is still a residue of hydrogen in the core and a steep gradient in hydrogen abundance outside that, so the evolution of the star slows down again in the early subgiant branch. Overall, this evolution track is a hybrid, with the early core hydrogen burning evolution similar to low-mass stars, and the late evolution similar to more massive stars.

The extent and lifetime of the convective core in stars slightly more massive than WOCS 11028 A is imprinted on the shape of the turnoff in the CMD and in the number of stars present. The noticeable gap seen at $12.7\lesssim G\lesssim 13$ in M67 is a demonstration of the rapid evolution that results from the disappearance of the convective core shortly before central hydrogen exhaustion. Similarly, the fairly large density of stars present just brighter than the gap shows that gas with large hydrogen abundance is residing just outside the exhausted center of the star, and the burning of this hydrogen allows the evolution to slow down again.

To define the shape of the turnoff, we derived bolometric fluxes and surface temperatures from SED fits for probable single stars at critical points in the CMD. We identified WOCS 7006 / Sanders 1003 as a likely single cluster member at the bright end of the gap based on CMD position and lack of radial velocity variation (Geller et al., 2015), and WOCS 8011 / Sanders 1219 as a likely single member at the faint end. We also examined WOCS 3003 / Sanders 1034 as a representative of the end of the most heavily populated part of the early subgiant branch. The stars are shown in CMDs in Figure 3, and their derived properties are given in Table 3.

As was shown in Figure 14, several model physics parameters play a role in setting the shape of the turnoff and subgiant branch, and most affect core convection. As a result, there is some degeneracy in the parameters that can be used to model the turnoff. For our purposes here, we will focus on convective core overshooting. Overshooting has the most substantial effects on evolution tracks shortly before the time of core hydrogen exhaustion. Core hydrogen exhaustion is delayed by overshooting, which allows the redward evolution to extend to lower temperature, giving isochrones a deeper redward kink (see bottom middle panel of Figure 14). In addition, the morphology of the upper turnoff changes in response to overshooting, with increasing amounts resulting in a more horizontal early subgiant branch. (See Figure 13, and compare the BaSTI-IAC isochrones with and without overshooting, or younger PARSEC isochrones to older ones that have less overshoot.) Both the morphology of the upper turnoff and the luminosity and temperature positions of the stars measured near the gap in M67 are indicative of weaker convective overshooting than is used in the MIST isochrones. This is a means of calibrating overshoot for stars with masses around the gap (approximately $1.32M_{\odot}$, based on models that match WOCS 11028 A). The exact amount depends somewhat on other physics, like diffusion rates, CNO cycle reaction rates, and CNO element abundances.

It has long been known that the turnoff of M67 is not well matched by isochrones (VandenBerg et al., 2007; Magic et al., 2010). Currently published theoretical isochrones should not be used to calibrate the overshooting or read an age because they do not reproduce the characteristics of WOCS 11028 A (see Section 4.2.3), and would therefore have clear systematic errors. There are unfortunately several composition and physics factors that can affect the size (and even the presence) of a convective core in models, and relatively small changes in the extent of this core can strongly affect the amount of hydrogen fuel the core gets, and therefore, the length of the star’s life. But by using models that are at least able to match WOCS 11028 A, we should be able to reduce the size of the errors.

Figure 20 shows isochrones from MESA models having $Z=0.0195$ and $Y=0.2703$, with convective core overshooting set to half the default value, or $f=0.008$. Larger amounts of overshooting move the gap to lower temperatures, and appear to be ruled out. With this set of input parameters, we are able to approximately reproduce the HR diagram positions of the three stars (described above) that reside at identifiable points in the cluster CMD using an age around 3.75-4 Gyr. In addition, Stello et al. (2016) found an average giant star mass of $1.36\pm 0.01M_{\odot}$ in their asteroseismic analysis of M67 giant stars using K2 data. This measurement is also consistent with an age near 4 Gyr for the MESA models used in Figure 20. The agreement with all of these observations implies that the models are approximately reproducing the hydrogen abundance profile in the stars at core exhaustion.