BINARY FREQUENCIES IN A SAMPLE OF GLOBULAR CLUSTERS. I. METHODOLOGY AND INITIAL RESULTS

To perform artificial star tests and later CMD analysis, we need to define the main-sequence ridge line (MSRL) of the observed CMD first. We constructed the MSRL from the data instead of fitting theoretical isochrones to the CMD, as the latter method is not sufficiently accurate to define the ridge line. To define the MSRL from the CMD, we used a moving box method. In each step, we defined a box with a height of 0.1 mag in the F606W axis, and a width of the entire range in the color (F606W–F814W) axis. Then we searched the peak value of the color histograms in that box with a color bin-size of 0.0015 mag. The color value at the peak and the middle point of the box in the F606W axis direction are the defined MSRL for stars in that box. The process was repeated in increments along the main sequence and the resulting line was then smoothed with the moving average method with a period of 50. In Figure 3, the green line is the defined MSRL.

With the defined MSRL, we performed artificial star tests in Dolphot. The F606W magnitudes of those fake stars were randomly generated while the F814W magnitudes were derived from the MSRL, so that all those fake stars are on the MSRL. The (X, Y) coordinates of those fake stars were randomly chosen from the area of the drizzled reference frame, but were not at the empty non-data CCD area. Each fake star in this list was added to each FLT frame one at a time, and was recovered using the same photometry process by Dolphot. The final photometric output file was screened using the same criteria as used to obtain the CMD, as described above; we added 100,000 fake stars for each image. Figure 4 shows two examples of the artificial star tests, NGC 5053 (low stellar density cluster) and NGC 1851 (high stellar density cluster). The green line is the input fake stars, which are all on the MSRL. The scattered black dots are the recovered fake stars (69,995 stars for NGC 5053 and 52,205 stars for NGC 1851).

Zoom In Zoom Out Reset image size

From the input and recovered fake star lists, we studied the photometric accuracy and uncertainty. We compared the photometry for two clusters, NGC 5053 (low stellar density cluster) and NGC 1851 (high stellar density cluster), as their HST/ACS observations are very similar in exposure times: NGC 5053 (F606W: 340 s × 5, F814W: 350 s × 5) and NGC 1851 (F606W: 350 s × 5, F814W: 350 s × 5).

In Figure 5, we plot the magnitude differences for input and recovered stars for F606W filter (first row) and F814W filter (second row). These plots indicate photometric discrepancies, and all of them showing differences within ±0.2 mag. At fainter magnitudes, the magnitudes of the recovered stars tend to be fainter than the input ones. Rows 3 and 4 in Figure 5 show the photometric uncertainties for F606W and F814W filters, respectively, while row 5 in Figure 5 shows the uncertainties in color (F606W–F814W). From those plots we can see that most stars in NGC 1851 have similar photometric accuracy and uncertainties as in NGC 5053, except that in NGC 1851 there are more stars with large errors, even for bright stars.

Zoom In Zoom Out Reset image size

Row 6 in Figure 5 shows the completeness curve (recovery percentage), which is the ratio of the total number of recovered stars (N(out)) to the total number of input stars (N(In)) in a particular magnitude interval. The completeness curve for NGC 5053 is quite flat down to the 26.5 mag, indicating an even star recovery percentage over a broad magnitude range. However, for NGC 1851, the star recovery percentage is not evenly distributed, with a relatively higher recovery percentage for brighter stars than for fainter ones. This can be seen from the fake CMDs in Figure 4, where there are fewer faint stars in NGC 1851 than in NGC 5053, even though they have similar exposure times. This indicates that faint stars are difficult to recover in dense stellar region due to high and non-uniform background caused by the multitude of stars.

The most important source of contamination in identifying binaries arises from the superposition of two unrelated stars that occurs by chance when the star cluster is projected from three dimensions to two dimensions. Two or more stars along the line of sight and within the minimum angular resolution will be measured as one star, and are indistinguishable on the CMD from real binaries. It is very difficult to screen for blended stars, but statistically we can determine their distribution through Monte Carlo simulations. The blending fraction is proportional to the radial stellar number density of a cluster, decreasing from the core to outside region. The general way to measure the blending fraction is through artificial star tests, in which fake stars from MSRL are added to real images and are recovered with the same photometric processes as real stars. Any fake stars overlapping with real stars will represent a blended source, and the recovered magnitudes will be their sum. As long as the added fake stars follow the radial light distribution of a globular cluster, it represents the similar blending fraction as that cluster. The drawback for this method is that it is computationally time-consuming. Consequently, we performed the Monte Carlo simulations for an appropriate set of conditions and then modeled the results analytically, using Poisson statistics (see the Appendix).

In each simulation, we added fake stars to the real star map, where the fake stars have the same luminosity function and light distribution as the real ones, so it should have a similar blending fraction as the observed cluster. To generate the fake star list, we chose fake stars with their V magnitudes sampled by the observed luminosity function, and their I magnitudes calculated from the MSRL (i.e., all the fake stars are on the MSRL and have the same luminosity function as the real stars). The number of total fake stars added is equal to that of the total observed stars to produce a similar condition. Each fake star was assigned a radius r from the center of the cluster, which was sampled by the radial light distribution of real stars, so that the radial light distribution of fake stars is similar to the real ones. Then the azimuth angle is randomly assigned to the fake star to calculate the coordinates (x, y). To that position (x, y), any real star with distance less than 1 resolution element (the minimum resolution size Dolphot used) will be considered as a blended star, and the total magnitude is added. To find those blended stars, we used our k-dimensional tree algorithm by comparing lists of the observed stars and fake stars. To those blended sources we added photometric errors from the observed stars (see Figure 5). We repeated this simulation 30 times to obtain an average number or an average residual color distribution of blended stars.

Field stars, both faint foreground stars and bright background stars, can also contaminate the binary population on CMDs and thus affect the accuracy of the measurement of the binary fraction, especially for low Galactic latitude clusters. They affect the number of binaries in the binary region on the CMD, as well as the number of single stars on the main-sequence. High galactic latitude clusters do not have many contaminating field stars. For example, NGC 4590 (b = 3605) only has 24 field stars (simulated from the model of Robin et al. 2003) in the ACS field of view, or 0.06% of the total observed stars in that field. Low Galactic latitude clusters, however, have many contaminating field stars, such as NGC 6624 (b = −791), which has 51,325 field stars in the ACS field of view, about half of the total observed stars in that field. The best way to select cluster members is by proper motion, which, however, requires at least two epochs of HST observations separated by years. An alternative way is to construct the field star model from the theoretical model of the Galaxy to statistically account for those field stars. We used the Stellar population Synthesis model of the Galaxy (Robin et al. 2003) to simulate field stars at the cluster position, with the size of ACS CCD chips, 202'' by 202''. The V and I Johnson–Cousin magnitudes of the generated field stars are corrected for the reddening first and are converted into ACS F606W and F814W magnitudes by the transformations of Sirianni et al. (2005), then are randomly added to the original images along with Poisson noise. Then we adopt the same photometry method with Dolphot to recover those field stars. The recovered field stars will have similar photometry errors as the cluster stars, as well as the completeness and blending effect.

To estimate the high mass-ratio binary fraction, we divided the MSRL-subtracted CMD into different regions to count stars (see Figure 6, left panel). In Figure 6, the green line is the equal-mass binary population. The red line on the right side of the main sequence is the binary population with a binary mass ratio of 0.5. The binary mass ratio of 0.5 is chosen because in most cases, this binary population is beyond the 3σ photometric errors of the main-sequence stars. The dashed blue lines are the upper and lower limits of usable stars for both the F606W magnitude and the residual color. The main-sequence star (or single star) region S is defined between the red lines. The binary region B is where the residual color is greater than the 0.5 binary mass-ratio curve and less than 0.2 (i.e., on the right side of the red lines but on the left side of dash blue line). The residual region R is where the residual colors are beyond the blue side of the symmetric line for the 0.5 binary mass-ratio curve but greater than −0.2.

We count stars separately in those regions for three types of star populations, observed stars, simulated blended stars, and simulated field stars. The simulated blended stars and field stars are generated with methods mentioned in Sections 5.1 and 5.2. Region S contains mainly single stars (main-sequence stars), with some contaminating field stars. Region R contains main-sequence stars with large photometric errors, and with some contaminating field stars. Region B contains mainly binaries with mass ratios greater than 0.5, with some blended stars and field stars. It also has some main-sequence stars with large photometric errors, and can be estimated with the number in region R. So the high mass-ratio binary fraction fb(high q) can be calculated by the expression:

$$\begin{equation} fb({\rm high}\ q)=\frac{n^{{\rm obs}}_B-n^{{\rm blend}}_B-n^{{\rm field}}_B- \left(n^{{\rm obs}}_R-n^{{\rm field}}_R \right)}{n^{{\rm obs}}_t-n^{{\rm field}}_t} \end{equation} \tag{ 1 }$$

where $n^{{\rm obs}}_B$, $n^{{\rm blend}}_B$, and $n^{{\rm field}}_B$ are the star numbers in region B for the observed stars, blended stars, and field stars, respectively. $n^{{\rm obs}}_R$ and $n^{{\rm field}}_R$ are the star numbers in region R for observed stars and field stars, respectively. $(n^{{\rm obs}}_R-n^{{\rm field}}_R)$ represents any residual main-sequence stars with large photometric errors. This number should be also subtracted from region B assuming a symmetric distribution. The quantities $n^{{\rm obs}}_t$ and $n^{{\rm field}}_t$ are the total star numbers in the dashed blue line box for observed stars and field stars, which is the whole restricted region we use during the analysis. The 1σ error estimate for the binary fractions is estimated using Poisson errors and error propagation.

To avoid contamination from photometric errors, we followed three procedures. First, we only select stars 3–4 mag below the turn-off point to exclude faint stars with larger photometric errors. Second, we add two Gaussian models to represent the photometric errors during fitting, one for the main component with similar small errors, and the other for stars with larger errors. Tests show this to be a suitable strategy. Third, we introduce the parameter q_cut in the binary model, which represents the minimum binary mass ratio that can be extracted from the data.

To estimate this, we first fit the residual color distribution with two Gaussian models, which represent the photometric errors of the main-sequence stars. The fit is only applied from 0.02 on the red side to all the blue side, since the blue side is not affected by both blending stars and real binaries, the broadening is only due to photometric errors. The fit is quite good (with χ² close to 1). The positive residual of the fitting on the red side is due to binaries, with contamination from blending stars and field stars.

After subtracting the photometric errors, field stars, and blended stars, the residual color distribution on the red site is only from the contribution of physical binaries. We summed all the residuals on the red side and divided it by the total number of stars without the field stars, which gives the binary fraction including low mass-ratio binaries.

There is another way to account for low mass-ratio binaries, and one can even constrain the binary mass-ratio distribution function. Here, we constructed the binary population by fitting to an additional binary mass-ratio distribution function.

We assume the binary mass-ratio distribution function has the following power-law form: f(q)∝q^x, where q ≡ M_s/M_p, and with a minimum value of q_min, and a maximum value of 1 (equal mass binaries). The minimum mass for the secondary star is set to 0.2 M_☉ (the observed lowest mass from luminosity function of F606W band), thus this will set a minimum value of q_min not to be 0. First we assume here three different cases for the power of x, x = 0, a flat mass-ratio distribution; x = −1, which leads to a peak at low mass ratios; x = 1, which leads to a peak at high mass ratios. Whenever there are enough stars, we can fit for x rather than assign a value.

To construct the physical binary population, we first assume a binary fraction of f_b, or a total of N*f_b stars to be in the binary systems, where N is the total number of the observed stars. We then extract V magnitudes with the number N*f_b from the observed V magnitude luminosity function to be the V magnitudes of the primary stars. The V magnitudes of the secondary stars are calculated using a mass ratio q extracted from an assumed mass-ratio distribution and an assumed mass–luminosity relationship: L∝M^3.5. Their I magnitudes are derived from the MSRL, and the combined binary magnitudes are calculated. Each binary system has added to the photometric errors at their V and I magnitudes to simulate the photometric spreads. Finally, we can obtain the residual color distribution for the constructed binary population by subtracting the color of the MSRL. We repeat this process 30 times for each binary fraction f_b, and construct the average profile for the physical binary population as the way we applied it to blended stars.

Now, we have models for single stars (two Gaussian models), field stars (see Section 5.2), the superposition of stars (see Section 5.1), and binary model with the known fraction f_b (see Section 6.2). The total sum of all these models will produce the final model that is compared to the observed residual color distribution profile with a χ² test. For a given cluster, the models for single stars, field stars, and the blending stars are fixed as they depend only on the observed luminosity function, the Galactic positions, and the observed light distribution. The only model that changes during the fitting process is the binary model, as it varies with the binary fraction f_b, and the power law index x from the binary mass-ratio distribution function.

We developed the bisection method to search the best-fit f_b (i.e., fits with the minimum χ² value) at each assumed x. In this method, we first chose a wide initial range of binary fraction values (initial range is from 0 to 1). The χ² values (model comparing to observed counts) were calculated for three binary fraction values, left most, right most, and the middle. Then the χ² values at the left and right were compared to the middle one, and the left most or right most binary fraction value will be replaced by the middle one if its χ² value is greater than the middle one. So the search range for the binary fraction now is reduced by half, and we calculated the χ² value at the middle for the new range, and repeated the process again until the difference of χ² values or the binary fraction range approached limiting value. The final middle value of the binary fraction range is the best-fit binary fraction with the minimum χ² value. For clusters with enough bins, we also fit the power of x.

A fitting example is shown in Table 3, for NGC 4590. The first three rows in the table show the fitting results with the power x, and the minimum binary mass ratio q_min fixed, only with the binary fraction as a free parameter in the fit. The error range in the table is estimated by changing the parameter so that the χ² value changes by 1.0 (or 1σ confidence level). The best-fit value favors the model with the power x = −1, which also gives a higher value of binary fraction (10.8 ± 0.4)%. With the power x, and q_min free to fit, the fit improves significantly (χ²/dof = 88.4/82), with a binary fraction of (10.8 ± 0.3)%. Figure 7 shows the residual color distribution fitted by using only the Gaussian model (upper) and with the best-fit model (lower) for NGC 4590. The symmetric spread of the main-sequence is due to photometric errors, which can be fitted by the Gaussian model fairly well, as there are no large systematic residuals on the blue side. The asymmetric spread of the main-sequence is due to binary populations and blending of stars, which is shown as positive residuals on the red side on the upper panel. In the lower panel, models for binaries and blending of stars are included, and this best-fit model fits well to the observed data. Figure 8 shows the model components for the fitting (upper) and the enlarged view (lower) for NGC 4590. For this cluster, we estimate about 50 field stars in the field of view of HST, so they are a negligible contaminant. In Table 4, we show the binary fractions estimated with the three methods discussed above, the high mass-ratio method (the second column), the star counting method (the third column), and the χ² fitting method (the fourth column), and we expect the binary fraction obtained with the global methods (the star counting method and the χ² fitting method) to be higher than the high mass-ratio method, as we approach lower mass-ratio values.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this analysis, we divide the whole ACS field of view into three annular bins, with their centers at the cluster center obtained from Harris (1996; 2010 edition). The bin sizes were chosen (iteratively) so that each bin contains roughly one-third of the total stars recovered from the whole field, which leads to similar error bars for the binary fraction in each region.

For high density clusters, the CMD quality for the central bin is the worst among the three, which shows larger photometric spread and a lower faint star recovery rate than the other two. This is understandable, as the higher star not only increases the background level, making fainter stars more difficult to detect, but also increases the blending probability, making the PSF determination poorer. For low density clusters, we do not observe large variations in the CMD qualities.

In Table 5 and Figure 9, we show the results of the radial analysis on the binary fraction for NGC 6981 as an example. In Table 5, we list the sizes of the annular bins, binary fractions obtained by the high mass-ratio method, the counting method, and the χ² fitting method, and the dof/χ² for the fitting method. In Figure 9, we plot the binary fractions obtained from different methods and from different annular bins against their positions relative to the cluster center, in units of the half mass radius. We clearly see a decreasing trend of the binary fraction toward the outside of this cluster, an effect discussed for the sample of 35 globular clusters in Paper II.

Zoom In Zoom Out Reset image size

The key aspect to estimating the binary fraction using the CMD method is the highly accurate photometry for the CMDs. As the photometric errors become larger, the spread of the main-sequence becomes larger, which will smear out the signals from binaries with small mass ratio. In order to decrease the photometric errors, we need to increase the exposure time. This is true for low density clusters, in which stars are quite isolated, and the PSF can be determined quite well. For most low density clusters in our sample, the photometric errors are approaching their theoretical limit. High density clusters, however, have very crowded central regions, where the high and varied background make the PSF determination uncertain. Most of the errors are not from low S/N ratios but from uncertainties in the PSFs and the backgrounds. Increasing the observing time for those crowded clusters would not help lower the photometric errors. Instead, developing more sophisticated PSF photometry algorithm for crowded region is needed, which is still not fully mature.

Theoretical modeling shows that the dispersion of the metallicity of a globular cluster can also cause a spread in color on the main sequence. Figure 10 shows that a metallicity (Z) difference of 0.002 (from red to blue line at certain V magnitude), equivalent to δ[Fe/H] = 0.30, can cause a spread of about 0.054 mag in color, and 0.284 in the V magnitude. So given the uncertainty in [Fe/H] of 0.03, the spread in color will be 0.005, and shift in V magnitude will be 0.028. The observed intrinsic color spread is smaller compared to the typical width in color of the main-sequence with the HST observations (about 0.012 for NGC 5053 in our sample), and the shift in V magnitude will not affect the color distribution. Thus the intrinsic metallicity dispersion can be negligible. For multi-population systems (such as NGC 2808), however, this will not be the case (Piotto et al. 2007; Pasquini et al. 2011).

Zoom In Zoom Out Reset image size

Observations on globular clusters located near the Galactic center can be affected by the existence of large and differential extinction of the foreground dust. Alonso-Garcia et al. (2011) discuss a technique to correct this effect. From our sample, most clusters are well above the Galactic plane and with the reddening E(B − V) less than 0.1, which are not important to spread the color of the main sequence comparing to their photometric errors. For clusters near the Galactic plane, however, the reddening can be very large. Along with the heavily contaminated field stars, the determination of binary fractions in those clusters are quite uncertain.

At the same time as this work was being carried out, another group was working toward a similar goal and recently published their comprehensive work (Milone et al. 2012). Although both efforts follow established approaches that use the CMD, there are some differences, which we identify and compare the values derived from the two independent approaches. One important difference is the software used to obtain photometry in crowded fields, which makes extensive use of the PSF libraries. We used Dolphot (V1.2) while Milone et al. (2012) used the proprietary algorithms described by Anderson et al. (2008), developed specifically for crowded field photometry. We used the stellar field model of Robin et al. (2003) while Milone et al. (2012) used the model of Girardi et al. (2005).

Both we and Milone et al. (2012) performed extensive artificial star tests, although with slight differences in how completeness was defined and how finely the globular cluster stellar density was subdivided. Most other procedures were essentially identical, including spline ridge-line fitting to define the main sequence or the magnitude range of the main sequenced used for analysis. The two clusters described here have small reddening, so we did not need the sophisticated corrections applied by (Milone et al. 2012), nor were there multiple epochs of data to be considered.

The two clusters that we discuss here were also analyzed by Milone et al. (2012) and we find general agreement, although the results are expressed slightly differently. For NGC 4590, our binary fractions for q > 0.5 and within the half mass radius was (6.2 ± 0.3)% while (Milone et al. 2012) obtain a somewhat lower value of (5.3 ± 0.7)%. For the total binary fraction, (Milone et al. 2012) doubles the f(q > 0.5) value, which assumes a flat distribution in q, obtaining (10.6 ± 1.4)%. We fit a flat functional form to the CMD distribution and obtain a binary fraction of (9.4 ± 0.7)%. Since these are the same data sets, the differences, which are comparable to the uncertainties, most likely reflect systematic differences between the approaches.

The other comparison that can be made is the radial distribution of NGC 6981, which we find drops by about a factor of four to five from a bin within r_h to one that extends to 1.9 r_h (from (9.6 ± 0.6)% to (2.1 ± 0.3)%). This decrease is similar to Milones mean result for their sample (Milone et al. 2012) but for this particular object, they find a smaller decline, from about 5% to 3%, with error bars of about 1% for each value. The reason for this difference is not obvious to us. In a companion paper, we will compare our sample to theirs, which will provide the statistical power to identify significant differences and systematic effects.