© The Authors 2023

1. Introduction

A key issue in cool star science is to understand the extent to which the rotation rate of a star can be used to infer its age. While empirical data are available in a series of young open clusters, they are lacking for cooler spectral types in older clusters. Here we confirm related prior results and extend the empirical knowledge of this extent to K- and M-type stars of age 4 Gyr by presenting measured rotation periods in the open cluster M 67.

The fact that the rotation rate of a star depends on its age was originally proposed by Skumanich (1972) following v sin i measurements for stars of solar mass in open clusters. Older stars were found to rotate slower than younger ones, that means the stars spin down as they age, approximately as a power law. This relationship was later extended to other late-type stars by Barnes (2003) who showed the existence of a mass-dependence of the spindown rate and suggested that the age of a star can be inferred from its measured rotation period and mass (or a suitable proxy of the latter, such as color). This marks the onset of what is now known as gyrochronology.

Using stellar rotation for age estimates of late-type main sequence stars opens up new possibilities in age dating for stars whose classical parameters change only very little over the course of billions of years (e.g., temperature and luminosity), which can only be used reliably in young stars (e.g., Li abundance), or which change as long-term averages but vary significantly on short timescales (e.g., activity). However, to be able to use gyrochronology reliably, it needs to be first established if and how stellar rotation changes systematically with age for stars of different masses. Now, if the rotation periods of stars change systematically with age, then regardless of the details of the dependencies on stellar mass and age, appropriate measurements of open cluster stars of known ages can be used to set up a series of “mileposts” at suitable ages to derive stellar ages from their measured rotation periods. This has led to the exploration of accessible open clusters of differing ages to construct the sequence of spindown empirically as a function of mass and age. Open clusters are the preferred means to explore the spindown as they provide sets of stars of different masses with well-established ages and metallicities.

Fortunately, the relationship has so far been found to be single-valued as a function of mass at least for clusters older than the Hyades (∼600 Myr). Unfortunately, only a small number of such clusters are readily available. Outstanding among them is M 67; with its sun-like age of 4 Gyr, it is the oldest cluster explored in this regard today. However, only two distinct groups have been measured reliably to date: F- and G-type stars by Barnes et al. (2016) and mid to late M-type stars by Dungee et al. (2022). Here we present new measurements for stars spanning from early-G to early-M types, overlapping with the prior work and bridging the gap between them. Crucially, we find that the rotation periods of single cluster members in M 67 continue to define a single-valued relationship with stellar mass at this age. Consequently, we claim that rotation is usable as an age indicator even though its mass and age dependencies are likely more complex than was thought earlier.

The systematic exploration of rotation in stellar open clusters using timeseries photometry can be traced back to van Leeuwen & Alphenaar (1982; see also van Leeuwen et al. 1987) in the Pleiades (∼125 Myr) and Radick et al. (1987) in the Hyades. However, v sin i measurements were typical at that time and, despite the ambiguity introduced by the unknown inclination i, showed traces of the patterns we observe today (e.g., Stauffer et al. 1984; Soderblom et al. 1993; Queloz et al. 1998). Subsequently, CCD photometry from ground-based facilities and, later, the onset of space-based CCD photometry (mostly thanks to the Kepler mission) has allowed large-scale studies of various readily available open clusters. Photometric rotation periods, that means measuring the periodic brightness variations induced by activity features traversing the stellar disk, remove the ambiguity of unknown inclination.

Numerous ground-based, and more recently, space-based photometric studies have together constructed a set of mileposts with open clusters, including the Pleiades (125 Myr; van Leeuwen et al. 1987, revisited by Rebull et al. 2016), the Hyades (650 Myr; Radick et al. 1987, revisited by Douglas et al. 2019), Praesepe (700 Myr; Agüeros et al. 2011), NGC 6811 (1 Gyr; Meibom et al. 2011a; Curtis et al. 2019), NGC 6819 (2.5 Gyr; Meibom et al. 2015), and Ruprecht 147 (2.7 Gyr; Gruner & Barnes 2020; Curtis et al. 2020). Fritzewski et al. (2020) have also recently shown that the distribution of stellar rotation periods in the 125 Myr-old clusters NGC 2516, Pleiades, M 35 (Meibom et al. 2009), M 50 (Irwin et al. 2009), and Blanco 1 (Cargile et al. 2014) are indistinguishable. For additional work, and especially for rotation in the wider context of stellar activity, readers may refer to Chap. 5 in Basri (2021).

Recent work has advanced to large field star samples (e.g., Reinhold & Hekker 2020; David et al. 2022; Distefano et al. 2023) which provide a much more extensive parameter space than the clusters alone, especially with respect to metallicity but lack the definitive nature of clusters to allow their usage for calibration. Furthermore, especially concerning slower rotating stars, the detection and identification of accurate stellar rotation periods is somewhat problematical. As Tan & Basri (2020, see also Basri & Nguyen 2018) have pointed out, a significant fraction (≈10%) of the widely used McQuillan et al. (2014) sample is likely listed with only half the actual period due to double-dipping. Field stars require more rigorous and extensive vetting than cluster stars, which have uniform metallicity, age, and distances among them. If not, unaccounted non-solar metallicities, binarity, or an evolutionary state even marginally past the main sequence are prone to lead to flawed conclusions. Given those shortcomings, open clusters remain the primary calibrators for the evolution of stellar rotation, and provide the main motivation for this work.

M 67, also known as NGC 2682, has long been of great interest because it is one of the oldest open clusters known, has an age near that of the sun and is located close enough to allow extensive study. In Table 1, we give an overview of some of the fundamental parameters of M 67. Apart from the work on stellar rotation in M 67 by Barnes et al. (2016) and Dungee et al. (2022) there have been the studies of Gonzalez (2016a,b) and Esselstein (2018, hereafter E18) whose findings challenge not only the prior work on M 67 but also those for other clusters and one bedrock principle of gyrochronology: the single-valued nature of the spindown relation. We also investigate and make some sense of their findings below.

This paper is structured in the following way. In Sect. 2, we describe the observations that this work is based on, as well as the stellar content of the Kepler K2 superstamp. This section includes a membership analysis based on GDR3 data which is detailed in Appendix A. Section 3 addresses the creation of light curves from the K2 data. The section gives a short overview of the steps carried out while Appendix B describes all the details regarding the difficulties of K2 data and our ways of overcoming them. The light curves thus constructed are subjected to a period analysis which is described in Sect. 3.2. The results of the analysis are presented in Sect. 4. There, we construct a Color-period diagram (CPD) for M 67 after scrutinizing each and every star in the sample. We then go on to discuss our findings and their implications in the contexts of prior work on open cluster and the modeling of rotational evolution in Sect. 5 before we finish with some conclusions in Sect. 6.

2. Observational data

In this section we describe the photometric data we use throughout this work: the Kepler K2 Campaign 05 superstamp. We obtain the stellar content of the field of view from the Gaia catalogs and carry out a membership analysis. The cluster sample constructed is then used to reconfirm the age of M 67 and estimate reddening parameters in different colors.

2.1. The K2 superstamp

Our goal is to explore the cluster center of M 67 to its fullest. Therefore, we ignore the few Kepler K2 light curves supplied by the mission itself, and construct our own based on a new, more extensive target list. This allows us to use the richness of the cluster center without being bound by the preselection of targets in the K2 program. We show in Appendix B.9 the meager extent of the original K2 program, especially when it comes to the superstamp. However, this also means that we now have to solve by ourselves the problems addressed by the mission itself when it provided the K2 light curves.

After the conclusion of its primary mission, the Kepler telescope started its K2 program. Between 2014 and 2019, the telescope was pointed at 19 different fields along the ecliptic for approximately 80 d each (denoted as campaigns, C## hereafter). Prior to and during this long observation run several parts of the telescope ceased to function. However, Kepler was able to continue its mission until it ultimately ran out of coolant in 2019, long after its originally designated lifetime. Some of the defects make K2 data problematic to use if not properly addressed.

M 67 was in the field of view during three of those campaigns, namely C05, C16, and C18. Special attention was paid to the cluster during C05 and data for an extended region around the cluster center were collected. This superstamp covers the central part of the cluster (cf. Fig. 1) from April to July in 2015 and is the basis for the present work. The coverage during C16 (Dec. 2017 to Feb. 2018) and C18 (May 2018 to Jul. 2018) may provide additional opportunities for future explorations.

Fig. 1. Overview of the K2 C05 superstamp coverage of M 67. The upper panel shows one superstamp image on a logarithmic gray-scale. Stars recognized by GDR3 are overplotted, color-coded according to our membership evaluation (see Sect. 2.2 and Appendix A for details). The lower panel shows a DSS2 (red channel) image of the same region. The extent of the superstamp is indicated in red. The small inset in the corner shows the location of the superstamp in the Kepler FOV C05. Both panels are shown such that their orientations are as similar as possible (see coordinates in panel b).

During C05, Kepler obtained 3663 images in its slow-cadence mode. These span a baseline of approximately 70 d. Because of data transmission limitations, Kepler did not transmit full frame images, but only certain preselected regions. These are stored and made available as target pixel files (TPF). For C05, TPFs were selected such that they cover a continuous region around the center of M 67. Cody et al. (2018) used those TPFs to create a continous region, the superstamp, around the cluster, measuring 400 by 400 pixels in size (≈0.5° ×0.5°). The coverage of the cluster region in the superstamp and its location in the Kepler FOV are illustrated in Fig. 1. Additionally, they created a new, time dependent astrometric calibration for each individual image, overwriting the default, time independent K2 pipeline results to account for the jitter of the targets. We note that this astrometric solution turns out to be both very good and useful, and it is of fundamental importance for this work. We note further that Cody et al. (2018) omitted a small number of cadences, resulting in only 3620 superstamp images for M 67.

2.2. Stellar content of the superstamp

The superstamp covers a region of about 0.25 deg², aligned with the cluster center. The most consistent and extensive catalogs for any random region of the sky are from the Gaia mission and thus we use Gaia DR3 (Gaia Collaboration 2022, GDR3 hereafter) as a foundation for our work. We obtain a subset of GDR3 data for the region covered by the superstamp. This subset contains about 2000 stars, but the K2 mission only provided light curves for 96 of these. Due to the motion of the telescope, individual targets on the edges may shift in and out of the field of view. Here, we did not limit the stars by brightness or any other parameter. However, we note that, as expected, the parameter coverage of the individual stars varies widely in the sample, generally becoming more sparse toward the faint end. We note further that the spatial resolution of Gaia greatly exceeds that of Kepler¹, while on the other hand Gaia omits certain significantly brighter sources². This requires us to carefully assess which star’s light we are analyzing.

M 67 is obvious in the field-of-view (FOV) sample in both photometry and astrometry (see Fig. 2). In a color-dependent sense, the brightest stars in the field belong to M 67, with the exception of a few very bright foreground stars. We complement the GDR3 sample with data from the Two-micron All Sky Survey (2MASS; Cutri et al. 2003; Skrutskie et al. 2006), US Naval Observatory B catalog (USNO; Monet et al. 2003), Guide Star Catalog 2.4.2 (GSC; Lasker et al. 2008), and the Ecliptic Plane Input Catalog (EPIC; Huber & Bryson 2017, via their 2MASS IDs, where available) and the SIMBAD database³. This crossmatch provides us with multiband photometry and radial velocity information for a large number of stars. Most of the radial velocity measurements are from groundbased surveys by Geller et al. (2015, 2021) and Donor et al. (2018). We note that these surveys prioritized cluster stars and, as such, the availability of radial velocities is strongly biased toward member stars and is brightness limited, with the cutoff being approximately G = 17 mag. We supplement those radial velocities with additional ones provided in GDR3. Based on this, admittedly heterogeneously available information, we designate each star in the field either as a member, or as field star. Details about this designation are described in Appendix A. Figure 3 depicts the resulting cluster sample which amounts to 971 members and 1042 field stars. Around 80% of the members can be assumed to be main sequence stars, spanning a brightness range from G = 13 mag down to G = 21 mag.

Fig. 2. Overview of GDR3 stars in the field of view and their M 67 membership status (following Sect. 2.2 and Appendix A). In all panels, blue stars denote members while non-members are gray. Panel a shows a color-magnitude diagram in Gaia G − GRP color, panel b the proper motions, panel c a histogram of radial velocities, and panel d a histogram of the stellar parallax. Dashed red lines indicate the position of M 67.

We also include a designation for stars displaying signs of binarity. Those are derived either from the extensive information present in the SIMBAD archive or from a star’s position in the CMD indicating a photometric binary. The former leads to a star being assigned binary status if it is listed in SIMBAD as an eclipsing binary (EB*), spectroscopic binary (SB*), or cataclysmic variable (RSCVn or CataclyV). Their designations are based on various catalogs and studies, too numerous to be listed here, do not include all actual binaries in the sample, or could misidentify a single star as a binary. However, judging from the overall picture that emerges, we find this designation to be relatively reliable and consistent. An identification of binaries is crucial because multiplicity can lead to a number of issues, for example, tidal interactions or mass exchange may change the stellar angular momentum and as such disqualify a star for our purposes. Those stars need to be discarded from our final sample.

Henceforth, we will use G − G_RP as the color parameter for the stars. This has the advantage that it is essentially independent of Gaia G_BP, which is either very uncertain or unavailable for the faint and red stars. The disadvantage of this choice is the unavailability of a measured reddening parameter in this color. We will address this issue below. Plots in other colors (such as B − V, V − K, or G_BP − G_RP) that are potentially useful to the community are provided either parallel to the main plot or as supplementary plots in Appendix C.3. We will not use the individual parallaxes of the stars to obtain absolute magnitudes. Stellar photometry is much better constrained and more widely available. As such, for a subsequent isochrone fit to the cluster, we will apply the cluster parallax (and reddening, cf. Table 1) to the isochrone to match the cluster data instead.

2.3. Color magnitude diagram for M 67

To reconfirm the age of the cluster and to determine the E(G − G_RP) reddening and extinction A_G, we proceed with an isochrone fit to the cluster sample (cf. Fig. 3). We opt for an empirical determination of E(G − G_RP) for M 67 based on the offset observed between a distance-corrected isochrone and the observed cluster sample. We work with the isochrones from the Padova and Trieste Stellar Evolutionary Code (PARSEC⁴; Bressan et al. 2012; Chen et al. 2014, 2015). In principle, the isochrone fit has four free parameters: the cluster metallicity, cluster age, reddening, and extinction. However, we can constrain some of them. The cluster metallicity has been measured repeatedly in the past, with values ranging from [Fe/H]= − 0.01 to [Fe/H] = 0.03. Thus M 67 is essentially of solar metalicity. We adopt [Fe/H] = 0.03 (Netopil et al. 2016) and note that changes within the range described above do not change the result below in a significant way.

Fig. 3. Isochrone fit to the identified cluster members in two different Gaia colors. The PARSEC isochrone G-band magnitude is adjusted for the parallax ϖ and extinction AG using Eq. (2). The PARSEC GBP − GRP color is reddened by E(GBP − GRP) from Eq. (1). The dashed line indicates the nominal position of the equal-mass binary sequence (main sequence stars only). An alternative version using B − V and V − K is shown in Fig. C.1.

Following the method described in Gruner & Barnes (2020), we use the coefficients obtained by Casagrande & VandenBerg (2018) for the relation between the extinction in Gaia G_BP − G_RP and Johnson B − V colors. We adopt E(B − V) = 0.04 mag⁵ (Taylor 2007) and calculate

(1)

and combine it with the relation that we already used for Ruprecht 147 (Gruner & Barnes 2020) for extinction

(2)

With A_G constrained, we adjust E(G − G_RP) to best reproduce the observed cluster sequence and find

(3)

This fit uses the distance modulus listed in Table 1. As a byproduct of this calculation, we have obtained E(G_BP − G_RP). We note that this value also produces a consistent fit between the isochrone and the cluster in G_BP − G_RP (see panel b of Fig. 3). Below, we also provide V − K as a Gaia-independent color that omits the issues of blue bands. For the corresponding reddening, we use the prescription by Martin & Whittet (1990) to obtain

(4)

M 67 is generally believed to be about 4.0 Gyr old, somewhat younger than the sun. In the process of fitting isochrones to the cluster sequence, we notice that the isochrones suggest a slightly younger age for the cluster, namely t = 3.6 Gyr. This is mostly caused by the shape of the turn-off point and the subgiant sequence (see zoomed plots in Fig. 3). However, a 4 Gyr PARSEC isochrone still provides a reasonably good fit to the cluster stars. After all, isochrone fitting is not the purpose of this paper. We note that in no color do the isochrones reproduce the faint red end of the main sequence completely satisfactorily. We note further that the observed red giant branch is slightly bluer than the isochrone prediction. This seems to point to an underlying problem in the isochrones (or the input physics) rather than the photometry as it is also visible in other, non-Gaia colors (cf. Fig. C.1).

3. Light curves and rotation periods from the K2 superstamp

In this section we describe the creation and correction of stellar light curves from the K2 superstamp. We lay out briefly the problems inherent to K2 data and our approach to dealing with them. A more thorough explanation of all the technical details is provided in Appendix B.

During the K2 mission, some vital parts of the Kepler telescope became dysfunctional. The loss of certain parts of the detector (see missing CCDs in panel b inset in Fig. 1) does not affect us. However, the well known pointing problems of the telescope do. Essentially, the telescope was in a constant state of drift throughout the observations, causing the stars to move slowly across the detector. This drift was periodically corrected (i.e., at ≈6 h intervals) by firing the telescope’s thrusters (Van Cleve et al. 2016). Consequently, stars move across the detector during the run. The movement is small for individual exposures and as a result there are no noticeable star trails on the images. Between the thruster firings, however, a star moves up to two pixels (≲8″) across the detector. There are also significant sensitivity variations between the pixels and within the individual pixel. This means that due to the changed positions of an otherwise constant star on the detector, the recorded flux varies from image to image. Those changes, while systematic, are unique for each pixel and its environment. As such, there is no general way to correct for those systematics in a simple, wholesale manner. However, those variations are rather fast and there have been approaches to correct those systematic effects (e.g., K2SC and EVEREST programs by Aigrain et al. 2016; Luger et al. 2016, respectively) with varying levels of success and usability of the resulting corrected light curves for various purposes.

Unfortunately, none of these prior light curves are very well suited for our purpose. The main reason is their availability for only a limited number of targets. K2SC and EVEREST (and essentially all other works in this regard) operate on the EPIC catalog and the sample of stars observed by Kepler⁶. The list of observed stars originates from the original proposals that shaped the Kepler K2 mission. However, for the superstamp there are only 96 individually designated targets with light curves whereas the field contains more than 2000 stellar sources. Therefore, we create our own light curves, directly based on the superstamp FFI and a list of sources obtained from Gaia DR3 rather than the EPIC catalog. This means that we have to perform the photometry and corrections thereof from scratch. Below, we lay out the principal ideas and steps that we employ to create the light curves we subsequently investigate for rotation signals.

3.1. Light curve extraction and correction

We have constructed a procedure aiming to extract light curves and to remove the artificial variations introduced by the effects mentioned above. The plan is to create an empirical model that captures the systematic position dependence of the flux in order to remove it. The position dependence turns out to be of a very complex nature and requires special attention to the individual stars. As such, it is rather labor intensive and requires iteration between different steps until the best result is achieved. Ultimately, we are able to create light curves that are sufficiently free of systematics for a large number of stars.

In Appendix B, we provide a detailed explanation of the extraction and correction process applied to the data in order to create the light curves. There we include all the technical details and illustrate them on a sample star. Here, we only describe the general steps taken.

First, we introduce a naming convention that we will use throughout this work. The total flux variation of an extracted light curve for a star is the combination of three things. First there is the intrinsic astrophysical variability of the star, which is the variation we are actually seeking. This flux is modified by the motion of the star across the detector in combination with the detailed pixel sensitivity. We will refer to the changes caused by this effect as the instrumental systematics. Superimposed on those are quasi-systematic, long-term patterns in the data which are common between similar stars, and which we will refer to as trending or trends. We illustrate the correction steps for an example star in Fig. 4.

Fig. 4. Example of the light curve reduction process using Gaia DR3 604971466769552128. Panel a shows the raw light curve (blue dots), panel b shows an intermediate light curve after the systematics correction (blue dots) together with the PCA correction (red) that is applied to the create the final reduced light curve that is shown in panel c. The orange line in panel c shows the reduced light curve with 1.5 d binning to highlight the now visible intrinsic stellar variability. The insets in a and b zoom in on a 2 d span of the light curve. The drift-induced flux changes and the jumps caused by the realignment of the telescope are readily visible in panel a, while their disappearance can be seen in panel b.

We define an aperture mask for each star on the superstamp and create the raw light curve by simply summing up the enclosed flux. The aperture mask is defined manually and adjusted individually for each star to provide the best (i.e., cleanest, low noise, best systematics removal) light curve possible. Depending on the star being considered, the instrumental systematics appear to introduce flux variations up to a few tenths of a magnitude (cf. panel a of Fig. 4). We correlate the flux with positional changes of the star on the detector. The time-dependent position is obtained from the world coordinate system (WCS⁷; Greisen & Calabretta 2002) which is part of the superstamp FFI data for each individual image and was created by Cody et al. (2018). We proceed to model this behavior (flux as a function of position) with a fifth-order polynomial in order to remove it. Here, we also need to handle jumps in the data that do not allow us to process the light curve as a whole but only in chunks as described in Appendix B. This is probably the most intricate step and provides us with a set of light curves largely free of the instrumental systematics (cf. Fig. 4b).

The processed light curves still contain trending at this point, which we then correct using Principal Component Analysis (PCA). For this, we select a subsample which excludes all stars with obvious intrinsic variability. For any given star that is to be corrected we select all similar stars (in terms of brightness and position) from this subset as the basis for the PCA. This procedure needs such selectivity as there are significant differences between the stars in the trends with respect to those parameters. And with an 80 d light curve, expected rotation periods of 20−30 d and trends of similar length, there is a limited time baseline available for the light curve to express those variabilities. Degeneracy between them is a real issue and it again requires attention to the individual stars to identify. From the individually tailored PCA basis we calculate a correction based on two components and apply it to the light curves. This results in the final light curve product (cf. panel c of Fig. 4) that we provide and use in the subsequent period analysis.

This process is not successful for all stars; faint stars and those in particularly crowded areas cannot be processed satisfactorily. In fact, most light curves still show remnants of systematics and trends. However, those are now typically small compared with the stellar signal itself and as such more of a nuisance than a real stumbling block. We remove all stars for which we cannot obtain light curves of reasonable quality. This includes stars that are located in problematic areas on the detector, that is, those where large systematic effects impact groups or rows of pixels. Some of these regions are visible in the upper panel of Fig. 1 (horizontal and vertical structures, e.g., around y = 245). Our usage of aperture photometry, does not allow us to separate heavily blended sources. Furthermore, the correlation between location and flux cannot be adequately reproduced and removed for every star. This is especially true for stars which show very rapid intrinsic variability (i.e., variability similar to the systematics). We are also forced to remove the first 2 d of the light curves and a central part (172−178 d) as those simply defy any correction attempt with our approach.

The resulting light curves show a wide variety of signals in addition to the rotation signals we seek. We also identify pulsations and eclipsing binaries (some of them with secondary eclipses). Traces of the instrumental systematics remain for a number of stars; however, those are often minor compared with the observed intrinsic signal. We note that we do not (and cannot) attempt to extract and fine tune each and every star in the field of view. Each target apparently requires individual attention, from the design of the pixelmask to the evaluation of the resulting light curve. Given our science goal, we limit ourselves only to M 67 members, together with a sufficient number of non-members to provide a good basis for the PCA.

We also emphasize that the light curves we produce with this method are not intended to rival large scale correction endeavors such as EVEREST or K2SC. Our light curves are purpose-built, with assumptions made in the creation that are invalid for other purposes (e.g., astroseismology). However, we include a comparison between our results and light curves provided in the larger endeavors in Appendix B.9 to validate our results. We will provide all light curves of stars for which a signal was identified as part of the auxiliary data to this publication. Figure C.4 displays the light curves for those 47 stars that are single M 67 main-sequence stars and where we were able to identify a periodic rotational signal.

3.2. Period analysis

For the period analysis, we continue with a hands-on approach for the individual light curves. Unlike the case of most ground-based data, space-based data are both well-sampled and, for the most part, equally sampled. Owing to that, a long-period signal can usually be identified by eye easily when present (cf. panel a in Fig. 5). Manual inspection also permits the identification of a periodic signal when spot evolution or remnants of systematics or trends are present. Automated algorithms often struggle or fail outright in such cases. However, we subject all light curves to an array of algorithms to verify and quantify a generally obvious signal (panel b in the figure). We run period finding algorithms employing a Lomb-Scargle (LSC; Scargle 1982) periodogram⁸, a Clean analysis (CA; Roberts et al. 1987), and phase dispersion minimization (PDM; Stellingwerf 1978).

Fig. 5. Result of period analysis on the light curve for Gaia DR3 604971466769552128. Panel a shows the light curve produced, one which exhibits a clear periodic signal. The power spectra in panel b are obtained using multiple methods as indicated. The Clean (CA) and Combined power spectra are normalized to their maxima and the PDM spectrum is shown as 1 − θ. Panel c shows the phase-folded light curve, folded with P = 27.9 d, corresponding to the peak in the combined spectrum.

Whenever the results of the individual algorithms are inconsistent, a manually determined period supersedes the algorithms. Typically, the algorithms agree among each other, with the one most prone to failure being LSC, and the most consistently reliable being PDM. This is to be expected as the spot-induced features in the light curve appear to deviate strongly from a sinusoidal shape for the stars we are most interested in. We provide the derived power spectra together with the light curve data in the auxiliary files. Figure C.4 includes plots of the phase folded light curves.

Period errors are derived in much the same way as the periods themselves. The presence of data systematics and spot evolution makes the usage of an automated algorithm for error determination problematical. Therefore, we decided to determine errors manually based on the phase folded light curve. We folded the light curve with different periods and set the error range such that for every period covered a phase folded light curve provides a reasonable result in matching the spot-induced features. The uncertainty arises from a number of factors: (1) the number of periods covered, that is, how many repetitions of the features we observe in the K2 C5 baseline, (2) the general signal-to-noise of the variation, and (3) the amount of spot evolution, together with the remaining systematics. Generally, this results in errors that are in the order of ±5−10%. This procedure is supported during the comparison between our derived periods and those from the literature (see below in Sect. 5). Furthermore, ignoring the three outliers, 37 of the 44 remaining stars of our sample overlap with the trend line (see Fig. 9 below), suggesting that our uncertainties are on the order of 1.5σ.

4. Results

In this section we present the stars for which periodic signals were identified. We investigate outliers, construct a Color-Period diagram (CPD) for M 67 of single MS-stars to be used in the subsequent discussion, and present the final CPD for M 67.

4.1. The raw Color-Period Diagram

We identify 136 stars in the FOV that exhibit periodic signals of which 96 are M 67 members which exhibit signals that can be attributed to stellar rotation and 83 of which are MS stars. Apart from the few stars that do not have a convective envelope (namely the blue stragglers with spectral types earlier than mid-F), the sample covers spectral types from early-G to mid-M. Figure 6 gives an overview of the stellar rotation period distribution found. The detected periods range from a few hours up to 38 d, the latter abutting the detection limit of K2 C05. A significant fraction of the cluster sample (53 stars, 43 on the MS) are binaries; those are typically both fast rotating and among the bluer stars of the sample.

Fig. 6. Stars for which we have identified rotation periods. Panel a shows a color-period diagram (CPD) and panel b the corresponding color-magnitude diagram (CMD). Stars belonging to M 67 according to our membership analysis are shown as blue circles (the main sequence stars) and red squares (post-MS stars). Open symbols (in all of the above) mark those that show signs of binarity. Error bars are suppressed here for visibility reasons.

To enable the usage of this sample as a calibration set for gyrochronology, we also remove all stars that detract from this purpose. We define these as any that violate one of the following statements: (1) The stars are M 67 cluster members. (2) Their colors are a valid proxy for their mass, which means they are unevolved. (3) They are not likely to have experienced angular momentum transfer that has caused externally induced changes in their rotation rate, and have spun down over time based on the processes described above.

Adherence to the first statement is simply addressed through our membership analysis. Figure C.2 shows that the field stars (red symbols there) are scattered all across the CPD, corresponding to a wide range of (unknown) ages. Accordingly, we remove them from our sample.

The second statement is violated by evolved stars. Stellar evolution causes stars to change their colors as well. Changes on the main sequence are small enough to be a nonissue here. However, when stars leave the main sequence they become significant. This means that the color ceases to be a valid proxy for mass for all stars that have left the main sequence and have undergone drastic changes in their colors (typically turning redder toward the red giant branch). We will address those stars separately below in Sect. 5.5 with respect to what happens when we forgo color and revert to assess the stars directly via their masses. For now, we omit them from our sample by removing all stars brighter than G = 13.2 mag (red symbols in Fig. 6).

The third statement is violated by stars that have experienced angular momentum transfer in the past. Prime candidates for this are stars in close binary systems. Angular momentum is typically transferred from orbit to rotation, causing a spin up (rejuvenation of the rotation rate). Consequently, removal of binaries from studies like ours is performed almost habitually. Here we see the effects of binarity clearly thanks to our sample size. A significant fraction of the cluster sample (53 stars, 43 on the MS) are binaries; and those are both fast rotating and typically among the bluer stars of the sample. Examination of Fig. 6 shows that whereas the single cluster members are concentrated at the long-period end of the distribution, already suggesting a sequence against color, the binaries are scattered over the entire rotation period range, with almost all having shorter periods than the single stars.

To our knowledge, this phenomenon has not been observed as starkly before. After an extensive radial velocity and photometric survey of the young open clusters M 35 (∼150 Myr) and M 34 (∼200 Myr; Meibom et al. 2009), Meibom et al. (2011b) noticed that the rotational distribution of binaries beyond the influence of tides was marginally skewed toward faster rotation. Other work in young open clusters (e.g., in NGC 2516 and NGC 3532 by Fritzewski et al. 2020, respectively) found that while certain (presumably very close) binaries were anomalous, the vast majority of the photometric binaries were rotationally indistinguishable from the single stars. In the far older (2.5 Gyr) but much sparser cluster Ruprecht 147 (Gruner & Barnes 2020) noted that three of the four binaries therein were off the single star sequence. However, that is too few stars to be able to draw conclusions. In contrast, the situation here is obvious 34 out of the 43 main sequence stars with signs of binarity exhibit slightly to considerably faster rotation than the single star sequence.

In fact, the scatter is so great that it seems unlikely that photometric or tighter binaries of this age are suitable for gyrochronology. Their ages will have to be determined by other means. Although we only have indications for photometric binarity for many of our stars (increasingly toward the fainter end), this already appears to be enough to identify a star as unsuitable. We stress that binarity does not mean that a star does not agree with the rotation of a similar single star per se; however, it is obvious in our results that it is considerably more likely than not that rotation is affected. Consequently, we remove all stars with signs of binarity from our sample, including the nine binaries whose positions agree with the single star sequence.

The above pruning leaves us with the 47 stars that are displayed in Fig. 7, enabling a closer look at the emerging period distribution and its features. We find a group of slower rotating stars with rotation periods between 15 and 35 d, spanning from early-G to early-M. They follow a somewhat sinusoidal shape, with a maximum at (G − G_RP)₀ ∼ 0.65 and a local minimum at (G − G_RP)₀ ∼ 0.85. We identify this with the classical slow rotator sequence. We have indicated this group with a simple trend line (cf. solid line in the figure). It is created from a simple cubic interpolation to points listed in Table C.1. We stress that this line is drawn solely to guide the eye and indicate the cluster sequence. We will employ similar indications in color-period diagrams for other clusters below. Another group of stars apparently follows the same distribution but at half the period. This half-period sequence is formed by double dipping stars which feature two spots whose signals in the light curve are indistinguishable and as such appear to have only half the actual period.

Fig. 7. Color-period diagram for our sample, now including rotation period uncertainties. Only member stars on the main sequence, and with no indications of binarity are displayed. An approximation of the emerging cluster distribution (solid gray line) is overplotted, together with its half-period counterpart (dashed line).

4.2. Double dipping stars

Stars can exhibit more than one significant star spot (or group). In fact, as Basri & Nguyen (2018; see also Tan & Basri 2020) have shown, it is more likely for stars with longer rotation periods to exhibit more than one spot at a given time. We recall that we use the term spot as a simple handle for any coherent activity (combination of) features on the stellar surface that is visible as a significant modulation in our light curves. This explicitly includes faculae, which are assumed to be the dominant flux-altering feature at the age of M 67 (e.g. Reinhold et al. 2019). Following the results of Basri & Nguyen (2018), at the rotation rates we find for M 67 it is up to ten times more likely. We indeed find that about half of the stars in our sample exhibit signs of more than one surface feature. This may create an additional problem. If those features are sufficiently similar in shape and close to a phase shift of 0.5, they may become indistinguishable. This means that we are likely to find a number of stars whose determined rotation periods will be half of their actual periods. Sufficient numbers of such stars will form a sequence at half the actual period, a half-period sequence.

Indeed, we do observe this behavior (cf. Fig. 8) 43 stars define a long-periodic sequence with periods between 15 and 38 d, while four stars follow the half-period sequence. Following Basri & Nguyen (2018) and the distribution of other cluster stars, it is a reasonable assumption that those four stars have rotation periods that are actually twice the identified ones. Therefore, we adopt final rotation periods for those four stars that are twice the measured ones. We note that such light curves are also visible in the prior rotation period work on warmer stars in the M 67 work of Barnes et al. (2016).

Fig. 8. CPD for our M 67 sample emphasizing stars on the half-period sequence. The cluster sequence, approximated by an interpolated line, is overplotted in gray. The dashed line shows its half-period counterpart. Stars for which we double the measured periods (red) are connected to their new positions by dashed lines.

4.3. Color-period diagram for M 67

With the half-period sequence stars accounted for and with binaries and evolved stars eliminated, we can construct a color-period diagram for M 67 that can be compared with other clusters and with spindown models. The final emerging CPD contains 47 stars and is shown in Fig. 9. Table C.2 lists these stars.

Fig. 9. Final CPD of M 67 based on our study. The corresponding data are displayed in Table C.2.

The emerging distribution of rotation periods, spanning from early-G to mid-M stars, shows a clear sinusoidal sequence. It rises from around 18 d periods for G0 stars, the earliest type stars at the age of M 67 that is still on the MS, to periods around 30 d for K3 stars. Following the distribution redward, the periods decrease again to around 25 d at M0 only to rise again to 30 d and longer for mid-M. This distribution will be compared with prior work and discussed in more detail in the next section, following a discussion of the outliers immediately below.

4.4. Outliers

Our sample contains three fast rotating stars with periods P ≤ 10 d that deviate strongly from the slow rotator sequence occupied by all other stars. We suspect that the likely origin of this deviation is undiscovered binarity impacting the stellar angular momentum evolution. The stars in question are actually inconspicuous among the sample as regards their renormalized unit weight error (RUWE) in GDR3. (A larger than average RUWE value can indicate an underlying issue with the astrometric solution, which itself may originate in undiscovered multiplicity).

Gaia DR3 604911204083987584 (P = 7.4 d, (G − G_RP)₀ = 0.55 mag) sits right within the area of the CPD that is typically occupied by binaries (cf. Fig. 6). It does have two very close companions (0.5 and 4.3 mag fainter) within 5″. Those cannot be separated during light curve extraction; however, it is clear from an investigation with different pixelmasks that the star itself is the source of the observed variability. Geller et al. (2015) actually designate the star as a non-member based on radial velocity work but Gaia astrometry indicates a high probability member. Consequently, we retain it in the CPD.

Gaia DR3 604922229264424448 (P = 9.7 d, (G − G_RP)₀ = 0.84 mag) could potentially be part of the half-period sequence. However, closer inspection shows that it lies significantly below, a deviation even more pronounced if one were to double its period. Furthermore, its light curve shows clear signs of continuous spot evolution which would mean both spots evolving identically, if it were a double-dipping star. This seems unlikely.

For those two objects, angular momentum exchange (e.g., in a binary system) appears to be the most likely explanation. The third outlier is the reddest and simultaneously fastest rotating star of our sample Gaia DR3 604969061592133376 (P = 1.8 d, (G − G_RP)₀ = 1.11 mag). It could still be part of the fast rotator sequence, similar to what can be observed in other, younger clusters. It would be the first known star that is part of the fast rotator (denoted convective in Barnes 2003) sequence at a confirmed age older than 1 Gyr. However, it may also simply be a star with hidden binarity.

5. Discussion

In this section our results for M 67 are placed in the larger context of the rotational evolution of stars. First, we compare our results with findings in past studies and investigate deviations and inconsistencies. In the process, we address certain issues with prior work. We then build a sample of combined knowledge for rotation data from M 67 and compare it with the few old open clusters studied to date and also with the predictions of rotational spindown models.

5.1. Comparison with prior work on M 67

There have been three prior studies of M 67 using K2 light curves from C05: Barnes et al. (2016), Gonzalez (2016a, with a followup in Gonzalez 2016b, G16 hereafter), and E18 (see also Esselstein et al. 2018). Those studies were performed on the presearch data conditioning (PDC) light curves provided in the K2 archive⁹. These light curves were created from aperture photometry and include a correction that, to a certain degree, accounts for systematics and trending. However, as pointed out by Barnes et al. (2016), the light curves are not free of either problem. E18 made additional use of self-extracted light curves from the C05 superstamp and C16 light curves. All studies operated on the low spatial resolution given by the EPIC catalog. For this comparison we match their samples with GDR3. In addition to the work mentioned above, a recent study by Dungee et al. (2022) created light curves from 3 yr worth of ground-based observation. This latter work is limited to M-dwarfs.

Barnes et al. (2016, 20 stars) and Gonzalez (2016a,b, 98 stars) operated on the light curves provided by the Kepler mission itself; only a few of those were from the superstamp region (cf. Fig. B.14). Thus, the overlap with our sample is small. E18 (30 stars) includes parts of the superstamp, but the overlap is still small. We also have some overlap with Dungee et al. (2022). Table 2 shows the stars that are common between the samples. The agreement in derived rotation periods is very good, with period differences typically ≲2 d. The only star with a significant difference is EPIC 211397501¹⁰, for which Gonzalez (2016b) reported the half-period. (It is identified as a half-period sequence star by us). In the comparison with Dungee et al. (2022) we find that seven rotation periods (including three that were nominally rejected from their sample) are in good agreement with our periods.

5.1.1. Barnes et al. (2016)

Given that instrumental systematics and trending are still present in the data, Barnes et al. (2016) subjected their sample to a PCA to remove trending and, partially, systematics. They identified 20 stars with a rotational signal, shown with red squares in panel a of Fig. 10, together with our sample. Generally speaking, they subjected each light curve to the same vetting that we performed, and for the common stars, found periods in agreement with ours to within the uncertainties.

Fig. 10. Comparison of our results with those of Barnes et al. (2016). Panel a shows a CPD as in Fig. 9 and overplotted in red are the measurements by Barnes et al. (2016). The gray lines, solid and dashed, are our indication of the cluster and half-period sequences, respectively. Panel b shows a CMD of the same stars and the inset (panel c) the proper motions centered on M 67. Gray dots in panels b and c are members without period determinations.

The Barnes et al. (2016) sample forms a sequence spanning from early G stars with rotation rates around 18 d to mid-K stars with rotation rates bset with our findings. This new subset reduces the Gonzalez (2016b) sample from an area to a sequence in the CPD. A particular improvement is seen for the bluest stars of the G16 sample, where we now exclude almost all of the long period G-type stars. This improvement is also visible in the comparison with our data shown in Fig. 11. The agreement is now much better than with the entire sample. The largest part of the subset occupies the downward sloping section of our sequence, starting at the period maximum at early-K stars and going toward the minimum at early-M, thereby matching our distribution well. However, the downward sloping section appears slightly red-shifted (Δ(G − G_RP)₀ ≈ 0.1 mag). It also connects well to the Barnes et al. (2016) sample and extends it toward redder stars with a small range of overlap. Interestingly, all five stars of the Gonzalez (2016b) sample that lie on the half-period sequence made it into our vetted subset.

Fig. 11. Same as Fig. 10 but with a comparison to results of Gonzalez (2016b; red), showing both retained and discarded stars from that sample.

Finally, we note that G16 adopted 998 light curves and identified rotational signals in 441 (≈44%) for their M 67 study. This is a very high fraction, even for a much younger, more active cluster, where stellar rotation periods are far easier to identify. For comparison, we identified rotation signals in < 10% of the sample and Gruner & Barnes (2020) identified rotation periods in 21 out of 102 (≈20%) stars for the 2.7 Gyr-old open cluster Ruprecht 147. We therefore conclude that the results of Gonzalez (2016a) represent an overly inclusive set of stars (unfortunately including misidentifications) in addition to those with true rotational periodicity. Our sample is far more exclusive in comparison. When systematics and instrumental trends are excluded, the two distributions can be seen to be substantially similar.

5.1.2. Esselstein (2018)

The study of E18 (see also Esselstein et al. 2018) had a slightly different focus than the ones before. Their main objective was to assess the detectability of rotation signals from light curves of stars as old and inactive as in M 67. They devised a sophisticated methodology around the artificial injection of variability signals into real Kepler data. Mostly as a byproduct, they identified rotation periods from archived light curves (SAP) and self-extracted data.

E18 provides rotation periods for 30 stars (their Table A.13) that span from the turn-off point to early K-stars (cf. panel a in Fig. 12). Their sample includes a large number of fast rotators while the slow rotators show significant scatter. As before, we scrutinize the sample by identifying binaries and also subject their light curves to manual inspection. E18 lists identifications for binaries and probable binaries which we complement with a photometric binary designation from Gaia photometry. We come up with a list of 11 binaries. Additionally, three of their stars are clearly evolved past the main-sequence. All fast rotators in the E18 data either show signs of binarity or are post-MS.

Fig. 12. Same as Fig. 10 but with a comparison to results of E18 (red), again identifying retained, discarded, and binary stars.

The E18 slow rotators agree reasonably well with our distribution (as verified with Barnes et al. 2016 and G16). However, theirs shows significantly greater scatter than the other studies. Similar to our investigation of the G16 sample, we inspect the light curves to look for inconsistencies. This, however, is more difficult than before as their light curves are largely only available to us in plotted form (see their Appendix B). We identify five stars for which we cannot verify the rotation signal. However, their elimination from the sample does not change the overall picture. As before, we adopt their period as is even though we may disagree about the derived value for individual stars. One star of the sample (EPIC 211404554) sits arguably on the half-period sequence.

This scrutiny suggests that E18 does agree reasonably well with the overall picture of M 67 constructed here. Inconsistencies can largely be explained through binarity and post-MS evolution. However, the slow rotator sequence has greater scatter than other results. We argue that this originates in the detection method used. E18 used a Lomb-Scargle analysis. While that is a good choice for fast rotating stars, it ceases to be for slow rotators. This is because the sine wave that is assumed in the LSC becomes less and less suitable for slow rotators due to the occurrence of multiple spots and significant spot evolution between consecutive rotations. Phase dispersion and autocorrelation are much better suited for this period regime. As such, we treat the results of E18 with caution.

5.1.3. Dungee et al. (2022)

A recent study by Dungee et al. (2022) has used ground-based data acquired over a 3 yr baseline to investigate the rotation of M dwarfs in M 67. They have obtained periods for 383 stars and applied iterative outlier rejection leading them to adopt 64 of the M 67 M-dwarfs as their cluster sample. The full sample spans a huge range in rotation periods, from a few to several hundred days. After outlier rejection the remaining sample (dubbed as converged), a clear sequence appears to emerge, starting around 28 d periods for late-K stars and showing an increase in period at mid-M to almost 60 d. Figure 13 shows their periods together with ours.

Fig. 13. Same as Fig. 10 but comparing with the results of Dungee et al. (2022; red). Dungee et al. (2022) assume a flat period error of 10%.

We find very good agreement between our results and those of Dungee et al. (2022). Their sample has its blue boundary at late K-stars and extends toward mid-M. This connects directly to a well populated area in our sample and extends it redward where our data is thinning out but is still in direct agreement with theirs. They have included a photometric cut to their data, eliminating photometric binaries and consequently, we do not find any problematic binary systems in the sample.

We do not have access to the light curve data from Dungee et al. (2022; only plotted images) to inspect them in similar fashion as we did with the other samples. However, we also see no reason to do so as the agreement between our results and theirs is excellent and their long-baseline ground-based data is not subject to the problems specific to Kepler data.

5.1.4. A lesson from prior studies

Based on the above comparisons, we emphasize the need for a very careful approach in the identification of rotation periods from space-based data, especially in the present period regime. Trending and systematics could mask or mimic spot-like variability in a light curve and spot evolution may further complicate period identification. Despite concerns about objectivity, it appears to be essential to supplement all algorithmic period determinations with manual inspection and judgment. We suggest far more exclusivity with similar rotation period derivations generally, but especially in the presence of systematics and when operating with a rather limited observation baseline.

5.2. Assembly of final sample

Before comparisons with models and to work in other clusters can be performed, we need to assemble a final sample. This consists of the following steps:

1. We accept all stars of our sample that are both on the MS and that are not identified as binaries.

2. We assume that the stars occupying the half-period sequence are in fact double-spotted stars and as such we double those observed periods. (This assumption is valid to the extent that at the rotation rates we observe, stars tend to be in the double spotted domain. Furthermore, we have several stars among our sample that are clearly double spotted, but with the spot shape and phase shift such that they cannot be mistaken as being single spotted).

3. We add the stars from the Barnes et al. (2016) sample, but omit the binaries therein.

4. We omit the Gonzalez (2016a) and E18 samples entirely.

5. We include the 64 stars marked as converged from the Dungee et al. (2022) sample.

6. We do not remove stars that are common between the samples.

The addition of stars from Barnes et al. (2016) significantly strengthens the emergent distribution and bolsters the region around mid-G type stars, while adding the D22 stars completely defines the M-dwarf region.

We match the sample to the 2MASS, USNO, and GSC catalogs to obtain a large set of measured stellar colors for each star. The availability of particularly blue bands is understandably sparse toward fainter stars. We also extend this search to non-M 67 stars, namely those included in rotational studies for other clusters (cf. overview in Sect. 1). Those will allow a more reliable comparison between the different cluster sequences that emerge later. The final assembled sample, offering as complete a view of the rotational distribution of M 67 is shown in Fig. 14.

Fig. 14. Color-period diagram for our final assembled sample (combined with the results of Barnes et al. 2016), now including only main sequence stars with no indications of binarity. Only the slow rotator sequence is shown.

5.3. Empirical comparison with observations of other clusters

We now compare our rotational distribution for M 67 with those of certain other open clusters. For this comparison we use data for NGC 6819 (Meibom et al. 2015; 2.5 Gyr), Ruprecht 147 (Gruner & Barnes 2020; Curtis et al. 2020; 2.7 Gyr), and NGC 6811 (Meibom et al. 2011a; Curtis et al. 2019; 1.0 Gyr) as they are the only open clusters equal or older than 1 Gyr for which rotation periods have been determined. We crossmatch their individual catalogs analogously to what was done in Sect. 5.2 to obtain a broader and more equally described sample. This requires accounting for the individual reddenings. We obtain E(B − V) for each cluster and calculate E(G_BP − G_RP) from that. We use the following prescription to obtain E(G − G_RP): As the transformation from E(B − V) to E(G_BP − G_RP) can be obtained with a simple multiplicative factor (as an approximation) it is reasonable to assume the same works for the transition from E(G_BP − G_RP) to E(G − G_RP). To determine the relevant factor we use our empirical result for E(G − G_RP) in Eq. (3) (namely E(G − G_RP) = 0.03 ± 0.005) for M 67 and the calculated E(G_BP − G_RP). The results are listed in Table 3.

(5)

We use this to calculate E(G − G_RP) for all clusters.

M 67 is the oldest open cluster with reliable rotation periods available. It is therefore unsurprising, indeed expected, that its sequence lies above all other clusters in a color dependent sense (cf. Fig. 15). Similar to Ruprecht 147 and NGC 6811, M 67 shows a flattening of its sequence around early-K spectral type, although this feature appears to move to the red with increasing age. The cooler stars (later than mid-K) in Ruprecht 147 indicated a slight downturn of the sequence, with a minimum around late-K/early-M and a steep increase in period with redder color. However, given the small sample size in that range, the existence of the downturn was uncertain. Our new data for M 67 now confirms the existence of this downturn and validates the distribution seen for Ruprecht 147.

Fig. 15. CPD of the three open clusters 1 Gyr or older. Shown are NGC 6811 (Meibom et al. 2011a; Curtis et al. 2019) at 1 Gyr, NGC 6819 (Meibom et al. 2015) and Ruprecht 147 (Gruner & Barnes 2020; Curtis et al. 2020) at 2.5 Gyr, and M 67 (Barnes et al. 2016; Dungee et al. 2022; and this work) at 4 Gyr. The color coding groups cluster data by age. Overplotted is a simplified representation of the emerging cluster sequences by age (solid lines color coded by age). The sun is shown with its usual symbol. Figure C.3 shows equivalent plots in the colors (GBP − GRP)0, (B − V)0, and (V − K)0.

It is visually apparent that M 67 displays the greatest scatter in rotation period. This fact may be surprising in that the stars are expected to have all converged to the slow rotator sequence. The scatter is, however, a consequence of the difficulty of the period measurement and not an intrinsic property of the distribution. Additionally, the relatively short baseline of the space-based observations, combined with differential rotation results in significant uncertainties on the periods. To help visualize the cluster distributions and unclutter subsequent plots, we have created rough approximations of the cluster sequences (based on cubic interpolations to the points listed in Table C.1) as shown with colored lines in Fig. 15.

5.4. Comparison with rotation models

The long-term goal of studies like ours is to understand the rotational evolution of stars and the physics of magnetic braking, to describe the evolution in an empirically constrained model, and to use it to derive stellar ages via gyrochronology. Numerous stellar spindown models have been created over the years with varying levels of detail in their physical descriptions. As our new data extends the knowledge of stellar rotation rates into a parameter space that has not been explored before (lower mass and higher ages), we go on to examine how well the models perform against our new findings.

We know from prior work in Ruprecht 147 (e.g., Curtis et al. 2020; Gruner & Barnes 2020) that models of rotational evolution do not reproduce well the observed distributions of cluster stars at higher ages. Any lingering doubts about discrepancies between the models and the data are removed now. The deviating shape shown by Ruprecht 147 is also present in M 67, and is more pronounced. We know now that rotational evolution deviates significantly from model predictions. As such, we limit this section to comparison up to this conclusion and omit a longer discussion on the intricate model details. For a more detailed discussion see the comparison to Ruprecht 147 in Sect. 6 of Gruner & Barnes (2020).

We also limit ourselves to a comparison with the spindown descriptions of Barnes (2010) and Spada & Lanzafame (2020). These are chosen among the numerous models available as they are related in their descriptions and exemplify the difference based on the inclusion of an additional parameter. This comparison is essentially similar to the one performed in Gruner & Barnes (2020), but the observed deviations there are even more pronounced here, and the larger sample allows a stronger conclusion. We note that the models are only available in certain colors and for particular ages. Spada & Lanzafame (2020) published a few time steps (notably 2.5, 4.0, and 4.57 Gyr in the relevant age range) and are available only in (B − V)₀ color. Barnes (2010) offers a mathematical prescription that we use similarly to Gruner & Barnes (2020) (see their Sect. 6.1) which intrinsically allows any given age and a variety of Johnson and 2MASS colors. B − V is especially problematic with respect to reliability and availability for stars as red and faint as we deal with here. As such, an examination based on a more reliable and consistent set of colors is preferred. We opt for Gaia G − G_RP as it omits the problems of Gaia G_BP − G_RP in the blue and is accessible from one consistent source (unlike e.g., V − K). However, we do include those colors below to offer a complete picture.

To perform the comparison in G − G_RP, we do need the models in said color. We employ a color transformation that is constructed from stellar models in the PARSEC isochrones by comparing the different synthetic colors for stars of the same mass, age, and metallicity. This approach has its weaknesses, for example it is problematic toward the very red, but it is sufficient for our comparison. The comparison that is shown in Fig. 16 is now based on the following: (1) all stellar colors are measured, that is, obtained from one of the above mentioned catalogs and dereddened according to the description in Sect. 2.3, (2) both models in B − V and Barnes (2010) in V − K are obtained directly from the sources, and (3) V − K for Spada & Lanzafame (2020) and both models in Gaia colors are obtained via the transformations. We use the same color transformation on our empirical cluster sequence which was originally constructed in G − G_RP (see Table C.1).

Fig. 16. Comparison between our final sample of measured rotation periods (this work combined with the results of Barnes et al. 2016; Dungee et al. 2022) and the stellar rotation models of Barnes (2010) and Spada & Lanzafame (2020). Each panel performs the comparison in a different color system: (a) G − GRP, (b) GBP − GRP, (c) B − V, and (d) V − K. The errors on B − V are suppressed for visibility reasons.

5.4.1. The Barnes (2010) model

We find good agreement between our results and the Barnes (2010) model prediction for stars earlier than K5 (in the overlap region with the Barnes et al. 2016 data). Redward of K5, the model predicts a steady increase in rotation period with color which is not seen in the data. Instead, the observed stars turn significantly downward toward shorter rotation periods. After this point, the model diverges away from the observed stellar distribution. Similar behavior was observed in Ruprecht 147 by Gruner & Barnes (2020) but the effect is stronger here, and thanks to the larger sample, can now be stated with much more certainty. It is to be noted that the period increase for the reddest stars in the sample is much steeper than predicted in the model. This leads to an apparent reconvergence between model and observations for mid-M stars. As has now been seen with multiple clusters, the Barnes (2010) model fails to predict the stellar rotation rates for older stars later than mid-K. Additional physics is apparently required to describe the spindown adequately.

5.4.2. The Spada & Lanzafame (2020) model

The model of Spada & Lanzafame (2020) incorporates the Barnes (2010) description for magnetic braking and adds a two-zone model where an additional parameter describes the coupling between those two zones and the angular momentum exchange between them. Generally speaking, angular momentum from the radiative core is transferred to the convective envelope, reducing the stellar spindown to the point where a star may even spin up slightly. This effect is mass-dependent. It should be noted that the Spada & Lanzafame (2020) model only describes stars of the slow rotator sequence, in contrast to the Barnes (2010) model. However, this is not directly relevant here, since at the age of M 67, all stars have converged to the slow rotator sequence.

Similar to the Barnes (2010) model, the Spada & Lanzafame (2020) model provides good agreement for stars bluer than mid-K. For later spectral types, it is closer to the data than the Barnes (2010) model, but still does not appear to be able to reproduce the observations in detail. Both the model and observed cluster sequence show first a flattening of the distribution followed by a subsequent downturn in the CPD and eventually a steep increase from blue to red. However, those features appear not to be aligned. The Spada & Lanzafame (2020) model predicts a flattening of the distribution and a potential downturn around M0, whereas the observations show said downturn already at mid-K. Furthermore, the amplitude of this downturn in the model is much shallower than observed. Given that the model incorporates a mass-dependent coupling that is responsible for the observed downturn, it is possible that a recalibrated model that incorporates our findings and the recent results for Ruprecht 147 could improve the level of reproduction of the data and increase the predictive power of the model.

5.4.3. Implications for modeling of rotating stars

The new M 67 rotation periods thus continue to challenge theoretical rotational evolution models. No rotational evolution model to date is capable of reproducing the observed behavior for old cool stars later than mid-K. This is unsuprising, as the calibration of the models for a long time could only rely on warmer stars. And it is only recently that we have started to see that the behavior for cool stars differs strongly from that anticipated. This discovery began with the work of Curtis et al. (2019) on NGC 6811, whose behavior in the red sparked the idea of a (temporarily) stalled spindown, found further indication in the works of Gruner & Barnes (2020) and Curtis et al. (2020) for Ruprecht 147, and is now confirmed by the new data on M 67. The best results are achieved by the Spada & Lanzafame (2020) parameterized two-zone model and small adjustments to the model could potentially provide a better match to the observations.

5.5. Stars at and around the turn-off point

In contrast to younger clusters where all stars redward of the Kraft break are still on the main sequence, M 67 is old enough that this becomes an issue. For those evolved stars, their colors (and temperatures) are no longer a valid proxy for the stellar mass. Consequently, we need to return to the mass as the relevant indicator, and examine the stars omitted above. We estimate the stellar masses of our sample stars from an isochrone fit with a PARSEC isochrone of 3.6 Gyr age (cf. panel a of Fig. 17). The uncertainty for this mass estimate is typically on the order of ΔM ≲ ±0.05 M_⊙.

Fig. 17. CMD and three CPDs highlighting M 67 sample stars evolved beyond the Main sequence. Panel a shows a CMD of our sample with the stars color coded according to their masses. Differing symbols distinguish different types of binaries. Panel b shows a corresponding CPD in Gaia colors and identical symbols. Panel c plots the rotation periods of the stars against their estimated masses. An approximate indication of the position where stars cease to have an outer convective zones (i.e., the Kraft break) is overplotted with a dashed gray line. The solid black line corresponds to the scaling on the right and denotes the ratio between current (R(t)) and ZAMS stellar radius (R0). Panel d shows an inferred CPD with stars assigned colors they had when on the ZAMS.

Plotting the rotation periods against mass (panel c) instead of color (panel b) unclutters the region around the turn-off point. It also demonstrates that those stars appear to populate a direct extension of the cluster sequence toward the blue. Panel d of the same figure projects the stars back to their zero-age main sequence (ZAMS) color, estimated from a 100 Myr isochrone and the derived stellar masses. For this comparison, we have included binaries whose type does not suggest a tight system, that means, excluding eclipsing and cataclysmic systems, and those whose CMD position makes a mass estimate difficult, that is, photometric binaries. As expected, eclipsing and cataclysmic systems show much faster rotation rates than comparable stars of the same mass/ZAMS color.

The impact of post-main sequence expansion of the star’s envelope is likely to be a small effect on their current rotation periods as their radii have not increased significantly (≲2%, cf. panel c of the figure). This is estimated from a comparison of the radii of ZAMS stars and those at the age of M 67.

6. Conclusion

We have studied space-based photometric data from Kepler’s K2 mission for the 4 Gyr open cluster M 67 to examine this key object in its role as a gyrochronology calibrator and to extend our knowledge of the evolution of stellar rotation. We used the K2 superstamp created from Campaign 5 target pixel files together with a stellar sample based on Gaia DR3 and performed membership analysis. We constructed light curves of stars in the superstamp region based on aperture photometry and devised a correction algorithm to deal with the well known K2 systematics. We identified periodic signals in 128 of these light curves, and created a color-period diagram using the 47 stars which are believed to be effectively single main-sequence stars. Those span spectral types from early-G type stars to mid-M and as such reach a region hitherto unexplored.

Our data connect well to prior studies on M 67, especially Barnes et al. (2016). Most of the data of Dungee et al. (2022) lies redward of ours but their results agree with ours in the overlap region and extend it smoothly toward the very red. The work of Gonzalez (2016a) and E18 is also shown to be consistent with our work if they are suitably pruned of binaries and light curves with remaining instrumental effects. Despite the current extension toward significantly later-type stars, open cluster rotation period data remains compatible with the idea that effectively single stars populate a unique surface in rotation period-mass-age space. However, the suggestion of a diminishment of spindown for K-type stars by Curtis et al. (2019) based on NGC 6811 data, and subsequently supported by Gruner & Barnes (2020) and Curtis et al. (2020) for Ruprecht 147 now appears to be secure. The shape of the rotation period surface is more complex than originally envisaged, and deviates strongly from the classical Skumanich-style, mass-dependent predictions. A comparison with models of rotational evolution shows that the models appear to be inadequate for stars redder than mid-K. We find that the model of Spada & Lanzafame (2020) currently provides the closest description of the spindown by invoking a third mass-dependent timescale that parameterizes internal angular momentum transport. We conclude that future models likely need to include three distinct physical processes to account for slow, fast, and low-mass rotators if they are to accurately describe stellar spindown.

One fact about stellar rotation appears to continue to hold in the aftermath of the newly-obtained data. Single stars continue to occupy a unique surface in rotation period-mass-age space. However, its shape is now found to be more complex than was indicated schematically in Meibom et al. (2015). The warmer part of the sequence appears to behave mostly Skumanich-like, albeit with a strong mass dependence. Current models already reproduce this behavior reasonably, and it is likely that a good description can be found using only one or two mass-dependent timescales. But it is not true for the cooler part of the surface where the behavior becomes more complex. It is likely, and the transition from the Barnes (2010) model to the Spada & Lanzafame (2020) model strengthens this thought, that the description of the cooler part requires the inclusion of at least a third parameter, an additional distinct mass-dependence of stellar spindown. Finally, it appears that even photometric binary stars of M 67 age are so diverse in their rotational behavior as compared with single stars that they seem to be unsuitable for gyrochronology at present.

Acknowledgments

We are grateful to an anonymous referee for insightful suggestions that helped to improve the quality of this paper. S.A.B. gratefully acknowledges support from the Deutsche Forschungs Gemeinschaft (DFG) through project number STR645/7-1. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts. This work has made use of data from the European Space Agency (ESA) mission Gaia (cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This paper includes data collected by the Kepler mission and obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the Kepler mission is provided by the NASA Science Mission Directorate. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. Based on photographic data obtained using The UK Schmidt Telescope. The UK Schmidt Telescope was operated by the Royal Observatory Edinburgh, with funding from the UK Science and Engineering Research Council, until 1988 June, and thereafter by the Anglo-Australian Observatory. Original plate material is copyright © of the Royal Observatory Edinburgh and the Anglo-Australian Observatory. The plates were processed into the present compressed digital form with their permission. The Digitized Sky Survey was produced at the Space Telescope Science Institute under US Government grant NAG W-2166. This work made extensive use of the python3 implementations of the numpy (Harris et al. 2020), astropy (Astropy Collaboration 2013, 2018, 2022), scipy (Virtanen et al. 2020), numba (Lam et al. 2015), and scikit-learn (Pedregosa et al. 2011) libraries. Indications of dwarf stars’ spectral types are based on “A Modern Mean Dwarf Stellar Color and Effective Temperature Sequence” (Pecaut & Mamajek 2013, and continuously updated afterwards, Version 2022.04.16; https://www.pas.rochester.edu/~emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt).

References

Appendix A: M67 cluster membership

In this section we describe the process that we used to determine cluster membership. For this, we use the stellar photometry (G and G − G_RP), proper motions (μ_Ra and μ_Dec), parallax ϖ, and radial velocity v_rad to the extent available in the sample. Based on a comparison with the cluster parameters (cf. Table 1), each star is then classified as either a member, candidate, or a field star in each of these parameter categories. Stars that lack a particular parameter are designated as unknown. Those individual designations are then combined into a single membership assessment.

For parallax and radial velocity, we define a range around the cluster value in which each star is designated a member when the measured value overlaps with that range within its error. A second, larger range is defined in similar fashion for the designation of candidate status (see below). For proper motions, we use a circular area around the cluster that is compared with the error ellipse spanned by the proper motion errors in the same way. The ranges around the cluster are defined as

(A.1)

(A.2)

(A.3)

(A.4)

(A.5)

(A.6)

An additional constraint for the parallax criterion is that the parallax needs to exceed its error by a factor > 3. Furthermore, we limit the proper motion offset to ±8 mas/yr. Table A.1 provides an overview of the sample sizes in each category determined by this process.

Because a cluster follows a complex sequence in the color-magnitude diagram (instead of being approximately represented by a single point), we need a more sophisticated approach. Instead, for each star we calculate the χ² distance to a cluster isochrone. We use a PARSEC isochrone corresponding to the age and metallicity of M 67 and apply dilution, extinction, and reddening according to Table 1. Similar to all other parameters, this χ² is then used to evaluate the photometric cluster membership.

Following the assessment as described above, we have evaluated four individual membership criteria. Figure A.1 shows the assigned category to the sample for each parameter. To combine them into a single, general membership designation, we assign a numeric value to each category. Here, member counts as 1.0, candidate as 0.5, and field star as 0.0. Those values are then added up for each star and divided by the number of available criteria for each star. If the resulting fraction for a star evaluates to a value ≥0.5, we designate that star as a member, wheras stars with lower values are designated as field stars. The result of this combination is displayed in the final column named M67 in Table A.1 and plotted accordingly in Fig. 2. Stars with fewer than two criteria available are treated as field stars. This is only relevant for stars with G > 20 mag, for which photometry is often the only available parameter.

Fig. A.1. M 67 cluster membership evaluation. The panels (a) to (d) show the individual parameters used for our membership analysis: photometry, proper motions, radial velocity and parallax, respectively. Each point in panels (a) and (b) is a Gaia DR3 target in the field of view. Panels (c) and (d) show histograms (gray) of the sample. The parameter numbering in the legend refers to the number of parameters that are available per star. The color coding indicates the assigned membership status. We note that this refers solely to the membership on the respective parameter of panel, not on the total membership. For an overview of the total membership, see Fig. 2. The dashed lines in panels (b), (c), and (d) indicate the parameter ranges for member and candidate designations based on Eqs. (A.1) to (A.6) (dark and light blue, respectively). The red circle in panel (b) indicates the 8 mas/yr cutoff.

Appendix B: K2 systematics correction

In this section, we describe the steps that were taken to extract and correct light curves from the K2 superstamp. In Sect. 3 have we already introduced the notation regarding intrinsic variability, instrumental systematics, and trending which we will continue to use here. The method for correcting the instrumental systematics was iterative, following the understanding of their nature and an exploration of different ways of approaching the problems present. Below we describe the approach that provided the best results. The method involves a detailed understanding of the systematics present in the data. We illustrate the process on the example of the sample star Gaia DR3 604895948360165888¹¹.

B.1. General description of the systematics

We begin with a description of the systematics because their precise characteristics arise from multiple sources, not apparent immediately. However, those details are crucial to obtaining good results.

As noted before, during the K2 mission several parts of the Kepler telescope ceased to function. Among those were the reaction wheels, crucial for the stability of the telescope’s pointing. Their malfunction caused the telescope to drift and to lose its pointing on the sky. This drift was then regularly corrected with the telescope thruster, pointing the telescope back to its intended location. The drift was small and firings were frequent to prevent targets from moving out of the field of view (aside from objects on the edges).

However, this constant motion away and back caused an object to meander across the detector. While the drift itself is slow and the stars do not leave visible trails on the images, they still fall on slightly different clusters of pixels on the detector. Figure B.1 shows the visible difference in the stellar positions on the detector between two different cadences. This shift is small; a little more than two pixels in a diagonal motion (at least for the K2 superstamp data). With the large pixels of Kepler (4″), however, this motion spans ≈10″ a significant amount in crowded regions, especially when one considers the overlap of the individual point spread functions (PSF).

Fig. B.1. Visualization of the two most extreme positions of the telescope drift during C05. Upper and lower panel both show superstamp images (Δt ≈ 47 d) indicating the peak-to-peak positional variation, that is, the minimal and maximal x coordinate of a star. The CCD pixel range is identical in both images. The red dots represent the mean position of a star throughout the observing run, and the filled red crosses connected to them mark the position at the time of the image. Red outlines crosses indicate the extreme opposite positions of a star. The white crosses indicate identical pixel coordinates in both panels to help to visualize the magnitude of the positional changes.

While motion alone could probably be corrected relatively easily, an additional effect introduces further complexity; the sensitivity of the detector is not constant across its surface. To be more precise, even on a subpixel scale, two areas of the CCD record different fluxes despite being illuminated by the same amount.

Considering both the constant drift and varying sensitivity, this means that the recorded flux of an intrinsically constant star is different depending on the epoch of observation. And since the individual regions of the detector are largely independent of each other, those changes in flux are different from position to position (and with that from star to star) on the CCD. Their systematic changes can be small and barely noticeable in a particular star while causing flux differences exceeding 15% for others. Applications that seek to identify variations in the mmag range are therefore difficult.

We are fortunately seeking stars that are not very variable and which change only over the course of tens of days, whereas the systematics occur of time scales of a few hours. As such, we can assume that our targets of interest are constant during a single drift, with only a small change in flux occurring from one drift to another.

However, the motion is more complicated than the description so far; it changes over the course of the campaign, drifts are unequal in length and amplitude, and there are jumps in the pointing. Therefore, we cannot adopt a mean flux from each drift; more precautions have to be taken to obtain reasonable light curves. Fortunately, the motions do allow us to identify common patterns which are then transferred into chunk-wise processing of the light curves. Images taken during the thruster firing (where the motion during the exposure is significant) exist, but those are fortunately rare (and only one or two at a time) and are consequently simply omitted by us. We also note that the sensitivity changes across the CCD are sufficiently well behaved to allow us to model them with a polynomial function, simplifying the problem.

We use aperture photometry for this work. Photometry that models the point spread function (PSF) generally allows disentangling stars in crowded regions to a certain extent, but there are two problems with the application of PSF photomtery to K2 data. The low spatial sampling and large FOV introduce complicated PSF shapes that vary across the detector, and even vary strongly depending on where a star is located on a subpixel scale. This alone can be solved with a sufficiently detailed empirical PSF, for example, see the PATHOS project (Nardiello et al. 2019; Nardiello 2020; Nardiello et al. 2021; Messina et al. 2022) that operates this way on TESS data.

However, given the special problems of K2 as described above we also face the issue that the PSF changes for a star from image to image, while also being essentially unique for each star in a particular image. The nominal strength of PSF fitting, the assumption that stars have similar PSF shapes despite differing in their fluxes, is not valid for these data.

B.2. Details of telescope drift

To obtain the information about the motion of the star across the detector, we utilize the world coordinate system (WCS) for each image provided by Cody et al. (2018). This allows us to calculate the position p₀(t) = (x₀(t),y₀(t)) of a star in each image. While we use the GDR3 ICRS coordinates for this, it would, however, not be a problem to use J2000 (or other) coordinates instead since we are only interested in relative changes.

For simplicity going forward, we will always refer to the telescope drift from the viewpoint of the star moving across the detector. As already noted, the motion of stars is not as simple as a straightforward back-and-forth motion. It is shown in Fig. B.2 for one example. The majority of the motion is in a diagonal direction, with an additional slow drift perpendicular to that, the latter jumping with some regularity. These jumps force us to divide the light curve into individual chunks of similar behavior. Since these are all on the same detector, the emerging patterns are the same for each star, allowing us to construct one mask that works for all.

Fig. B.2. Motion of the sample star Gaia DR3 604895948360165888 across the detector during C05. Each point represents the central position in one image color coded with the time of observation. Panel (a) gives the original detector coordinates and panel (b) those after the rotation applied. The extent of 1″ on-sky is indicated. We note the difference in scaling between the x and y coordinate axes, and the change from (a) to (b).

The first step is to simplify this motion by transforming the underlying coordinate system (x, y) to one where the majority of the motion is along one coordinate axis. For that we determine the gradients m_i between two consecutive cadences (p(t_i − 1) and p(t_i)) as

and take the respective median . We note that we purposefully do not take the mean value or use a linear fit to the data. Both of those approaches fail to provide good approximations of the motion due to the frequent jumps in the data. We then rotate the coordinate system around the angle such that

creates a new set of coordinates p′=(x′(t),y′(t)). Figure B.2 shows the new coordinates for same target. Note the differing scaling for y′ and x′ in the plot. The motions clearly follow patterns within certain segments of the data in both coordinates. This becomes even more obvious when we plot x′ and y′ over time as shown in Fig. B.3.

Fig. B.3. Detector coordinates of the sample star Gaia DR3 604895948360165888 over the course of C05. Panel (a) shows the x and y coordinates while panel (b) shows x′ and y′ after the rotation. Note that significantly smaller scale for y′ compared with x′.

It turns out that these individual segments of similar behavior are impossible to process together because their individual motion patterns are still too distinctive. Following this, we mask the segments and identify regions of similar behavior for collective processing. Based on the patterns, we mark 12 different groups of pixel behavior which we will henceforth refer to as segments. As can be seen in Fig. B.4, the individual segments are neither identical in size nor are they distinct in time. The latter fact will be advantageous for us below.

Fig. B.4. Motion of the sample star Gaia DR3 604895948360165888 across the detector during C05. The correspondence of each point with the segments is color coded. The gray shaded area is inspected more closely in Fig. B.7 below.

This pattern of behavior is indeed real and traceable in the image and not, for example, an artifact of the creation of the WCS. One can easily show this by combining the superstamp images into a movie, with the position of a star according to the WCS indicated and then observe the changes in the observable PSF of the star. These match one another, indicating that the WCS solution is an accurate representation of the positions (and therefore of the motion) of the stars.

There is one potential assumption we have not talked about yet which, if valid, would allow for an additional simplification or validation of our results. It is reasonable to assume that when a star falls two times on a very similar position, the detector response should be identical. With that, all observed flux variations between those two individual cadences should be free of the instrumental systematics, that means only composed of trending and the intrinsic signal. A closer inspection of Fig. B.4 shows that the distribution of positions of a star would allow this for at least some cadences (upper right region of the figure). However, upon further investigation, we notice that this assumption is not valid in most cases. This can be shown by selecting a star where the photometric noise is negligible and that exhibits a clear signal of variation (e.g., a M67 giant branch star). When we now compare the fluxes of two similar locations and compare those to the overall behavior of the raw light curve, we find that they do not match. We can identify the reason for that. Figure B.4 shows that most the of cadence pairs that we can use for this are made of cadences from different segments. If we limit ourselves to cadence pairs that are only from one segment, the assumption does indeed hold. However, there are too few of those available to be a significant help in the reduction process. At this point we have to look back at the origin of the detector coordinates. They come from the superstamp, which was created from the individual TPFs. Thus they are not necessarily connected to the actual physical pixel coordinates of the detector and may be shifted (by integer multiples of pixels). Judging from what we see, it is reasonable to assume that the distinction into the segments we see finds its origin in a jump in actual pixel coordinates and is not reflected in the superstamp coordinates. It may be possible to verify this assumption from the telescope telemetry itself or the individual TPFs. However, this verification (or rebuttal thereof) does not provide any additional value given that the segments are obvious in the data and therefore, we do not pursue this line of inquiry any further.

Table B.1 shows the extent of the individual segments. Within those segments, we can now express the flux f as a function of the coordinate x′ short of one more technical detail we will address below. However, we first need to describe the flux extraction from the superstamp images, because it also needs to be performed segment-wise.

B.3. Flux extraction from the full frame images

The FFI underwent background subtraction in the K2 pipeline. We see no reason to revisit this, and accept it as is. However, we will introduce an additional background subtraction to reduce the impact of the light from bright stars surrounding a particular target (as detailed below).

As described in Sect. 3, we use aperture photometry with an individually defined pixelmask rather than a fixed aperture or PSF photometry. A fixed circular or elliptical aperture differs too much from the very non-Gaussian PSF present in the FFI. Furthermore, taking only fractions of the flux stored in a given pixel based on the degree of aperture coverage introduces additional artifacts in the recorded flux.

To include all the relevant flux of the star, we need to account for the motion of the star as well. The peak to peak motion covers more than two pixels. In a situation where the PSF of a star of intermediate brightness roughly covers a cluster of 5 × 5 pixel, this is a significant effect and cannot be neglected. The brightest stars have PSFs about three times this size, whereas the faintest stars can be adequately covered by 2 × 2 pixel. Thus we need to account for the motion in the design of the aperture mask. For stars that are isolated, this does not present a problem as the aperture can just be designed large enough to cover all eventualities, though that would accumulate unnecessary noise. For stars in crowded regions, this is not a valid approach since the apertures begin to overlap and flux from neighboring stars is recorded.

Thus, we need an aperture that moves with the star. As already stated above, taking only fractions of the flux from a pixel introduces additional unwanted effects. Thus, we cannot have the pixelmask gradually moving with the star. Another approach that adds and removes entire pixel to and from the mask based on the current position of the star fails due to similar problems. Given the originally diagonal motion of the star and the fact that our determined segments separate mostly perpendicular to this motion, we use the following approach: we design different pixelmasks for different segments. It is, however, not necessary to design different pixelmasks for every segment. The coordinates allow different segments to share the same pixelmask. We refer to this grouping as a block, denoted with the letters A to E, for the five groups we found. Table B.1 includes the block designation in the column of the same name. We note that the blocks D and E only contain regions that will be cut due to the peculiar motion during the involved segments. Therefore, we focus on the blocks A – C from here on. The raw flux is extracted in D and E with a copy of the pixelmask for A.

We go on to design a pixelmask for each star, based on its visible shape in an image that averages all superstamp FFIs belonging to a certain block. This design is performed manually and iteratively. We endeavor to create a pixelmask for each star in each block that includes most of the flux of a star while simultaneously excluding flux from surrounding stars. This is, however, not always possible and so there are stars for which it is impossible to define a pixelmask that allows good extraction of the flux. We discard those from the sample. The design is iterative because the work showed that an incomplete pixelmask may leave strong imprints on the extracted flux that cannot be corrected. Thus, we iterate between pixelmask design and flux correction until we achieve light curves for the individual stars that are appropriate for our science case. In Figure B.5, we show the pixelmask for the example star in block A together with a comparison to a high resolution image from the Digitized Sky Survey 2 (DSS2¹²).

Fig. B.5. Pixelmask and background selection for Gaia DR3 604895948360165888 in block A. Panel (a) shows an averaged K2 superstamp image around the target. Panel (b) shows a comparison with a DSS2 (red channel) image of the same region. K2 pixels outlined in red mark the chosen pixelmasks pixels selected for the background estimate are marked with an X. The DSS2 image has the K2 pixelgrid indicated by the dotted lines and the solid line marks the edge of the superstamp.

Given the extent of manual work on each target in the FOV, we limit ourselves to targets of interest and a certain number of additional stars to verify the method. At this point we also discard all stars which cannot be separated well-enough from neighboring stars.

B.4. Extraction of the light curve and background subtraction

We can now extract the raw flux f_raw(t) based on the pixelmasks found by simply adding the flux from all included pixels for each cadence. Before doing so, we create a mean image for each block and determine all pixels fainter than the median in each. We assume that those pixels constitute the background. However, we do not use all those pixels for each star. From this selection, we use the ones closest to the extracted star. Here we sort the pixels by their distances to the mean position of the star in the block, exclude a potential overlap with the pixelmask itself, and include all pixels until we have at least three times as many background pixels as we have in the pixelmask. We then use the median of those pixels as the background and subtract it, weighted by the number of pixels included. This procedure allows the removal of stray light from stars in the vicinity. The number of pixels included and the use of the median are based on extensive testing of different variations and the final procedure presented appears to provide the best result overall. This provides us with the raw light curve f_raw(t) which is shown in Fig. B.6 for our example star. The typical systematic patterns due to the telescope motion are obvious. We note that due to differently shaped pixelmasks, the individual segments are generally not as well aligned as seen in Fig. B.6 but show offsets (cf. Fig. B.15). Those are eliminated when the light curves are realigned (see below).

Fig. B.6. Raw light curve fraw(t) of Gaia DR3 604895948360165888.

B.5. Correction of instrumental systematics

With the raw, but background subtracted light curves f_raw(t) at hand, we can now proceed to correct the instrumental systematics. This process continues to be on a star-by-star basis and is also carried out segment-by-segment. When we investigate the position dependency of the flux f_raw(x′) (see Fig. B.8) we see that the correlations are well-behaved enough to be modeled with a low order polynomial. To a large degree, the trends are almost linear but even a slight curvature can have a significant impact. We find that a 5th-order polynomial s the best compromise for the fit and use it to reproduce f_raw(x′) for each segment. However, one more difficulty has to be accounted for before doing so. When we directly fit a polynomial to f_raw(x′), it is very dependent on the range covered in x′ and the number of points responsible for said coverage.

At this point we introduce the a third (and last) masking variable which we will refer to as a slice. A slice is a short section of the light curve which covers one continuous episode of motion between two jumps. In Figure B.3, each apparent upward (or downward) stripe for x′ is distinguished as a slice. (There are jumps and gaps that are not apparent in the plot, causing the number of slices to be higher than obvious). By virtue of the design of the process, each slice in its entirety is part of only one segment. Figure B.7 provides a zoomed-in view of the light curve in Fig. B.4, with the slices distinguished by the color coding. Slices vary strongly in their extent for all involved parameters (x′, t, number of cadences). They typically cover 5 to 15 cadences, with the smallest containing only one individual cadence. Table B.1 includes a column that lists the number of slices in each segment.

Fig. B.7. Designation of the slices in a short light curve segment. Panels (a) and (b) show the time dependency of x′(t) and y′(t), respectively. The alternating color coding distinguishes individual slices and different symbols refer to different segments as indicated in panel (b). Only a small cutout from the full light curve is shown, corresponding to the shaded area in Fig. B.4. It contains 40 slices for the sample star Gaia DR3 604895948360165888.

We now utilize the slices that were introduced above. All slices of one segment exhibit the same systematic behavior (cf. Fig. B.8). However, given the fact that they are spread over several days (the segments cover typically about half the baseline of the light curve), even slow variations in the light curve leave their imprints in the mean flux of a slice. Furthermore, the individual slices have a different degree of coverage in x′. This complicates a polynomial fit to the data. We need an additional technical step to proceed.

Fig. B.8. Position dependence of the recorded raw flux fraw(x′) and the correction thereof. Panel (a) shows the raw flux from Gaia DR3 604895948360165888 as a function of the detector coordinate and clearly shows the linear dependence. Flux values are plotted for the same slices as in Fig. B.7. Each line represents one slice. Panel (b) shows the same slices but with the flux after the correction is applied. Panels (c) + (d) are similar to (a) + (b), only for the full light curve of Gaia DR3 604895948360165888 and the color coding representing the time as in Fig. B.2. Panels (e) + (f) are the same as (c) + (d) but for Gaia DR3 604917629355039360 for which the position dependency becomes nonlinear.

We designate the one slice which has the most extensive coverage in x′ in each segment as the prime slice. Each slice is then rescaled to this prime slice. We proceed to fit a 5th-order polynomial p₅(x′) to the resulting distribution of fluxes, and apply it as a correction to the individual slice (their unscaled fluxes!) as to create the first step in the correction process f_cor1. We multiply with the average flux of segment because the polynomial normalizes the fluxes to the segment mean. The result can be seen for a few selected slices in Fig. B.8. This processes is carried out for all slices in a segment and for all segments. We note that this process conserves long term flux changes and differences between individual slices. The resulting light curve for our example star Gaia DR3 604895948360165888 is shown in Fig. B.9.

Fig. B.9. Light curve of Gaia DR3 604895948360165888 after the application of the first correction for the instrumental systematics. Both panels show the same light curve. The color coding in the upper panel indicates the block as in Fig. B.5. In the lower panel, the colors indicate the segments as in Fig. B.4.

This process reveals that two regions of the light curve cannot be processed in this way. The first ≈3 d and a central part (t ≈ 172 − 177 d) correspond to exceedingly long slices that exhibit behaviors that cannot be found anywhere else and therefore those regions cannot be corrected this way (cf. the same regions in Fig. B.4). We mask those regions and cut them from the final light curves. Furthermore, there are individual cadences in the light curve from times of fast telescope motion (thruster firing). Those cannot be processed as well. We mask them in the correction step and later, when all other segments are processed, recreate those points from a linear interpolation in the corrected light curve. Those regions are entirely part of the blocks D and E (cf. the comment column Tab. B.1).

Thanks to the processing described above, the relative fluxes between the individual slices is conserved and the (fragmented) light curve of a segment may act as a valid light curve all by itself. However, this is not true when we compare the fluxes of the individual segments. They are not aligned. This is, however, not a result of our processing but a consequence of the instrumental systematics and different pixel masks that gives each segment a slightly different recorded average flux (cf. Fig. B.6). To align those fluxes, we merge the segments successively with a scaling to the chronologically first sector. The scaling factor is determined by minimizing the average standard deviation in the merged light curve for a series of windows 15 cadences wide along the light curve. This rescaling provides us with the next level of correction f_cor2 which is shown for example star in Fig. B.10. For the chronologically first segment we set f_cor2 = f_cor1.

Fig. B.10. Light curve of Gaia DR3 604895948360165888 after the application of the second correction, the realignment of the segments. Both panels show the same light curve. The color coding in the upper panel indicates the block. In the lower panel, the colors indicate the segments as in Fig. B.4. We emphasize that the flux axis is identical between the Figs. B.6, B.9, and this one.

The corrected and realigned flux f_cor2 now represents the light curve in which we have removed all instrumental systematics. However, this is not fully true for all light curves. A polynomial fit can only do so much and an insufficient pixelmask may introduce effects that cannot be corrected. As such, the process fails for stars in very crowded regions and on the edges of the superstamp. The slices generally cover a few hours. Intrinsic variations that occur on the same timescale may be misidentified by the polynomial. However, such rapidly varying targets are not of interest to us here, allowing us to ignore this problem. Despite this, we note that rapid signals are generally still visible reasonably well in the final light curves, thanks to their large amplitudes.

B.6. Cleaning the data

From this point onward, there is no further separation into blocks, segments, or slices, and the light curves are always treated as a whole. In the next step we apply σ-clipping to the data to clean it of outliers. This is done only at this stage, because the instrumental systematics may create or hide actual flux outliers and a σ estimate is dominated by it. For each point along the light curve we calculate the mean and standard deviation in a 1 d window around the point and replace it with a linear interpolation between its neighbors if it exceeds a 3σ deviation from the mean. We note that this may create artifacts for fast transients like transits or flares. However, since we are not interested in such phenomena, we elect to ignore those issues.

At this point we have constructed light curves for target stars that are corrected for instrumental systematics and cleaned of outliers. The light curves exhibit a gap of a few days around the middle and still include long term trending effects. The removal of the latter will be performed using Principal Component Analysis (PCA).

B.7. Principal Component Analysis

PCA is performed in a very similar manner to the light curve processing for Ruprecht 147 (from K2 C07, see Gruner & Barnes 2020, Appendix B). However, we modify the process slightly and describe the deviations below.

Because we do not have a large number of (∼10 k) light curves from the entire campaign for this dataset, we select stars for the PCA basis from our processed sample. We identify 769 stars for the said basis, covering a large range of brightnesses, colors, and locations on the detector (cf. Fig. B.11). From this basis, it becomes obvious that the trending signal is not universal for all stars. It shows a strong brightness dependence and a weak location dependence. We therefore limit the PCA basis applied for each target to stars of similar brightness and location. This means we take only the 125 nearest stars on the FFI within ΔG ± 1 mag. Both of these choices are the result of testing different parameter combinations; these particular ones provided the best results throughout the sample. We note that, depending on the brightness, there could be fewer than 125 stars overall in the brightness range and the PCA basis for such stars is correspondingly smaller.

Fig. B.11. Star selection for the Principal Component Analysis, indicating those used and those not used.

We smooth the light curves before entering them into the PCA to remove the impact of the noise from the analysis. This smoothing is done by replacing each flux value by the mean of a 1 d wide window around this point. We note that this smoothing is only in place to create the PCA correction and is not used beyond that. Furthermore, we normalize all light curves to their respective means before entering them into the PCA. This smoothing removes rapid-trending signals from the PCA input. This means that such signals are not accounted for in the PCA correction and thus remain in the final light curve. However, their impact is small with respect to our purposes and as such we do not address them any further. Traces of this can be seen in the light curve plots in Fig. C.4 and are generally of two different types:

1. A very rapid fluctuation, as seen for instance, between t = 150 d and t = 160 d for Gaia DR3 604914880575997056 and

2. A distinct short term signal, as seen for instance, between t = 150 d and t = 170 d for Gaia DR3 604915773929008384.

The (long-term) trending signals are dominant and the PCA can identify their different shapes relatively easily. We use only two components for the reconstruction; this is enough to remove the trending signals adequately. This is a consequence of the limited PCA basis (brightness and position), which preremoves very different trending signals from the basis. Figure B.12 shows the application of the PCA correction on the example star and on an essentially constant star.

Fig. B.12. Application of the PCA correction to the example star Gaia DR3 604895948360165888 (panels a and b), and the fainter and essentially constant star Gaia DR3 604909310002693632 (panels c and d). Panels (a) and (c) each show the result of the instrumental systematics correction (blue) and the PCA correction (red) whereas panels (b) and (d) show the final light curve after the PCA correction was applied.

At this point we are in possession of the final light curves f_fin which are corrected for instrumental systematics, cleaned for outliers, and are free of trending. Figure 4 summarizes the individual steps for the example star. This process is by no means perfect, and does not work for all targets. We only use a carefully selected subsample of all potential targets, and remove all light curves where residuals of the systematics or trending cause uncertainties. Most removed light curves suffer from an insufficient pixelmasking (mostly due to being in a crowded area in the superstamp FFI) or because the corresponding star coincides with a region of bad pixels. Both problems make corrections nearly impossible. Only for a few stars does the correction itself directly fail. For those, the flux does not behave well enough to be adequately modeled with a polynomial. However, their numbers are small enough not to merit the effort required for an adjusted correction procedure.

The above-described process is specifically designed for the K2 C05 superstamp around M67. However, it should be possible to adapt it for other parts of the K2 survey as well. The principal idea should hold and, assuming that the detector variability stays similar in its magnitude, provide good results without modifying the core components. The parts that would need to be modified are those that are fine-tuned to the superstamp. The pixelmasks for the individual targets are the first among these. Secondly, the masking (block, segment, slice) would also need to be adjusted to a different campaign. However, we can reasonably assume that our masks work for other regions of C05 without modification. Caveats such as the difficulties with rapidly varying stars and crowded regions do, of course, remain. A precondition for such work is the existence of a reference for the telescope motion, and with that of motion of the stars on the detector. The K2 TPFs do not provide this, but some of the other superstamps created by Cody et al. (2018) do. We note that the SPICE kernels for the K2 mission only include position and velocity of the telescope but not its orientation, a necessity to be able to leverage those products.

B.8. Performance of the reduction

The correction process described above works well for stars that are isolated and whose PSF does not exceed the FFI range, that is, stars for which we can define a pixelmask that includes all of the stellar flux and none from surrounding sources. This condition is not fulfilled on the edges of the superstamp, and especially in the central region of the cluster. Figure B.13 highlights some details that show the performance and limitations of our corrections.

Fig. B.13. Example light curves illustrating the capabilities and limitations of our correction process. Panel (a) shows the light curve of an eclipsing binary with a significant primary eclipse and a secondary eclipse that only becomes really apparent after the correction was carried out. Panel (b) shows a star where the instrumental systematics completely obscure the intrinsic signal despite their being of comparable amplitude. Panel (c) shows the light curve of a star located at the edge of the superstamp and which suffers from artifacts created by the correction process as a consequence. Each panel includes highlighted regions to show the details of the relevant effects.

Whenever the pixelmask cannot, for whatever reason, be created in a way that includes all the stellar flux the recorded raw flux f_raw includes a second systematic effect. This effect invalidates the original assumption we had to make, that the raw flux only changes because of variations in detector sensitivity, and is otherwise somewhat constant on short time scales. With that, it very much depends on the individual star’s location on the detector whether our correction still works or whether we under- or over-correct the apparent systematics. The star we used as an example for the correction process is located on the edge of the superstamp (cf. Fig. B.5) but is well-behaved enough to allow a good correction process. Meanwhile, Gaia DR3 604896429396491136, whose light curve is depicted in panel (c) of Fig. B.13, shows remnants of an over-correction in the highlighted area due to proximity to the superstamp edges and subsequent loss of flux.

B.9. Comparing our results with K2SC, SAP, EVEREST

We validate our light curves by comparing them with those produced as part of the original survey and from two commonly used works that implemented a systemtatics correction, K2SC, and EVEREST. This comparison again highlights one important fact that supports our approach with creating our own light curves – not many light curves are available for targets within the superstamp. Figure B.14 shows the distribution of stars with available light curves.

Fig. B.14. Spatial distribution of stars with available light curves in the archives for different data products. Panel (a) shows a map centered on M67 with the extent of the superstamp indicated (black box). Targets with archived light curves within (red and blue) and outside (black) of superstamp are overplotted. Panel (b) shows the same stars in a CMD.

Only 96 are within the field of view of the superstamp, as compared with GDR3 listing ∼2000 sources. This most likely originates from the differing nature of the superstamp as compared with the normal TPFs. The latter are typically associated with a proposed star and therefore are automatically processed. This is not true for the superstamp. As far as we are aware, all other projects that dealt with the K2 systematics only operate on the same sample. Therefore, other extant work does not venture outside this 96-star sample.

We extracted the 96 light curves in the FOV of the superstamp from the Kepler archive¹³. These include the light curves based on simple aperture photometry (SAP, f_SAP) on the TPFs and corrected light curves that resulted from the Presearch Data Conditioning (PDC, f_PDCSAP) module of the Kepler pipeline (Smith et al. 2020). The PDC is employed to remove data systematics and trending while retaining astrophysical signals. This generally works reasonably well but leaves traces of both artifacts. In some cases the trend correction removes astrophysical signal. For our comparison we employ the SAP flux as well as the PDCSAP with the aim of comparing the most similar products from each. For the former we compare with our raw extracted flux. However, we generally do not expect a large agreement in the shape of the raw light curves on short term scales due to the differently shaped aperture mask used. Long term effects are expected to be the same. The PDCSAP flux f_PDCSAP is compared with our final light curve f_fin (lower panel of Fig. B.15) as both have the same level of processing (systematics and trend removal), whereas the SAP flux f_SAP is equivalent to our raw flux f_raw (upper panel of Fig. B.15).

Fig. B.15. Comparison between the Kepler PDCSAP light curves and the ones created by us. In both panels, our light curves are shown in blue and the archival light curves are shown in red. Panel (a) shows the uncorrected aperture photometry fluxes (flux_raw from ours and SAP_FLUX for theirs). Panel (b) shows the final light curve product – flux_fin from ours and PDCSAP_FLUX for theirs.

We can see that, as expected, the raw light curves generally do not agree with respect to the short-term systematics while having the same long-term behavior. However, after the correction is applied, we find very good agreement between the Kepler and our pipeline results. Despite our relatively simple approach for an empirical correction, we not only do match the quality of the official product but partially exceed it. Generally, our approach is slightly better in the removal of the instrumental systematics.

The K2 Systematics Correction (K2SC) pipeline implements Gaussian processes to correct for the telescope jitter (Aigrain et al. 2015, 2016). It is geared toward exoplanet detection; as such the pipeline is set up in a way that also removes the astrophysical signal from the final product f_k2sc. However, they provide the removed systematics as part of their data product, differentiated in a position-dependent f_posi and a time-dependent part f_time. Long term variability is retained in the time-dependent part, together with trending. Recombining their final light curve with both those parts yields back the PDCSAP light curve (f_PDCSAP = f_k2sc ⋅ f_time ⋅ f_posi) which was the starting point for the K2SC pipeline. Thus, the K2SC light curves benefit from the long-term trend correction in the Kepler PDC pipeline. Technically, with their procedure they also implement a correction for instrumental systematics twice. To obtain a light curve that is corrected to a similar degree as our final product, we recombine the final flux and the time-dependent trending part. Figure B.16 shows the comparison. The additional systematics removal pays off, as the light curves are improved as compared with PDCSAP, and reach a fidelity similar to, or perhaps are even superior to ours. All 96 PDCSAP light curves are processed with K2SC. The final K2SC light curve (f_k2sc) is not comparable to any other, as it blatantly removes any long term signal, including physical ones (cf. upper panel of Fig. B.16). A comparison makes sense between our light curves before the PCA (f_cor2) and the K2SC one with the time-dependent signal included (f_k2sc ⋅ f_time, cf. lower panel of Fig. B.16).

Fig. B.16. Same as Fig. B.15 but for the comparison between the K2SC light curve of Gaia DR3 604917491916095872 and the ones created by us.

The EPIC Variability Extraction and Removal for Exoplanet Science Targets (EVEREST) light curves (Luger et al. 2016, 2018) use a pixel-level decorrelation (PLD) method that operates on pixel level light curves and, since its update to v2.0, includes the telescope motion for the correction. They are, as far as we know, the best available light curves based on K2 data and the study of variability longer than a few days. Only long-term trending is a remaining issue; such trends are generally more pronounced than in other light curves. However, with those retained, there is also no problem with accidentally removed astrophysical signal. We have worked successfully with them in the past (Gruner & Barnes 2020). However, they are even more limited in their availability regarding the superstamp region. Only three stars in the FOV of the superstamp have an EVEREST light curve. And all three of those stars are very faint. Two are white dwarfs, while the third is a very faint red dwarf at the brightness limit of Gaia (cf. panel (d) of Fig. B.14). All three stars are nearly constant in their light curves, and dominated by noise and some weak long-term trending. We use the red dwarf for the comparison in Fig. B.17. The extent to which we can correct for instrumental systematics is similar in both products, however, EVEREST reaches a slightly higher photometric precision. PDCSAP and K2SC fail to fully remove either trends or systematics, and their photometric precision is lower compared with ours or EVEREST for the faint stars.

Fig. B.17. Same as Fig. B.15 but for the comparison between the EVEREST, PDCSAP, and K2SC light curves of Gaia DR3 604898490980773888 (panels a, b, and c, respectively) and the ones created by us.

Appendix C: Supplements

C.1. Supplementary tables

C.2. The data provided

Table C.2 gives an overview of the most relevant sample of the stars for which we found periodic signals. As a supplement to this work we provide the extend version containing all 136 stars with additional columns in machine-readable form. The columns included are described in Table C.4. Additionally, we provide the light curves, including, raw fluxes, the intermediate steps, the PCA, and all masking.

C.3. Supplementary figures

Fig. C.1. Same as Fig. 3 but for B − V and V − K colors.

Fig. C.2. Color-period diagram (CPD, panel a) and color-magnitude diagram (CMD, panel b) of the complete superstamp sample. Blue symbols indicate the M 67 members as discussed above. Overplotted in red are the field stars for which we found a periodic signal. Note that the colors here are not reddening corrected. Open symbols denote binaries.

Fig. C.3. Same as Fig. 15 but for (GBP − GRP)0, (B − V)0, and (V − K)0 colors.

C.4. Sample light curves

Fig. C.4. Light curves for our sample stars. Each row in the figure is a star, with the name indicated. The left panel shows the light curve, the right panel their phase folded counterparts. The period used for the phase folding is shown in the upper right. The light curve panel displays the Gaia DR3 id in the lower left corner. Additional labels indicate whether a star’s light curve exhibits multiple spot features, i.e., is double dipping (DD, upper right corner in the left panel) or is a member of the half-period sequence (HPS, lower left corner of the right panel). The orange line indicates a 0.01 binning. Note that double dipping refers to all stars whose light curves exhibit signs of more than one spot. Stars are sorted as in Table C.2, by their Gaia Id.

Fig. C.4. continued.

Fig. C.4. continued.

Fig. C.4. continued.

Fig. C.4. continued.

All Tables

All Figures

Fig. 1. Overview of the K2 C05 superstamp coverage of M 67. The upper panel shows one superstamp image on a logarithmic gray-scale. Stars recognized by GDR3 are overplotted, color-coded according to our membership evaluation (see Sect. 2.2 and Appendix A for details). The lower panel shows a DSS2 (red channel) image of the same region. The extent of the superstamp is indicated in red. The small inset in the corner shows the location of the superstamp in the Kepler FOV C05. Both panels are shown such that their orientations are as similar as possible (see coordinates in panel b).

In the text

Fig. 2. Overview of GDR3 stars in the field of view and their M 67 membership status (following Sect. 2.2 and Appendix A). In all panels, blue stars denote members while non-members are gray. Panel a shows a color-magnitude diagram in Gaia G − GRP color, panel b the proper motions, panel c a histogram of radial velocities, and panel d a histogram of the stellar parallax. Dashed red lines indicate the position of M 67.

In the text

Fig. 3. Isochrone fit to the identified cluster members in two different Gaia colors. The PARSEC isochrone G-band magnitude is adjusted for the parallax ϖ and extinction AG using Eq. (2). The PARSEC GBP − GRP color is reddened by E(GBP − GRP) from Eq. (1). The dashed line indicates the nominal position of the equal-mass binary sequence (main sequence stars only). An alternative version using B − V and V − K is shown in Fig. C.1.

In the text

Fig. 4. Example of the light curve reduction process using Gaia DR3 604971466769552128. Panel a shows the raw light curve (blue dots), panel b shows an intermediate light curve after the systematics correction (blue dots) together with the PCA correction (red) that is applied to the create the final reduced light curve that is shown in panel c. The orange line in panel c shows the reduced light curve with 1.5 d binning to highlight the now visible intrinsic stellar variability. The insets in a and b zoom in on a 2 d span of the light curve. The drift-induced flux changes and the jumps caused by the realignment of the telescope are readily visible in panel a, while their disappearance can be seen in panel b.

In the text

Fig. 5. Result of period analysis on the light curve for Gaia DR3 604971466769552128. Panel a shows the light curve produced, one which exhibits a clear periodic signal. The power spectra in panel b are obtained using multiple methods as indicated. The Clean (CA) and Combined power spectra are normalized to their maxima and the PDM spectrum is shown as 1 − θ. Panel c shows the phase-folded light curve, folded with P = 27.9 d, corresponding to the peak in the combined spectrum.

In the text

Fig. 6. Stars for which we have identified rotation periods. Panel a shows a color-period diagram (CPD) and panel b the corresponding color-magnitude diagram (CMD). Stars belonging to M 67 according to our membership analysis are shown as blue circles (the main sequence stars) and red squares (post-MS stars). Open symbols (in all of the above) mark those that show signs of binarity. Error bars are suppressed here for visibility reasons.

In the text

	Fig. 7. Color-period diagram for our sample, now including rotation period uncertainties. Only member stars on the main sequence, and with no indications of binarity are displayed. An approximation of the emerging cluster distribution (solid gray line) is overplotted, together with its half-period counterpart (dashed line).
In the text

	Fig. 8. CPD for our M 67 sample emphasizing stars on the half-period sequence. The cluster sequence, approximated by an interpolated line, is overplotted in gray. The dashed line shows its half-period counterpart. Stars for which we double the measured periods (red) are connected to their new positions by dashed lines.
In the text

	Fig. 9. Final CPD of M 67 based on our study. The corresponding data are displayed in Table C.2.
In the text

Fig. 10. Comparison of our results with those of Barnes et al. (2016). Panel a shows a CPD as in Fig. 9 and overplotted in red are the measurements by Barnes et al. (2016). The gray lines, solid and dashed, are our indication of the cluster and half-period sequences, respectively. Panel b shows a CMD of the same stars and the inset (panel c) the proper motions centered on M 67. Gray dots in panels b and c are members without period determinations.

In the text

	Fig. 11. Same as Fig. 10 but with a comparison to results of Gonzalez (2016b; red), showing both retained and discarded stars from that sample.
In the text

	Fig. 12. Same as Fig. 10 but with a comparison to results of E18 (red), again identifying retained, discarded, and binary stars.
In the text

	Fig. 13. Same as Fig. 10 but comparing with the results of Dungee et al. (2022; red). Dungee et al. (2022) assume a flat period error of 10%.
In the text

	Fig. 14. Color-period diagram for our final assembled sample (combined with the results of Barnes et al. 2016), now including only main sequence stars with no indications of binarity. Only the slow rotator sequence is shown.
In the text

Fig. 15. CPD of the three open clusters 1 Gyr or older. Shown are NGC 6811 (Meibom et al. 2011a; Curtis et al. 2019) at 1 Gyr, NGC 6819 (Meibom et al. 2015) and Ruprecht 147 (Gruner & Barnes 2020; Curtis et al. 2020) at 2.5 Gyr, and M 67 (Barnes et al. 2016; Dungee et al. 2022; and this work) at 4 Gyr. The color coding groups cluster data by age. Overplotted is a simplified representation of the emerging cluster sequences by age (solid lines color coded by age). The sun is shown with its usual symbol. Figure C.3 shows equivalent plots in the colors (GBP − GRP)0, (B − V)0, and (V − K)0.

In the text

Fig. 16. Comparison between our final sample of measured rotation periods (this work combined with the results of Barnes et al. 2016; Dungee et al. 2022) and the stellar rotation models of Barnes (2010) and Spada & Lanzafame (2020). Each panel performs the comparison in a different color system: (a) G − GRP, (b) GBP − GRP, (c) B − V, and (d) V − K. The errors on B − V are suppressed for visibility reasons.

In the text

Fig. 17. CMD and three CPDs highlighting M 67 sample stars evolved beyond the Main sequence. Panel a shows a CMD of our sample with the stars color coded according to their masses. Differing symbols distinguish different types of binaries. Panel b shows a corresponding CPD in Gaia colors and identical symbols. Panel c plots the rotation periods of the stars against their estimated masses. An approximate indication of the position where stars cease to have an outer convective zones (i.e., the Kraft break) is overplotted with a dashed gray line. The solid black line corresponds to the scaling on the right and denotes the ratio between current (R(t)) and ZAMS stellar radius (R0). Panel d shows an inferred CPD with stars assigned colors they had when on the ZAMS.

In the text

Fig. A.1. M 67 cluster membership evaluation. The panels (a) to (d) show the individual parameters used for our membership analysis: photometry, proper motions, radial velocity and parallax, respectively. Each point in panels (a) and (b) is a Gaia DR3 target in the field of view. Panels (c) and (d) show histograms (gray) of the sample. The parameter numbering in the legend refers to the number of parameters that are available per star. The color coding indicates the assigned membership status. We note that this refers solely to the membership on the respective parameter of panel, not on the total membership. For an overview of the total membership, see Fig. 2. The dashed lines in panels (b), (c), and (d) indicate the parameter ranges for member and candidate designations based on Eqs. (A.1) to (A.6) (dark and light blue, respectively). The red circle in panel (b) indicates the 8 mas/yr cutoff.

In the text

Fig. B.1. Visualization of the two most extreme positions of the telescope drift during C05. Upper and lower panel both show superstamp images (Δt ≈ 47 d) indicating the peak-to-peak positional variation, that is, the minimal and maximal x coordinate of a star. The CCD pixel range is identical in both images. The red dots represent the mean position of a star throughout the observing run, and the filled red crosses connected to them mark the position at the time of the image. Red outlines crosses indicate the extreme opposite positions of a star. The white crosses indicate identical pixel coordinates in both panels to help to visualize the magnitude of the positional changes.

In the text

Fig. B.2. Motion of the sample star Gaia DR3 604895948360165888 across the detector during C05. Each point represents the central position in one image color coded with the time of observation. Panel (a) gives the original detector coordinates and panel (b) those after the rotation applied. The extent of 1″ on-sky is indicated. We note the difference in scaling between the x and y coordinate axes, and the change from (a) to (b).

In the text

	Fig. B.3. Detector coordinates of the sample star Gaia DR3 604895948360165888 over the course of C05. Panel (a) shows the x and y coordinates while panel (b) shows x′ and y′ after the rotation. Note that significantly smaller scale for y′ compared with x′.
In the text

	Fig. B.4. Motion of the sample star Gaia DR3 604895948360165888 across the detector during C05. The correspondence of each point with the segments is color coded. The gray shaded area is inspected more closely in Fig. B.7 below.
In the text

Fig. B.5. Pixelmask and background selection for Gaia DR3 604895948360165888 in block A. Panel (a) shows an averaged K2 superstamp image around the target. Panel (b) shows a comparison with a DSS2 (red channel) image of the same region. K2 pixels outlined in red mark the chosen pixelmasks pixels selected for the background estimate are marked with an X. The DSS2 image has the K2 pixelgrid indicated by the dotted lines and the solid line marks the edge of the superstamp.

In the text

	Fig. B.6. Raw light curve fraw(t) of Gaia DR3 604895948360165888.
In the text

Fig. B.7. Designation of the slices in a short light curve segment. Panels (a) and (b) show the time dependency of x′(t) and y′(t), respectively. The alternating color coding distinguishes individual slices and different symbols refer to different segments as indicated in panel (b). Only a small cutout from the full light curve is shown, corresponding to the shaded area in Fig. B.4. It contains 40 slices for the sample star Gaia DR3 604895948360165888.

In the text

Fig. B.8. Position dependence of the recorded raw flux fraw(x′) and the correction thereof. Panel (a) shows the raw flux from Gaia DR3 604895948360165888 as a function of the detector coordinate and clearly shows the linear dependence. Flux values are plotted for the same slices as in Fig. B.7. Each line represents one slice. Panel (b) shows the same slices but with the flux after the correction is applied. Panels (c) + (d) are similar to (a) + (b), only for the full light curve of Gaia DR3 604895948360165888 and the color coding representing the time as in Fig. B.2. Panels (e) + (f) are the same as (c) + (d) but for Gaia DR3 604917629355039360 for which the position dependency becomes nonlinear.

In the text

	Fig. B.9. Light curve of Gaia DR3 604895948360165888 after the application of the first correction for the instrumental systematics. Both panels show the same light curve. The color coding in the upper panel indicates the block as in Fig. B.5. In the lower panel, the colors indicate the segments as in Fig. B.4.
In the text

	Fig. B.10. Light curve of Gaia DR3 604895948360165888 after the application of the second correction, the realignment of the segments. Both panels show the same light curve. The color coding in the upper panel indicates the block. In the lower panel, the colors indicate the segments as in Fig. B.4. We emphasize that the flux axis is identical between the Figs. B.6, B.9, and this one.
In the text

	Fig. B.11. Star selection for the Principal Component Analysis, indicating those used and those not used.
In the text

Fig. B.12. Application of the PCA correction to the example star Gaia DR3 604895948360165888 (panels a and b), and the fainter and essentially constant star Gaia DR3 604909310002693632 (panels c and d). Panels (a) and (c) each show the result of the instrumental systematics correction (blue) and the PCA correction (red) whereas panels (b) and (d) show the final light curve after the PCA correction was applied.

In the text

Fig. B.13. Example light curves illustrating the capabilities and limitations of our correction process. Panel (a) shows the light curve of an eclipsing binary with a significant primary eclipse and a secondary eclipse that only becomes really apparent after the correction was carried out. Panel (b) shows a star where the instrumental systematics completely obscure the intrinsic signal despite their being of comparable amplitude. Panel (c) shows the light curve of a star located at the edge of the superstamp and which suffers from artifacts created by the correction process as a consequence. Each panel includes highlighted regions to show the details of the relevant effects.

In the text

	Fig. B.14. Spatial distribution of stars with available light curves in the archives for different data products. Panel (a) shows a map centered on M67 with the extent of the superstamp indicated (black box). Targets with archived light curves within (red and blue) and outside (black) of superstamp are overplotted. Panel (b) shows the same stars in a CMD.
In the text

	Fig. B.15. Comparison between the Kepler PDCSAP light curves and the ones created by us. In both panels, our light curves are shown in blue and the archival light curves are shown in red. Panel (a) shows the uncorrected aperture photometry fluxes (flux_raw from ours and SAP_FLUX for theirs). Panel (b) shows the final light curve product – flux_fin from ours and PDCSAP_FLUX for theirs.
In the text

	Fig. B.16. Same as Fig. B.15 but for the comparison between the K2SC light curve of Gaia DR3 604917491916095872 and the ones created by us.
In the text

	Fig. B.17. Same as Fig. B.15 but for the comparison between the EVEREST, PDCSAP, and K2SC light curves of Gaia DR3 604898490980773888 (panels a, b, and c, respectively) and the ones created by us.
In the text

	Fig. C.1. Same as Fig. 3 but for B − V and V − K colors.
In the text

	Fig. C.2. Color-period diagram (CPD, panel a) and color-magnitude diagram (CMD, panel b) of the complete superstamp sample. Blue symbols indicate the M 67 members as discussed above. Overplotted in red are the field stars for which we found a periodic signal. Note that the colors here are not reddening corrected. Open symbols denote binaries.
In the text

	Fig. C.3. Same as Fig. 15 but for (GBP − GRP)0, (B − V)0, and (V − K)0 colors.
In the text

Fig. C.4. Light curves for our sample stars. Each row in the figure is a star, with the name indicated. The left panel shows the light curve, the right panel their phase folded counterparts. The period used for the phase folding is shown in the upper right. The light curve panel displays the Gaia DR3 id in the lower left corner. Additional labels indicate whether a star’s light curve exhibits multiple spot features, i.e., is double dipping (DD, upper right corner in the left panel) or is a member of the half-period sequence (HPS, lower left corner of the right panel). The orange line indicates a 0.01 binning. Note that double dipping refers to all stars whose light curves exhibit signs of more than one spot. Stars are sorted as in Table C.2, by their Gaia Id.

In the text

	Fig. C.4. continued.
In the text

	Fig. C.4. continued.
In the text

	Fig. C.4. continued.
In the text

	Fig. C.4. continued.
In the text