The K2 M67 Study: A Curiously Young Star in an Eclipsing Binary in an Old Open Cluster^*

M67 was observed during Campaign 5 of the K2 mission. For the central region of the cluster, all pixel data were recorded to form "superstamp" images for each long-cadence (30 minutes) exposure. One primary eclipse and two secondary eclipses were observed for WOCS 12009 during the nearly 75 day campaign, enabling the immediate confirmation of the spectroscopic orbital period. Once the eclipsing binary had been identified, we were able to identify observations of earlier eclipses in other ground-based data sets, as discussed below. On its own, the K2 light curve is particularly important to model for its higher-precision photometry.

Because of the incomplete gyroscopic stabilization during the K2 mission, systematic effects on the light curve are substantial. We used several different light curves for WOCS 12009 as part of our investigation of systematic error sources. Our primary light-curve choice is the one from Nardiello et al. (2016b), which employed an effective point-spread function and high-resolution imagery of the field to identify and subtract neighbors (hereafter "PSF-based approach"; Libralato et al. 2016). Use of the Nardiello et al. detrended light curve or small photometric apertures generally involved larger numbers of discrepant data points or more complicated light variation around the primary eclipse. We chose to use their raw 4.5 pixel aperture photometry and further detrend the light curve using a 60-point median moving boxcar filter, corresponding to a 30 hr time span. This allows a small amount of contaminating light from other stars in the field into the aperture, but we allow for this contamination.

We experimented with two additional K2 light curves. The first is the K2SFF light curve derived using the method of Vanderburg & Johnson (2014) and 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 this was removed by fitting the out-of-eclipse points with a low-order polynomial and subtracting the fit. A few eclipse measurements (and others out of eclipse) were discarded from this light curve because of thruster firings and other events that produced discrepant points. The second additional light curve came from the K2SC pipeline (Aigrain et al. 2016), which models and removes position-dependent systematics and time-dependent trends from K2 light curves themselves. We used a 2 pixel aperture in this light curve.

Because the K2 light curves have high signal-to-noise, they have a substantial influence on the best-fit binary star parameter values (see Section 3.3). 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 take the scatter in the binary model parameter values as an indicator of systematic error resulting from the light curves.

In order to obtain color information for the binary, we obtained eclipse photometry in BVR_cI_c using several facilities. The earliest eclipse observations we identified in archival data were made of part of a secondary eclipse on 2003 April 21 (JD 2452750) at the Mount Laguna Observatory (MLO) 1 m telescope in the I_c band by one of us (ELS). After the discovery in K2 data, new observations were made at the MLO 1 m in BVR_cI_c filters. The observational details are given in Table 1. These images were processed using standard IRAF¹⁴ routines to apply overscan corrections from each image, subtract a master bias frame, and divide a master flat-field frame.

Archival observations of a secondary eclipse in V were found within data from the multi-telescope campaign on M67 reported by Stello et al. (2006). The only eclipse observations were made at the Danish 1.5 m telescope at the La Silla Observatory. To incorporate these observations into our analysis, we matched the out-of-eclipse level in the light curve with those of our other ground-based data sets, which involved a −0.004 mag shift.

The majority of our new ground-based eclipse observations were made with the telescopes of the Las Cumbres Observatory Global Telescope (LCOGT) Network. We mostly used 1 m telescopes, but a 0.4 m was used for one night of observation. The imagers on the 1 m telescopes came in two varieties: the SBIG (156 × 156 field; 4096 × 4096 pixels; 2 × 2 pixel binning) and Sinistro (26' × 26' field; 4096 × 4096 pixels) cameras. The SBIG camera on the 0.4 m had a 29' × 19' field and used binned 2 × 2 pixels. We subsequently used images that had been processed through the standard LCOGT pipeline in our photometric analysis below.

The brightness measurements were derived from multi-aperture photometry using DAOPHOT (Stetson 1987). We undertook a curve-of-growth analysis of the 12 apertures photometered per star in order to correct all measured stars to a uniform large aperture (Stetson 1990). We then corrected the data for each filter to a consistent zero point by using ensemble photometry (Honeycutt 1992; Sandquist et al. 2003). This essentially uses all measured nonvariable stars on the frame to determine magnitude offsets resulting from differences in exposure time, airmass, atmospheric transparency, and the like. Our implementation iteratively fits for position-dependent corrections that result from variations in point-spread function across the frame. These steps each brought noticeable reductions in the amount of scatter in the light curves. We found that an additional step of subtracting the light curve of a nearby bright nonvariable star also helped improve the agreement of zero points from night to night. In the case of WOCS 12009, we used the star WOCS 6008/Sanders 1242, which is about 1.3 mag brighter than WOCS 12009 and about 30'' away on the sky.

After the initial reductions, we found that there were significant color-dependent residuals in the light curves from different stars that resulted from the use of different telescope/camera/filter combinations. These were identified and characterized as linear relations using nonvariable stars in the field but were only judged to be significant in the B and I_c filters. The correction relations were used to correct all light curves to the system of the MLO filter set as part of the final ensemble photometry analysis. The specific corrections that were applied to the WOCS 12009 light curves were less than 0.008 mag for all observational setups except one in I_c (0.014 mag). The eclipse observations are shown in Figure 1, and the light curves are provided in tables in the electronic version.

Zoom In Zoom Out Reset image size

Radial velocities of M67 stars have been monitored over a very long period, and there is a large database of WOCS 12009 observations from different telescope/spectrograph combinations spanning almost 29 yr. The earliest 23 observations we use here employed the CfA Digital Speedometers (DS; Latham 1985, 1992), which covered 516.7–521 nm around the Mg i b triplet. Additional observations were taken on the 1.5 m Tillinghast reflector at Whipple Observatory using the Tillinghast Reflector Echelle Spectrograph (TRES; Fűrész 2008). The TRES spectra contained 51 echelle orders covering the whole optical range and into the near-infrared (from about 386 to 910 nm).

Ten spectroscopic observations were taken as part of the WIYN Open Cluster Survey (WOCS; Mathieu 2000) using the WIYN¹⁵ 3.5 m telescope on Kitt Peak and the Hydra multi-object spectrograph instrument (MOS), which is a fiber-fed spectrograph capable of obtaining about 80 spectra simultaneously (∼10 fibers set aside for sky measurements and ∼70 for stellar spectra). Observations used the echelle grating with a spectral resolution of 20,000 (velocity resolution ∼15 km s⁻¹). The spectra are centered at 513 nm with a 25 nm range that covers an array of narrow absorption lines around the Mg i b triplet. Spectroscopic observations for WOCS 12009 were completed using total integrations of either 1 or 2 hr per visit that were split into three subintegrations to allow for rejection of cosmic rays.

Image processing was done using standard spectroscopic IRAF procedures. First, science images were bias- and sky-subtracted, then the extracted spectra were flat-fielded, throughput-corrected, and dispersion-corrected. The spectra were calibrated using one 100 s flat field and two bracketing 300 s ThAr emission lamp spectra for separate integration times. A more detailed explanation of the acquisition and reduction of the CfA DS, TRES, and WIYN Hydra observations can be found in Geller et al. (2015).

We also collected two new spectra from 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 resolution R = 115,000 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 obtained archival spectra from two additional sources. Eight archived H-band infrared (1.51–1.70 μm) spectra with a spectral resolution R ∼ 22,500 for WOCS 12009 were available from the APOGEE project (Holtzman et al. 2015) on the 2.5 m Sloan Foundation Telescope. Before making our radial velocity measurements, we used APOGEE flags to mask out portions of the spectrum that were strongly affected by sky features, and we continuum-normalized the spectra using a median filter. Pasquini et al. (2008) tabulated three radial velocity measurements for the primary star as part of their survey of solar analogs in M67, based on observations from 2007 February using the FLAMES/GIRAFFE spectrograph at the VLT. Their spectra covered 644–682 nm, including the Hα line, and we obtained the spectra through the ESO Phase 3 Data Archive. We continuum-normalized the spectra that had been processed and calibrated through the ESO pipeline.

The radial velocities were derived using a spectral disentangling code following the algorithm described in González & Levato (2006). In each iteration step, the spectrum for each component is isolated by aligning the observed spectra using the measured radial velocities for that component and then averaging. This immediately deemphasizes the lines of the nonaligned stellar component, and after the first determination of the average spectrum for each star, the contribution of the nonaligned star is subtracted before the averaging to improve the spectral separation. The radial velocities can also be remeasured using the spectra with one component subtracted. This procedure is repeated for both components and continued until a convergence criterion is met. We use the broadening function (BF) formalism (Rucinski 2002) to measure the radial velocities using synthetic spectra from the grid of Coelho et al. (2005) as templates. The BFs can significantly improve the accuracy of radial velocity measurements in cases when the lines from the two stars are moderately blended (Rucinski 2002). In the case of the WOCS 12009 binary, the spectra of the two stars have significant differences, so that lines of each star can usually be identified even when the velocity separation is less than 10 km s⁻¹. As a result, we can derive reasonably reliable velocities for the two stars at small velocity separations when we couple the disentangling technique (especially the subtraction of the averaged spectrum of each star) with BFs.

The disentangling procedure requires spectra covering the same wavelength range taken at very different orbital phases, and so we had to work with different subsets of our spectra. The CfA DS spectra covered a wide range of orbital phases but had a limited wavelength range, so we disentangled these spectra separately. Many of the spectra had comparatively low signal-to-noise ratios that made it impossible to identify the secondary component. We used all spectra in determining the average spectrum for the primary but only used the eight spectra with the most clearly detected secondary component to determine its average spectrum. We were only able to measure reliable velocities for the secondary star from these spectra and two others.

We disentangled the TRES and HARPS-N spectra together because they covered similar wavelength ranges at similar resolution, and it would have been difficult to disentangle the two HARPS-N spectra on their own. Because the TRES and HARPS-N spectra covered the entire optical range, we decided to disentangle them in the wavelength ranges of individual TRES echelle orders. We selected 14 orders with a large number of strong lines, avoiding the bluest and reddest ends of the optical range (due to low signal on the blue end and telluric lines on the red end). This had the additional benefits that we could determine the stellar velocities as means of the velocities derived from the different orders and also get uncertainty estimates from the standard error of the mean. These uncertainties were typically around 0.08 km s⁻¹ for the primary star and 0.4 km s⁻¹ for the secondary.

The WIYN spectra were disentangled on their own because we were able to obtain spectra covering a sufficient set of different phases, although we had to take some care in the selection of spectra used to calculate the averaged spectra for the two stars. We used eight of the 10 spectra for the brighter primary star and six of the best for the secondary. Uncertainties were initially derived from the BF fits and scaled to produce a reduced χ² value of 1 based on a best Keplerian fit to the WIYN radial velocities alone. Typical uncertainties were around 0.15 km s⁻¹ for the primary star and 1 km s⁻¹ for the secondary star.

The APOGEE spectra also had to be disentangled on their own because their spectral range did not overlap with the other data sets. 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.2 km s⁻¹ for the primary star and 0.7 km s⁻¹ for the secondary star. The velocity was corrected to the barycentric system using values calculated by the APOGEE pipeline. We compared our measured velocities for the primary star with those of the APOGEE pipeline, and we found a phase-dependent difference with a sign consistent with the idea that the tabulated APOGEE velocities were pulled slightly (0.5–1.0 km s⁻¹) by the presence of the unmodeled secondary star.

We also ended up disentangling the three FLAMES/GIRAFFE spectra (Pasquini et al. 2008) separately. For two of the three spectra, we clearly detected and resolved BF peaks for the secondary star on opposite sides of the primary star peak. For the third spectrum, the BF peaks were moderately blended. However, because the two unblended spectra were at very different orbital phases, we were successfully able to disentangle the spectra in this set and were able to derive velocities for both stars in all three spectra. Reliable uncertainties could not be derived from scatter around a best-fit solution for these measurements on their own, however. As a result, we derived uncertainty estimates from their scatter around the best fit to the entire velocity data set.

Finally, we employed one radial velocity measurement from Jones et al. (1999) made in 1996 March using the HIRES spectrograph at the Keck 10 m telescope. Because the secondary star was not subtracted from the spectrum before their measurement, it is possible that there is a small bias in the direction of the system velocity for that measurement.

Corrections were applied to the velocities from different spectroscopic data sets. For the CfA DS velocities, we used run-by-run corrections derived from observations of velocity standards. Run-by-run offsets measured for the TRES spectra were no more than 17 m s⁻¹, and although we did apply them, they were much smaller than the measurement uncertainties. Velocity offsets for individual Hydra fibers in the WIYN data set were determined and applied to the WOCS 12009 spectra, although these were of a different sign and no more than about ±0.2 km s⁻¹ in general.

Because the synthetic spectra used in disentangling do not account for the gravitational redshifts of the stars, our velocities should have an offset relative to velocities derived using an observed solar template. For consistency with the large survey of Geller et al. (2015), we subtracted the combined gravitational redshift of the Sun and gravitational blueshift due to the Earth (0.62 km s⁻¹) from the velocities we tabulated. We did not apply the gravitational redshift correction to the Jones et al. Keck HIRES velocity because it was derived using a solar template. 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⁻¹ versus 0.64 km s⁻¹ for the Sun), but we did not correct for these small differences. We did, however, allow for different system velocities for the two stars in our fits to allow for effects like this and convective blueshifting that could otherwise affect the measured radial velocity amplitudes K₁ and K₂ (and the stellar masses) if not modeled. The measured radial velocities are presented in Table 2.

Even with these corrections to the velocities, we are attentive to the possibility that systematic differences might exist from data set to data set. The phased radial velocity measurements are plotted in Figure 2, showing the best fit to the orbit using all measurements. Among the measured primary star velocities, mean residuals (observed minus computed) were fairly small: 0.18 km s⁻¹ for WIYN, 0.05 km s⁻¹ for TRES, −0.07 km s⁻¹ for HARPS-N, −0.23 km s⁻¹ for APOGEE and FLAMES, and −0.45 km s⁻¹ for CfA DS. Because we do not see clear evidence of systematic shifts between spectroscopic data sets, we have not applied any corrections of this sort. We have opted not to scale the uncertainties on the data sets individually but instead scale uncertainties on all velocity measurements for each star to separately produce a reduced χ² measurement of 1. This preserves information on the relative precision of measurement in each data set (which affects their weighting for the final fits) that was most often determined a posteriori. The primary star uncertainties were scaled upward by a factor of 2.1, and the secondary star uncertainties were scaled by a factor of 1.38.

Zoom In Zoom Out Reset image size

For the purpose of later discussions of the reliability of the stellar mass and radius measurements, we summarize here that we have consistent radial velocity information on both stars from five different spectrographs in a variety of optical and infrared wavelength ranges. This gives us confidence that systematic errors in the radial velocities do not significantly affect the masses and radii we will derive.

There have been many studies of M67 that have included proper-motion measurements and estimations of membership probabilities. All of the studies agree that WOCS 12009 is a highly probable member. The lowest published membership probability (71%) comes from Zhao et al. (1993), but all others are greater than 93% (Sanders 1977; Girard et al. 1989; Yadav et al. 2008; Krone-Martins et al. 2010; Nardiello et al. 2016a). In our later use of proper-motion membership probabilities to produce a CMD cleaned of field stars, we made use of all of these sources, computing a straight average of all of the probabilities published for each star to assign a final probability.

Geller et al. (2015) also published membership probabilities based on their radial velocity survey, and again WOCS 12009 was assigned a high probability of membership (98%). Our modeling of the radial velocities presented below indicates that both stars in the binary have high membership probabilities based on the system velocities we derive: γ₁ = 33.81 and γ₂ = 33.75 km s⁻¹. Geller et al. found the mean cluster velocity to be 33.64 km s⁻¹ with a radial velocity dispersion of 0.59^+0.07_−0.06 km s⁻¹.

Based on all of the kinematic evidence, the binary is a high-probability cluster member. Additional evidence comes from the distance moduli derived from the radius, T_eff, and V magnitude of the stars (see Section 4.2) and from the position of the components of the binary in the CMD once the light from the two stars has been disentangled (see Section 4.3).

Surface temperature is somewhat important for constraining the limb-darkening coefficients to be used in the eclipse modeling (Section 3.3), but it is even more important for the distance modulus determination in Section 4.2.

Photometric determinations of T_eff using the total system photometry (see Table 3) will be affected by the contributions of the secondary star, but the calculation will be a lower limit to the temperature of the primary star and closer to the actual temperature the larger the fraction of the system light the primary contributes. Jones et al. (1999) found 5889 K from BV photometry. Using the Stromgren photometry of Balaguer-Núñez et al. (2007) and the T_eff calibration of Önehag et al. (2009), we calculate 5920 K from the (b − y) color and a reddening E(b − y) = 0.030. The APOGEE pipeline provides a temperature of 5950 ± 70 K from its best fit to its infrared H-band spectra using synthetic templates, but again, this will be somewhat affected by the light of the secondary star, which is more prominent in the infrared.

With luminosity ratios determined from the binary star analysis below, we are able to separate the light of the two stars and estimate temperatures for them. For this purpose, we used photometry from Sandquist (2004) in the optical along with the color–T_eff relations of Casagrande et al. (2010). Assuming a reddening E(B − V) = 0.041 ± 0.004 (Taylor 2007), the weighted averages of values computed from (B − V) and (V − I_C) were 6110 and 4760 K for the primary and secondary stars, respectively. These values are in good agreement with temperatures in the synthetic spectra that optimized our detection of BF signal for radial velocity measurement (6000 and 6250 K for the primary and 4500 and 4750 K for the secondary).

To simultaneously model the ground-based radial velocities and photometry and Kepler K2 photometry, we used the ELC code (Orosz & Hauschildt 2000). We assumed the stars are spherical, and we did not attempt to model spots. The core of the model fit involved a set of 14 parameters: orbital period P, time of conjunction (primary eclipse) t_c, velocity semi-amplitude of the primary star K₁, mass ratio q = M₂/M₁ = K₁/K₂, systematic radial velocities¹⁶ γ₁ and γ₂, the combinations of eccentricity (e) and argument of periastron (ω) $e\cos \omega$ and $e\sin \omega$, inclination i, ratio of the primary radius to average orbital separation R₁/a, ratio of radii R₁/R₂, temperature of the primary T₁, temperature ratio T₂/T₁, and contamination of the K2 light curve by other sources (independent of the BVR_CI_C light curves). The K2 contamination parameter is used to account for possible dilution of the light curves by stars that are not physically associated with the binary. Kepler pixels are large by the standards of most ground-based cameras, and in the M67 cluster field, light contamination is a greater concern than in the majority of the field of view. The two stars of greatest significance here are II-216 (V = 19.45) and WOCS 6008/Sanders 1242 (V = 12.33), but their distances (23'' and 286, respectively) will mean that their contributions are relatively small. Because the K2 light curves use long-cadence data with 30 minute exposures, they also effectively measure the average flux during a 30 minute exposure. To account for this, the ELC code integrates the computed light curves in each observed exposure window.

Limb darkening is a continuing uncertainty in binary star fits, and to model this, we fit for 18 additional parameters using the Kipping (2013) algorithm: two coefficients in the quadratic limb-darkening law for each star in the Kepler and BVI_C filter bands (and two coefficients for the secondary star in R_C because to date, we only have secondary eclipse observations). Briefly, Kipping's algorithm has the advantages that it reduces the correlation between limb-darkening-law parameters and allows for a straightforward sampling of physically realistic values. We modified the Kipping (2013) algorithm by restricting the second parameter u₂ to values less than 0.5 in order to force the limb-darkening law to be concave down as it approaches the limb. By systematically exploring the parameter space, we can include our limb-darkening uncertainties in our uncertainties of the stellar parameters.

The quality of the model fit was quantified by an overall χ². We applied the temperature of the primary star T₁ =6100 ± 100 K (see Section 3.2) as an observed constraint, allowing the value to vary but applying a χ² penalty based on deviation from the determined value. Because the photometric and spectroscopic data sets were pulled together from a variety of sources, the relative weights of different data influence the determination of χ². For radial velocity measurements, we utilized empirical estimates of measurement uncertainties as much as possible, as described in Section 2.3 for each spectroscopic data set. These empirical measurements were then scaled separately to give a reduced χ² of 1 for the primary and secondary star data sets. Although some binary star parameters are derived almost exclusively from the velocities, others (like eccentricity e and argument of periastron ω) are affected by the photometry as well, so it is important to try to weight the velocity measurements properly with respect to the photometry. For the photometry, we scaled all of the a priori measurement uncertainties upward by the same factor in order to produce a reduced χ² value of 1 based on separate best fits to the data for each filter. Once this was done, we computed simultaneous fits to the photometry and velocity data sets.

The model parameter space was sampled using an implementation of the parallel emcee algorithm (Foreman-Mackey et al. 2013) for Markov chain Monte Carlo. This algorithm allows for efficient sampling of alternate parameter sets using a group of "walkers" (in our case, 1000) that move in parameter space around the values providing the best fit to the data. A proposal for a jump in the position of one walker is taken along a line connecting to but extending beyond another randomly chosen walker. The proposed position is accepted with a probability that decreases as χ² increases. We found that the walkers converged to stable distributions within 2000 steps, and these steps were discarded before forming posterior probability distributions from chains of 6000 step length, skipping every 15 steps. Approximate 1σ parameter uncertainties were derived from the parts of the posterior distributions containing 68.2% of the remaining models from the Markov chains. The results of the parameter fits are provided in Table 4. In the last column, we provide the consensus parameter values for the most important quantities, including a statistical error estimate derived from the results of models in and around the best fit, as well as a systematic error estimate derived from the scatter among the runs with the three different K2 light curves. When there was disagreement, the K2SFF results were generally the outliers. This was particularly noticeable in the luminosity ratios and in the secondary star radius and gravity. We attribute this mostly to the slightly deeper eclipses in the K2SFF light curve and, to a lesser extent, the clearly discrepant points during the one observed primary eclipse that had to be deleted. For the radii and gravities, the combined values were derived from an average of the results for all three K2 light curves, while we averaged the results for the PSF-based approach and K2SC light curves for the luminosity ratios. The systematic error quotes reflect these differences.

Table 4.  Best-fit Model Parameters for WOCS 12009

Parameter	K2 Light Curves Law			Combined
	PSF-based Approach LC	K2SC LC	K2SFF LC
T1 (K)	6100 ± 100 (constraint)
γ1 (km s−1)	33.815 ± 0.009	33.814 ± 0.009	33.815 ± 0.009	33.815 ± 0.009
γ2 (km s−1)	33.745 ± 0.010	33.746 ± 0.009	33.745 ± 0.010	33.745 ± 0.009
P (day)	69.728087	69.728107	69.728117	⋯
σP (day)	0.000014	0.000013	0.000014	⋯
tC − 2450000	7180.14126	7180.14068	7180.14063	⋯
${\sigma }_{{t}_{C}}$	0.00012	0.00012	0.00016	⋯
i (deg)	89.565 ± 0.019	89.568 ± 0.016	89.637 ± 0.013	⋯
q	0.674 ± 0.004	0.674 ± 0.004	0.674 ± 0.003	0.674 ± 0.004
$e\cos \omega$	−0.036166 ± 0.000009	−0.036126 ± 0.000010	−0.036134 ± 0.000005	⋯
$e\sin \omega$	0.0377 ± 0.0016	0.0376 ± 0.0016	0.0361 ± 0.0016	⋯
e	0.0522 ± 0.0012	0.0522 ± 0.0011	0.0510 ± 0.0011	⋯
ω (deg)	133.7 ± 1.1	133.8 ± 1.2	135.1 ± 1.3	⋯
K1 (km s−1)	25.64 ± 0.04	25.65 ± 0.04	25.63 ± 0.04	25.64 ± 0.04
K2 (km s−1)	38.06 ± 0.19	38.06 ± 0.20	38.05 ± 0.22	38.06 ± 0.20
R1/R2	1.447 ± 0.042	1.460 ± 0.035	1.537 ± 0.020	⋯
R1/a	0.01221 ± 0.00007	0.01221 ± 0.00008	0.01225 ± 0.000014	⋯
R2/a	0.00844 ± 0.00025	0.00836 ± 0.00021	0.00797 ± 0.00016	⋯
T2/T1	0.770 ± 0.005	0.772 ± 0.006	0.772 ± 0.006	⋯
K2 contam.	0.031 ± 0.004	0.0027 ± 0.0029	0 + 0.0003	⋯
L2/L1(B)	0.0694 ± 0.0049	0.0686 ± 0.0043	0.0655 ± 0.0028	0.0690 ± 0.0049 ± 0.0004(sys)
L2/L1(V)	0.1081 ± 0.0067	0.1076 ± 0.0055	0.0940 ± 0.0025	0.1079 ± 0.0067 ± 0.0003(sys)
L2/L1(R)	0.1439 ± 0.0090	0.1420 ± 0.0070	0.1209 ± 0.0033	0.1429 ± 0.0090 ± 0.0010(sys)
L2/L1(I)	0.1724 ± 0.0099	0.1734 ± 0.0096	0.1554 ± 0.0052	0.1729 ± 0.0099 ± 0.0005(sys)
M1/M⊙	1.109 ± 0.015	1.112 ± 0.013	1.111 ± 0.014	1.111 ± 0.015
M2/M⊙	0.749 ± 0.005	0.749 ± 0.005	0.747 ± 0.006	0.748 ± 0.005
R1/R⊙	1.070 ± 0.008	1.070 ± 0.007	1.073 ± 0.006	1.071 ± 0.008 ± 0.003(sys)
R2/R⊙	0.739 ± 0.019	0.732 ± 0.016	0.699 ± 0.010	0.713 ± 0.019 ± 0.026(sys)
$\mathrm{log}{g}_{1}$ (cgs)	4.425 ± 0.006	4.426 ± 0.006	4.422 ± 0.005	4.424 ± 0.002
$\mathrm{log}{g}_{2}$ (cgs)	4.575 ± 0.021	4.582 ± 0.019	4.624 ± 0.011	4.607 ± 0.021 ± 0.032(sys)

Download table as:  ASCIITypeset image

The masses and radii of the eclipsing stars are important because they are high-precision, distance- and reddening-independent measurable characteristics. In the case of M67, these kinds of measurements can potentially provide strong constraints on the cluster age.

In Figure 3, we compare the results of our binary star model fits with theoretical expectations from isochrones, and in Figure 4, we compare against precisely measured members of eclipsing binaries (Southworth 2015).¹⁷ For the comparisons in the figures, the error ellipses for WOCS 12009 include contributions from statistical and systematic errors added in quadrature, and we have employed the chemical compositions used by the modelers in calibrating their isochrones to the Sun. The properties of the secondary star (WOCS 12009 B) are in rough agreement with isochrone predictions, although there is some dependence on the chemical composition and the implementation of the stellar physics in the models. Also of note, the secondary star is slightly closer to the isochrones than any other star in the range from about 0.7 to 0.8M_⊙. This is probably because WOCS 12009 B is in a much longer period binary than other measured main-sequence stars. For example, V404 CMa (Rozyczka et al. 2009) has a period of just 0.452 day, and ASAS J082552–1622.8 (Hełminiak & Konacki 2011) has a period of 1.528 days. Tidal interactions and the resulting stellar activity are likely to be responsible for the relatively large radii in systems like these (see Torres et al. 2010 and references therein). WOCS 12009 B is probably a better approximation to a star that has evolved in isolation and is therefore valuable in its own right.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On the other hand, while the primary star (WOCS 12009 A) has characteristics putting it among other well-measured eclipsing stars in the field, it is much smaller than predicted for its mass and reasonable ages for M67. Isochrone fits to the CMD have returned ages from about 3.5 to 4.0 Gyr (Sarajedini et al. 2009) and 3.6 to 4.6 Gyr (VandenBerg & Stetson 2004), with a large part of these ranges resulting from model physics and chemical composition uncertainties. Stello et al. (2016) quoted an age of about 3.5 Gyr from their asteroseismic analysis of pulsating giants in M67. The disagreement between WOCS 12009 A's radius and the predictions from models with ages at the low end of the range quoted for M67 is more than 15%. The star's radius is larger than that of a zero-age main-sequence star and is consistent with a star of approximately 1.4 Gyr age.

To illustrate the difficulties in explaining the mass and radius of this star consistently with the cluster age, we explored one-dimensional stellar evolution models (Figure 5). Stellar models were calculated using the 1D Yale Rotating Evolution Code (YREC; Pinsonneault et al. 1989) using the nuclear cross-sections, opacities, equations of state, and model atmospheres described in Section 2.3 of Somers & Pinsonneault (2016). For the base model, we adopt the Asplund et al. (2009) present-day solar photospheric abundance ratio Z/X = 0.0181, assume helium and heavy metal diffusion as calculated by Thoul et al. (1994), calibrate the solar mixing length and helium abundance to reproduce the solar luminosity and radius at 4.57 Gyr, and assume no rotation. As parameter variations, we test the alternate solar mixture given by Grevesse & Sauval (1998), assume diffusionless envelopes, adopt different proto-solar helium abundances, vary the convective efficiency by changing the mixing length (e.g., Chabrier et al. 2007), and permit modest core overshooting of ∼0.2 pressure scale heights.

Zoom In Zoom Out Reset image size

The secondary star's characteristics are less sensitive to uncertainties in age, stellar physics, and composition than the primary star's, so the overall agreement of models and observations of the secondary is a useful validation. Reasonable changes to the chemical composition of the stars are unable to explain a radius discrepancy of this magnitude, although higher metallicity or lower helium abundance would push the models in the right direction. In the case of helium alone, the abundance would have to be moved in a direction opposite what is observed for galactic chemical evolution and reduced below Big Bang nucleosynthesis values (Y = 0.245) to come close to the observed radius. The metallicity of M67 stars has been repeatedly shown to be very close to solar, so to explain the radius discrepancy with metallicity would require a change several times any reasonable uncertainty in M67's value.

Stellar physics uncertainties also seem to fall short of explaining the small radius. The inclusion of diffusion makes a model star marginally larger, so this cannot alleviate the issue. More efficient convection (through a larger mixing-length parameter α_ML, for example) can produce a smaller star, but again, an explanation of the radius issue would require a large change. For a star with a thinner surface convection zone than the Sun, convection is expected to be less efficient due to radiative losses. Convective core overshoot tends to make the star slightly smaller, but even for overshoot several times larger than inferred from cluster comparisons (0.2 pressure scale heights), this effect is also too small by a factor of at least three.

To summarize, we find no chemical composition or stellar physics uncertainties that can explain the small size of the star WOCS 12009 A.

We can further test the membership of the binary by deriving the luminosity of each star and comparing the absolute with the apparent magnitudes. The photometry of the two stars was deconvolved using luminosity ratios from the multi-filter light-curve analysis, and photometric temperatures were derived from the colors and cluster reddening (see Section 3.2). The temperatures, along with radii derived from the binary analysis and theoretical bolometric corrections from Casagrande & VandenBerg (2014), allow us to derive the absolute magnitudes M_V. We find (m − M)_V = 9.74 ± 0.03 and 9.82 ± 0.09 for the primary and secondary stars, respectively. The quoted uncertainties include contributions from the measurement of the apparent magnitudes, the calculation of the luminosity, and the calculation of the bolometric correction. These values are in excellent agreement with each other and with previous measurements of M67's distance modulus. For comparison, some notable recent measurements are 9.76 ± 0.06 ± 0.05 (statistical and systematic uncertainties) from 10 solar analogs (Pasquini et al. 2008), 9.72 ± 0.05 using metal-rich field dwarfs with Hipparcos parallaxes (Sandquist 2004), and 9.70 ± 0.04 from K2 asteroseismology of more than 30 cluster red giants (Stello et al. 2016).

With well-studied binary stars, we can bring important and rarely available new mass information to the study of the CMD and the comparison of observational data with theoretical models. To that end, we used luminosity ratios derived from the binary star modeling to deconvolve the light in different wavelength bands for the two WOCS 12009 stars and determine their CMD positions. In the discussions below, we use the high-precision photometry of Sandquist (2004), along with the identifications for high-probability single-star members from that study. The original list was selected based on being at the blue edge of the band of main-sequence and binary stars in the CMD, but we have verified that the majority of the list were classified as single members with no detectable radial velocity variability by Geller et al. (2015). Those that were identified as binaries were eliminated. So, to the greatest degree we can, these stars should trace out where the main sequence and turnoff should be for the age of M67.

Because the statistical uncertainties in the system photometry are well below 0.01 mag and the system photometry has been calibrated in the same way as the rest of the cluster stars we use in the CMDs, the relative photometry should be reliable. In the CMDs (see Figure 6), the secondary star falls among probable single main-sequence stars at V ≈ 16.5. The primary star is slightly bluer than the main sequence at V ≈ 14.1, as we would expect for a star that is less evolved than the other cluster stars.

Zoom In Zoom Out Reset image size

When we shift isochrones to force a fit to the mass and photometry of the secondary star (see Figures 7 and 8), there is decent agreement between the cluster main-sequence stars and the theoretical predictions for the Victoria–Regina and MIST isochrone sets. The landscape is similar for the (V − I_C, V) CMD (see Figure 8), although the MIST isochrones do not match the main sequence as well.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The mass of the primary star would be a more valuable constraint on the masses of evolved stars in the cluster and on the cluster's age if it had a conventional evolution history. However, if we utilize the model isochrones in a relative sense, we can infer with some confidence that a single cluster star with the mass of WOCS 12009 A (about 1.11M_⊙) and an uneventful evolution history should reside at V ≈ 13.8, about 0.35 mag brighter than WOCS 12009 A and just faintward of the "hook" feature in the CMD. Models tied to WOCS 12009 B also predict that a solar-mass star in M67 should have V ≈ 14.57. This is quite close to the V magnitude of the star Sanders 1095 (1787 in Yadav et al. 2008) that Pasquini et al. (2008) identified as one of their two best solar analogs.

In the CMD, the turnoff of M67 potentially gives us a lot of insight into the physics of the stellar cores. Our group is currently working on additional binary systems that will allow us to identify the precise masses of stars at the turnoff, so we will not undertake a full discussion of the topic here. The CMDs in Figures 7 and 8 illustrate the substantial differences in the isochrones produced by different groups, and, without clarity on the "correct" physics and chemical composition inputs to include in the models, systematic errors in the cluster age will be unavoidable.

The noticeable gap seen at V ≈ 12.9 is an indicator of rapid evolution through that part of the CMD, probably due to the presence of a small convective core maintained by CNO cycle hydrogen-burning reactions and their strong temperature dependence. However, there are several composition and physics factors that can affect the size (and even the presence) of a convective core in models. VandenBerg et al. (2007) and Magic et al. (2010) studied extensive sets of stellar models to try to identify the important factors, but Magic et al. in particular found that there are degeneracies in how the CMD is affected. Composition factors like the metallicity of M67 relative to the Sun and the actual solar mixture affect the abundance of CNO cycle catalyst nuclei, while physics like diffusion, the prescription for convective core overshoot, and CNO cycle reaction rates can affect the strength of the energy release in the stellar cores. While some of the models approximately reproduce the brightness of the gap (which occurs in the straight-line segment immediately faintward of the bluest point on the isochrone), models generally are not as successful in matching the blue protrusion among the M67 main-sequence stars at $V\approx 13.7$. To the extent that we are able to check with current models, the cluster CMDs are consistent with typically quoted ages of 3.8–4.2 Gyr.

Every indication we have is that the binary WOCS 12009 is a member of the M67 cluster, but the primary star is significantly smaller than expected for reasonable cluster ages, and it is also hotter and less luminous than expected for its mass. Before discussing possible explanations, we assemble other observations that can comment on the puzzle.

Using our highest-resolution spectroscopic observations (with the TRES and HARPS-N spectrographs), we looked at the rotational motions of the two stars in the binary but were unable to detect rotation in the BFs of the two stars. The upper limits are ${v}_{\mathrm{rot}}\sin i\lt 5$ km s⁻¹ from the TRES observations and <4 km s⁻¹ from the HARPS-N observations. Although slow rotation is expected for single solar-type stars in an old cluster like M67, the low rotation speed of the primary has some impact on our understanding of its genesis, as discussed later.

The mass of WOCS 12009 A should place it at the faint end of the cluster's Li plateau with abundances A(Li) ≈ 2.5 (Pace et al. 2012). Stars in the plateau have outer convective zones that do not mix Li atoms 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 have been detectable if it was present with the plateau abundance. Our spectral disentangling allows us to mostly remove the secondary star's contribution to the spectrum and to produce an average spectrum for the primary star as a by-product. The averaged FLAMES spectrum for the primary is shown in Figure 9, as it has the highest signal-to-noise ratio in that wavelength range. There is still no sign of the Li lines for the primary star. Jones et al. (1999) provided an upper limit for the system (A(Li) < 1.32) that puts the star's abundance at least about 0.8 dex lower than stars at the same absolute magnitude, but because the star is probably underluminous compared to normal cluster stars of the same mass (see Section 4.3), the discrepancy is more than an order of magnitude (Figure 10).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Li observations could be explained by a stellar union: stars of about 0.95 M_⊙ or less would have depleted their surface Li to undetectability well before the present age of the cluster because the surface convection zone easily reaches down to Li burning temperatures.¹⁸ In the union of two such stars, a more massive star with depleted Li would result.

This could also explain the smaller-than-expected radius of the primary star. In hydrodynamic simulations of stellar collisions, it is found that the material from the two stars gets sorted radially in the merger remnant by entropy to good approximation (Lombardi et al. 2002). Because lower-mass stars have lower temperatures and higher densities in their cores, the core of the lower-mass star ends up at the core of the remnant. This effectively restarts the remnant with a higher core hydrogen content and makes it look like a younger main-sequence star once a Kelvin–Helmholtz timescale has passed and the remnant has thermally adjusted to its new structure. As a result, if the star's size makes it look like it is 1.4 Gyr old now, that is likely to be fairly close to the time elapsed since a merger took place.

Two stars could combine as part of the merger of a close binary pair or the collision of two stars during a resonant multi-star interaction. Both are likely to produce a remnant that is depleted of surface Li, but the route for combining the stars could leave signatures on their motions. Is there a reasonably probable mechanism that could produce a slowly rotating remnant star in a wide low-eccentricity (e = 0.05) binary?

Unless the combining stars were in a nearly head-on collision, the remnant star is likely to have started with rapid rotation. The primary star in WOCS 12009 probably has a surface convection zone, though, meaning that some spindown via a magnetic wind is likely. If the remnant was made approximately 1 Gyr ago, we can ask whether there was sufficient time for it to spin down as an effectively single star. The rotation of stars in the 1 Gyr old cluster NGC 6811 was studied during the Kepler mission (Meibom et al. 2011). In that cluster, stars with (B − V) color similar to that of WOCS 12009 A have spun down from their zero-age main-sequence values to a period of around 10 days. For a star the size of WOCS 12009 A, a 10 day period roughly corresponds to rotation near our ${v}_{\mathrm{rot}}\sin i$ detection limit of 5 km s⁻¹. We computed rotational stellar models using YREC (see Section 4.1) to check this idea further, as shown in Figure 11. A 1.1 M_⊙ star has enough of a surface convection zone that spindown via a magnetic wind could bring a rapidly rotating star (modeled here with a period P_rot = 0.5 day) down to our detection limit in 1.5 Gyr. A detected rotation period for the primary, as opposed to our current upper limit, could, in principle, be used to infer the time since merger, for the same reason that a rotation period for a single noninteracting star can be used as an age indicator. This would be one interesting test of the merger hypothesis.

Zoom In Zoom Out Reset image size

It is plausible that a rapidly rotating merger remnant could have spun down to the levels we see today, but other observations may test this. For example, the Ca ii H and K emission index (an activity indicator) was measured for WOCS 12009 by Giampapa et al. (2006), and they found that it has a relatively high value relative to other M67 stars but is not out of the range of the best-vetted single-star candidates (Curtis 2017). Mamajek & Hillenbrand (2008) computed $\mathrm{log}R{{\prime} }_{\mathrm{HK}}=-4.70$ from the Giampapa et al. data, which probably makes it more active than the 95% range for the Sun. An examination of spot modulation on the primary star might also be able to reveal whether it is rotating significantly more rapidly than single cluster stars of similar color. Such measurements are currently underway under the umbrella of the K2 M67 Study.

The other major question is whether the orbital characteristics of the current binary could be natural products of any potential formation mechanisms. Below, we consider stellar collisions and the evolution of hierarchical triple systems.

A stellar collision during a multibody dynamical interaction is a relatively common outcome if the incoming binary or binaries are tightly bound. Fregeau et al. (2004) examined binary–single and binary–binary scattering experiments with collisions allowed and found, as expected, that low eccentricities are strongly disfavored in binaries formed after stellar collisions: the conditions required for the outgoing stars to wind up on a nearly circular orbit are too narrow. Even with a maximal amount of expansion by a merger remnant before it thermally relaxes, tidal interactions are unlikely to be able to help circularize the orbit of WOCS 12009, which is currently relatively wide (a/R₁ ≈ 82).

Formation involving a hierarchical triple system that dynamically formed as part of a multibody interaction along with the subsequent merger of the close pair faces some of the same difficulties. For guidance, we examined the scattering study of Antognini & Thompson (2016). They found that when triple systems are dynamically formed, a compact configuration (with the ratio of the outer orbit periastron to the inner orbit semimajor axis r_P,outer/a_inner ∼ 10) is a preferred outcome. If the current orbit of WOCS 12009 corresponds to the outer orbit of such a compact triple, the inner orbit could have had a period of around 2 days. This is long enough that the inner binary would not have merged rapidly via tidal interactions or magnetic braking. However, an eccentricity as low as 0.05 for the outer orbit is again strongly disfavored as part of the dynamical formation of the triple.

A primordial triple system with an inner binary that merges after about 3 Gyr is also a possibility. First, considering the low eccentricity of the current orbit, we note that low-eccentricity binaries with periods above 10⁴ days are known among solar-type binaries in the field (Raghavan et al. 2010).¹⁹ A few percent of binaries with orbital solutions in the 2.5 Gyr old cluster NGC 6819 (Milliman et al. 2014), M67 (M. Pollack et al. 2017, in preparation), and the 7 Gyr old cluster NGC 188 (Geller et al. 2009) are known to have e < 0.1 and P > 20 days. (For NGC 6819 and M67, there are known long-period systems with e < 0.06). For known triple systems, the outer orbit is usually eccentric, but there are known examples with e < 0.05. So it is not out of the question to think that the low eccentricity can result from the star formation process.

The inner binary might have formed with an initial period near 2 days as part of the star formation process, allowing angular momentum loss mechanisms like tides and magnetic braking (e.g., Stȩpień 2011) to produce a merger. The majority of the time (∼2–3 Gyr) would be needed to have the binary lose enough energy to come into contact and ultimately merge. The M67 systems AH Cnc, EV Cnc, and HS Cnc are contact or near-contact binaries near the cluster turnoff (Stassun et al. 2002; Sandquist & Shetrone 2003; Pribulla et al. 2008; Yakut et al. 2009) that could be in this stage of their orbital evolution, and theoretical studies indicate that sub-turnoff mergers could be of order 3%–4% of the upper main-sequence population for an old cluster like M67 (Andronov et al. 2006).

A primordial triple with an inner binary having a period substantially longer than 3 days might also be capable of producing a system like WOCS 12009, but stability against disruption is a concern if the inner binary has too long a period. As a guide, Mardling & Aarseth (2001) used a criterion for stability against escape of a component star that implies r_P,outer/a_inner ∼ 3.4 for the WOCS 12009 system, which translates to an inner orbit with a_inner ≈ 0.12 au or P_inner ≈ 14 days. For an inner binary formed with an initial period closer to 14 days, Kozai–Lidov oscillations due to the third star would probably be needed to drive the binary's eccentricity up before tidal dissipation could produce a shorter period (Fabrycky & Tremaine 2007). This would require misalignment of the inner and outer orbital planes, but Kozai–Lidov cycles can have a very short timescale (comparable to but longer than the period of the outer orbit), so there would be no difficulty with producing a close binary within the age of M67. Once close binary interactions become dominant, Kozai cycles are suppressed, and the binary is expected to evolve as if it is isolated from the third star.

For a stable primordial triple formed with an outer orbital period near 70 days, either of these paths would likely result in a circularized inner binary with a period of a few days, which will later merge. If a Kozai mechanism was responsible for the initial stage of decreasing the separation of an inner binary system, we might be able to observe the signature of an inclination difference between the orbit of the inner binary and the orbit of the tertiary star (Fabrycky & Tremaine 2007). The rotation axis of the merged star would probably retain the direction of the angular momentum vector for the close binary. Because the current system reflects the orbital plane of the outer binary in the original triple, and because the binary is currently eclipsing, the rotation axis of the primary star could be inclined to our line of sight by about 40°. Thus, depending on its axial orientation, it might be possible to observe muted rotational effects, such as the Rossiter–McLaughlin effect, on the radial velocities or rotational modulation of the light curve.