The statistical properties of early-type stars from LAMOST DR8

Massive binary stars play a crucial role in many astrophysical fields. Investigating the statistical properties of massive binary stars is essential to trace the formation of massive stars and constrain the evolution of stellar populations. However, no consensus has been achieved on the statistical properties of massive binary stars, mainly due to the lack of a large and homogeneous sample of spectroscopic observations. We study the intrinsic binary fraction $f_{\rm b}^{\rm in}$ and distributions of mass ratio $f(q)$ and orbital period $f(P)$ of early-type stars (comprised of O-, B-, and A-type stars) and investigate their dependences on effective temperature $T_{\rm eff}$, stellar metallicity [M/H], and the projection velocity $v\sin{i}$, based on the homogeneous spectroscopic sample from the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST) Data Release Eight (DR8). We found that $f_{\rm b}^{\rm in}$ increases with increasing $T_\mathrm{eff}$. The binary fraction is positively correlated with metallicity for spectra in the sample. Over all the $v\sin{i}$ values we considered, the $f_{\rm b}^{\rm in}$ have constant values of $\sim$50\%. It seems that the binary population is relatively evenly distributed over a wide range of $v\sin{i}$ values, while the whole sample shows that most of the stars are concentrated at low values of $v\sin{i}$ (probably from strong wind and magnetic braking of single massive stars) and at high values of $v\sin{i}$ (likely from the merging of binary stars). Stellar evolution and binary interaction may be partly responsible for this.There are no correlations found between $\pi$($\gamma$) and $T_{\rm eff}$, nor for $\pi$($\gamma$) and [M/H]. The uncertainties of the distribution decrease toward a larger sample size with higher observational cadence.


Introduction
Many studies show that most stars are in binary systems (Heintz 1969;Abt & Levy 1976;Duquennoy & Mayor 1991;Machida et al. 2008), especially early-type stars (Sana et al. 2012;Chini et al. 2012;Duchêne & Kraus 2013;Dunstall et al. 2015;Moe & Di Stefano 2017). Statistical analyses of binary stars, including intrinsic binary fractions, distributions of orbital period, and mass ratios, are important tracers of stellar formation and basic physical inputs for binary population synthesis Liu 2019;Han et al. 2020). In particular, massive binaries and their statistical properties are essential for understanding compact objects' formation. These include double black holes, double neutron stars, and neutron star-black hole binaries, as they are the dominant gravitational wave sources for the Laser Interferometer Gravitational-wave Observatory (LIGO), Virgo and KAGRA detectors (Abbott et al. 2016a,b;Chen et al. 2018;Han et al. 2020;Langer et al. 2020).
Among the statistical properties, the binary fraction f b is one of the most critical parameters, and it has been widely explored over the past few decades (Carney 1983;Henry & McCarthy 1990;Duquennoy & Mayor 1991;Fischer & Marcy 1992;Mason et al. 1998;Latham et al. 2002;Kratter et al. 2010;Rastegaev 2010;Raghavan et al. 2010;Duchêne & Kraus 2013;Sana et al. 2013;Tanaka & Omukai 2014;Moe & Di Stefano 2017;Liu 2019). These studies show that the binary fraction may depend on the mass, metallicity, and population age of the sample stars. The results from the various studies differ, and are contro-Article number, page 1 of 10 arXiv:2209.09272v3 [astro-ph.SR] 7 Nov 2022 versial in some cases. For example, Sana et al. (2012) reported that Galactic O-type stars have an estimated binary fraction of 69%, while Kobulnicky et al. (2014) suggested a binary fraction of 51% for O-type stars using a different sample of stars. Latham et al. (2002) reported no correlation between the f b and metallicity for stars in their sample. Later works by Raghavan et al. (2010), Tian et al. (2018), and Liu (2019) suggest that an anticorrelation of the binary fraction and metallicity appears in their studies. Alternatively, a positive correlation between these two quantities was reported from Carney (1983) and Hettinger et al. (2015).
For the studies mentioned above, many works were based on small sample sizes, and on heterogeneous samples, especially for massive stars or early-type stars. The ongoing spectroscopic surveys, such as the Sloan Digital Sky Survey (SDSS) (York et al. 2000) and the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST) (Cui et al. 2012;Liu et al. 2020), provide good opportunities to study stellar statistical properties with huge homogeneous samples. Yuan et al. (2015) found an average binary fraction of 41% ± 2% for field FGK-type stars using spectroscopic sample from the SDSS. The fractions (at least 10%) of carbon-enhanced metal-poor (CEMP) stars and doublelined spectroscopic binaries, based on metal-poor stars selected from SDSS, are investigated by Aoki et al. (2015). Based on low-resolution observation from the LAMOST DR5, Luo et al. (2021) reported a binary fraction of 40% for O-and B-type stars in the sample using spectra with more than three observations. Later work by Guo et al. (2022) utilizing the medium-resolution spectra with a cadence of more than two from the LAMOST DR7 suggests a binary fraction of 68% for the early-type OB stars.
In the work of Guo et al. (2022) (hereafter Paper I), the relative radial velocity is obtained by a maximum likelihood estimation, and a star is recognized as a binary according to the maximum variation of radial velocity. The intrinsic binary fraction can be obtained by correcting for any observational biases appearing in the sample using a series of Monte Carlo simulations from Sana et al. (2013). Since most of the early-type star spectra from Guo et al. (2022) only have two observations, significant errors would exist, although the large number of stars in the sample may reduce the errors to a certain degree. This influence has not been examined and corrected due to the small cadence number of spectra from LAMOST DR7. For the same reason, the dependence of the statistical properties on metallicity is not investigated in Paper I.
Motivated by the recent release of more than six million medium-resolution spectra from LAMOST DR8, in this paper we examine the uncertainties of the method used in Paper I and Sana et al. (2013), using both the sample size and the number of observations. We aim to improve the investigation of the binary fraction for early-type stars in the LAMOST DR8 database using an increased sample size (886 stars) with a higher observational cadence (≥ 6) to determine the dependences of the statistical properties of these stars on their metallicity.
The structure of the paper is as follows. We introduce the data of LAMOST-MRS in Section 2. In Sect. 3 we describe our work of dividing the sample into different groups based on their parameters and radial velocity (RV) measurements, and the criterion for identifying binaries. We briefly describe the Monte Carlo method used to correct the observational bias, the consistency check, and the details of testing the suitability of this method in Section 4. We report the results of the relationships between the statistical properties and the stellar spectral types, as well as the metallicity for early-type belonging to various subgroups from

LAMOST data and sample
LAMOST is a four-meter quasi-meridian reflecting Schmidt telescope located at the Xinlong station of the National Astronomical Observatory, implemented with 4000 fibers. Both medium-resolution spectrographs with a resolving power of R ∼ 7500 and low-resolution spectrographs with a resolving power of R ∼ 1800 were installed on the telescope (Cui et al. 2012;Zhao et al. 2012;Deng et al. 2012). In October 2018, LAMOST began a new medium-resolution survey (MRS) to obtain multi-epoch observations. The spectra of MRS observations made from the blue arm have a wavelength range of 495 ∼ 535 nm, and cover a wavelength range of 630 ∼ 680 nm for the red arm (Liu et al. 2020). Guo et al. (2021) utilized the spectra from the LAMOST MRS database to derive the atmospheric parameters of 9,382 early-type stars via a data-driven technique called stellar label machine (SLAM) (Zhang et al. 2020a,b). In this study, we adopt the early-type stars from Guo et al. (2021) to investigate their statistical properties.

Sample selection and RV measurements
3.1. Classifying the sample As mentioned in Sect. 1, both the sample size and observation frequency affect the uncertainty of binary fraction and the distributions for mass ratio and orbital periods. In the work of Sana et al. (2013) about 93% of the samples have more than six observations. Luo et al. (2021) found that an increase in the observational frequency of a sample may lead to smaller uncertainties. We thus selected the sample with more than six observations from the catalog of Guo et al. (2021) for our study. There are a total of 886 stars in the sample. We also did some tests on the effect of observational frequencies and sample size to examine the uncertainties of the method, and the results are shown in Sect. 5.
Based on the atmospheric parameters given by Guo et al. (2021), we divided the samples into different groups according to three observables: effective temperature (T eff ), metallicity ([M/H]), and projected rotational velocity (v sin i). We simply divided the sample stars into low, medium, and high groups based upon their spectral types, metallicity, and v sin i because the aim of the grouping is to explore trends between those atmospheric parameters and the binarity of the samples. For T eff , we divided

Radial velocity measurements
The data reduction was carried out via the standard LAMOST 2D pipeline, and the wavelength calibration was accomplished using the Sc and Th-Ar lamps for MRS spectra. The typical accuracy of wavelength calibration is 0.005nm pixel −1 for LAMOST-MRS (Luo et al. 2015;Guo et al. 2022;Ren et al. 2021). However, Liu et al. (2019) and Zhang et al. (2021) concluded that significant systematic errors exist among the RVs obtained from spectra collected by the different spectrographs, the exposures of LAMOST MRS surveys, and the temporal variations of the RV zero-points. Therefore, Zhang et al. (2021) applied a robust self-consistent method using Gaia DR2 RVs to deter-mine the RV zero-points (RVZPs) from the exposures for each spectrograph of LAMOST after measuring the RVs based on the cross-correlation function method.
We cross-matched our 886 early-type stars with the updated RV catalog 2 for LAMOST DR8 from Zhang et al. (2021) and adopted their corrected RV measurements. In Table 2, for each star we list the observation ID, coordinates, MJD date, S/N, T eff , [M/H], v sin i, RV, and the uncertainty (σ i ).

Criterion for the binary
To identify the binaries in our sample, we adopted Equation 4 from Sana et al. (2013), which states that a star is a binary if its RVs satisfy where v i( j) is the RV measured from the spectrum at epoch i (j) and σ i( j) is the associated uncertainty. This criterion has been applied in several similar studies, for example Sana et al. (2012) According to statistical standards, the confidence threshold of 3σ suggests a probability of finding 30 false detections in a sample of 1000 data entries, and this threshold is sufficient to filter out any outliers. We adopted a stricter criterion of 4σ, as suggested by Sana et al. (2013). This threshold criterion indicates a probability of finding one false detection in every 1000 entries.
The threshold C is used here to filter stars with pulsations, which may also lead to significant RV variation. Sana et al. (2013) adopted C=20 km s −1 for O-type stars, and Dunstall et al. (2015) adopted C=16 km s −1 for B-type stars, based on the kinks in the RV distributions of their sample. We did not find the appearance of a kink feature in our sample, hence simply adopted C= 16 km s −1 as that in Dunstall et al. (2015) because our sample is dominated by B-type stars.
The observed binary fraction f obs b could be underestimated by the constraint of C because of the binary stars, even with large RV amplitudes, but small phase differences of the spectroscopic observations might be missed, and similarly binaries with relatively small RV variations might also be missed. Moreover, a large value of C would lead to a small value of f obs b from the same sample. However, this biases could be corrected by the Monte Carlo simulation described in Sect. 4. Here we just show (see Fig. 2) the results of the corrections for the O-B3 group. In the figure, the circles represent the observed binary fractions ( f obs b ) of binaries in the sample (i.e., without applying any corrections), and these distributions are heavily affected by the arbitrary choice of C values. The squares stand for the intrinsic binary fractions ( f in b ) after correcting any observational biases using a Monte Carlo simulation (discussed in Sect. 4.1). As shown in the figure, the intrinsic binary fractions are constantly distributed over the C values, suggesting that it is safe to adopt the C = 16 km s −1 threshold criterion Mahy et al. (2022).

Monte Carlo method
Here we use the approach described in Sana et al. (2013), where several Monte Carlo (MC) simulations are run to assess the intrinsic binary fraction ( f in b ) based on f obs b . To perform the simulations, we need to construct two synthetic cumulative distributions (CDF) of RV variation (∆RV), the peak-to-peak RV variation of each star, and the minimum timescale between the exposures (∆MJD) by assuming the orbital configurations. We use a power law to describe the distribution for both the orbital period P and mass ratio q: f (P) ∝ P π and f (q) ∝ q γ . We do not use (log P) ∝ (log P) π , as is done in the literature, because the linear form is more sensitive to short-period binaries (Guo et al. 2022). Table 3 lists the variable ranges for the above parameters and power index. We adopt the initial mass function (IMF) from Salpeter (1955) with α = −2.35 3 for the primary masses ).
In general, two-body kinetic systems in a Keplerian orbit are described by the orbital parameters, which are inclination (i), semimajor axis (a), the argument of periastron (ω), eccentricity (e), and the epoch of periastron (τ). The inclination is randomly drawn over an interval from 0 to π/2 and satisfies a probability distribution of sin(i). The semimajor axis is correlated to the orbital period (P). The ω satisfies a uniform distribution randomly drawn from it from 0 to 2π. We use e η to describe the distribution in which η is set to -0.5 from Sana et al. (2013), and τ is selected randomly in units of days.
With the above assumptions we can obtain the simulated RV to further simulate the CDF distribution of ∆RV and ∆MJD. The example for the simulated CDF of ∆RV and ∆MJD is shown in Fig. 3. The constructed simulated CDFs are compared with the observed ones by using Kolmogorov-Smirnov (KS) test, and the simulated and observed fractions of detected binaries ( f sim b and f obs b ) are compared using binomial distribution . We then use the optimum values from the global merit function (GMF) projection, following Sana et al. (2013), to estimate the final results. The GMF consists of the KS probabilities of ∆RV and ∆MJD distributions and the binomial probability where N b is the number of binaries in the observations, while N is the sample size.
To examine the applicability of the GMF for stars (more than six observations) from the LAMOST-MRS, we constructed 24600 synthetic CDFs with the known input sets of f in b , π, and γ. Figure 4 shows the choices of the known input sets of the parameters. We then compare the input parameters with the simulated parameters and show the projections of residuals for the comparisons in Fig. 5. The residuals are the differences between the output parameters predicted by the simulations and the input known parameters. The projection is the distribution of the residuals in which the histogram is the marginal distribution for each parameter. For each parameter, we use the 50th percentile of marginal distribution from the residual distributions to estimate offset and mean values of the 16th percentile (p16) and the 84th percentile (p84) to estimate errors. We estimated the uncertainty (unc) using the offset and error of each parameter through Carlo simulation without systematic bias (p50=0), which verifies the reliability of our method.
We also reanalyzed the 360 O-type stars in the sample of Sana et al. (2013). We obtained statistical properties of the sample consistent with those of Sana et al. (2013) (see Paper I for details).

Impact of observational frequency and sample size
In order to investigate the effect of observational frequencies and the sample sizes on the final results, we performed the three sets of tests. First, we examined the impact of observation frequency by using the same sample of stars (N = 154, the same sample size as O-B3), but with different observational frequencies n (i.e., n = 2, 3, 4, 5 and 6). The result is shown in the left panel of Fig. 6, where the black open and filled squares (plus signs) represent the error and offset (uncertainties). More detailed calculations are described in Sect. 4.1. We see that the uncertainty of each parameter becomes smaller with increasing observational frequency, as expected (Luo et al. 2021). In particular, the offset for all three parameters becomes zero when n = 6.
We then investigated the effect of the sample size. We used three samples with the number of stars N=100, 200, and 300, and the observational frequency n set to 5 (the result of n=6 is similar to n=5; here we take 5 as an example). The middle panel of Fig. 6 shows the results. It indicates that the uncertainty decreases with the increasing number of stars, especially for π and γ.
We find that both big sample size and a high observational frequency would reduce the uncertainty of the method. However, for the real data the sample size generally decreases with increasing observational frequency. We therefore make the uncertainty estimates using the number of observations and the sample size of the real sample. We consider four samples here, with n ≥ 3, 4, 5, and 6. The number of stars in each sample is N=459, 300, 240, and 154. The results are shown in right panel of Fig. 6. We see that the uncertainties of binary fraction f in b are very similar (around 0.1) for all the cases, while the uncertainties of π and γ are slightly smaller for the cases of n ≥ 3, 4 than that for the cases of n ≥ 5, 6. The sample size could be the cause for this. This test indicates that the number of stars in the sample may indeed reduce the uncertainty of the results induced by the observational frequency. In the following, we show all the results from the samples with observational frequencies greater than 3, 4, 5, and 6, and selected the case of n ≥ 6 as the standard to compare with previous studies.

Results and discussions
Figures 7 to 9 show the intrinsic statistical properties of f in b , π, and γ for various samples, with observation frequency n ≥ 3, 4, Table 3. Range of different parameters and power indexes used in MC simulation. The column Power Index shows π and γ, the power index of P (orbital period) and q (mass ratio), respectively. f in b is the binary fraction for simulation. The ranges of P and q are shown in the column Parameter Range, while the ranges of the power index are shown in the column Index Range. The last row gives the step of each power index.

Parameter Power law Parameter Range Power Index Index Range
Step --0.20 -1.00 0.04 Notes. As for stars in group O-B3, the range of γ is -2.5 -2.50 since earlier types are assumed to have slightly larger values Dunstall et al. 2015).  5, and 6, and the dependences of these parameters on the spectral type, metallicity, and projection velocity v sin i. Figure 7 shows the dependence of f in b , π, and γ on the spectral type. We note that the intrinsic binary fraction f in b is positively correlated with T eff (i.e., the intrinsic binary fractions increased toward the early-type stars), which is consistent with the result of Raghavan et al. (2010), Duchêne & Kraus (2013), Moe & Di Stefano (2017), and Guo et al. (2022). Marks & Kroupa (2011) also argued that binaries born preferentially efficiently with larger primary mass (Liu 2019).

Dependence on T eff
For the cases of n ≥ 5 and 6, the binary fraction is 76% ± 10%, 60% ± 10%, and 48% ± 10% for stars of O-B3, B4-B7, and B8-A types, respectively. The binary fraction for stars of types B4-B7 and B8-A in the cases of n ≥ 3 and 4 becomes smaller ( f b = 48%±10%) which is expected. For a large observation frequency (i.e., n ≥ 5 or 6), there is a high possibility of obtaining the maximum RV difference of a binary, resulting in more binaries to be verified in comparison to the cases of relatively low observation frequencies. The results for n ≥ 5 and 6 are very similar, indicating that 5 and 6 times are the optimal numbers of the observational frequencies for this method. Figure 10 shows the f in b for the sample of n ≥ 6 in comparison to that in the literature. We see that the result of O-B3-type stars in our study agrees well with that of Galactic O-type stars ( f in b = 69% ± 9%, filled squares) from Sana et al. (2012), and that f in b in our sample is significantly higher than that given by Luo et al. (2021) for the OB stars from LAMOST DR5 (about 40%, with filled circles). In the study of Luo et al. (2021), most of stars (80%) in the sample have three observations. As seen in Fig. 7, the binary fraction is about 48% for B4-A stars for the sample of n ≥ 3, which is consistent with the result of Luo et al. (2021) within the uncertainties. The binary fraction in Guo et al. (2022) is systematically lower than that of this study for each spectral type since most of the samples in Guo et al. (2022) have only two observations, leading to many binaries being missed by the method, as mentioned in Sect. 4.
The values of π and γ remain almost constant with the spectral type of stars, possibly indicating the same origins of the binary star formation with various temperatures in this range. Our study gives π = −0.9 ± 0.35, −0.9 ± 0.35, and −0.9 ± 0.35, and γ = −1.9±0.9, −1.1±0.9, and −2±0.9 for stars of O-B3, B4-B7, and B8-A type, respectively.  (Hettinger et al. 2015). Machida et al. (2009) suggested that a metal-poor cluster is more likely to form binaries via cloud fragmentation, and this may explain some observations with high f in b for metalpoor stars. However, Korntreff et al. (2012) argued that many short-period systems would merge shortly after formation since the density of gas for a newly formed cluster increases due to dynamical evolution of the cluster. Blue stragglers, formed by binary evolution or dynamical interaction of clusters, provide evidence for these processes (Yanny et al. 2000;Latham et al. 2002;Lu et al. 2010).

Dependence on [M/H]
Similarly, we do not find obvious trends for the π and γ val- Intrinsic statistical properties of f in b , π, and γ for stars with observational frequencies greater than 3, 4, 5, and 6 in groups of classification based on T eff . The x-axis in the left panel is based on spectra type, while in the right panel it is based on T eff . The top panel represents the trend of f in b , and the middle (bottom) represents the trend of π (γ). The blue dot-dashed and red solid lines represent the results of stars classified based upon T eff in groups of B8-A, B4-B7, and O-B3 with observational frequencies greater than five and six, respectively, while the gray and yellow dashed lines represent frequencies greater than three and four.

Dependence on the projection velocity of v sin i
The results for the dependence of statistical properties on the projection velocity are shown in Fig. 9. There is no evidence showing the correlation between the statistical parameters and v sin i. The intrinsic binary fraction f in b is around 50% for all cases, as shown in the figure. For the sample of n ≥ 6, f b = 56% ± 10%, 48% ± 10%, and 48% ± 10%; π = −0.9 ± 0.35, −0.7 ± 0.35, and −0.9 ± 0.35; γ = −3.3 ± 0.9, −3.6 ± 0.9, and −3.9 ± 0.9 for v sin i <35 km s −1 , 35 ≤ v sin i<70 km s −1 , and v sin i ≥70 km s −1 , respectively.
The absorption lines of fast rotators (v sin i over 200 km s −1 ) become very flat in medium-resolution spectra, which can only be detected by high-resolution spectra. These objects were therefore filtered out when we selected early-type stars based on equivalent widths of H and HeI lines (Guo et al. 2022). The number of the fast rotators is very small and does not have a significant impact on the results. On the other hand, the flat absorption lines would cause large RV errors, leading to an overestimate of f in b (since Equation (1) could be satisfied more easily). In any case, the stars in groups of v sin i < 70 km s −1 are not affected. For stars with v sin i > 70 km s −1 , we display the ∆RV distribution in Fig. 11 in comparison with the other two groups. We found that the three distributions are similar, indicating that the measurement of v sin i in our sample is not signif-icantly affected by the RV measurement. So the dependence of f in b on v sin i would not change due to these errors.
Although the values of v sin i estimated by the machine learning are systematically lower than those derived from highresolution spectra (see Fig. 8 of Guo et al. 2021), the result that f in b has no correlation with v sin i does not change. This is consistent with previous studies; for example, Fig. 1 at the website 4 (Głȩbocki & Gnaciński 2005) shows that O, B, and A stars rotate faster than the stars later than A, while there is no obvious dependence between spectral type and v sin i among O, B, and A stars. Our study is for O, B, and A stars, and is in close agreement.
Since the value of v sin i may give hints to star formation and binary evolution, we further compare the frequencies of binary population and the whole sample on the value of v sin i. . We see obvious differences in the CDFs between the binary populations and the whole samples.
In Fig. 12, for stars of O-B3 type (in red), the CDF of the whole sample, in comparison to that of the binary population, is larger when v sin i < 40 km s −1 and gradually becomes smaller after that. This variation indicates that the likely single stars in the sample have a low-vsini group and a high-vsini group, while the binaries have a relatively flat distribution of vsini. This is possibly related to stellar evolution and binary interaction. The lowvsini group of the likely single stars is the result of strong wind and magnetic braking of massive single stars (Matt & Pudritz 2005;Ekström et al. 2008Ekström et al. , 2012, while the high-vsini group is likely from the merger of binary stars . Binary interaction gives the binary population a relatively wide range of vsini (or orbital periods if tidally locked). The case of stars with of types B4-B7 (in blue) is similar, but the high-vsini group of the likely single stars seems to have a relatively low vsini value (in the range of 50-100 km s −1 ) in comparison to that of O-B3 stars. It is a little different for B8-A-type stars (in black). The CDF of the binary population is larger than that of the whole sample when v sin i < 100 km s −1 , and becomes smaller after that. The CDF of the binary population also indicates a relatively wide range of binaries. In particular, there are some binaries having vsini slightly higher (v sin i > 95 km s −1 ) than that of likely single stars.
The differences of the CDFs for various metallicity groups are more obvious than those grouped by spectral type. v sin i < 40 km s −1 . About 70% of both samples have v sin i less than 40 km s −1 . After v sin i >∼ 40 km s −1 , the CDF of the whole sample increases faster than that of the binary population and is close to 1.0 when v sin i =∼ 240 km s −1 , indicating the majority of the remaining 30% of the stars have v sin i in 40-240 km s −1 in the whole sample. Similarly, the CDFs of the binary population span a very wide range of v sin i values. About 1.5% of the binaries have v sin i that can be larger than 250 km s −1 , as shown by the CDF. We also see that a small part of the binary population have higher values of v sin i in comparison to that of the whole sample for the case of low [M/H].
Several factors affect the values of v sin i, from the fragmentation of clouds and the stellar wind of massive stars to the binary or dynamical interaction. It is hard to assess all the effects in quantity presently, and it is beyond the scope of this manuscript. It is clear to us, as shown in Figs. 12 and 13, that the binary population is quite evenly distributed over a wide range of v sin i values. This property probably comes from binary evolution or interaction (i.e., the components of a binary are tidally locked with the orbit).

Conclusion
We collected 886 early-type stars with more than six observations from LAMOST DR8, divided the sample in three ways We first used the variations of radial velocities of each sample to identify the binaries in the sample and to obtain the observed binary fraction. Based on the Monte Carlo simulations in statistics, we corrected observational biases and estimated the intrinsic statistical properties of the binary fraction, the distributions of orbital period, and the mass ratio.
Our study shows that the intrinsic binary fraction increases with increasing T eff , consistent with what is found in the literature. We also find that the binary fraction is positively correlated with metallicity in our sample, consistent with Carney (1983) and Hettinger et al. (2015), but opposite to what is found in Tian et al. (2018) and Liu (2019). We do not find correlations between binary fraction and projection velocity v sin i. However, it seems that the binary population is relatively evenly distributed over a wide range of v sin i values, while the whole sample shows that most of the stars are concentrated at low values of v sin i, and at high values of v sin i in some cases. Binary evolution may partly account for this.
We examined the uncertainties in our method induced by the sample size and observation frequency systematically in the paper and found that the uncertainties of the statistical properties decrease with increasing sample size and the observation frequency. For the real sample of the LAMOST DR8, the uncertainty of f in b is similar (around 0.1) to that of the samples with observation frequency n ≥ 3, 4, 5, and 6. For the cases of n ≥ 5 and 6, we obtain the binary fraction of 76% ± 10%, 60% ± 10%, and 48%±10% for stars of O-B3, B4-B7, and B8-A type, respectively. The binary fraction becomes smaller ( f in b = 48%) for the samples of n ≥ 3 and 4 since it has a low possibility to obtain the maximum RV difference of a binary, resulting in fewer binaries to be verified in these cases, in comparison to that of relatively high observation frequencies.