The Frequency and Mass-Ratio Distribution of Binaries in Clusters I: Description of the method and application to M67.

Standard Gaia EDR3 photometry is available in three photometric passbands, $G_{BP}$ and $G_{RP}$ integrated from low resolution spectra, and $G$ measured from the astrometric field CCDs (Evans et al., 2018; Riello et al., 2021). Hereafter we refer to these as $B$, $R$ and $G$. We first made the flux-excess corrections to the $B$ and $R$ magnitudes as recommended by Riello et al. (2021). For our data, these corrections were almost negligible. There is a known problem with overestimated Gaia $B$ flux for faint red stars (see Section 8.1 in Riello et al. (2021)) so for our analysis we use colour-magnitude diagrams in $(G-R,G)$.

Initially we selected all Gaia EDR3 data from a cone of radius 1 deg centered on M67 (NGC 2682). This data included the cluster as well as many foreground and background stars. We next made a coarse cut on proper motion, selecting those stars that had a proper motion $(\mu_{\alpha},\mu_{\delta})$ that lay within a circle of radius 2 mas yr^-1 centred on (-11, -3) mas yr^-1 (Fig. 1). We next filtered out stars with a parallax uncertainty greater than 20$\%$ and fitted a 5-D gaussian in galactic latitude ($l$), galactic longitude ($b$), $\mu_{\alpha}$, $\mu_{\delta}$, and parallax ($\varpi$) to the remainder. By trial and error, we imposed a cut on this 5-D distribution at a level for which we were satisfied with a moderately clean colour-magnitude diagram, Fig 2. The parameter boundaries for the adopted cuts are listed in Table 1. This process has resulted in a very clean CMD, with less contamination than that of Yadav et al. (2008), albeit not as deep as that study. It appears similar in quality to the M67 Gaia CMD derived by Gao (2018) using machine-learning methods.

The final CMD shows a main sequence with an obvious population of binary stars distributed up to 0.75 mag above the main-sequence ridge line, plus a few outliers that may be triple systems, or evolved binaries. Binary systems with two equal-mass main-sequence stars appear 0.75 mag above the main sequence. Binaries with lower mass ratios appear redder and fainter than the equal-mass case, see for instance the tracks in Fig. 3 of Elson et al. (1998).

There is some evidence for mass segregation in M67. Fan et al. (1996) found a half-mass radius of 7 arcmin for blue stragglers, 10 arcmin for upper main sequence stars, and 12 arcmin for the lower main sequence. This is confirmed from radial velocity studies (Mathieu & Latham, 1986; Geller et al., 2021). We might therefore expect main-sequence binaries to be more centrally concentrated than single stars. Our analysis covers stars drawn from a radius of $\sim 30$ arcmin, so effectively averages over any mass segregation of binaries. The tidal radius of M67 is estimated to be 42 arcmin (Fan et al., 1996; Francic, 1989), and its core radius is 8.24 arcmin (Davenport & Sandquist, 2010).

For the analysis presented below, we restricted ourselves to considering stars that are clearly on the main sequence or are binaries consisting of two main-sequence stars. We have manually removed the few stars that are clearly above the equal-mass-binary main-sequence. Additionally we have retained only stars where the primary (or only) component lies within $13.5<G<18$. These cuts are shown as magenta curves in Fig. 2, and were found necessary in order that our binary star distribution can be modelled below with the same primary-star mass function as the single stars.

Each 2-dimensional data point on the CMD, $(G-R,G)$ is formed from G and R measurements with independent uncertainties, but the $G-R$ combination is not independent of $G$. Consequently a $(G-R,G)$ datum has an associated covariance matrix

where $\sigma_{G-R}^{2}=\sigma_{G}^{2}+\sigma_{R}^{2}$ and it can be shown (see Appendix A) that the off-diagonal elements, $\sigma_{G-R,G}=\sigma_{G}^{2}$. As viewed in the $(G-R,G)$ plane, each data point can be considered as having a tilted bivariate gaussian probability distribution.

To model the stellar population of the cluster, we rely on having an isochrone that accurately describes the main sequence magnitude and colour as a function of mass. We used isochrones from the MESA Isochrones and Stellar Tracks project¹¹1http://waps.cfa.harvard.edu/MIST/ (Dotter, 2016; Choi et al., 2016; Paxton et al., 2011; Paxton et al., 2013, 2015). After some trial and error, we adopted a 5 Gyr isochrone with $[{\rm Fe/H}]=+0.06$ as being the best fit to the upper main sequence and turnoff. This is consistent with the metallicty determinations of Hobbs & Thorburn (1991) ($[{\rm Fe/H}]=-0.04\pm 0.12$) and Önehag et al. (2014) ($[{\rm Fe/H}]=0.06$). The age, however, is long compared with most recent determinations, which are generally close to 4 Gyr. The age inferred for M67, however, depends on the adopted theoretical isochrone calculations (see Fig. 10 of Yadav et al. (2008)), in particular the treatment of convective overshooting. In Fig. 2 we shown the M67 Gaia CMD with MESA isochrones for 4, 5 and 6 Gyr and $[{\rm Fe/H}]=+0.06$. We based our fit on the main-sequence turnoff region, and found that only the 5 Gyr isochrone provided a correct width for the subgiant branch and matched the "hook" at the top of the main sequence. The specific choice of isochrone is inconsequential for the purpose of this paper since we will eventually marginalise over all parameters associated with the mass distribution.

The MESA isochrones deviate from the observed main-sequence ridge line near the bottom of the main sequence. We introduce a colour correction to the isochrone below $G=16.5$, implicitly adopting the $G$ vs mass relation from the isochrone.

Depending on the quality of the photometry and the chosen passbands, a given CMD will have a threshold in $q$ below which binaries are practically indistinguishable from single stars. In Fig. 3 we show sections of our final trimmed CMD data, along with the colour-corrected isochrone for the main sequence and tracks representing binary stars with a fixed primary mass and different values of $q$. We indicate the locations of binary stars with $q=$ (0.4, 0.5, 0.6). For our data, $q=0.5$ is a sensible threshold to adopt, as we perceive that such binary stars are sufficiently displaced from the main sequence to be recognised as such. The model we will develop in the following section does not depend on this threshold, but it is important for the interpretation of the model outputs.

It is common in observational studies for measurement uncertainties to be underestimated (or overestimated). In practice, this is equivalent to assuming that the underlying true distribution has an intrinsic width. Such a width could arise from effects such as differential reddening or small systematic errors in the photometry. We allow for the possibility of an intrinsic width or incorrect error bars in our model by adopting a free parameter, $h$, that scales all data error bars. In general this manifests as $\textbf{{S}}_{G,G-R}\rightarrow h^{2}\textbf{{S}}_{G,G-R}$ in the probability calculations that follow. From examination of Fig. 2, it appears that the data may become more noisy towards the lower main sequence, so we allow $h$ to vary linearly along the main sequence as $h=h_{0}+h_{1}(G-G_{0})$, where we adopt $G_{0}=16.0$. Error bar scaling thus introduces two parameters to our model.

We begin with the distribution of outliers, since that is the simplest of the three individual likelihoods. In what follows, we adopt the notation

for a normalised bivariate gaussian centred at $(0,0)$ with covariance matrix S, evaluated at $x$. We also note that a probability density in $x$ should strictly be denoted by the derivative, $dP/dx$. For brevity, we adopt the looser notation, $P(x)$.

From the standard rules of probability we can expand the outlier likelihood function

The first term in the integral is the (gaussian) probability of the given data point ${\bm{D}_{k}}$ given a "true" CMD location ${\bm{D}^{\prime}}$. The second term is the probability distribution on the CMD for outliers, for which we adopt a very broad bivariate gaussian. Consequently, the integral becomes a convolution

and the likelihood is thus a bivariate gaussian. In principle, ${\bm{D}_{O}}$ and ${\bm{S}_{O}}$ could be regarded as elements of ${\bm{\theta}}$, but in practice it is sufficient to set them as constants. Here we adopt ${\bm{D}_{O}}=(0.75,16)$ and ${\bm{S}_{O}}={\rm diag}((0.75)^{2},(4.0)^{2})$. Generally $\textbf{{S}}_{k}$ is negligible compared to ${\bm{S}_{O}}$ and can be removed from the final covariance with no effect.

The above description is certainly an imprecise representation of the outlier distribution but, as noted by Hogg et al. (2010), much of the power of a mixture model such as this one comes from simply including an outlier distribution, even if it is inaccurate in detail.

We assume that the true locations of single stars in the CMD lie along the isochrone line, and are drawn from some mass distribution. This mass distribution is truncated at the top of the main sequence, and fades away more gradually at the bottom of the main sequence due to declining Gaia detection efficiency with increasing G magnitude. We model this observational mass distribution as

Here, the first term is a power law in $M$. The second term is a logistic function that acts as a smoothed step, centred at $M_{0}$ and with characteristic width $1/k$, to control the bottom of the mass distribution. The top of the distribution is subject to a sharp cutoff by the the third term, a Heaviside step function at mass $M_{\rm max}=1.186\,{\rm M}_{\odot}$, corresponding to G=13.5. The normalisation parameter $C_{M}$ is computed numerically so that $\int_{-\infty}^{\infty}P(M|\gamma,k,M_{\rm 0},M_{\rm max})dM=1$. For our M67 data, we finally opted to apply a hard cut near the bottom of the main sequence, so could have replaced the logistic function with a second step function, however we retain this form for generality of the model.

Since there is a one-to-one mapping, $D^{\prime}(M)$, from a given $M$ to a location on the CMD the single-star likelihood function can be expanded as

We can perform this integration by breaking the mass distribution function into a sum of discrete mass steps of width $\Delta M$ and height $P(M_{l}|{\bm{\theta}})$ located at masses $M_{i}$,

In practice, results from using this expression are equivalent to employing the expression we will develop in the following section for binary stars, with a mass ratio set to be close to zero.

We assume that the more massive star, $M_{1}$, in each binary system is drawn from the same distribution as the single-star mass function, Eqn 10. Each secondary star has a mass $M_{2}=qM_{1}$, where $q$ is drawn from a parameterised probability density. After testing various forms for this distribution, we adopted a general cubic polynomial,

where $\tilde{P}_{k}(q)$ are the shifted Legendre polynomials, an orthogonal basis set over the range $0<q<1$. This general form allows the distribution to bend up or down at either or both ends. The form of the distribution as written requires no normalisation since

By including the derivative parameters $\dot{a}_{k}\equiv\partial a/\partial M$, the shape of the distribution is allowed to vary with mass, around some reference value $M_{\rm ref}$, which we choose as the centre of the mass range covered by the CMD main sequence.

We consider this general model for the distribution as described, and also the more restricted versions with $\dot{a}_{k}$ set to zero.

To compute the likelihood, we represent the probability densities for $M_{1}$ and $q$ as linear combinations of 50 predefined gaussian basis functions,

where the centroids, $M_{1,i}$ and $q_{j}$ are equally spaced across the mass and mass-fraction ranges, $0.1<=M_{1}/M_{\odot}<=1.1$ and $0<q<=1$. The gaussian widths are set to $\sigma_{M_{1}}=0.01$ and $\sigma_{q}=0.01$, which we found by trial-and-error to produce a good representation of the underlying functions. The choice of gaussians for these basis functions will become apparent below.

It is a requirement of the formalism that follows that all $a_{i}$ and $b_{j}$ are zero or positive. The coefficients $a_{i}$ and $b_{j}$ can be computed analytically as linear fits to numerical representations of Eqns 10 and 15. However, the optimal linear fit coefficients are not guaranteed to be positive. Instead we use the non-negative least squares algorithm from Lawson & Hanson (1974), which iterates to a solution.

A model magnitude for a binary star is obtained by adding the fluxes of the two components in the appropriate bandpass. This is done by interpolating magnitudes for masses $M_{1}$ and $qM_{1}$ from the isochrones, converting magnitudes to fluxes, adding the fluxes, and reconverting to magnitudes. This procedure is performed separately for $G$ and $R$, then the resulting values are combined to form $(G-R,G)$.

Each combination of $(M_{i},q_{j})$ forms a normalised bivariate gaussian in $(M,q)$ space with covariance matrix, $\textbf{{S}}_{M,q}={\rm diag}(\sigma_{M}^{2},\sigma_{q}^{2})$. This transforms to a covariance matrix in $(G-R,G)$ space, $\textbf{{S}}_{G-R,G}=\textbf{{J}}\textbf{{S}}_{M,q}\textbf{{J}}^{T}$, where the Jacobian

The partial derivatives in the Jacobian are not dependent on the data or model parameters and can be pre-computed for any given $(M,q)$ by interpolating colours and magnitudes from the isochrone. The combination of all the $(M,q)$ bivariate gaussians map to a set of tilted overlapping bivariate gaussians in $(G-R,G)$ space.

Similar to what we have done previously, we expand the binary-star likelihood

This can be further decomposed as

where ${\bm{D}_{ij}}\equiv{\bm{D}}(M_{i},q_{j})$ and $\textbf{{S}}_{ij}\equiv\textbf{{J}}\textbf{{S}}_{M_{i},q_{j}}\textbf{{J}}^{T}$. Here, the choice of gaussian bases for $M$ and $q$ becomes apparent. The integral is again a convolution that can be evaluated analytically, so that

Expanding Eqn. 4,

Inserting the individual likelihoods developed above, we can write this in a compact nested form,

where LSE is the log sum exp function,

and

As described in the preceding sections, our general model has 13 parameters, ${\bm{\theta}}=(\gamma,k,M_{0},a_{1},a_{2},a_{3},\dot{a}_{1},\dot{a}_{2},\dot{a}_{3},f_{\rm B},f_{\rm O},h_{0},h_{1})$. To compute the posterior probability distribution for ${\bm{\theta}}$, we must define a prior distribution, $P({\bm{\theta}})$. We assume that this distribution is separable for most parameters, except $a_{k}$ and $\dot{a}_{k}$, which we couple to limit the total change in shape of $P(q)$ along the main sequence. Additionally, we require $h_{1}$ to be positive, and scale its potential distribution with $h_{0}$. In the absence of any compelling previous knowledge, and after some trial and error, we generally adopt sensible uniform, normal or truncated-normal distributions as follows:

where $M_{\rm min}$ is the minimum main-sequence mass, $\Delta M$ is half of the mass range of the CMD main sequence, $U(a,b)$ is the uniform distribution between $a$ and $b$, and $\mathcal{N_{T}}(a,b,c,d)$ is the normal distribution centred at $a$ with standard deviation $b$ truncated below $c$ and above $d$.