Smoking Gun of the Dynamical Processing of the Solar-type Field Binary Stars

TL;DR

This study analyzes the binary star properties in the solar neighborhood, revealing how dynamical processing in star clusters influences binary fractions and mass-ratio distributions across different stellar populations.

Contribution

It introduces a hierarchical Bayesian model to connect binary properties with metallicity and primary mass, highlighting the impact of early dynamical processing on binary star evolution.

Findings

High-$\gamma$ populations have lower binary fractions and are more affected by dynamical processing.

Binary fractions correlate with metallicity in high-$\gamma$ groups and with primary mass in low-$\gamma$ groups.

Dynamical processing preferentially destroys binaries with smaller mass and wider separation.

Abstract

We investigate the binarity properties in field stars using more than 50\,000 main-sequence stars with stellar mass from 0.4 to 0.85\,$M_\odot$ observed by LAMOST and {\emph Gaia} in the solar neighborhood. By adopting a power-law shape of the mass-ratio distribution with power index of $\gamma$, we conduct a hierarchical Bayesian model to derive the binary fraction ($f_{b}$) and $\gamma$ for stellar populations with different metallicities and primary masses ($m_1$). We find that $f_b$ is tightly anti-correlated with $\gamma$, i.e. the populations with smaller binary fraction contains more binaries with larger mass-ratio and vice versa.The high-$\gamma$ populations with $\gamma>1.2$ have lower stellar mass and higher metallicity, while the low-$\gamma$ populations with $\gamma<1.2$ have larger mass or lower metallicity. The $f_b$ of the high-$\gamma$ group is anti-correlated with [Fe/H] but flat with $m_1$.Meanwhile, the $f_b$ of the low-$\gamma$ group displays clear correlation with $m_1$ but quite flat with [Fe/H]. The substantial differences are likely due to the dynamical processing when the binaries were in the embedded star clusters in their early days. The dynamical processing tends to destroy binaries with smaller primary mass, smaller mass-ratio, and wider separation. Consequently, the high-$\gamma$ group containing smaller $m_1$ is more effectively influenced and hence contains less binaries, many of which have larger mass-ratio and shorter period. However, the low-$\gamma$ group is less affected by the dynamical processing due to their larger $m_1$. These are evident that the dynamical processing does effectively work and significantly reshape the present-day binary properties of field stars.

Tables

Table 1. Table 1: Best-fit f b subscript 𝑓 𝑏 f_{b}  and γ 𝛾 \gamma values at different [Fe/H] and m 1 subscript 𝑚 1 m_{1}  bins.

[Fe/H]	$m_1$			$m_1$$=0.65$			$m_1$
dex	N	$f_b$	$\gamma$	N	$f_b$	$\gamma$	N	$f_b$	$\gamma$
(-0.65,-0.40)	2109	$0.27\pm 0.03$	$0.60\pm 0.27$	3074	$0.35\pm 0.03$	$0.31\pm 0.22$	2903	$0.31\pm 0.03$	$0.63\pm 0.25$
(-0.50,-0.25)	2874	$0.25\pm 0.02$	$0.94\pm 0.28$	6234	$0.33\pm 0.02$	$0.56\pm 0.16$	6172	$0.39\pm 0.02$	$0.48\pm 0.16$
(-0.35,-0.15)				7801	$0.16\pm 0.01$	$2.74\pm 0.29$	9401	$0.22\pm 0.01$	$1.77\pm 0.21$
(-0.25,-0.05)				7984	$0.13\pm 0.01$	$2.98\pm 0.42$	12583	$0.16\pm 0.01$	$2.21\pm 0.29$
(-0.15,0.05)				5211	$0.15\pm 0.02$	$1.30\pm 0.49$	10578	$0.15\pm 0.02$	$1.61\pm 0.40$
(-0.05,0.15)				2138	$0.20\pm 0.07$	$0.38\pm 0.76$	5953	$0.13\pm 0.02$	$1.92\pm 0.59$
(+0.05,0.40)							2212	$0.22\pm 0.12$	$0.02\pm 0.92$
[Fe/H]	$m_1$			$m_1$			$m_1$
dex	N	$f_b$	$\gamma$	N	$f_b$	$\gamma$	N	$f_b$	$\gamma$
(-0.65,-0.40)	2273	$0.27\pm 0.03$	$0.85\pm 0.32$
(-0.50,-0.25)	5085	$0.44\pm 0.03$	$0.30\pm 0.16$	3906	$0.43\pm 0.04$	$0.44\pm 0.20$
(-0.35,-0.15)	6036	$0.37\pm 0.03$	$0.66\pm 0.19$	5263	$0.60\pm 0.06$	$-0.08\pm 0.17$	3761	$0.53\pm 0.05$	$0.24\pm 0.21$
(-0.25,-0.05)	8910	$0.20\pm 0.01$	$1.95\pm 0.29$	6680	$0.40\pm 0.04$	$0.35\pm 0.21$	6485	$0.78\pm 0.08$	$-0.45\pm 0.15$
(-0.15,0.05)	10184	$0.13\pm 0.01$	$2.98\pm 0.48$	7505	$0.18\pm 0.02$	$1.88\pm 0.38$	8259	$0.57\pm 0.06$	$-0.25\pm 0.18$
(-0.05,0.15)	7501	$0.10\pm 0.01$	$3.15\pm 0.75$	5881	$0.11\pm 0.05$	$3.18\pm 1.36$	8277	$0.40\pm 0.05$	$-0.06\pm 0.22$
(+0.05,0.40)	4167	$0.10\pm 0.02$	$2.98\pm 1.09$	3660	$0.09\pm 0.02$	$3.39\pm 1.46$	5953	$0.45\pm 0.11$	$-0.47\pm 0.31$

Table 2. Table 2: Summary of the relationships in γ 𝛾 \gamma -low and high groups.

Condition		$\gamma >1.2$
$f_b$ vs. [Fe/H]	no relation	anti-correlation
$f_b$ vs. $m_1$	correlation	no relation
$\gamma$ vs. [Fe/H]	anti-correlation	marginal correlation
$\gamma$ vs. $m_1$	anti-correlation	marginal correlation
$m_1$	larger $m_1$	smaller $m_1$
[Fe/H]	lower [Fe/H]	higher [Fe/H]

Table 3. Table 3: The MS star catalogue with probability of binary and mass-ratio. The full catalogue is at https://github.com/liuchaonaoc/Binarity/blob/master/table3.txt .

LM obsid	Gaia source id	ra	dec	$M_G$	$M_G$ error	$A_G$	$\Delta M_G$	$m_1$	$m_1$ error	$\log P_s/P_b$	$q$	$q$ 15%	$q$ 85%
		h:m:s	d:m:s	mag	mag	mag	mag	$M_{\odot }$	$M_{\odot }$
601119	1893391222402009728	22:11:51.26	+28:37:44.8	7.217	0.083	0.313	0.067	0.57	0.03	1.96	0.46	0.29	0.56
603081	1896846231534740352	22:04:18.06	+30:22:39.0	5.866	0.082	0.261	-0.962	0.65	0.04	-26.62	1.0	0.98	1.0
603241	1896831594286160384	22:04:38.07	+30:12:40.5	6.282	0.083	0.251	0.072	0.71	0.04	1.4	0.51	0.29	0.61
604003	1895407486209553664	22:10:48.25	+30:45:29.4	6.692	0.081	0.236	-0.115	0.68	0.03	1.04	0.72	0.54	0.76

Equations

$$p(M_G^{\prime }|\varpi ,G)\propto \exp \left(-\frac{(10^{-(G-M_G^{\prime }+5)/5}-\varpi {)}^2}{2{\sigma }_{\varpi }^2}\right),$$

$$p(M_G|[Fe/H],(G-K{)}_0)=\sum _{i=1}^{n}K(\frac{M_G-M_{G,i}}{{\sigma }_{M_{G,i}}}),$$

$$p(q|\gamma )\propto q^{\gamma },$$

$$p(f_b,\gamma |$$

$$\int \sum _{c}p(\Delta M_i|c,q,Z_i,m_{1,i})p(c|f_b,Z_i,m_{1,i})\times$$

$$p(q|\gamma ,Z_i,m_{1,i})dq,$$

$$\int \sum _c$$

$$(1-f_b)p_s(\Delta M_i)\int p(q|\gamma ,Z_i,m_{1,i})dq+$$

$$f_b\int p_b(\Delta M_i|q,Z_i,m_{1,i})p(q|\gamma ,Z_i,m_{1,i})dq,$$

$$p_s(\Delta M_i)=\frac{1}{\sqrt{2\pi }{\sigma }_i}\exp \left(-\frac{\Delta M_i^2}{2{\sigma }_i^2}\right)$$

$$p_b$$

$$\exp \left(-\frac{(\Delta M_i-\Delta M_{model}(q,Z_i,m_{1,i}){)}^2}{2{\sigma }_i^2}\right),$$

$$p(f_b,\gamma |\{\Delta M_i\},Z,m_1)=\prod _{i}p(f_b,\gamma |\Delta M_i,Z,m_1).$$

$$E_b\propto \frac{m_1^{2}q}{a},$$

$$P_{s,i}=(1-f_b)p_s(\Delta M_i)\int p(q|\gamma )dq.$$

$$P_{b,i}=f_b\int p_b(\Delta M|q)p(q|\gamma )dq.$$

$$p(q|\Delta M)=p(\Delta M|q)p(q|\gamma ).$$

$$\ln \mathcal{L}\propto$$

$$-\frac{(G_{obs}-G_{model}(M_{ini},age){)}^2}{2{\sigma }_G^2}$$

$$-\frac{(Z_{obs}-Z_{model}(M_{ini},age){)}^2}{2{\sigma }_Z^2},$$

$$(G_{BP}-G_{RP}{)}_0=$$

$$+4.23(\pm 0.55)\times 10^{-7}{T}_{\mathrm{eff}}^2$$

$$-3.01(\pm 0.29)\times 10^{-3}{T}_{\mathrm{eff}}$$

$$+8.25(\pm 0.50).$$

$$k_X=$$

$$+c_5{A}_0+c_6{A}_0^2+c_7(G_{BP}-G_{RP}{)}_0{A}_0,$$

Full text

††pubyear: 2019††pagerange: Smoking Gun of the Dynamical Processing of the Solar-type Field Binary Stars–B

1 Introduction

Binary stars are very common in stellar systems with different scales, from star cluster to the whole galaxy. Some literatures even thought that the most of stars live in binary systems (Heintz, 1969; Abt & Levy, 1976; Duquennoy & Mayor, 1991). The statistical analysis of binary systems play important role to constrain star formation in different environments. The statistical properties that can well describe binary systems include the fraction of binary, the distribution of mass-ratio, the distribution of orbital period, and the distribution of the orbital eccentricity. These properties are usually compared between populations with different age, metallicity, and primary mass (Duquennoy & Mayor, 1991; Henry & McCarthy, 1990; Fischer & Marcy, 1992; Mason et al., 1998a, b).

The association of binary fraction with metallicity has been broadly studied by lots of observational (Carney, 1983; Latham et al., 2002; Raghavan et al., 2010; Rastegaev, 2010; Moe & Di Stefano, 2013; Gao et al., 2014; Hettinger et al., 2015; Gao et al., 2017; Tian et al., 2018; Badenes et al., 2018; Moe, Kratter & Badenes, 2019, etc.) and theoretical works (Machida et al., 2009; Kratter et al., 2010; Myers et al., 2011; Bate, 2014; Tanaka & Omukai, 2014, etc.). Recently, lots of works, especially depending on large volume of datasets, found that binary fraction, at least for close binaries, is anti-correlated with metallicity (e.g. Raghavan et al., 2010; Gao et al., 2014, 2017; Tian et al., 2018; Badenes et al., 2018; Moe et al., 2019), while a few others claimed that $f_{b}$ is positively correlated (e.g. Carney, 1983; Abt & Willmarth, 1987; Hettinger et al., 2015) or not correlated with metallicity (Latham et al., 2002; Moe & Di Stefano, 2013).

Among these works, Moe et al. (2019) reanalyzed five different datasets of close binaries with sample size from several hundreds to several ten thousands from literatures and concluded that the binary fraction is strongly anti-correlated with metallicity. The datasets used in their studies contain close binaries with separation smaller than $10$ AU. As complementary, El-Badry & Rix (2019) discussed about the variation of the binary fraction as a function of metallicity for wide binaries with separation of $50$–$50\,000$ AU. They found the anti-correlation between the binary fraction and metallicity exists only when the separation of the two companions is between $50$ and $100$ AU. It becomes flat for binaries with larger separations as $100<a<200$ AU.

The binary fraction of main sequence (MS) stars is also found positively correlated with primary mass in many works (see the review by Duchêne & Kraus, 2013, and references therein). For very low-mass field binaries (m_{1}$$<0.1 $M_{\odot}$), the binary fraction is only about $22$% (Allen, 2007). For low-mass binaries ($0.1<m_{1}<0.5$ $M_{\odot}$), the fraction increases to around $26$% (Delfosse et al., 2004; Dieterich et al., 2012). For solar-type stars ($0.5<m_{1}<1.4$ $M_{\odot}$, which is also defined as the MS stars with luminosity from $0.1$ to $10$ $L_{\odot}$ by Raghavan et al. (2010)), the fraction further increases to around $41$–$50$% (Raghavan et al., 2010). And for the intermediate-mass ($1.5<m_{1}<5$ $M_{\odot}$) and massive stars ($m_{1}>16$ $M_{\odot}$), the fraction is larger than $50$% (Chini et al., 2012; Mason et al., 2009; Sana et al., 2013).

For the well studied solar-type binaries, the arguments about mass-ratio distribution are contradicted.  Duquennoy & Mayor (1991) derived a mass-ratio distribution biased to lower mass-ratio for solar-type stars, which makes it very similar to the canonical initial mass function (IMF) (Miller & Scalo, 1979; Kroupa et al., 1990). However, later studies claimed different shape. Raghavan et al. (2010) found that the mass-ratio distribution is rather flat with a peak at $q>0.9$, which is not like any of the canonical IMF (Bastian, Covey & Meyer, 2010). Goldberg, Mazeh, & Latham (2003), on the other hand, showed a double peaked mass-ratio distribution. Duchêne & Kraus (2013) compiled many references and found that if mass-ratio distribution is described with a power law, then the power index decreases with primary mass. Again, they demonstrated that the mass-ratio distribution, which is equivalent with the present-day mass function of the secondary, is significantly biased from any of the canonical IMF.

A large fraction of previous studies on the binary fraction and mass-ratio distribution were based on small dataset with number of stars from a dozen to a few hundreds, which suffers from substantially large Poisson noise. In recent years, some of the works (Gao et al., 2014, 2017; Hettinger et al., 2015; Badenes et al., 2018; Moe et al., 2019; El-Badry et al., 2018b; El-Badry & Rix, 2019) have made use of large survey data.

In light of these studies, we make the measurement of the binary fraction and mass-ratio distribution for stars with different metallicity and primary mass using the LAMOST data (Zhao et al., 2012) combined with Gaia DR2 (Gaia Collaboration et al., 2016, 2018). LAMOST DR5 have collected more than 8 million stellar spectra with spectral resolution of $R=1800$ and limiting magnitude of $r=18$ mag (Deng et al., 2012). It provides stellar effective temperature, surface gravity, and [Fe/H], with precision of about 110 K, 0.2 dex, and 0.15 dex (Luo et al., 2015; Gao et al., 2015). These parameters enable us to derive the stellar mass for solar-type stars. On the other hand, Gaia DR2 provides very accurate photometries in $G$, $BP$, and $RP$ bands as well as astrometry (Gaia Collaboration et al., 2018; Luri et al., 2018; Lindegren et al., 2018). With accurate luminosity derived from parallax, field stars can be mapped in the Hertzprung-Russell (HR) diagram with clear binary sequence located around 0.75 mag above the single main-sequence (see Gaia Collaboration et al., 2018). Hence, one can study the statistical properties of binaries directly in HR diagram, similar to the binary studies based on the photometries of star clusters or OB associations (e.g. Kouwenhoven et al., 2007; Milone et al., 2012; Li, de Grijs & Deng, 2013; Yang et al., 2018).

The paper is organized as follows. Section 2 describe how we select the solar-type stars and how to estimate their parameters. Section 3 develops the hierarchical Bayesian model which can simultaneously derive binary fraction and mass-ratio distribution. Section 4 shows the results for stellar populations with various metallicity and primary mass. From the Bayesian model, we also provide probabilities to be binary of the sample stars with the mass-ratio estimates. Section 5 raises discussions about the caveats of the method and comparisons with other works. Finally, brief conclusions are drawn in section 6.

2 Data

2.1 Data selection

We select common stars cross-identified from LAMOST DR5 (http://dr5.lamost.org), Gaia DR2 (Gaia Collaboration et al., 2018), and 2MASS (Skrutskie et al., 2006). The cross-identification is conducted by comparing the right ascensions and declinations in equinox of J2000.0 among the three catalogs. Because the diameter of the LAMOST fiber is $3.3$ arcsec and the mean seeing during LAMOST observations is around $3$  arcsec (Shi, J-R. private communication), the spatial resolution of LAMOST should be around $5$ arcsec, which is used as the radius of circle in the cross-identification. Only the star pairs with nearest positions within the circle are selected.

First, we derive absolute magnitude in $G$-band of the stars from Gaia parallax based on the likelihood distribution. which can be written as

[TABLE]

where $M_{G}^{\prime}$ is the Gaia $G$-band absolute magnitude without extinction correction, $\varpi$ and $\sigma_{\varpi}$ are the parallax and its uncertainty, respectively. $M_{G}^{\prime}$ equals to $M_{G}+A_{G}$, where $M_{G}$ and $A_{G}$ are the Gaia $G$-band absolute magnitude with extinction correction and the interstellar extinction in $G$-band, respectively. We estimate the maximum likely $M_{G}^{\prime}$ and is uncertainty for each star by arbitrarily drawing 10 000 points from the distribution defined by Eq. (1).

We then select samples using the following criteria:

1. 1

$M_{G}^{\prime}>4$ mag; 2. 2

\varpi$$>3 so that most of the stars are located within $\sim 333$ pc; 3. 3

\varpi$$/$$\sigma_{\varpi}$$>5; 4. 4

phot_bp_mean_flux$/$phot_bp_mean_flux_error$>20$, phot_rp_mean_flux$/$phot_rp_mean_flux_error$>20$, and

phot_g_mean_flux$/$phot_g_mean_flux_error$>20$; 5. 5

$\sqrt{\chi^{2}/(\nu-5)}<1.2\times\max(1,\exp(-0.2(G-19.5)))$, where $\chi^{2}$ and $\nu$ represent for astrometric_chi2_al and astrometric_n_good_obs_al, respectively; 6. 6

K_{s}$$<14.3 mag and \sigma_{K_{s}}$$<0.2 mag; 7. 7

$T_{\rm eff}$ derived from LAMOST is in between $3800$ and $<6500$ K; 8. 8

signal-to-noise ratio at $g$ band of LAMOST spectra is larger than 20.

The first criterion ensures that the hot MS stars are removed.The second criterion selects the stars located within $\sim 333$ pc so that the data comply with the requirement of the completeness. The faintest apparent magnitude for a star with mass of $0.08$ $M_{\odot}$ at distance of $333$ pc is between $14.5$ and $16.1$ mag depending on different metallicity. This is bright enough for the companion in a binary to contribute luminosity that is able to be detected by LAMOST. Therefore, the volume cut enables the detection of the faintest secondaries. The criteria 3–6 is set to select stars with accurate photometry and astrometry. Among them, criterion 5 is adopted from El-Badry & Rix (2018). Criteria 7 and 8 help to select high quality spectra from LAMOST so that the stellar parameters derived from the spectra are reliable. Meanwhile, the two criteria also comply with the criterion 2. After applying these selection criteria, we obtain 101 061 MS stars in total.

2.2 Stellar mass

We adopt the present-day mass as the initial mass for the primary stars. First, because in this work we only focus on the main-sequence stars with stellar mass of 0.4-1 $M_{\odot}$, the mass loss due to stellar wind is very small and ignorable. Second, the fraction of mass transferring binary system is less than one percent (Rucinski, 1994; Chen et al., 2016) for the solar-mass main-sequence stars, hence the change of the stellar mass due to mass transfer is also ignored.

The initial stellar mass of each star is determined by comparing their $T_{\rm eff}$, $\log g$, and [Fe/H] with PARSEC isochrones (Bressan et al., 2012) by means of likelihood (see the details in Appendix A).

Although the stellar parameters obtained from the LAMOST spectra are dominated by the primary stars in a binary system, the secondary stars may averagely reduce the effective temperature by $\sim 50$ K when $q>0.6$ and induce additional 0.1 dex and 0.05 dex dispersions in $\log g$ and [Fe/H], respectively, for the case of LAMOST spectra (El-Badry et al., 2018a). Considering these systematics and uncertainties, the total uncertainty of the likelihood approach is about $0.05$ $M_{\odot}$.

2.3 Interstellar extinction correction

The determination of the interstellar extinction of the stars is described in Appendix B. The uncertainty of the extinction in $G$ band is around $0.05$ mag.

We further select the stars with $1.5<(G-K_{s})_{0}<2.6$, where $(G-K_{s})_{0}$ is the dereddened color index, so that the data can cover stellar mass from $0.4$ to $1.0$ $M_{\odot}$. A few stars without initial mass estimates due to their unusual photometries are excluded. After this cut, the samples are reduced to 76 507 stars.

2.4 Spatial distribution

We adopt the Bayesian distance derived from the parallax of Gaia by Bailer-Jones et al. (2018). The top panels of Figure 1 demonstrate the distributions of the stars in $M_{G}$ and mass vs. distance planes. It shows that $M_{G}$ covering from 4 to 8 mag does not show correlation with distance from 50 to 333 pc. Meanwhile, stellar mass in the range from 0.4 to 1 $M_{\odot}$ is also independent on distance. The bottom-left panel of the figure shows that [Fe/H] is uncorrelated with distance. These mean that $M_{G}$, stellar mass, and [Fe/H] of the samples are unbiased within the volume from 50 to 333 pc.

Although LAMOST survey only focuses on northern sky, it basically continuously cover the sky of about 20 000 square degrees, from low to high Galactic latitude, as seen in the bottom-right panel of Figure 1. Therefore, the samples selected in this work is representative to the stellar populations in the solar neighborhood.

2.5 Color-absolute magnitude diagrams of the stars of different [Fe/H]

Because different [Fe/H] may shift the main-sequence up and down in HR diagram, they can blur out the binary sequence. Therefore, we separate the stars in different [Fe/H] bins. The uncertainty of the LAMOST derived [Fe/H] is around 0.15 dex (Gao et al., 2015). Thus, we set the bin size of [Fe/H] lager than 0.2 dex. On the other hand, we need to keep sufficient number of stars in each [Fe/H] bin for further separation in primary mass in section 4. Therefore, we adopt 8 [Fe/H] bins with ranges of $-1.25<$[Fe/H]$<-0.55$, $-0.65<$[Fe/H]$<-0.40$, $-0.50<$[Fe/H]$<-0.25$, $-0.35<$[Fe/H]$<-0.15$, $-0.25<$[Fe/H]$<-0.05$, $-0.15<$[Fe/H]$<+0.05$, $-0.05<$[Fe/H]$<+0.15$, and $+0.05<$[Fe/H]$<+0.40$. Note that the neighboring bins are partly overlapped with each other such that the number of stars in each bin is as large as possible. The overlapping would not harm the final results, but only smooth the results to some extent. The density distributions of the sample stars in $M_{G}$ vs. $(G-K_{s})_{0}$ plane with different [Fe/H] bins is displayed in Figure 2. The density map is calculated in the grid of $(G-K_{s})_{0}$ and $M_{G}$ with size of $0.05\times 0.025$ mag. The colors in the figure codes the logarithmic density of stars. It is clearly seen two sequences in the color-absolute magnitude diagrams such that the unresolved binary sequence is located at upper side of the single star sequence due to their larger luminosity contributed by both companions.

The shift of the absolute magnitude of an unresolved binary star from the single sequence reflects the mass-ratio of the binary system, which is denoted as q$$=m_{2}/m_{1}, where $m_{1}$ and $m_{2}$ are the masses of the primary and secondary stars, respectively. When $q$ is small ($\lesssim 0.7$), the contribution to the luminosity from the secondary star is only around 0.1 mag. Therefore, the composed luminosity for such binaries is similar to the single stars. When $q$ becomes larger ($\gtrsim 0.7$), the contribution from the secondary becomes obvious, hence the luminosity of the whole system increases substantially. When $q$ is close to 1, the composed luminosity is doubled, which is equivalent to 0.75 mag brighter.

To quantify the difference between the single and binary stars in $M_{G}$, we firstly address the ridge lines of the single MS stars in a finer [Fe/H] grid with bin size of $0.05$ dex. For each [Fe/H] slice, we further separate the stars into $(G-K_{s})_{0}$ bins with size of $0.05$ mag. At each color index bin, we apply a kernel density estimation (KDE) to produce the density distribution of $M_{G}$. Specifically, the density of $M_{G}$ is

[TABLE]

where $K(t)=\exp(-t^{2}/2)$ is the kernel function; $M_{G,i}$ and $\sigma_{M_{G,i}}$ are the absolute magnitude and its uncertainty, respectively, of the $i$th star. We then locate the intrinsic $M_{G}$ of single stars at the peak of $p(M_{G}|[Fe/H],(G-K)_{0})$. Figure 3 shows the color-coded intrinsic $G$-band absolute magnitude, denoted as $M_{G,single}([Fe/H],(G-K)_{0})$, of the single star in [Fe/H] vs. $(G-K_{s})_{0}$ plane. Note that there five points marked as black frames in the figure are failed to derive $M_{G,single}([Fe/H],(G-K)_{0})$ due to the missing data. We then provide the linear interpolated values for the five points based on their neighboring points.

For a given star with fixed [Fe/H] and $(G-K_{s})_{0}$, we then calculate the differential absolute magnitude $\Delta M_{G}=M_{G}-M_{G,single}([Fe/H],(G-K)_{0})$, in which $M_{G,single}$ is selected from Figure 3 with nearest [Fe/H] and $(G-K_{s})_{0}$. The uncertainty of $\Delta M_{G}$ is contributed from several ways. First, it is mostly from the astrometric and photometric uncertainties. Then, the de-reddening process also contributes an additional uncertainty. Finally, $M_{G,single}$ also brings some uncertainty. Combining all these together, we obtain the typical uncertainty of $\Delta M_{G}$ as $0.08$ mag.

Figure 4 shows the distributions of the stars in $\Delta M_{G}$ vs. $m_{1}$ plane for different [Fe/H] populations. It is seen that a large fraction of the stars are aligned with \Delta M_{G}$$=0, which is the location of the single stars. The vertical dispersion of these stars is roughly consistent with the uncertainty of $\Delta M_{G}$. Another notable groups of stars are located at around \Delta M_{G}$$\sim-0.75 mag, which are mostly the binary stars with mass-ratio close to 1. For metal-poor populations ([Fe/H]$<-0.55$), the coverage of $m_{1}$ is from around 0.45 to 0.7 $M_{\odot}$, while for metal-rich ones ([Fe/H]$>-0.15$), the range of $m_{1}$ is from 0.6 to 1.0 $M_{\odot}$.

Note that, for single stars, $m_{1}$ stands for the stellar mass of the single star, while for binaries, it represents for the primary stellar mass. For convenience, we do not distinguish them in notation but simply use $m_{1}$ for both cases.

It is also seen that for stars with super-solar metallicity, as shown in the two bottom-right panels of Figure 4, the dispersion of $\Delta M_{G}$ is larger when m_{1}$$>0.9 $M_{\odot}$. It implies that there are more stars located below the single main-sequence in the HR diagram with mass larger than 0.9 $M_{\odot}$. This larger dispersion, which could be associated with the so called blue sequence discussed by Yang et al. (2018), induces larger uncertainties in binary fraction and mass-ratio distribution estimates. Therefore, we cut off the stars with mass larger than $0.85$ $M_{\odot}$ to avoid it.

3 Hierarchical Bayesian model

A broadly used approach to calculate the fraction of binary stars especially in star clusters is to define an empirical threshold in color-absolute magnitude diagram to separate single and binary stars. Then, a following Monte Carlo simulations, which adopt some presumed model of binaries, e.g. binary fraction and mass-ratio distribution, are conducted to compare with the observations. Subsequently, the intrinsic binary fraction and mass-ratio distribution can be obtained from the simulations which converge their results to the observations (e.g. Li et al., 2013; Sana et al., 2013).

In this work, instead to follow the previous method, we propose a hierarchical Bayesian model to simultaneously solve out the fraction of binaries and the distribution of mass-ratio of the binary systems. In Bayesian model, we concentrate to the posterior probability distribution of the binary fraction and parameters determining the shape of the mass-ratio distribution, rather than derive the best-fit binary fraction and mass-ratio distribution.

Now assume that the distribution of mass-ratio in MS binary system follows a power law, that is

[TABLE]

where q$$=m_{2}/m_{1} represents for the mass ratio of a binary and $\gamma$ is the power index. It is noted that the power-law assumption is not a precise description of the distribution of mass-ratio (Duchêne & Kraus, 2013). However, it is simple and can essentially provide the information about the relative number of binaries with low and high mass-ratio.

For each MS star, the observed absolute magnitude is not only determined by its stellar mass, metallicity, and age, but also affected by two other factors: 1) the luminosity of the unresolved companion if it is a binary and 2) the random error propagated from astrometry and photometry. To quantify the two effects, we firstly need a prior knowledge about whether a star is an unresolved binary. The probability of a star to be a binary is equivalent with the fraction of binary ($f_{b}$). If the star of interest is a binary, we also need to know in a prior about the mass-ratio of the binary so that we can derive the additional luminosity contributed by the unresolved companion. The prior knowledge about the mass-ratio is provided by the power law of the mass-ratio distribution with the power index, $\gamma$, in Eq. (3). Therefore, there are two unknown parameters in the prior distributions, $f_{b}$ and $\gamma$. Then, the posterior distribution of $f_{b}$ and $\gamma$ given the observed differential absolute magnitude $\Delta M_{G_{i}}$ for the $i$th star can be written as

[TABLE]

where $Z_{i}$ and $m_{1,i}$ are metallicity and primary mass for the $i$th star, respectively. $\Delta M_{i}$ stands for the differential absolute magnitude. $c$ represents for the class of the star, either a single or a binary. Obviously, $p(c|f_{b})$ follows Bernoulli distribution and thus the term within the integral in the right-hand side of Eq. (3) can be rewritten as

[TABLE]

where $p_{s}$ and $p_{b}$ stand for the likelihood distribution of $\Delta M_{i}$ for single and binary star, respectively. We adopt that both of them follow normal distribution and then we have

[TABLE]

for single star and for binary it becomes

[TABLE]

where $\Delta M_{model}(q,m_{1,i},Z_{i})$ is the model differential absolute magnitude considering that the star has a companion with mass-ratio of $q$ and $\sigma_{i}$ is the uncertainty of $\Delta M_{i}$.

The posterior distribution of $f_{b}$ and $\gamma$ for a group of stars with same (or very similar) [Fe/H] and $m_{1}$ can be obtained by multiplying all the single star posterior distributions, such as

[TABLE]

4 Results

We separate the stars into different metallicity and stellar mass groups, as shown in Figure 5. The stars are grouped into [Fe/H] and $m_{1}$ bins with various sizes so that we can find balance between the number of stars in each bin and the resolution of metallicity and stellar mass. Along $m_{1}$, the center positions of the bins are at 0.525, 0.6, 0.65, 0.7, 0.75, and 0.8 $M_{\odot}$. The width of the bins are 0.15 for the first and 0.1 $M_{\odot}$ for the other bins. The bin size of $m_{1}$ are a factor of 2 larger than the uncertainty of $m_{1}$ for individual stars. Along [Fe/H], the centers of the bins are at -0.9, -0.525, -0.375, -0.15, -0.05, 0.1, and 0.225 dex, with bin widths of 0.7, 0.25, 0.25, 0.2, 0.2, 0.2, 0.2, and 0.35 dex, respectively. The bin size is also larger than the uncertainty of [Fe/H] for individual stars so that the errors of [Fe/H] would not significantly affect the results. Note that these bins are partly overlapped with their adjacent bins such that more stars are fell in the bins and therefore can stabilize the results by smoothing the results.

Table 1 and Figure 5 show the numbers of stars fell into these bins. According to the Monte Carlo simulation conducted in section 5.2, we find that, in principle, the results are stable when the number of stars is larger than 2000 stars in a bin. Therefore, in the rest of the analysis we only consider the bins containing larger than 2000 stars, which are filled up with colors in Figure 5 . This cut will lose the stars with [Fe/H]$<-0.6$.

Then, we apply the hierarchical Bayesian model described in section 3 to the [Fe/H]–$m_{1}$ bins with more than 2000 stars. We run Markov chain Monte Carlo simulations using EMCEE software package (Foreman-Mackey et al., 2013). We adopt the median values and 15% (85%) percentiles of the random draws in MCMCs as the best fit values and the uncertainties, respectively, for $f_{b}$ and $\gamma$.

In the hierarchical Bayesian model, we set the lower limit of the mass of the secondary at 0.08 $M_{\odot}$. We do not take into account the brown dwarfs in this work.

The resulting $f_{b}$ and $\gamma$ for different [Fe/H] and $m_{1}$ bins are shown in Figure 6 and Table 1.

4.1 Overview of the results

The left panel of Figure 6 shows that the binary fraction is roughly larger at lower [Fe/H]. Meanwhile, $f_{b}$ is essentially larger when $m_{1}$ is larger at most [Fe/H] values except $-0.5$ dex. A ridge-like valley with f_{b}$$<0.22 is seen from [Fe/H]$\sim-0.25$ dex and m_{1}$$\sim 0.6 $M_{\odot}$ to [Fe/H]$\sim+0.25$ dex and m_{1}$$\sim 0.75 $M_{\odot}$. Within the ridge, $f_{b}$ increases with declining [Fe/H]. In the other regions outside this ridge, the values of $f_{b}$ are larger than $0.25$. Most of them are larger than $0.40$.

The right panel of Figure 6 shows the resulting $\gamma$ as a function of [Fe/H] and $m_{1}$. Similar to the distribution of $f_{b}$, the location of the bins with $\gamma>1.2$ also concentrate in the same ridge as the low $f_{b}$ bins, i.e. from [Fe/H]$\sim-0.25$ dex and m_{1}$$\sim 0.6 $M_{\odot}$ to [Fe/H]$\sim+0.25$ dex and m_{1}$$\sim 0.75 $M_{\odot}$. Furthermore, $\gamma$ roughly increases from $1.77$–$2.74$ at bottom-left to $3.18$–$3.39$ at top-right in the ridge. On the other hand, the values of $\gamma$ in other bins are around zero. These bins shows opposite trend to those in the ridge such that their values increase when [Fe/H] decreases.

Then, we draw the relationship between $f_{b}$ and $\gamma$ in Figure 7. Surprisingly, we find that they are strongly anti-correlated with each other. This implies that the less the binary fraction, the more binaries with high mass-ratio and vice versa. To our knowledge, it is for the first time that the anti-correlation is unveiled in the field stars.

The right-side panel of Figure 7 displays the distribution of $\gamma$ using KDE smoothing technique. It shows a gap at around $\gamma\sim 1.2$, separating the values of $\gamma$ into two groups: the low-$\gamma$ group with $\gamma<1.2$, which contains larger binary fractions, and the high-$\gamma$ group with $\gamma>1.2$, which contains lower binary fractions. The high-$\gamma$ group is marked with black frames in Figure 6, which is exactly overlapped with the ridge regions with low $f_{b}$ and high $\gamma$.

The gap is a real feature rather than an artifact. The bin sizes are quite large than the measurement errors of [Fe/H] and $m_{1}$ and each bin contains sufficient number of samples. Therefore, the gap would not be due to the random spike dominated by noise. Moreover, the bins have been smoothed by partly overlapping with the adjacent bins. This makes the gap more robust. However, the nature of the gap is not clear.

Since the $f_{b}$–$\gamma$ pairs are naturally separated into two groups, in next sections, we separately discuss the features of the two groups.

4.2 Low-$\gamma$ group

First, although the bins of the low-$\gamma$ group cover a wide range of [Fe/H] and $m_{1}$, as shown in Figure 6, only with [Fe/H]$<-0.3$ dex does $m_{1}$ continuously cover from 0.5 to 0.75 $M_{\odot}$. In the regime with [Fe/H]$>-0.3$ dex, most of the bins of the low-$\gamma$ group have $m_{1}$ larger than $0.7$ $M_{\odot}$. Two bins at ([Fe/H], $m_{1}$)=($+0.1$ dex, $0.6$ $M_{\odot}$) and ($+0.2$ dex, $0.65$ $M_{\odot}$) are isolated from other low-$\gamma$ bins.

Then, we look at how $f_{b}$ is associated with [Fe/H] and $m_{1}$ for the low-$\gamma$ group in the top panels of Figure 8. In all four panels of the figure, the resulting $f_{b}$ and $\gamma$ of low-$\gamma$ group are drawn in blue. We fit the relationship of $f_{b}$ ($\gamma$) with [Fe/H] and $m_{1}$ respectively using a linear model. We find that $f_{b}$ is not related to [Fe/H], because the slope $\partial f_{b}/\partial{\rm[Fe/H]}=0.08\pm 0.14$ dex*-1* (see the blue line in the top-left panel) is roughly zero. However, $f_{b}$ is clearly correlated with $m_{1}$ with substantially positive slope of $\partial f_{b}/\partial m_{1}=1.15\pm 0.26$ $M_{\odot}$-1 (see the blue line in the top-right panel).

Finally, as shown in the bottom panels of the figure, $\gamma$ is tightly anti-correlated with both [Fe/H] and $m_{1}$ with negative slopes of $\partial\gamma/\partial{\rm[Fe/H]}=-1.16\pm 0.29$ dex*-1* and $\partial\gamma/\partial m_{1}=-3.04\pm 0.77$ $M_{\odot}$-1, respectively.

4.3 High-$\gamma$ group

As mentioned in Section 4.1, the high-$\gamma$ group bins concentrate in the ridge-like region in the $m_{1}$ vs. [Fe/H] plane (see the black frames marked in all panels of Figure 6). This means that the high-$\gamma$ group only contains data with [Fe/H]$>-0.3$ dex and $m_{1}$ smaller than 0.75 $M_{\odot}$. It is also noted that at the most metal-rich regions, the $m_{1}$values of the high-$\gamma$ bins show lower limits, which is 0.6 $M_{\odot}$ at [Fe/H]$=+0.1$ dex and 0.65 $M_{\odot}$at [Fe/H]$=+0.25$ dex.

Then, we investigate the correlation of $f_{b}$ with [Fe/H] and $m_{1}$ in the top panels of Figure 8. The red stars indicate the $f_{b}$ values of the high-$\gamma$ group and the red lines display the best linear fits of the relationship of $f_{b}$ with [Fe/H] and $m_{1}$. It shows that $f_{b}$ is tightly anti-correlated with [Fe/H] based on the substantially negative slope of $\partial f_{b}/\partial{\rm[Fe/H]}=-0.21\pm 0.04$ dex*-1*. In the mean time, $f_{b}$ is not related with $m_{1}$, since the slope $\partial f_{b}/\partial m_{1}=-0.18\pm 0.19$ $M_{\odot}$-1 is not significantly different from zero.

Finally, we look at the correlation of $\gamma$ with [Fe/H] and $m_{1}$ for the high-$\gamma$ group in the bottom panels of Figure 8. We find that $\gamma$ is weakly correlated with both [Fe/H] and $m_{1}$. Although the slopes $\partial\gamma/\partial{\rm[Fe/H]}=2.05\pm 1.20$ dex*-1* and $\partial\gamma/\partial m_{1}=5.05\pm 3.35$ $M_{\odot}$-1 are larger than zero, the significances are less than $2$-$\sigma$. The red stars located in the bottom panels also show larger dispersion than the low-$\gamma$ group (blue rectangles). Therefore, we infer that the correlations of $\gamma$ with either [Fe/H] or $m_{1}$ is quite marginal.

To summarize these complicated relationships, we list them in Table 2.

4.4 Smoking gun of the dynamical processing

In this section, we explain the results described in above sections based on the current theories of binary formation and evolution. We firstly give four assumptions.

First, we assume that the mass function of the secondary is universal, i.e. no matter [Fe/H] and $m_{1}$, the secondary mass function are all the same. It is not necessarily to be the canonical IMF given by Chabrier (2005) or Kroupa & Bouvier (2003) or others. Note that although Myers et al. (2011) has discussed theoretically about the invariability of the IMF with metallicity, it is not clear whether the IMF of secondary is same as that of the primary.

Second, according to Moe et al. (2019), the close binaries are more likely formed from disk fragmentation, which is associated with the metallicity of the gas (also c.f. Kratter et al., 2008, 2010; Machida et al., 2009; Tanaka & Omukai, 2014). However, for the wide binaries, as found by El-Badry & Rix (2019), they are not related to metallicity, which is consistent with the scenario that the wide binaries are formed from gas fragmentation (Tohline, 2002; Offner et al., 2010; Myers et al., 2011; Bate, 2012, 2014). Consequently, we assume that the binaries are formed both from disk fragmentation, which contributes mostly to the close binaries depending on metallicity, and gas fragmentation, which produces most of the wide binaries not depending on metallicity.

Third, we assume that the binary fraction increases with primary stellar mass based on the simulation by Tanaka & Omukai (2014). This assumption is also supported by many observations (Delfosse et al., 2004; Raghavan et al., 2010; Sana et al., 2012, 2013).

Finally, it is known that most stars are formed in embedded clusters (Lada & Lada, 2003). Marks & Kroupa (2011) argued that the dynamical processing in the embedded clusters, in which binaries were born and spent their early times, prefers to efficiently disrupting binaries with larger separations, smaller primary mass, and smaller mass-ratios. This is because the binding energy of a binary system can be expressed as

[TABLE]

where $a$ is the separation of the two companions (Marks, Kroupa & Oh, 2011). The dynamical processing occurred in embedded clusters has a very short time scale of a few million years (Belloni et al., 2018). We assume that the dynamical processing is generally very efficient and can change the present-day field binary properties to different extents depending on the primary mass, separation, mass-ratio etc.

Now, we can apply all above assumptions to the two groups of stars and explain why they show the observed $f_{b}$ and mass-ratio distribution features as listed in Table 2.

Firstly, because the high-$\gamma$ group contains smaller $m_{1}$, they are easier to be affected by the dynamical evolution. The dynamical processing destroys many binaries with lower $q$ and wider separation since they have smaller binding energy. Consequently, this leads to smaller binary fraction for high-$\gamma$ populations. In the mean time, the loss of lower mass-ratio binaries also results in larger $\gamma$ since the existing binaries are more biased to larger $q$. Moreover, the loss of wide binaries tends to increase the fraction of close binaries in the present-day samples. According to the second assumption that the formation of close binaries is metallicity-dependent, the close binary dominated populations exhibit the tight anti-correlation between [Fe/H] and $f_{b}$.

Meanwhile, the high-$\gamma$ group is not only biased to lower $m_{1}$, but also biased to higher metallicity. This means that the metal-rich stars may have experienced more effective dynamical processing. It is not clear how metallicity plays the role in the formation of the embedded clusters and how it affect the dynamical processing. This observational evident may be helpful in clarifying how metallicity works in this issue.

On the other hand, because the low-$\gamma$ group contains larger primary masses, it is less affected by the dynamical processing. Therefore, the populations in this group have relatively larger present-day binary fractions. The flat $f_{b}$–[Fe/H] relationship is evident that the populations with smaller $\gamma$ contain more metallicity-independent wide binaries, which smear out the $f_{b}$–[Fe/H] anti-correlation contributed by close binaries. These support that the low-$\gamma$ group may be less or even not affected by the dynamical processing during the stage of the embedded clusters.

The anti-correlations between $\gamma$ and [Fe/H] and between $\gamma$ and $m_{1}$ in low-$\gamma$ group may reflect the various star (and also binary) formations influenced by various characteristics, e.g. cluster mass, IMF, metallicity, of the embedded clusters. Future theoretical modeling making use of these observational evidences will be very useful to investigate the nature of star formation.

We should emphasize that the explanations of high- and low-$\gamma$ groups are based on a number of assumptions compiled from various theoretical and observational works. Because the theoretical works are usually limited by their ingredients and techniques, they do not exactly reflect the realities. Therefore, these explanations are not exclusive. Although it seems that they are well consistent with the assumption and the observations, other mechanisms are not ruled out.

4.5 Identification of binaries

In each [Fe/H]–$m_{1}$ bin, we are able to identify which star is likely a binary, since the resulting $f_{b}$ and $\gamma$ are the hyper-parameters of the prior in Eq. (3). We can then use Eq. (3) to estimate the probabilities that a star to be either a single or a binary. Specifically, the probability of the $i$th star to be a single can be written as

[TABLE]

And the probability to be a binary reads

[TABLE]

Then, the odd of the star to be a single can be written as $P_{s}/P_{b}$. When $P_{s}/P_{b}$ is smaller than 1, the star is more likely to be a binary, while it is larger than 1, it is likely to be a single star. Figure 9 shows that $\log P_{s}/P_{b}$ is correlated with $\Delta M_{G}$ as expected. Given a star, the more likely to be a binary, the smaller the $\Delta M_{G}$.

When a star is identified as a binary, the mass-ratio can also be roughly approximated from $\Delta M_{G}$. We consider a posterior distribution for $q$ such that

[TABLE]

We adopt the best-fit $q$ to be at the maximum of the posterior distribution and the 15% and 85% percentiles of the distribution as the uncertainties. Figure 9 also indicates the mass-ratios based on Eq. (12).

Because only the filled regions in Figure 5 have resulting $f_{b}$ and $\gamma$, only the stars located in these regions can be classified as single and binaries. We finally calculate $\log P_{s}/P_{b}$ and $q$ for 57 978 stars with [Fe/H]$\gtrsim-0.6$ dex. Table 3 lists the whole 57 978 samples with $\log P_{s}/P_{b}$ and the best-fit mass-ratio. Note that $q$ estimated from Eq. (12) is only related to $M_{G}$ and $\gamma$, no matter whether a star is classified as binary or single, we provide $q$ value and its uncertainty for it. We suggest 5 410 of these samples with $\log P_{s}/P_{b}<0.1$ are very likely binaries.

5 Discussions

5.1 Comparison with other works

5.1.1 Binary fraction

We compare the binary fraction from our results with Moe et al. (2019) (black dashed line in the top-left panel of Figure 8). The slope of the anti-correlation between $f_{b}$ and [Fe/H] is $\partial f_{b}/\partial{\rm[Fe/H]}\sim-0.12$ dex*-1* from Moe et al. (2019), which is similar to $-0.21\pm 0.04$ dex*-1* from the high-$\gamma$ group. The binary fraction from Moe et al. (2019) is larger by about 0.1 than that of the high-$\gamma$ group. Compared with the binary fraction of low-$\gamma$ group, however, we find that the $f_{b}$ from Moe et al. (2019) is significantly smaller. Because Moe et al. (2019) studied the close binaries ($a<10$ AU), it is not surprising that their $f_{b}$ is more similar to high-$\gamma$ group, which is likely affected by dynamical evolution and hence more dominated by close binaries, than to low-$\gamma$ group, which contains more wide binaries.

El-Badry & Rix (2019) shows that the anti-correlation is substantial in binaries with separation smaller than 100 AU but disappears in wider binaries. This also agrees with our results for high-$\gamma$ and -low groups.

We compare the $f_{b}$–$m_{1}$ correlation with previous works. First, Delfosse et al. (2004) and Leinert et al. (1997) found that the multiplicity fraction of M dwarf stars ($0.1$–$0.5$ $M_{\odot}$) is $0.26$. The closest mass range in our work is m_{1}$$\sim 0.525 $M_{\odot}$, which has $f_{b}$ of $0.25\pm 0.02$. This is consistent with Delfosse et al. (2004). Second, Raghavan et al. (2010) selected stars in the mass range slightly larger than us. The multiplicity fraction of F6-G2 stars, which are slightly more massive than $\sim 0.8$ $M_{\odot}$, in Raghavan et al. (2010) is $0.50\pm 0.04$. Only low-$\gamma$ group contains m_{1}$$\sim 0.8 $M_{\odot}$and it falls between $0.40$ and $0.78$ depending on different metallicity. The fraction of their G2-K3 samples, which covers 0.7<$$m_{1}$$<1 $M_{\odot}$, is $0.41\pm 0.03$, while our results with m_{1}$$=0.7 $M_{\odot}$ at low-$\gamma$ group is from $0.40$ to $0.60$, which are in the similar ranges.

5.1.2 Mass-ratio Distribution

We then compare the resulting $\gamma$ with previous literatures. Duquennoy & Mayor (1991) argued that the mass-ratio distribution of the G type stars is quite similar to the canonical IMF (Miller & Scalo, 1979; Kroupa, Tout & Gilmore, 1990), while Raghavan et al. (2010) suggested a flat distribution of mass-ratio with a sharp peak at $q\sim 1$. However, we do not find any literature has discussed about the anti-correlation of the mass-ratio distribution with metallicity as displayed in the bottom-left panel of Figure 8 for the low-$\gamma$ group.

Duchêne & Kraus (2013) compiled many literatures and demonstrated that the power index of a power-law fitted mass-ratio distribution substantially changes with the primary mass (see their Figure 2). At m_{1}$$\lesssim 1 $M_{\odot}$, the smaller the primary mass, the larger the power index. This means that for low-mass binaries, the secondaries tend to have similar mass as the primaries, while the more massive binaries tend to have more low-mass companions. We fit the $\gamma$ values for binaries with all range of period and m_{1}$$<1 $M_{\odot}$, which are displayed as the red diamond symbols in the Figure 2 of Duchêne & Kraus (2013), using a power law and show it as black solid line in the bottom-right panel of Figure 8. In the range of $m_{1}$ between 0.4 and 0.85 $M_{\odot}$, the $\gamma$s derived from previous works prefer to a quite flat relationship with $m_{1}$. Compared to the low-$\gamma$ group, at around m_{1}$$\sim 0.8 $M_{\odot}$, our result is similar to Duchêne & Kraus (2013). However, the $\gamma$ of the low-$\gamma$ group has stronger anti-correlation with $m_{1}$ than Duchêne & Kraus (2013). Duchêne & Kraus (2013) also provided $\gamma$ for populations of binaries with smaller period (the blue squares in their figure 2). We simply connect the two points at m_{1}$$\sim 0.2 and $1.0$ $M_{\odot}$ and show it with a black dashed line in the bottom-right panel of Figure 8. Their $\gamma$ is essentially similar to the values of the high-$\gamma$ group. This is also consistent with the hypothesis that high-$\gamma$ group is dominated by shorter period binaries due to the effective dynamical processing.

5.2 Effect of sample size

In this section, we investigate the effect of the sample size in the hierarchical Bayesian model.

For this purpose, we generate mock samples with various numbers of mock stars and different $f_{b}$ and $\gamma$. For all simulated data, we adopt the primary stellar mass m_{1}$$=0.7 $M_{\odot}$ and [Fe/H]$=-0.3$ dex. We first select 6 $f_{b}$ values at $0.1$, $0.2$, …, $0.6$ and at 5 $\gamma$ values at $-1$, [math], …, $3$. Given any combination of $f_{b}$ and $\gamma$ from the above values, we arbitrarily draw three mock datasets with numbers of $N=500$, $2000$, and $10000$, respectively. Each mock star has a probability of $f_{b}$ to be a binary. If it is a binary, the mass-ratio $q$ is arbitrarily drawn from the power law with the given $\gamma$. In total, we have $6\times 5\times 3=150$ groups of mock data with different $f_{b}$, $\gamma$, and samples size.

The differential absolute magnitude $\Delta M_{G}$ for the mock data is calculated from PARSEC stellar model given the adopted $m_{1}$, [Fe/H], and $q$. We then run the hierarchical Bayesian model described in section 3 via MCMC for each group of the mock data.

Figure 10 shows the results of the mock data grouped with sample size. The first row shows the results for the mock data with $N=500$. The top-left panel displays that the estimated $f_{b}$ is essentially consistent with the true values (the presetting binary fraction when generating the mock data, $x$-axis). However, in the top-middle panel, a few estimated $\gamma$ are substantially biased, while the others well reproduce the “true” $\gamma$. The top-right panel shows the offsets of the estimated (black dots with error bars) parameters from the presetting parameters (gray circles). It is seen that when the true $f_{b}$ is smaller than $0.25$, the $\gamma$ tends to be overestimated. Meanwhile, the binary fractions are slightly overestimated with larger error bars to around $0.2$ when the presetting fraction is larger than $0.25$ and $\gamma<1$. The overall systematic offset of derived $f_{b}$ is $-0.01$ with uncertainty of $0.10$. The systematic bias of the estimated $\gamma$ is $0.46$ with uncertainty of $1.38$ for the sample size of $N=500$.

The middle row of Figure 10 shows the results for the mock data with $N=2000$. Because the sample sizes increases, the performance of the estimations of $f_{b}$ and $\gamma$ are substantially improved. The left and middle columns show that the estimates of $f_{b}$ and $\gamma$ well reproduce the presetting values. The right column shows no significant overestimation for $\gamma$ at f_{b}$$<0.2, while the slightly overestimation and larger uncertainty of $f_{b}$ still exists. The $\gamma$ is also slightly underestimated in most of the datasets at f_{b}$$>0.2 and $\gamma<0$. At f_{b}$$<0.2, the uncertainties of both $\gamma$ and $f_{b}$ are large. The systematic bias of the estimated $f_{b}$ is $0.02$ with uncertainty of $0.06$. The bias of the deriving $\gamma$ is $-0.21$ with dispersion of $0.47$. The random errors are significantly smaller than the results of the $N=500$ samples and the systematic bias in $\gamma$ is reduced by a factor of 2.

The bottom row shows the results for samples with $N=10000$. The bias of the $f_{b}$ estimates is $0.03$ with dispersion of $0.05$. For $\gamma$ estimates, they are $-0.29$ and $0.34$. Compared to the results for $N=2000$, the systematics does not improve any more, although the random errors mildly improve.

To summarize, the exercises with mock stars are evident that when the sample size is a few hundreds, the binary fraction estimates is reliable but suffers from large random errors. The mass-ratio distribution estimates, however, may be significantly overestimated when binary fraction is low. When $N\gtrsim 2000$, the estimations for $f_{b}$ and $\gamma$ are already robust. Although increasing sample size can reduce the random error, a slight systematic bias does not disappear.

We therefore conclude that, first, we need to keep the sample size to be larger than 2000 so that the random errors and systematics can be well controlled. This is why we only derive $f_{b}$ and $\gamma$ for the [Fe/H]-$m_{1}$ bins with $N>2000$. Second, the readers have to be aware of that the results shown in section 4 may slightly overestimate the binary fraction by around $0.03$ and underestimate $\gamma$ by $0.2$–$0.3$. Finally, although the covariance between $f_{b}$ and $\gamma$ does exist and leads to larger uncertainties when either $f_{b}$ or $\gamma$ is small, there is no strong degeneracy between $f_{b}$ and $\gamma$, meaning that they can be determined with reasonable uncertainties and thus, the anti-correlation shown in Figure 7 is real.

5.3 More discussion about the $f_{b}$–$\gamma$ anti-correlation

To double-check the anti-correlation between $f_{b}$ and $\gamma$, we compare the mean distributions of $\Delta M_{G}$ for the two groups of stars in Figure 11. The two distributions are normalized so that the maximum values are at 1. The low-$\gamma$ group (blue line) shows more fraction of stars at $\Delta M_{G}$ between $-0.5$ and $-0.1$ mag than the high-$\gamma$ group, while the two groups have similar fraction of binaries at \Delta M_{G}$$\sim-0.75 mag. This means that the low-$\gamma$ group contains more binaries with moderate $q$ than the high-$\gamma$ group, although they have similar fractions of binaries with q$$\sim 1. This leads to a two-fold affect. First, the fraction of binaries for the low-$\gamma$ group is larger, since it contains more binaries with moderate $q$. Second, $\gamma$ must be smaller for the low-$\gamma$ group, because that its mass-ratio distribution must be flatter than the high-$\gamma$ group since it contains more binaries with moderate $q$. This again confirms that the anti-correlation between $f_{b}$ and $\gamma$ is indeed a real feature.

5.4 Incompleteness of period distribution

In this work, we identify unresolved binaries from HR diagram, which requires the binaries to appear as a single point. When a binary shows recognizable spatial separation between the two companions, it would be mistakenly identified as two single stars rather than one binary. Therefore, our method may miss wide binaries in the samples and lead to incompleteness to some extent. In this section we evaluate the effect of the incompleteness.

The incompleteness of period distribution is mainly subject to the angular resolution of Gaia,which is around $2$ arcsec (Arenou et al., 2018). The corresponding cut-off period varies with distance, as shown with black lines in Figure 12. The cut-off in period leads to incompleteness of binaries in distribution of period. We then approximate the fraction of missing binaries with large periods beyond the cut-off points by adopting the distribution of period from Raghavan et al. (2010). The fractions of missing binaries is displayed as a function of distance with green lines in Figure 12. It is seen that, in the worst case at 50 pc, about one fourth binaries are missed. However, at 300 pc, only 13% binaries, most of which have wide separations, are cut off. Weighted by the number distribution of the observed stars along distance (bottom panel in Figure 12), the averaged fraction of the missing binaries is around 15%.

The missing wide binaries with period beyond the cut-off points become two single stars to the field. The secondaries of the missing binaries may have opportunity to be observed by LAMOST as an independent single star. However, since the secondary is fainter, they have less probability to be targeted in LAMOST survey according to the selection function of LAMOST, which is roughly flat along apparent magnitude (Carlin et al., 2012; Yuan et al., 2015; Liu et al., 2017). Moreover, lots of the resolved secondaries have masses lower than 0.45 $M_{\odot}$, which have been already excluded from our samples. We can reasonably assume that the completeness of LAMOST survey is around 20%, i.e., only one fifth of the stars in a field are finally observed by LAMOST. And we can coarsely assume that about half of the resolved secondaries observed by LAMOST have masses larger than 0.45 $M_{\odot}$ and hence are considered as single stars. Combining these two factors, we obtain that only 10% of the resolved secondaries from the 15% missing wide binaries mistakenly contribute to the single stars. This means that the number of the single stars only increase by a few percents due to the incompleteness of the wide binary. Therefore, we ignore the systematic bias of the number of single stars in our samples.

Finally, considering that the resulting binary fraction is slightly overestimated (see section 5.2) by $\sim 0.02$, the systematic underestimation of the binary fraction due to the incompleteness of the distribution of period reduces from 15% to 13%.

5.5 Effect from high order multiple systems

Because we identify binaries from their larger luminosity, we cannot distinguish the binary and high-order multiple system. Therefore, $f_{b}$ used in the whole paper is actually equivalent with the frequency of multiplicity, i.e. how many observed stellar objects are multiple stars. This is different with the component frequency, which is defined as the averaged number of companions in each stellar object (Duchêne & Kraus, 2013).

Undiscriminating triple (or higher order multiple) system would not change the fraction of binaries, but may distort the mass-ratio distribution. The mass-ratio of the binary in this work is defined as the stellar mass of the secondary over the stellar mass of the primary, which is not the case for high-order multiplicity. In principle, without considering the tertiary contribution in luminosity, we may mistakenly attribute the total incremental luminosity of the secondary, tertiary and other high-order companions to a more massive secondary. Therefore, without discriminating high-order multiple system, the mass-ratio distribution may be systematically bias to larger $q$.

We investigate 27 triple systems with spectral types or mass studied by Raghavan et al. (2010) and find that the mass of tertiary is much smaller than the secondary in 21 of them. Consequently, they would not significantly change the total luminosity. For instance, if the primary has mass of $0.7$ $M_{\odot}$, the mass of the secondary is $0.5$ $M_{\odot}$, then a tertiary with mass of $0.2$ $M_{\odot}$ only makes $\Delta M_{G}$ brighter by $0.02$ mag. When the mass of the tertiary increases to $0.4$ $M_{\odot}$, the $\Delta M_{G}$ increases by 0.08 mag, roughly comparable to the uncertainty of $\Delta M_{G}$. Therefore, when the mass of tertiary is significantly smaller than that of the secondary, its contribution in the mass-ratio distribution can be neglected.

We also find 5 of the 27 triple systems, most of which contain brown dwarf companions, have the secondary and tertiary with similar masses, which can potentially affect the mass-ratio distribution. However, since these stars are only a small fraction in triple systems and triple systems occupy only a small fraction in multiplicities, their effect in the mass-ratio distribution is small and hence can be neglected.

6 Conclusions

In this work, we use more than 50 000 solar-type MS stars selected from LAMOST and Gaia data to derive the binary fraction and mass-ratio distribution as functions of [Fe/H] and $m_{1}$. We develop a hierarchical Bayesian model to simultaneously derive both the binary fraction and the mass-ratio distribution in each small metallicity and primary mass bin.

For the first time, we find that the binary fraction is tightly anti-correlated with $\gamma$, which is the power index of the power-law shape mass-ratio distribution. The data can be separated into two groups with different $\gamma$. We find that the $f_{b}$ of the low-$\gamma$ group with $\gamma<1.2$ displays clear correlation with $m_{1}$ but quite flat relationship with [Fe/H]. The $\gamma$ for this group, on the other hand, shows anti-correlations with both [Fe/H] and $m_{1}$. For the high-$\gamma$ group with $\gamma>1.2$, $f_{b}$ shows significant anti-correlation with [Fe/H] but is not correlated with $m_{1}$. Although $\gamma$ of the group is mildly correlated with [Fe/H] and $m_{1}$, the large uncertainty prefers to that these correlations are not statistically substantial.

Moreover, looking at the distributions of $f_{b}$ and $\gamma$ in the $m_{1}$ vs. [Fe/H] plane, the two groups are quite different. The high-$\gamma$ group favorites smaller $m_{1}$ and larger [Fe/H], resulting in a concentrated ridge-like regime in the $m_{1}$ vs. [Fe/H] plane, while the low-$\gamma$ group distributes in the regions with larger $m_{1}$ or lower metallicity.

These complicated features strongly hint that the field binaries may experience efficient dynamical processing during the time when they still stayed in the embedded clusters in which they were born (Kroupa, 1995a, b). As a result, the properties of the present-day binaries have been substantially reshaped, especially for those with lower primary mass and higher metallicity. Using the resulting $f_{b}$ and $\gamma$ at each [Fe/H] and $m_{1}$ bin as a prior, we are able to identify 5 410 binary star candidates with $P_{s}/P_{b}<0.1$ from the 57 978 samples. And we also derived the mass-ratio for all 57 978 stars.

We find that the number of samples is critical to derive reliable mass-ratio distribution. If the sample size is only a few hundreds, the derived mass-ratio distribution for the population with small binary fraction would systematically bias to high mass-ratio. However, when the sample size is large enough, such as larger than 2000, we can reliably work out both binary fraction and mass-ratio distribution based on the hierarchical Bayesian model.

Acknowledgements

CL thanks Changqing Luo, Pavel Kroupa, Zhanwen Han, Lu Li, Cheng-Yuan Li, Xiao-Bin Zhang, Licai Deng, and Haijun Tian for their helpful discussions and suggestions. This work is supported by the National Natural Science Foundation of China (NSFC) with grant No. 11835057. This project was developed in part at the 2018 NYC Gaia Sprint, hosted by the Center for Computational Astrophysics of the Flatiron Institute in New York City. It was also partly developed at 2018 Gaia-LAMOST sprint workshop supported by the NSFC under grants 11333003 and 11390372. Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.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.

Appendix A Initial stellar mass estimation

The stellar mass of stars is determined by comparing the effective temperature ($T_{\rm eff}$), surface gravity ($\log g$), and metallicity ([Fe/H]) of stars with the PARSEC isochrones (Bressan et al., 2012). For the MS stars in this work, they are mostly not located in the star forming region and hence should be older than 1 Gyr. Therefore, we remove the isochrones with age younger than 1 Gyr so that the young isochrones would not affect the mass determination.

For each star, the initial stellar mass is determined by looking for the maximum likelihood in $T_{\rm eff}$, $\log g$, and [Fe/H] space. The likelihood distribution is written as

[TABLE]

where $T$, $G$, and $Z$ are the effective temperature, surface gravity, and metallicity, respectively. The uncertainties of $T_{\rm eff}$, $\log g$, and [Fe/H] are $\sigma_{T}=110$ K, $\sigma_{G}=0.2$ dex, and $\sigma_{Z}=0.15$ dex, according to Gao et al. (2015). Based on El-Badry et al. (2018a), we consider the averaged contribution of the unresolved secondaries in stellar parameter estimation and add systematic bias of 50 K in $T_{\rm eff}$, uncertainty of 0.1 dex in $\log g$, and uncertainty of 0.05 dex in [Fe/H] to the total uncertainties and adopt $\sigma_{T}=120$ K, $\sigma_{G}=0.22$ dex, and $\sigma_{Z}=0.16$ dex.

Because that the age for the late type MS stars is very difficult to be determined, we marginalize the age in Eq. (A) to derive the likelihood distribution of $M_{ini}$. We choose the value of $M_{ini}$ at the maximum likelihood and determine the uncertainty from the 15% and 85% percentiles. The likelihood approach can reach typical precision of about 0.04–0.06 $M_{\odot}$ for the stars with initial stellar mass between 0.4 and 1.0 $M_{\odot}$. For binary system, the stellar parameters are dominated by the primary stars. Therefore, we consider the mass estimates for the binaries as the mass of the primary.

Figure 13 displays the locations of 10000 randomly selected MS stars from the samples in $M_{ini}$ vs. $T_{\rm eff}$ plane with the color coded [Fe/H].

Appendix B Interstellar extinction correction

We first derive the color excess of $G_{BP}-G_{RP}$ and then estimate the extinction of $G$ band from it with adopted extinction coefficients for $G$ band. Figure 14 shows the distribution of the sample stars in $G_{BP}-G_{RP}$ vs $T_{\rm eff}$ plane. According to Ducati et al. (2001), Wang & Jiang (2014), and Xue et al. (2016), the bluest edge at a given $T_{\rm eff}$ is contributed by the stars with nearly no reddening in $G_{BP}-G_{RP}$. Jian et al. (2017) applied this approach and well determined the reddening of the LAMOST data. We follow the same method and adopt that the 5% percentiles of $G_{BP}-G_{RP}$ at given $T_{\rm eff}$ bins (red stars displayed in Figure 14) represent for the intrinsic color indices corresponding to $T_{\rm eff}$. Therefore, we can estimate the relationship by fitting the intrinsic colors with a cubic polynomial. We find the best-fit intrinsic color index is

[TABLE]

The extinction coefficients for Gaia bands are provided by Danielski et al. (2018) based on Gaia DR1 such as

[TABLE]

where $k_{X}=A_{X}/A_{0}$. Given that $E(G_{BP}-G_{RP})=(k_{BP}-k_{RP})A_{0}$, we can solve Eq (B) for $A_{0}$ by using coefficients, $c_{1}$, $c_{2}$, …, $c_{7}$, corresponding to $G_{BP}$ and $G_{RP}$, respectively. And then $A_{0}$ is brought back to Eq (B) with coefficients corresponding to $G$-band to solve for $A_{G}$. We adopt the coefficients, $c_{1}$, $c_{2}$, …, $c_{7}$ for $G_{BP}$, $G_{RP}$, and $G$, respectively, provided in Table 1 of Gaia Collaboration et al. (2018). As a by-product, we can also obtain $A_{BP}$ and $A_{RP}$ by applying the same method. Adopting the extinction law of Cardelli, Clayton & Mathis (1989) with $R_{V}=3.1$, we can also estimate the color excess of $E(G-K_{s})$ and subsequently give the intrinsic color index $(G-K_{s})_{0}$.

It is noted that the estimation of $(G_{BP}-G_{RP})_{0}$ does not account for the effect induced by metallicity. We tried to derive the intrinsic color index for different metallicity and found that the difference of the intrinsic color index with metallicity is only a few hundredth magnitude. Therefore, we ignore the effect of metallicity in the intrinsic color index.

The uncertainty of the reddening correction is contributed from several channels. First, the cubic polynomial relation may affected by the uncertainties from the observed color index, which is usually as small as a few mill-magnitudes. Then, the uncertainty of $T_{\rm eff}$ from the LAMOST pipeline, which is around 120 K, is propagated to $E(G_{BP}-G_{RP})$ and contribute to about $0.04$ mag. Finally, total uncertainty of $A_{G}$ is about $0.05$ mag.