Binary fractions of G and K dwarf stars based on the Gaia EDR3 and LAMOST DR5: impacts of the chemical abundances

Basing on the large volume \textit{Gaia} Early Data Release 3 and LAMOST Data Release 5 data, we estimate the bias-corrected binary fractions of the field late G and early K dwarfs. A stellar locus outlier method is used in this work, which works well for binaries of various periods and inclination angles with single epoch data. With a well-selected, distance-limited sample of about 90 thousand GK dwarfs covering wide stellar chemical abundances, it enables us to explore the binary fraction variations with different stellar populations. The average binary fraction is 0.42$\pm$0.01 for the whole sample. Thin disk stars are found to have a binary fraction of 0.39$\pm$0.02, thick disk stars own a higher one of 0.49$\pm$0.02, while inner halo stars possibly own the highest binary fraction. For both the thin and thick disk stars, the binary fractions decrease toward higher [Fe/H], [$\alpha$/H], and [M/H] abundances. However, the suppressing impacts of the [Fe/H], [$\alpha$/H], and [M/H] are more significant for the thin disk stars than those for the thick disk stars. For a given [Fe/H], a positive correlation between [$\alpha$/Fe] and the binary fraction is found for the thin disk stars. However, this tendency disappears for the thick disk stars. We suspect that it is likely related to the different formation histories of the thin and thick disks. Our results provide new clues for theoretical works on binary formation.

However, works focusing on the variations of the binary fraction with α elements are rare. The reason includes not only the relatively small data volume but also the bias of the individual method in period, mass ratio, and distance limitations, etc (Moe & Di Stefano 2017). To overcome the biases and limitations of previous techniques, Yuan et al. (2015a) proposed a stellar locus outlier (SLOT) method to estimate model-free binary fraction for large numbers of stars of different populations in large survey volumes. Applying the method to about 10,000 stars from the Sloan Digital Sky Survey (SDSS) Stripe 82, Yuan et al. (2015a) have determined the binary fractions of different stellar populations and analyzed their dependency on spectral type and [Fe/H]. They find the highest binary fraction in the Galactic halo and comparable values in the thin and thick disks. The result is consistent with Moe et al. (2019), but against those from Carney (1983) and Latham et al. (2002). For α abundances, Tian et al. (2018) suggest a positive correlation between [α/Fe] and binary fraction for solar-type stars in the thin disk. Recently, Mazzola et al. (2020) perform a detailed study on the close binary fraction as a function of stellar parameters including α abundances and conclude that α elements suppress multiplicity at most values of [Fe/H]. More efforts are needed to explore the variations of the binary fraction with α elements.
Thanks to the unprecedented data volume and quality of Gaia (Gaia Collaboration et al. 2016) and LAM-OST ) surveys, we can apply the SLOT method to a more advanced sample in terms of larger sample size and better data quality. The LAMOST-Gaia stars serve as an adequate sample to explore the binary fractions with various factors, especially [α/Fe] and stellar populations.
The paper is organized as follows. Section 2 describes the data selection, mainly focusing on late G and early K dwarfs in the solar vicinity. The SLOT method is described in Section 3. Section 4 presents the results for the individual samples and explores the impacts of stellar parameters on binary factions. The results are discussed in Section 5 and compared with previous works. Section 6 concludes by summarizing the main findings and discussing future work.

DATA
In this work, we use the Gaia Early Data Release 3 (EDR3; Gaia Collaboration et al. 2020) and the fifth data release of LAMOST (LAMOST DR5; Zhao et al. 2012;Luo et al. 2012), accompanying with the α abundances provided by Xiang et al. (2019). Gaia EDR3 provides 1.8 billion photometric data of G band and 1.5 billion of G BP and G RP bands based on its 36 months' observation. LAMOST database has accumulated more than 8 million stellar spectra in DR5 with a spectral resolution of R ∼ 1800. Effective temperature T eff , surface gravity log g, and metallicity [Fe/H] are delivered from the spectra through the LAMOST Stellar Parameter pipeline (LASP; Wu et al. 2011) with the precision of about 110 K, 0.2 dex, and 0.1 dex, respectively (Luo et al. 2015). Xiang et al. (2019) use a data-driven Payne approach to provide abundances of 6 million stars for 16 elements as well as the [α/Fe] with the precision of about 0.03 dex, which is defined as a weighted mean of [Mg/Fe] To grantee the quality of the sample, the following constraints are required: For Gaia data: 1) phot bp mean flux over error > 100 2) phot g mean flux over error > 100 3) phot rp mean flux over error > 100 4) phot proc mode = 0

5) duplicated source = False
For LAMOST data: 6) the signal-to-noise ratios for the g band (S/N g ) are larger than 20 7) T eff > 4500 K For spatial location:  Niu et al. (2021a,b), which present precise corrections on colors using the spectroscopy-based stellar color regression method (Yuan et al. 2015b), achieving a precision of about 1 mmag. The latter is provided by Yang et al. (2021), where a machine learning technique is applied to train the observed U BV RI magnitudes of about 10,000  Because we have specified the spatial location, we assume that all sources are beyond the source of reddening. All colors referred to hereafter are dereddened with the SFD dust map and the empirically determined reddening coefficients R(G − G RP ) and R(G BP − G RP ), which are dependent on temperature and reddening, as provided in Niu et al. (2021b). We have tested other 2D dust maps like Lenz et al. (2017) and Planck Collaboration et al. (2016) for reddening correction, and found no big differences. The absolute Gaia magnitude in the G band is calculated by M G = G + 5 + 5log 10 ( /1000) − A G , where is the Gaia parallax and A G = 1.890 E(G BP − G RP ) according to Wang & Chen (2019).
We select the main sequence (MS) stars in the T eff and log g diagram using the same criterion in Niu et al. (2021a). As shown in Figure 1, main sequence turn-off (MSTO) stars are eliminated in the Hertzsprung-Russell diagram following the criteria in Appendix A, as their intrinsic magnitudes can not be precisely determined by their colors and metallicities alone.
To avoid Malmquist bias, i.e., multiple systems can be detected further away than single stars in a magnitudelimited sample, distances are limited to 0.95 kpc. Therefore the upper limit of the G mag is about 16.5 where the color and magnitude corrections are valid.
Since one of the main purpose in this paper is to investigate the relationship between the α abundances and the binary fractions, we limit our sample to be within a range from 0.55 to 0.8 solar mass. A wide α abundances range is engaged by this. Stellar mass is estimated by interpolating the color mass relation provided by the PARSEC (Bressan et al. 2012;Tang et al. 2014) model at every 0.1 dex from −1.5 to 0.5 dex of [Fe/H], as shown in Figure B2 (see Appendix B).
Finally, our sample contains 86,141 objects. The distributions of stellar parameters of our sample are shown in Figure 2. In panel a, the discontinuity of the magnitude distribution around G ≈ 14 mag is due to the selection functions of the LAMOST (Niu et al. 2021a) and has no impact on this work. Panel b shows the spatial distributions. Panel c gives the distributions of the intrinsic color and effective temperature. Our sample covers wide ranges of [Fe/H] of about 2 dex and [α/Fe] of 0.5 dex. As expected, a bimodal distribution of [α/Fe] is clearly seen in panel d, corresponding to the chemically defined Galactic thin and thick disks (Fuhrmann 1998(Fuhrmann , 2004

SLOT method overview
Proposed by Yuan et al. (2015a), the SLOT method is developed to provide model-free estimations of binary factions for large numbers of stars. They have applied the SLOT method to two samples of Stripe 82 stars by combing the recalibrated SDSS photometric data with the spectroscopic information from the SDSS and LAM-OST surveys. After fitting the u−g, g −r, r−i, and i−z colors as a function of the g − i color and [Fe/H], they find that the fitting residuals are asymmetric, pointing to the presence of a significant population of binaries. This method works well for not only close binaries but unresolved wide binaries as well. It is also insensitive to the assumed mass ratio distribution and does not require multiple epoch data.
In Niu et al. (2021b), we have fitted the metallicitydependent stellar loci for the Gaia colors (see Table 1) and found the asymmetric G BP − G residuals (defined as loci predicting colors − observed colors, hereafter the same) after correcting the magnitude-dependent systematic effects of colors. To be more clearly, Figure 3 shows the distributions of the G BP − G residuals at different G magnitudes from 12 to 16 mag. The observed G BP − G residuals, as shown in the blue histograms, are fitted with Gaussian models. The results are over-plotted in the red lines. There are significant discrepancies between the symmetric Gaussian distributions (red) and the observed G BP − G residuals (blue) in the first 3 panels. The discrepancies weakens in the last panel when G ∼ 15.9 mag, due to the larger photometric errors at fainter magnitudes 1 . Fortunately, almost all the sample stars in our study are brighter than 16.5 mag in the G band, thus they do not suffer from this problem.
As explained by Figure 1 of Yuan et al. (2015a), the SLOT method estimates the binary fraction of a given sample by modeling the distributions of the observed color residuals with respect to the metallicty-dependent stellar loci. Here we reproduce the illustration as Figure 4 for better reading. For a given metallicity, single stars lie on the stellar locus, while binary/multiple systems deviate systematically from the locus. Two sets of Monte Carlo simulations are performed, assuming that all stars in the sample are either single or binary. Then by adjusting the relative fraction of stars in the two sets of simulations to fit the observed residual distributions, the binary fraction of the sample is determined. On the other hand, due to the existence of binary stars in the sample when fitting the metallicity-dependent stellar loci, systematic offsets µ are involved (shown by the gray dotted line in Figure 4). It only causes the mean values of the distributions of the observed G BP − G residuals not strictly zero, shifting less than 1 mmag from Figure  3, while the mean values of the simulated distributions are zero. Note that the µ and binary fraction f b are tightly correlated. For convenience, both the µ and f b are set as free variables in this work. In Section 4, one can see that the results of µ and f b are strongly correlated as expected.
A minimum χ 2 technique is adopted in the fitting: Considering the maximum of the binary influences on the G BP − G residuals and various uncertainties (see the next subsection), when calculating the χ 2 , G BP − G residuals are equally binned in the range from −7 to 4 mmag with a bin size of 0.275 mmag, because the simulated residuals barely reach values beyond these ranges. N i obs , N i single , and N i binary are the star numbers that are observed, predicted by the single star simulation, and predicted by the binary star simulation in the ith bin, respectively. Both the single and binary sets of simulations are carried out 100 times to eliminate the random errors. For a given binary fraction f b and an offset µ, the residuals of the single and binary sets are combined according to Equation 2, and then shifted toward the bluer side by the offset µ. N i sim is the finally predicted star number by the simulation in the ith bin. N bin = 40 is the bin number of the histogram. We vary f b from 0.0 to 1.0 at 0.01 intervals and µ from 0.0 to 1.5 mmag at 0.01 mmag intervals. χ 2 is calculated for each set of f b and µ. By choosing the global minimum χ 2 , the best solution is determined. In order to explore the binary fraction of the field star as a function of metallicities, we divide the whole sample into 35 bins in the [Fe/H] and [α/Fe] diagram. Bin width for [Fe/H] is fixed to be 0.1 dex, and it is flexible for [α/Fe] depending on the number density.
We use the Bootstrap Sample method to evaluate formal errors of the binary fractions. Random sampling is done by 1,500 times for one bin containing the most objects of 4,801. We fit the distribution of 1,500 f b values by a Gaussian model and get σ=0.04 as the formal error of the binary fraction of that bin. Errors for other bins are then estimated assuming they are inversely proportional to the square of the ratio of their numbers to 4,801, e.g., about 0.1 for a sample size of 1000.

Descriptions of the simulation
Here we describe simulations of the binary and single sets step by step. A schematic description is shown in Figure 5. For each object, the observed G BP − G RP , [Fe/H], and absolute G magnitude are the inputs. Error in this work is assigned and added through a normal distribution. We use superscript to indicate values after adding errors. When simulating the single set, steps 1-6 are as below: Step 1 : Given the intrinsic G BP − G RP and [Fe/H] , the intrinsic G BP − G is calculated with the metallicitydependent loci in Table 1.
Step 2 : The absolute G BP and G RP magnitudes are obtained by combining the absolute G magnitude and the intrinsic colors.
Step 3 : Random errors in the three magnitudes are added.
Step 4 : Then G BP − G and G BP − G RP are derived.
Step 5 : When yielding the simulated (G BP − G) , σ ZP,sys (see the next subsection) is acquired according to the G magnitude and added to the G BP − G . (G BP − G) loci is computed using G BP − G RP and [Fe/H] by the loci at the same time.
Step 6 : Simulated G BP − G residual defined as the difference between (G BP − G) loci and (G BP − G) is the output.
When simulating the binary set, we assume that the observed G BP −G RP is belong to the primary as it is the brighter component. Then the primary mass is determined based on the PARSEC model. The binary mass ratio distribution is assumed to follow a power law with the index γ= 0.3 (Duchêne & Kraus 2013). Mass of the secondary is obtained with a mass ratio q. When M 2 < 0.08M , we treat the binary system as one single star, performing the simulation following steps 1-6. When M 2 ≥ 0.08M , the intrinsic G BP − G RP of the secondary is interpolated based on the same PARSEC model assuming the metallicity of the binary system is uniform. Next, the absolute magnitudes of the primary and the secondary in three bands are separately calculated as steps 1 and 2 do, except the absolute G magnitudes are interpolated from the PARSEC model (see Appendix). Then the combined absolute magnitudes of the binary system are derived. Finally steps 3-6 are applied to the combined magnitudes and colors, G BP − G residual of the binary system are yielded.
Following the above flow chart, we could estimate the G BP − G residual for binary systems of different combinations of G BP − G RP and mass. In Figure 6, the residuals are estimated assuming all errors equal zero and [Fe/H]=0. We can see that the maximum binary influence on the G BP − G residual is about 3 mmag, happening at the intermediate mass ratios. That helps to set the boundaries when calculating the χ 2 .

Error treatments
As described in the previous subsection, different kinds of errors are considered in the simulation, including random and systematic errors of LAMOST [Fe/H] and Gaia photometry. For each object, Gaia EDR3 catalog provides flux errors of three Gaia bands. They are converted into magnitude errors as follows: where F and σ F are respectively the published flux and flux error, ZP is the photometric zero point, σ m is the magnitude error.
The uncertainties of the photometric zero points (ZP) in the VEGAMAG, σ ZP = 2.8 mmag for all three bands, are provided in Riello et al. (2020). However, the uncertainties in the three bands are probably strongly correlated. Therefore, the uncertainties in terms of colors should be much smaller than involving σ ZP of two bands independently. Besides, σ ZP are very likely to be magnitude-dependent, considering the unique observation mode of Gaia and the magnitude-dependent systematic effects discovered previously (Niu et al. 2021a;Maíz Apellániz & Weiler 2018;Casagrande & Vanden-Berg 2018).
We are mostly concerned with the effect of σ ZP on the G BP − G residuals in this work. To get an empirical estimation of this effect, which is represented by σ ZP,sys in this work, red giant branch (RGB) stars are used as ideal single stars in terms of photometry. Because in RGB-MS systems, the RGB stars are much brighter than the MS stars. A study of double-lined spectroscopic binaries (SB2) in APOGEE DR16 and DR17 (Kounkel et al. 2021) shows that only a few percent of SB2 appear to be RGB-RGB system. Even in that case, the asymmetric binary effect is tiny since the colors/masses of the two components should be similar. Therefore, G BP − G residuals of RGB stars should be fully accounted by errors of [Fe/H], random errors of the three bands, and the σ ZP,sys .
A number of 56,871 RGB stars are selected in the same way as Niu et al. (2021a). Corrections of the Gaia magnitudes and colors mentioned in Section 2 are also performed. The matallicity-dependent stellar loci of the RGB stars are listed in Table 1 as well. In order to explore whether σ ZP,sys is magnitude-dependent, the RGB stars are divided into bins of 0.2 mag width in the range from 11 to 16.4 mag. For each bin, we use the single star model shown in Figure 5 to fit their G BP − G residuals. The σ ZP,sys is varied from 0.0 to 2 mmag at 0.01 mmag interval, the corresponding χ 2 values are calculated. By choosing the global minimum χ 2 , the best solutions of the σ ZP,sys against the G magnitudes are determined. Figure 8 plots the smoothed curve of the σ ZP,sys . We can see that σ ZP,sys is magnitude-dependent, showing the minimum at around G ∼ 13.5 mag where owns the best data quality (Riello et al. 2020 Figure 7. black line: Logarithmic cut-off period as a function of distance with a total mass of 1 M . red line: Fractions of missing resolved binaries. Normal distributed phases and inclination angles of the binary system is adopted. The median value is chosen at each distance. 0.75 mmag when G < 15.5 mag, and then increases quickly towards fainter magnitudes. Only a few RGB stars (dwarfs as well) are of G > 16.4 mag, therefore the value of σ ZP,sys at G = 16.4 mag is adopted for stars with G > 16.4 mag.

Correction for resolved wide binaries
In the above simulation, we have assumed that binary stars are unresolved. However, due to the angular resolution of telescopes, 2 for Gaia (Arenou et al. 2018) and 2.5 for LAMOST (Liu 2019), a binary star would be resolved and mistakenly identified as two single stars when its two companions have a recognizable spatial separation ρ > 2 . Given the total mass of the binary system and the separation between two components, the cut-off period P cut−off could be estimated by the Kepler's Laws. Binary fractions with periods larger than the P cut−off would be underestimated. We correct for the effect of the resolved binaries in this subsection.
According to the target selection process of LAM-OST (Yuan et al. 2015c), targets are required to have no neighbour within 5 radius that is brighter than (m+1) mag, where m is magnitude of a given target star. Therefore, the fainter secondary has a little chance of being targeted in the LAMOST survey. We assume that LAMOST only targets the primary stars of the resolved binaries. We use N s , N unresolved b , and N resolved b to represent the number of single stars, unresolved binary systems, and primaries of the resolved binary systems in our sample, respectively. The binary fractions in the above simulation can be described as: Given the orbit distribution of the binary systems, the fractions of the resolved ones could be estimated via Monte Carlo simulations. As suggested by Raghavan et al. (2010), we choose the period distribution to follow a log-normal Gaussian profile with a mean of logP = 5.03 and σ logP = 2.28, where P is in unit of day. A total mass of 1 M is also assumed. We have verified that the result changes slightly with the assumed total mass. We further adopt uniform distributions of orbital inclinations and phases. In addition, studies have shown that the average eccentricity increases with the binary period but in a complex way (Halbwachs et al. 2003;Tokovinin & Kiyaeva 2016;Moe & Di Stefano 2017). The resolved binaries only contribute a few percent, not the main concern in this paper. Hence circular orbits are applied. The cut-off periods at different distances is firstly calculated and plotted in black in Figure 7. Based on the above assumptions and the angular resolution of Gaia, we calculate the fraction of the resolved binary systems as a function of distance. It is plotted in red and can be explained as The corrected binary fraction for a given sample is: where k is the averaged value weighted by the distance distribution. Most stars in our sample have distances around 0.6 kpc (see Figure 2), so that the corrections of the resolved binary fractions is typically about 0.1 times the unresolved binary fractions.

RESULT
In this section, we demonstrate the results of the fitting and the bias-corrected binary fractions of the individual bin at first. The average binary fractions of the whole sample and in Galactic disks and halos are given. We also investigate the dependency of binary fraction on the chemical abundances.

Binary fraction of the individual bin
To present the result of the best-fit model, we take one bin as an example, as shown in Figure 9. The full version of the fitting result for each bin can be found in Appendix C. In the top panel of Figure 9, the blue histogram shows the distribution of the observed G BP − G residuals. The green and red lines plot the modeled residuals of the single and binary sets, respectively. The and effects are discussed in the next Section. The corresponding χ 2 distribution is shown in the bottom panel of Figure 9. To highlight the result, we only contour the model sets with χ 2 less than 1.4 times the minimum χ 2 . The unresolved binary fraction f b and µ are tightly correlated, the higher f b , the larger µ, as we expected. As a result, the global minimum could be simply figured out through an ellipsoid fitting, resulting in a reasonably good minimum χ 2 value of 2.84. The red star donates the best-fit f b and µ. The error bar of f b from the Bootstrap Sample method is over-plotted. f corrected b is also subsequently obtained through Equation 5.
We summarize the f b and the corresponding formal error, f corrected b , and µ in Figure 10 and Table 2. For each individual bin, its ranges of [Fe/H] and [α/Fe] are illustrated by a box in Figure 10 and listed in Table 2. In Figure 10, the f b , f corrected b , and µ of each bin are texted within the box as exampled by the inset plot. The typical resolved binary fractions are 0.03∼0.07 depending on the distance distributions. The 2D histograms imply the number distribution of [Fe/H] and [α/Fe], where the thin and thick disks could be easily identified. We use the colors of orange and blue to mark the thin and thick disks, respectively. The bin with [Fe/H] < −1 is counted as the inner halo. Naturally, we calculate the average binary fractions of the thin disk, thick disk, and halo as listed in Table 3. f corrected b of the thick disk is clearly higher than that of the thin disk, 0.49 ± 0.02 versus 0.39 ± 0.02. It is due to the large difference of distributions of the chemical abundances in the thin and thick disks. While the binary fraction of the halo has a great probability to be larger than others. Notice that there are few objects in the halo in our sample, therefore the binary fraction of the halo has larger errors. The average binary fraction of the late G and early K type stars in the solar vicinity is 0.42 ± 0.01.  Hurley et al. 2002b) have shown that both mass of the primary and the metallicity affect the binary fractions. As the mass range of the sample in this work is narrow, we suppose its impacts on the binary fractions are tiny. Different distributions of the chemical compositions in the Galactic disks and halo therefore lead to different average binary fractions as we have shown above.

We separately regroup the thin and thick disk stars into ten bins in [Fe/H], [α/H], and [M/H], with
approximately equal width. The bias-corrected binary fractions as well as formal errors are deduced via the same procedures. They are scattered in Figure 11, thin disk stars in blue and thick disk stars in red. Linear fittings to these data show that f corrected b roughly decreases with the increasing metallicities, consistent with the previous works (e.g., Yuan et al. 2015a;Moe et al. 2019;Mazzola et al. 2020). This phenomenon is consistent with numerical simulations (e.g., Machida 2008;Machida et al. 2009;Tanaka & Omukai 2014) that cloud fragmentation is suppressed at higher metallicity. The slopes of the linear fittings of thin and thick disk stars are clearly different, thin disk stars suffer a steeper slope than that of the thick disk stars, by about a factor of two. Take  Thin and thick disk stars are believed to have different formation histories. In particular, α-elements enrichments are mainly done by the core-collapse supernovae explosions, which happened earlier than the iron enrichments by the type Ia supernovae explosions. Thick disk stars were formed on a short timescale when little iron was contained in the interstellar medium. Whereas thin disk stars were constantly formed in the late time (e.g., Freeman & Bland-Hawthorn 2002 and references therein). Therefore the formation and evolution process of the thin and disk stars could be different, consequently the binary fractions.
Another difference between the thin and thick disk stars is how the binary fractions vary with

Redder excess of the color residual
In this paper, when fitting the binary effect on the G BP − G residuals, the lower and upper limits of the G BP − G residuals are respectively −7 and 4 mmag. Fractions of the objects beyond these boundaries are typically 0.3 per cent for G BP − G residuals less than -7 mmag and 5 per cent for G BP − G residuals greater than 4 mmag. Reasons of the excess of the observing residuals in Figure C3 include: (1) underestimated photometric errors for few stars in Gaia; (2) contamination from close foreground or background objects; and (3) flux excess due to the unresolved multiple systems. The first effect should be symmetric, it barely affects the fitting results since the ignoring fractions are tiny. The second and third effects are expected to be similar with the binary effect. LAMOST has applied elaborate target selections to avoid the second situation. We also restrict the sample to the high galactic latitude stars, there are little chances of the second situation. The third case is worthy to be discussed.
It is hard to give reliable corrections on the binary fractions caused by the effect from multiple systems. Firstly, it is unpractical to model the various structures of multiple systems. The mass ratio distributions of the short-and long-period subsystems are different (Raghavan et al. 2010;Tokovinin 2008;Bate 2012). The period distributions of the subsystems also vary with different kinds of hierarchies (Tokovinin 2014a). Secondly, the missing multiple fraction is uncertain. Although the fraction of multiple systems is reported to be 0.11-0.17 for the solar-mass stars in general (Raghavan et al. 2010;Duquennoy & Mayor 1991;Bate 2012;Tokovinin 2014b), it changes with the orbital period (Allen et al. 2012;Makarov et al. 2008) in a complex way. Notice that a substantial portion of the multiple systems contain a very low-mass third/fourth component, whose contribution to the color residuals can be ignored. Therefore, such systems have been well accounted in our results.

Comparison with previous works
Comparing to previous works based on other approaches, our method is weakly limited by the binary orbital period. Summarized illustrations of the ranges of binary periods and mass ratios that can be detected with different methods are shown in Moe & Di Stefano (2017) and El-Badry et al. (2018b). For instance, the radial velocity monitoring method and light curve analysis method (for eclipsing binaries) may only work for close binaries. The former is suitable for log 10 (P/day) ≤ 3, while the later works only at even a shorter period of log 10 (P/day) ≤ 2. The common proper motion and visual binaries are only sensitive to wide binaries. The spectral fitting method could achieve log 10 (P/day) ≤ 8 but only for the binaries with intermediate mass ratio.
In addition to periods, these methods also suffer from other specific limitations. Taking several examples, radial velocity monitoring technique requires time-domain data. Photometric technique is confined to the systems with high inclination angles. Astrometric techniques like from Hipparcos and Gaia are biased to the bright binaries having intermediate separations based on the angular resolution and time baseline of the surveys. As for the SLOT method, its realization is solely taking advantage of the difference between the color deduced by the metallicity-dependent stellar loci (of the single star) and the combined color (of binary system). It covers both close and wide binaries, supporting log 10 (P/day) less than 6-8 depending on the distance distribution and instrument spatial resolution. Besides, it depends weakly on the mass ratio distribution and can be carried out with single epoch data.
We stress that most previous studies on binary fraction are restricted to the thin disk stars due to the selection bias. It would be reasonable to compare their results with ours of the thin disk. For example, the G2-K3 sample in Raghavan et al. (2010) has a binary fraction of 0.41 ± 0.03, consistent with with our 0.39 ± 0.02. Our work basically uses the the same method as Yuan et al. (2015a) except for taking the resolved binaries into consideration. Our sample is also comparable with the 0.9 < g − i < 1.2 mag and −1 < [Fe/H] < 0 sample in Yuan et al. (2015a). Their binary fraction is 0.35, in good agreement with the average fraction of the unresolved binary 0.37 in this work. Comparing with Mazzola et al. (2020), where the close binary fraction decreases by a factor of ∼2.4 from [Fe/H] = −0.5 dex to [Fe/H] = 0.25 dex, our binary fraction decreases by a factor of ∼2.3 from [Fe/H] = −0.5 dex to [Fe/H] = 0.25 dex. The decreasing factors seem to be in good agreement. We adopt a universal period distribution for binaries with different metallicities. However, recent works (e.g., Tanaka & Omukai 2014, Jayasinghe et al. 2020 indicate that the binary period shifts toward smaller separations as the metallicity decreases. If it is true, there would be an over-correction of binary fraction for the resolved metal-poor binary. We make a simple estimation that a log-normal Gaussian distribution with a mean of logP = 4 and σ logP = 1.5 for [Fe/H] < −0.2 reduces the binary fraction by ∼0.045 2 . Consequently, it would lead to a slower slope of the metallicity dependency. Given that our sample also contains wide binaries in addition to the close binaries, this could be accounted for that the metallicity dependency of the binary fraction changes with the increasing period (e.g., El-Badry & Rix 2019; Hwang et al. 2021;Bate 2019). As the period increases, the anti-correlation between the close binary fraction and metallicity weakens, even becomes positive.
One motivation of this paper is exploring the binary fraction as a function of the α enrichment. Figure 12 suggests that they are most likely positively correlated. This result is not affected by the assumed period distribution, as it is under the condition of a given [Fe/H]. Mazzola et al. (2020) take the opposite view. The indicator used in their paper is [α/H], which can be easily converted to [α/Fe] for a given [Fe/H] value. Their result that the anti-correlation between α abundances and multiplicity is mainly concluded from the α-poor sample, like −0.2 < [α/H] < 0.1 at solar [Fe/H] (see their Figure 10). The result could be affected by the wider mass ranges (0.5 -2.0 M ) of their sample compared to our sample (0.55 -0.8 M ). On one hand, highmass stars are generally younger than low-mass stars, and subsequently have lower α abundance. On the other hand, the close binary fraction increases by a factor of 2 from 0.5 to 2 solar mass (Moe & Di Stefano 2017). Therefore, the higher multiplicity of the lower α abundance stars could be a joint result of the higher mass and lower α abundance instead of lower α abundance alone. They do find that the close binary fraction flattens for [α/Fe] > 0.05 dex. Moreover, in their Figure B3 In this work, the binary fraction of individual stellar population is derived using the SLOT method and corrected for the effect of resolved binaries. The SLOT method is described in detail in Section 3. Compared with other methods, such as radial velocity monitoring and light curve analysis, our method suffers from less limitation on the period, mass ratio, and inclination, etc.
We use a distance-limited sample of 90 thousand late G and early K dwarf stars selected from Gaia EDR3 and LAMOST DR5 data, in which Gaia photometry has been delicately re-calibrated by Niu et al. (2021b) and Yang et al. (2021). Our sample covers wide ranges in chemical abundances, especially for α abundances comparing with the previous studies. We divide the sample into different bins based on their [Fe/H] and [α/Fe] values. Results of the individual bins are demonstrated in Subsection 4.1. All bins are classified into the thin disk, thick disk, and halo when exploring the tendencies of the binary fractions. The average binary fraction of the whole sample is 0.42 ± 0.01. The thin disk stars are found to have a lower binary fraction (0.39 ± 0.02) than the one of the thick disk stars (0.49 ± 0.02). The halo stars have the highest binary fraction (0.91 ± 0.09). Our estimations of the binary fractions are consistent with previous studies (Raghavan et al. 2010;Yuan et al. 2015a).
As shown in Subsection 4.2, we find that the binary fraction is not just simply anti-correlated with [Fe/H], [α/H], and [M/H]. The impacts of chemical abundances are different for thin and thick disk stars. It is shown that (1) metallicities depress the binary fractions more significantly in the thin disk than in the thick disk; (2) for a given iron abundance, the binary fractions increase with the increasing α abundance in the thin disk, the trend disappears in the thick disk. In particular, the slope of the binary fraction as a function of [Fe/H] alone is −0.37 per dex for the thin disk stars and −0.2 per dex for the thick disk stars. The impacts of [α/H] and [M/H] exhibit similar trends. Also, given a [Fe/H], [α/Fe] increase by 0.05 dex, the thin disk binary fractions typically increase by 0.1. These might be the consequences of the different formation histories of the thin and thick disks. Our results provide new clues for theoretical works on binary formation.
A single-power-law is assumed for the binary mass ratio distribution in this work. In the future, by using distance information, we are going to determine the mass ratio distribution and the binary fraction at the same time.

ACKNOWLEDGMENTS
We acknowledge the anonymous referee for his/her valuable comments that improve the quality of this paper significantly. We acknowledge helpful discussions with Prof. Liu  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.  Table A1.

B. PARSEC MODEL
The PARSEC model is applied in this work when calculating the stellar mass and the absolute G magnitude. We extract isochrones of different sets of [Fe/H] with step of 0.1 dex from −1.5 to 0.5 dex and age with step of 1 Gyr from 1 to 10 Gyr. In Figure B1, grey points are the sample we used, colored tracks are the color-magnitude tracks from the PARSEC model. Red lines that laying on the PARSEC model are empirically fitted using the grey points, in which 2-sigma clipping is performed. Some grey points scatter above the red line as well as the tracks are binaries. Our sample is in good agreement with the PARSEC model. We can see that the width of tracks caused by ages among the color ranges of our sample is narrow. Hence we just interpolate the tracks to have the M G . Figure B2 is the color-mass model of the PARSEC model. Stellar mass is derived in the same way as the M G .
C. FITTING RESULT OF THE INDIVIDUAL BIN Same as Figure 9, the fitting results of all bins are shown in Figure C3 and C4. The typical minimum χ 2 values are between 1.0 -5.0. Almost all µ values are less than 1 mmag. Despite the existence of binary stars in the sample when fitting the loci, our stellar loci listed in Table 1 are reasonably accurate for almost all cases.