DEATH OF A CLUSTER: THE DESTRUCTION OF M67 AS SEEN BY THE SLOAN DIGITAL SKY SURVEY

Traditionally in the absence of more detailed information, such as proper motions, radial velocities, or spectra, cluster membership can be crudely estimated based on the positions of stars in a color versus magnitude diagram (CMD). To maximize the contrast of the cluster against the surrounding field stars we must characterize the probability functions from both the field and cluster in color–magnitude space. These two CMD probability functions are then used to determine the probability of membership for every star.

The M67 CMD was determined by selecting the stars within 025 of the cluster center. Figure 2 shows the M67 CMDs for g − i versus g and r − z versus r in our SDSS data. The main sequence is clearly visible in both panels, and a faint equal-mass binary sequence can be seen above it. We chose to use the colors g − i and r − z because they better separate the M67 main sequence from the significant field contamination. The galactic field stars were crudely filtered out by rejecting any stars not within a few photometric σ of a spline hand-fit to the main sequence and equal mass binary sequence. A two-dimensional Gaussian smoothing algorithm was then applied to these rough cluster CMDs, creating a continuous CMD distribution for M67.

Zoom In Zoom Out Reset image size

The number of cluster stars in a given solid angle is described by the equation $n_{\rm {cl}} = \alpha f_{\rm {cl}}$, where α is the number of stars in a region, and $f_{\rm {cl}}({\rm color,mag})$ is a normalized probability function for the cluster in the color–magnitude plane. By normalizing the continuous CMD distributions described above, we created the $f_{\rm {cl}}$ functions for M67. This was done for both (g − i, g) and (r − z, r) independently, and are shown in Figure 3. The non-uniform distribution shown in Figure 3 is a result of sampling $f_{\rm {cl}}$ from the stars in the core of M67 which are not evenly distributed along the main sequence. R02 investigated whether the mass function of the cluster core would detract from the matched filter's ability to detect tidal debris that would generally have a different mass function. Their conclusion is that the method is robust against such biases, but we do bear this in mind in our analysis later on.

Zoom In Zoom Out Reset image size

The same task must be carried out for the field to describe the background filter. It is clear from Figure 1 that the field does not have a uniform density of sources in our sample. Instead we see the density of stars rise at lower galactic latitudes, near the galactic plane. To investigate the effect the changes in the field star population would have on our analysis, we sampled the background population in several regions. Figure 4 shows the background (g − i, g) CMD in four quadrants of the field. Other than the number of stars in each quadrant, the differences between these four CMD samples are minor, and we believe the changes in the CMD will not greatly affect our analysis, as was also concluded in R02. Thus, we make the assumption that the background has a uniform CMD structure across the field studied, and that the stellar density may be described by a low-order surface fit, as described in Grillmair & Johnson (2006).

Zoom In Zoom Out Reset image size

The background filter is referred to as $n_{\rm {bg}}({\rm color,mag})$, which, given any area on the sky provides the number of background stars in each color versus magnitude bin within that solid angle. It is valid to use $n_{\rm {bg}}$ and not $f_{\rm {bg}}$ here because the background is a continuous and relatively smooth population, and the number of stars in any area can be estimated. By using all of the stars around M67 with a radius between 15 and 8° (1,011,776) our background population is very well sampled. This annulus was chosen to avoid both the cluster and any potential tidal features, and the large hole in the data set in the southeast corner. Applying a boxcar smoothing algorithm and dividing by the area used, we create the two $n_{\rm {bg}}$ filters shown in Figure 5.

Zoom In Zoom Out Reset image size

Because the field population does not have a uniform distance or age, the familiar features of an ideal CMD (e.g., main sequence, turnoff, red giant branch, etc.) become severely blurred. The resulting color versus apparent magnitude diagram is known as a Hess diagram (e.g., see Alcock et al. 2000). In the Hess diagrams in Figures 4 and 5, several distinct features are visible in the CMD contours. The peak at (g − i, g) = (0.6, 17) is from the galactic disk (the thick disk according to R02). This feature is strongest in Figure 4(d), which has the lowest galactic latitude and thus should contain the most disk stars. The smooth feature at g − i ∼ 0.5 is the main-sequence turnoff for the halo population, which is naturally most prominent in Figure 4(a), at the highest galactic latitude of our sample. The large red buildup of stars at g − i ∼ 2.5 is the local field K and M-dwarf population, most of which are likely part of the disk. This is also most prominent in Figure 4(d). We have enforced magnitude cuts at both the bright and faint limits of our SDSS sample, and thus the $f_{\rm {cl}}$ and $n_{\rm {bg}}$ distributions have steep declines near the limiting magnitudes.

To produce the optimal contrast, the normalized cluster filter is divided by the scaled background. This yields $h=\frac{f_{\rm {cl}}}{n_{\rm {bg}}}$, which is shown in Figure 6. This produces a CMD matched filter that describes the stars which stand out from the background the strongest. For any given star, its color and magnitude will yield a value of h related to its probability of membership in M67. Naturally this probability function will promote regions of the CMD where there are many M67 members but few field stars, and conversely punish stars with CMD positions having large field populations.

Zoom In Zoom Out Reset image size

The spatial distribution of the cluster is then mapped by summing the values of h inside small (01, 01) spatial bins. We complete this entire process independently for both the (g − i, g) and (r − z, r) CMDs, and co-add the results, as per Grillmair & Dionatos (2006a). While the probability distribution h produces the maximum likelihood indicator for the stars, corrections for the field response to the h filter must be done. To solve for the actual number of cluster stars found in each spatial bin, α, the following formula from R02 is used:

$$\begin{eqnarray} \alpha &=& \left\lbrace \sum \left[\frac{f_{\rm {cl}}}{n_{\rm {bg}}}\right] - \int f_{\rm {cl}} d({\rm color,mag})\right\rbrace \!\!\bigg/ \nonumber \\ && \left\lbrace \int \frac{f_{\rm {cl}}^2}{n_{\rm {bg}}} d({\rm color,mag})\right\rbrace \,, \end{eqnarray} \tag{ 1 }$$

where $\frac{f_{\rm {cl}}}{n_{\rm {bg}}} = h$ is summed for all stars in a spatial bin as mentioned above, the integral of $f_{\rm {cl}}$ accounts for the background response in a solid angle dΩ, and the integral of $\frac{f_{\rm {cl}}^2}{n_{\rm {bg}}}$ is the signal-to-noise response of the filtering. This equation is the basis for our analysis with both the SDSS and 2MASS data sets, as well as the Monte Carlo simulations described in the following section.

In order to test the robustness of the matched filter method, we ran our algorithm using artificial clusters of known stellar composition and structure. By characterizing the efficiency and reliability of the algorithm we were able to examine issues of biases and errors for our results in Section 4. To accurately measure the efficiency and determine the greatest sources of error, we inserted artificial clusters at random locations within our field. These tests were run 1000 times on our SDSS data set in order to determine the reliability and significance of low-contrast features seen in Section 4.

The simulated clusters were created with a large range of core densities. The number of members was randomly chosen between 450 ⩽ N_STARS ⩽ 1350, while the core radius was randomly selected between 09 ⩽ R_core ⩽ 22'. This wide range was used to fully explore the densities and spatial compositions that open clusters might have.

The spatial composition of each cluster was created by randomly generating radial positions for every member, each being drawn from a Gaussian profile with a standard deviation radius of R_core. The angular positions were then chosen randomly between 0° and 360°, creating a roughly uniform circular distribution of stars. No tidal elongation was explicitly included; however, the Gaussian distribution of stars did at times produce asymmetric halos at large radii consisting of up to a few dozen stars.

The griz magnitudes were formed using the handmade splines fit to M67 in Section 3.1. A luminosity function for the g band was used which had the form g = (g^0.77_*) × 8 + 15, where g_* is an array of N_STARS elements filled with random values between 0 and 1. This function was chosen to provide an increasing number of stars with decreasing luminosity, but to have a deficiency as compared to the field star population of low-mass stars, as observed in old open clusters. We also found that it approximately reproduces the observed luminosity function of M67. An equal mass binary sequence was created by increasing the magnitudes of a subset of stars by 0.75 mag in every band. The fraction of equal mass binaries was chosen at random between 10% and 60%.

Appropriate photometric scatter was created by sampling errors for the actual SDSS data. The mean error and the standard deviation of the errors were calculated for griz bands in bins of 0.1 mag. These errors were then mapped to the artificial stars by their corresponding griz magnitudes, choosing for each star a random number within the Gaussian error envelope. Adding these errors to the artificial photometry created a very realistic scatter about the M67 splines which increased with increasing photometric magnitude.

For every model run, the simulated data were placed randomly in the field around M67, but was required to be greater than 15 radially away from M67 to avoid cross-contamination, and less than 8° away to stay well within the bounds of the SDSS data sample. The full matched filter analysis code was then run, sampling $f_{\rm {cl}}$ from a 025 radius around the model cluster and smoothed as before with the M67 $f_{\rm {cl}}$. The background $n_{\rm {bg}}$ was determined from an annulus around this with a radius between 15 and 6°, encompassing ∼500,000 stars on average, but taking care to avoid M67.

For every run of the simulation we recorded the simulated cluster position, surrounding background density, the estimated number of cluster members recovered by our algorithm, the radial surface density profile for each run, as well as every plot used in the analysis of M67. In Figures 7 and 8, we show example results for a random subset of our models. The resulting statistical analysis of the model runs show that we are able to reliably recover both the number of members and their locations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 shows the mean recovered radial surface density profile for our 1000 model runs, marked by stars. The squares show the mean radial surface density profile of the 1000 input models themselves. We are on average overestimating the surface density in the central bins, but are recovering the core radius well. Triangles mark the mean residual between the input models and the recovered profile. The error bars are one standard deviation in each residual bin. The small overestimation of the radial profile in the nucleus is attributable to the spatial smoothing kernel which will indicate more stars in the central radial bin than are actually observed.

Zoom In Zoom Out Reset image size

In every run of our models, the radial surface density profile is fit with a least-squares Gaussian profile. From the 1000 model profiles we found that, on average, our algorithm underestimated the core radius by 11%, and after correcting this bias we found a standard deviation of 060 for R_c. The number of stars detected by our algorithm was larger than the input models by, on average, 17%, and we computed a standard deviation on the number of stars recovered to be 67 stars. We also determined that the input binary fraction was 1.27 times more than our recovered binary fraction, and found a standard deviation of 6% for the resulting binary fraction. None of the recorded statistics changed as functions of the background density, suggesting that our matched filter algorithm consistently was able to remove the field contamination. We are attempting to improve these biases for future implementation of our algorithm.

Since each model run was required to be at least 15 away from M67, in some cases the artificial cluster can be close enough to M67 for it to be present on the radial profile plot at its furthest extent. This can be seen in Figure 8(a) for instance. These examples provide a useful benchmark for the robustness of the matched filter. Despite low-order changes in the local background population surrounding each run, and the variations in the model $f_{\rm {cl}}$ compositions and model core densities, M67 is always nearly perfectly recovered. The elongated features of M67 seen in these examples match those found in the analysis of the following sections.

Visual inspection of a subset of the models, as demonstrated in Figure 8, showed no examples of spurious extra-tidal features which deviated from the Gaussian distribution caused by fluctuations in the background population. Any elongated features appeared to be caused by the input model. Further, our algorithm itself, while slightly overestimating the number of stars in the clusters' core, does not produce random cluster halo signatures. We therefore conclude that the matched filter method is effective in characterizing and removing the field star population and providing an accurate representation of the cluster within the range of our sample. Any extra-tidal structures or elongations found in the vicinity of M67 are therefore considered to be real and intrinsic to the cluster.

By analyzing SDSS data for M67 using the technique outlined in Section 3.1, we were able to map the cluster distribution to a larger radial extent than has been done before. Figure 10 shows our spatial map of the open cluster M67. A small (01) boxcar smoothing kernel has been applied to the summed α data, and a low-order surface fit to the background was subtracted to remove any residual large-scale variation of the field density. The core of M67 is strongly detected and shows a circular distribution. The over-plotted circle, 25' in radius, denotes the furthest radius that BB03 detected the cluster against the surrounding field population using 2MASS data. There is however a significant low-density asymmetric halo of stars well outside the core. The contours from light to dark are increasing levels of density. They are defined by the equation level_j = MED(α) + (σ_α × 2^j), where MED(α) is the median value of α over all spatial bins, σ_α is the standard deviation of α, and j ≡ {1, 2, 3, 4, 5, 6}. To ensure that α was well characterized, we created a histogram of α from every 01 × 01 spatial bin in our SDSS sample. This is shown in Figure 11, and appears Gaussian in shape and is centered around α = 0. Since our contours begin at two standard deviations above the peak shown in Figure 11, and noting the reliability determined in Section 3.2, we argue that these features, discussed in Section 4.2, represent real structure around M67 that has previously gone unseen.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The radial surface density profile of M67 is shown in Figure 12. The strength of detection in each radial bin is the sum of α values in that bin, divided by the area that the bin encompasses. This is directly analogous to the spatial map in Figure 10, and the units of (α arcmin⁻²) are equivalent to (stars arcmin⁻²). The extended halo seen in Figure 10 is visible from a radius of ∼25' to ∼50' where it falls to background levels.

We fit this surface density profile using a King-like profile (King 1962) to our data, which has the form $n(r) = n_{{\rm bkgd}} + \frac{n_0}{1+(r/R_{c})^2},$ where R_c is the core radius, n_bkgd is the background surface density (not to be confused with the background filter $n_{\rm {bg}}$), and n₀ is the central peak surface density. This surface density profile equation was used previously by BB03 to model M67. This fit is shown in Figure 12 as the solid lines. Our King model fit provides R_c = 824 after correcting for biases from Section 3.2, somewhat larger than the previous 2MASS based measurement of R_c = 486 by BB03. From our testing with the model cluster, we estimate the error to be 060. We believe the difference between our radial profile and that of BB03 is due primarily to the lower limiting mass our study probes. We do give a direct comparison to BB03 by using the matched filter algorithm on the M67 2MASS data and provide more discussion in Section 4.3. After fitting the first-order King model used in BB03, we attempted to determine the tidal radius for M67, R_t. Using the equation n(r) = n_bkgd + n₀{[1 + (r/R_c)²]^−1/2 − [1 + (R_t/R_c)²]^−1/2}² from King (1962) and the values we determined for the first-order fit above, we found the estimated tidal radius for M67 to be R_t = 64' ± 21'. Errors here are determined from standard deviation of the rms scatter on our fit. This corresponds to a radius of 16.8 pc at the distance of M67, much larger than the determination using bright stars in M67 by Piskunov et al. (2007). By summing the number of stars in each radial bin in Figure 12 we estimate the total number of visible M67 members to be 1385 ± 67, where our error is again adopted from our models in Section 3.2.

Zoom In Zoom Out Reset image size

Having found evidence for a large halo around M67, we investigated the differences between the halo and core populations. Variations in the types of stars in the inner and outer regions of the clusters might suggest an origin for the elongated halo. In dynamically unevolved systems, the initial stellar luminosity function rises toward lower masses. BB03 use this property to investigate the evidence for mass segregation in M67. They created a set of 2MASS J-band luminosity functions in three regions (core, inner-halo, and outer halo). A higher fraction of low-mass and faint J-band stars were seen in the outer halo, and thus mass segregation is implicated.

Binary star systems are composed of two stars which are often too close to be visually resolved. This single point source is observed to be brighter than a single star of the same apparent color. These systems also have higher mass than single-star systems by definition. The effects of mass segregation discussed above are therefore applicable to binary stars which act as more massive single stars. Previous cluster studies have observed a stronger concentration of binary stars in the cores than in the outer envelopes of older clusters (e.g., Fan et al. 1996).

The binary fraction and luminosity function for each region of our data can be probed in much the same way as the spatial structure. Rather than summing α, the number of potential member stars, in each spatial bin, we summed α in three concentric regions for every magnitude bin (for the luminosity function) or across the main sequence (for the binary fraction). These three regions were the core, halo, and a large background annulus centered around M67.

Calculating the binary fraction across the entire magnitude range requires properly collapsing the CMD into the color plane. We employed a technique used by Clark et al. (2004) to collapse the main sequence of the globular cluster Palomar 13. Since a spline has already been fit by hand along the M67 main sequence in Section 3, we subtracted the main sequence spline color from every star's color, thus centering the main sequence on ΔColor = 0. Because the equal-mass binary sequence is 0.75 mag brighter than the main sequence, its ΔColor_bin position was easily tracked. We then divided ΔColor by the ΔColor_bin of the binary sequence, creating a reduced color R ≡ ΔColor/ΔColor_bin. This placed the main sequence at R = 0 and the equal mass binary main sequence at R = 1. To calculate the binary fraction we summed all of the α values in each reduced color bin, with a range of −1 ⩽ R ⩽ 2. The luminosity function is found by summing all of the α values in every magnitude bin, again cutting out stars outside the reduced color range −1 ⩽ R ⩽ 2.

In Figures 13–15, we investigate the M67 core, halo, and background populations, respectively. The central panel shows the reduced CMD, where the darker pixels show increasing α sums in each CMD bin. The right panel shows the sum of α values along the magnitude axis, while the bottom shows the same along the reduced color axis. Since the equal-mass binary main sequence was created by translating the main-sequence spline 0.75 mag brighter, a discontinuity in the calculation of ΔColor_bin arises at the end points. Thus, we have reduced our magnitude range by 0.75 on both the bright and faint end to avoid these biases.

The core population in Figure 13 (radius <08) shows a clear double Gaussian peak along the reduced color R-axis, and a deficiency of faint (low-mass) members along the magnitude axis. The halo in Figure 14 (0.8 < radius <12) contains a more flattened luminosity function and a diminished secondary peak in the binary fraction. These two figures however do not properly account for the contribution the field population would have, especially in the halo which samples a larger area and has a lower expected density of cluster members. Thus, in Figure 15, we made the same measurements for a much larger sample of stars (15< radius <6°) to estimate the field contamination in our binary fraction and luminosity function. This shows the response to the h filter that the background has.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The background contribution was normalized and then scaled by the area in the core and halo regions, respectively. This bias contribution was then subtracted from the M67 core and halo samples. Figure 16 shows the background subtracted binary fraction estimations for the core and halo of M67. A Gaussian profile was fit to the primary peak for R < 0.25. The residual in the core was nicely fit by a second Gaussian profile, centered at R ≈ 1. Integrating the two Gaussian profiles, we measure a fraction of binary point sources to be 21% in the core of M67. This is comparable to the Fan et al. (1996) determination of 16% and 22% from Montgomery et al. (1993). Correcting for the artificial cluster tests in Section 3.2 above, our estimation of the binary fraction is increased to 26%. This is necessarily a lower limit on the fraction of binary sources as many stars may be binary systems without detectable flux excess. No significant Gaussian binary star residual was found in the halo population, although the signal is noisier. Both the core and halo samples go to α = 0 for R < 0 and R>1. The luminosity function is given the same treatment, and the background corrected functions are presented in Figure 17. The halo shows a more flattened low-mass contribution. The core however contains a significant lack of low-mass members, even after accounting for the background field contribution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The dichotomy between the core and halo populations is consistent with qualitative ideas about mass segregation. As discussed in Section 1, an old cluster such as M67 is expected to be losing low-mass members to the outer regions of the cluster, while the core becomes increasingly more concentrated with higher mass stellar systems. This simplification, of course, assumes a single static population of stars and also ignores the effects of classical stellar evolution, whereby higher mass stars die before lower mass, and skew the observed present day mass function. Assuming an age for M67 of ∼4 Gyr, it would be expected that late A or early F stars would be leaving the main sequence, and thus would be brighter than our limiting magnitudes with SDSS. We do not therefore expect stellar evolution to contribute significantly to the observed luminosity functions of our sample. An et al. (2008) investigated the effects of crowding and stellar density on the photometric completeness of the SDSS pipeline, including the M67 region. They found that the automated pipeline recovered a comparable number of photometric sources as compared to manual reductions of the imaging for the open clusters M67, NGC 2420, and NGC 6791. Thus, we conclude that while manual reductions of wide-field imaging such as in our sample would be preferable, the SDSS pipeline is able to reliably probe the regions around high galactic latitude open clusters.

To provide a complete comparison of our matched filtering method with previous results for M67, we ran our algorithm on J-, H-, and K_s-band photometry from 2MASS. Using the 2MASS interface available on the internet entitled Gator,⁶ we retrieved a 10° × 10° box surrounding M67. This produced 264,354 point sources. The cluster CMD probability distributions $f_{\rm {cl}}(J-H,J)$ and $f_{\rm {cl}}(J-K_s,J)$ were drawn from all stars within a radius of 025 around M67, and the background CMD from everything with a radius greater than 15. Our Interactive Data Language (IDL) code for the SDSS data set was modified to use the 2MASS JHK_s data, and the analysis was carried out to measure the spatial distribution of M67. We did not measure the binary fractions or luminosity functions for the halo and core populations with the 2MASS data.

M67 as seen by 2MASS contains a different range in spectral types than is observed by the SDSS, reaching significantly higher masses. This is shown in Figure 2 of BB03, with the M67 CMD reaching a spectral type K0 at the faint limit, which we adopted to be J = 17. In addition to the upper main sequence, the 2MASS CMD also contains the turnoff and red giant branch. The average stellar density in our 2MASS sample was 0.73 stars arcmin⁻², whereas the average in our SDSS sample was 1.4 stars arcmin⁻² over the same spatial area.

Figures 18 and 19 summarize our 2MASS results for M67. These are analogous to Figures 10 and 12, respectively. Asymmetric distributions of stars are seen in Figure 18 which extend beyond the detection of the BB03 study. The elongation seen in Figure 10 is not well reproduced by the 2MASS data, however, the mass ranges are considerably different.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Two significant differences between the SDSS and 2MASS results are apparent in Figure 19: the 2MASS data yield a significantly higher surface density in the core, and the core radius fit with a King profile is much smaller (R_c = 412). This core radius corresponds very well with determinations by BB03 and Fan et al. (1996). Because the 2MASS sample contains higher mass stars than the SDSS sample, the discrepancy with Section 4.1 in core radius is expected. The higher mass stars in the 2MASS data will be more centrally concentrated in the core. These are too bright to be observed in the SDSS data and thus the SDSS will probe a mass range with a larger core radius.

Mass segregation differentiates all of the stars in a cluster according to their mass. The SDSS and 2MASS core radius comparison provides two estimates with a large separation of mass. To test that mass segregation was still evident in the 2MASS data alone we split our 2MASS sample into a bright bin and a faint bin. The separation was chosen to be J = 12.5 which roughly corresponds to the main sequence turnoff. The full matched filter was re-run for each 2MASS subset. The bright and faint samples yielded core radii of 370 and 467, respectively, indicating that mass segregation among the high mass sample is clearly evident, as found in Section 4.2 above.

The asymmetric halo around M67 shown in Figure 10 appears elongated roughly in-line with the proper-motion vector. Fan et al. (1996) have previously suggested an elongation of the core in M67 which roughly matches what is seen in Figure 18, but does not resemble what is seen in Figure 10. A strong correlation cannot be made between the Fan et al. (1996) results and our work because Fan et al. (1996) considered asymmetries through eight angular bins covering 45° each. Their optical data were also limited to a photometric depth of V = 21 and 0.3 mag errors making comparison to our SDSS data tenuous. Qualitatively the core radius and tidal features of Fan et al. (1996) match our 2MASS analysis of M67. Our matched filtering has been able to reproduce known structure for M67 using 2MASS, and revealed a slightly increased radius of detection for the halo compared to BB03 due to our more complete subtraction of the contaminating field population.