The Binarity of Early-type Stars from LAMOST medium-resolution Spectroscopic Survey

We collected a sample of over 800,000 stars from LAMOST-MRS. In these objects, ∼500,000 have stellar parameters ⁸ derived from the LAMOST stellar parameter pipeline (LASP) (Luo et al. 2015, A. L. Luo et al. 2021, in preparation). LASP adopts computational algorithms of the Correlation Function Initial (Du et al. 2012) and the Universite de Lyon Spectroscopic analysis Software (Koleva et al. 2009), and obtains T_eff by minimizing the χ² between model spectra and observed spectra (Luo et al. 2015). Due to the limitation of the pipeline, only stellar parameters for A and later type stars, i.e., with T_eff in range of 3100–8500 K, are given (Wu et al. 2011). The typical uncertainty of T_eff for LRS from the LASP is 110 K (Wu et al. 2014; Gao et al. 2015), while that for MRS was not given yet. Zhang et al. (2020a) showed that the accuracy of simulated MRS spectra is similar to that of LRS for stars with solar abundances.

Since we were only concerned with the multiplicity of early-type stars in this paper, we then eliminated stars with spectral type later than A5, i.e, with estimated T_eff < 8000 K given by the LASP. We included ∼300,000 stars without estimated T_eff in our initial sample because many early-type stars are in this part. However, a broad range of stars of various spectral types is mixed in such a sample. We thus eliminated spurious late-type stars in our sample by the following selection criteria.

The classic method to classify stars with MK types is comparing their spectra manually with standard stars (Gray & Corbally 2009, 2014). This approach is hard to handle with huge samples of observations. The photometric color indices are also widely used in classifying stars, but such a method is heavily affected by interstellar reddening. It is difficult to calibrate the reddening in the deep sky (Liu et al. 2015). In this work, we thus employed the EWs of several lines to exclude late-type stars (later than A5) in our sample (Liu et al. 2015). The EW for a given spectral line is defined as

$$\begin{eqnarray}&&\mathrm{EW}=\int \left(1-\tfrac{{F}_{\mathrm{line}}(\lambda )}{{F}_{\mathrm{con}}(\lambda )}\right)d\lambda ,\end{eqnarray} \tag{ 1 }$$

where F_line(λ) is the flux of the spectral line at a given wavelength λ, and F_con(λ) is the flux of pseudo-continuum and is derived from linear interpolation of fluxes on both sides of a selected line bandpass (Worthey 1994). Due to high effective temperatures for early-type stars, not many line features are shown in the LAMOST-MRS spectra. Consequently, we measured the EW of the Hα λ6564 Å, He i λ6679 Å and Mgb series (Liu et al. 2015, 2017). Table 1 features the information on index bandpass and the two continua bands for these lines (Cohen et al. 1998; Liu et al. 2015). Figure 1 displays the distributions for our initial sample as the black points in (EW_Mgb, EW_Hα ) and (EW_{He I}, EW_Hα ) panels.

Zoom In Zoom Out Reset image size

As the first step, we removed late A- and F-type (LAF) stars from the sample. To do this, we chose LAF stars in LAMOST-MRS, i.e., with 6500 K ≤ T_eff ≤ 8000 K (Habets & Heintze 1981) and calculated the EWs of Hα, He i and Mgb for these stars. Figure 2 depicts the probability density distributions of the EWs in the (EW_Mgb, EW_Hα ) plane (left panel) and in (EW_{He I}, EW_Hα ) plane (right panel), respectively. The color bar on the right side indicates the probability density. The distribution of the LAF stars within the 1σ (68.27%, black), 2σ (95.45%, green) and 3σ (99.73%, blue) regions are shown with black, green and blue lines respectively.

Zoom In Zoom Out Reset image size

In principle, we can remove most LAF stars from the sample if we discard the objects located in the region enclosed by the 2σ line. However, if we do like this, we will lose a majority of early A-type stars in our sample, since the spectral features of F-type stars are similar to those of A-type stars, resulting in EW density distributions of F-type stars being indistinguishable from those of A-type stars. This can be clearly seen in Figure 1 where we overlapped LAF stars within the region of 1σ and 2σ. So, we only eliminated LAF stars within the region enclosed by the 1σ line in this step. ⁹ We will double-check our final sample by visually examining the individual spectrum to exclude any spurious late-type star spectrum after applying all the criteria.

We applied a similar approach to exclude G-type stars with T_eff in the range of 5200–6500 K (Habets & Heintze 1981). The probability density distributions for the G-type stars are depicted in Figure 3. Stars within 2σ, corresponding to the 95.45% confidence level regions of the distribution, were excluded from the sample.

Zoom In Zoom Out Reset image size

Because the LAMOST-MRS pipeline does not assign T_eff to stars cooler than 3100 K, we are not able to apply the above method to remove M- and K-type stars. We cross-matched the stars in our sample using their coordinate designation with the Two Micron All-Sky Survey (2MASS) database (Cutri et al. 2003) (within 3''), and obtained the reddening-free quantity Q (Q = (J − H) − 1.70(H − K_S )) for the common stars. Negueruela & Schurch (2007) suggest that for M- and K-type stars, the Q values are expected to be within a range of 0.4 to 0.5. Therefore, we selected stars satisfying this condition (displayed as red dots in Figure 1) to remove them from our sample.

Since the He i line profile appears prominent in spectra of early-type stars, the line profile gradually decreases in strength toward cooler stars and is mostly absent in the spectra of A-type stars. We thus selected the stars with measured EW_{He I} ≤ 0.1 to remove any spurious late-type stars further. We use the information on T_eff given by Gaia Data Release 2 (DR2) ¹⁰ to remove any missing late-type stars in the sample. We cross-matched the data set within 3'' with Gaia DR2 to find common stars and then keep stars with estimated T_eff ≥ 7300 K (F₀) to be our final sample of early-type stars (Habets & Heintze 1981). ¹¹

Lastly, through visual inspection of the individual spectrum in the sample, we found a residual contamination rate of late-type stars of about 1.6% in our sample. We then removed these stars, as well as spectra displaying emission profiles from the sample. In our final data set, we collected a sample of 9382 stars with signal-to-noise ratio (S/N) > 40. All the targets have a simple magnitude cut-off G ∼ 13 (Liu et al. 2020).

In order to investigate the relationship of the multiplicity properties of early-type stars as a function of T_eff, we first need to group the OBA stars in our sample based upon their T_eff. However, no information on the T_eff for these stars is available. Also, it is difficult to distinguish the spectra of late O-type stars from those of early B-type stars due to their color degeneracy (Maíz-Apellániz et al. 2004) (or late B-type from early A-type stars, Liu et al. 2015). Liu et al. (2019c) has demonstrated the applicability of using the spectral line indices to classify the spectral types of early-type stars in LAMOST-LRS. We thus adopted their approach to divide the stars in our sample into four catalogs in T_eff. We collect a sample of template stars with known spectral types from literature (Hiltner 1956; Nesterov et al. 1995; Negueruela & Marco 2003; Beerer et al. 2010; Comerón & Pasquali 2012; Hou et al. 2015; Liu et al. 2019c). We cross-matched these stars with the LAMOST-MRS database to collect spectra of common stars, including 136 O-type stars, 308 B-type stars and 898 A-type stars. We then measured their EWs of the Hα λ6564 Å and He i λ6679 Å line profiles. In Figure 4, on the left panel, we plot the measured EWs in the (EW_{He I}, EW_Hα ) plane, in which the triangles and circles indicate the EW measurements for O- and A- stars, respectively, and the star symbols represent the B-type stars. The vertical color bar on the right side signifies the spectral sub-type distribution of B-type stars. Based upon the EWs, we group the stars into four groups of T1, T2, T3 and T4. As depicted in the right panel of Figure 4, the stars in group T1 comprise mostly O − B4-type stars, and most B5-type stars reside in group T2, while the majority of B7-type stars are distributed in the region T3. Cooler stars of B8 − A4 are illustrated in region T4 on the plot (Guo et al. 2021). On the right panel, we adopted the criteria to group our observed OBA stars in the sample.

Zoom In Zoom Out Reset image size

Sota et al. (2014) presented a spectral classification for a catalog of 448 Galactic O-type stars from GOSSS. We cross-matched the LAMOST-MRS catalog with GOSSS, and only four non-emission stars were found. We also cross-matched the LAMOST catalog with the published 146 OB stars from Kiminki et al. (2007), and only 14 common sources were found. All the 18 (4+14) sources above with the spectral type from O5.5 to B1 are in the T1 group indicating that it is reasonable to adopt the method to obtain the four groups.

In Table 2, we list the observation ID for each of the stars, their coordinates, Modified Julian Date (MJD), S/N, the associated EW measurements for all three line profiles, their classified group index, star name and V magnitude (from SIMBAD). The number of stars within each group is listed in Table 3.

The RV measurement is a vital ingredient of this study since detecting significant RV variations relies on RV differences. We adopted the maximum likelihood method described in Xiong et al. (2021) to measure the relative RVs. We briefly describe the procedures as follows. For each star, a set of multiple exposures was obtained from the LAMOST-MRS, and we selected the exposure with the highest S/N as the template spectrum. The likelihood distribution for the ith star observed at time t of the relative radial velocity (RV_rel) is

$$\begin{eqnarray}&&p({\mathrm{RV}}_{\mathrm{rel}})={\prod }_{\lambda }\tfrac{\exp \left[-\tfrac{{\left({f}_{t,i}(v,\lambda )-{f}_{\mathrm{tem},i}(\lambda )\right)}^{2}}{2({\sigma }_{t,i}^{2}+{\sigma }_{\mathrm{tem},i}^{2})}\right]}{\sqrt{2\pi ({\sigma }_{t,i}^{2}+{\sigma }_{\mathrm{tem},i}^{2})}}\end{eqnarray} \tag{ 2 }$$

where f_tem,i (λ) and f_t,i (v, λ) represent the observed flux of a template spectrum and individual exposure, respectively, while σ_tem,i (σ_t,i ) is the observed flux error of a template spectrum (individual exposure). A Gaussian fit was applied to each of the likelihood distribution functions, and we estimated the peak value of the fit as our RV_rel measurement, and the uncertainty (σ_unc,i ) of the RV measurements was obtained from the standard deviation of the Gaussian fit. Such an approach yields a systemic error of 0.25 km s⁻¹ (σ_sys).

Systematic errors exist among the RVs obtained from spectra collected by different spectrographs and exposures of LAMOST-MRS surveys (Liu et al. 2019b; Zhang et al. 2021). Therefore, we matched our catalog, including 2911 early-type stars, with the one from Zhang et al. (2021) and found that the median values of the RVs before correction and after correction are ∼0.6 km s⁻¹. In the paper, we used a strict criterion to select the binary stars (see Equation (3)). This means that the influence of uncorrected RV on the derived binary fraction is little and can be ignored. We also compared the observed binary fraction based on the RVs from Zhang et al. (2021) and Xiong et al. (2021) and found that the difference between the observed binary fractions derived in the two ways is less than 1.5%. Therefore, we kept using the RV_rel from Xiong et al. (2021) in the following analysis.

In our final sample of 9382 candidates, only 2911 stars were observed more than twice (see column 3 in Table 3). The following analysis of this work is based on the sample of 2911 stars. The measured RV_rel and the uncertainties (σ_unc,i ) for each star are tabulated in Table 4. The number distribution of the stars as a function of observing baseline is depicted in Figure 5. In Figure 6, we display the number of observations for stars in each temperature group.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We adopted the method from Sana et al. (2013) (see their Equation (4)) to identify binary candidate systems in our sample (Sana et al. 2012; Dunstall et al. 2015; Banyard et al. 2021; Mahy et al. 2021). We considered a star as a binary if RV_rel satisfies the criterion

$$\begin{eqnarray}&&\tfrac{| {v}_{i}-{v}_{j}| }{\sqrt{{\delta }_{i}^{2}+{\delta }_{j}^{2}}}\gt 4\ \mathrm{and}\ | {v}_{i}-{v}_{j}| \gt C,\end{eqnarray} \tag{ 3 }$$

where v_i(j) and δ_i(j) are the RV_rel and the associated uncertainty measured for the spectrum at epoch i (j) respectively. Here we divided the uncertainty into two parts

$$\begin{eqnarray}&&{\delta }_{i}=\sqrt{{\sigma }_{\mathrm{unc},i}^{2}+{\sigma }_{\mathrm{sys}}^{2}},\end{eqnarray} \tag{ 4 }$$

where σ_unc,i is the uncertainty of the ith star observed at time t while σ_sys is the typical systematic error (see Section 4.1).

The threshold C is adopted to eliminate the stars with significant RV variation caused by pulsation of the photosphere or atmospheric activity. Sana et al. (2013) and Dunstall et al. (2015) report that this value is assigned from a kink on the RV distribution plot (see Figure 6 in Sana et al. 2013). We followed a similar approach of searching for a possible kink from the RV distribution to constrain the threshold C value. However, no such pattern is apparent in our RV plot shown in Figure 7. In order to eliminate the pulsational variable stars in our sample, we adopted the typical period and RV amplitude of δ Scuti type stars ¹² to perform a Monte-Carlo simulation to constrain the threshold C value, resulting in a value of 15.57 km s⁻¹ (Breger 1979, 2000; Fernie 1994).

Zoom In Zoom Out Reset image size

According to the criterion (3), the observed binary fractions (${f}_{{\rm{b}}}^{{\rm{o}}}$) of these four groups are 24.6% ± 0.5%, 20.8% ± 0.6%, 13.7% ± 0.3% and 7.6% ± 0.3%, respectively. The ${f}_{{\rm{b}}}^{{\rm{o}}}$ for the stars in each group is plotted as red dots and displayed in Figure 8. The error bars were estimated through bootstrap analysis from Raghavan et al. (2010). The result affirms that ${f}_{{\rm{b}}}^{{\rm{o}}}$ decreases as the spectral type moves from O-type to A-type stars.

Zoom In Zoom Out Reset image size

In fact, this way may underestimate the fraction of binary stars, since the binaries with relatively small RV variations in the measurements based on the LAMOST-MRS spectra are missed. In addition, if the phase differences of the spectroscopic observations of LAMOST are small, even the binary stars with large RV amplitudes might be missed. Besides, a different threshold C-value may affect the result of ${f}_{{\rm{b}}}^{{\rm{o}}}$, but all the bias above can be corrected when we investigate the intrinsic binary fraction through Monte-Carlo simulations in the next section.

Zoom In Zoom Out Reset image size

Following the approach as noted by Sana et al. (2013), we need to adopt the Monte-Carlo simulation to construct synthetic cumulative distribution function (CDF) of RV variance (ΔRV is the maximum RV variance among individual exposure) and ΔMJD (minimum timescale between the exposures). We assume the probability density distribution of orbital period (P) and mass ratio (q) of binary system satisfies the power-law f(P) ∝ P^π and f(q) ∝ q^κ , respectively. We modified the distribution of orbital period in a log-scale of (log P) ∝ (log P)^π to linear form to be more sensitive to short-period binaries. In Table 5, we list the variable ranges for the parameters, and we set the step size to be 0.1 for π and κ, and 0.04 for f_b. We adopt the initial mass function (Salpeter 1955) to describe the mass distribution of stars in our grouped catalogs of T1 ∼ T4.

Generally, six parameters are used to describe two-body kinetic systems in a Keplerian orbit: inclination (i), the argument of periastron (ω), true anomaly (ν), eccentricity (e), semimajor axis (a) and the epoch of periastron (τ). The inclination satisfies a probability distribution of sin(i) and is randomly drawn over an interval from 0 to π/2. The argument of periapsis satisfies a uniform distribution and is randomly drawn from 0 to 2π. We use e^η to describe the distribution and set η as −0.5, the same as in Sana et al. (2013). The semimajor axis depends on the orbital period (P). The τ is selected at random with the unit of days.

We followed the same method as discussed in Appendix C of Sana et al. (2013), and used our selected parameters with assigned range to perform the simulation to obtain cumulative distributions of ΔRV and ΔMJD. We followed a similar approach as Sana et al. (2013) to obtain the global merit function (GMF) using the simulated CDFs and our observations. The final results are given by the projection of the GMFs, i.e., π, κ and f_b .

In order to verify the robustness of applying the methodology from Sana et al. (2013) to our data set, we performed a self-consistency validation test, as well as testing for the suitability of adopting the GMF.

5.2.1. Self-consistency Test

5.2.2. The Applicability of the GMF

In order to test the applicability of applying the GMF to determine the probability density distribution of sample stars from the LAMOST-MRS, we first constructed 68 synthetic CDFs of ΔRV and ΔMJD using the input sets of π, κ and f_b. The choices of such values are listed in Table 6. Based upon the constructed input synthetic CDFs and results discussed in Section 5.1, we then obtained the 68 projections of GMFs. A comparison of the output synthetic GMFs to the input parameters of π, κ and f_b is displayed in Figure 10.

Zoom In Zoom Out Reset image size

Table 6. The Input Parameters of π, κ and f_b

Group	π	κ	fb
T1 (T2)	[−2.2, −1.5, −0.5, 0.5, 1.5, 2.2]	[−2.2, −1.5, −0.5, 0.5, 1.5, 2.2]	[0.28, 0.4, 0.6, 0.8, 0.92]
	κ = −0.5 fb = 0.6	π = −0.5 fb = 0.6	π = −0.5 κ = −0.5
T3 (T4)	[−2.2, −1.5, −0.5, 0.5, 1.5, 2.2]	[−3.2, −2.5, −1.5, −0.5, 0.5, 1.2]	[0.12, 0.2, 0.4, 0.6, 0.72]
	κ = −1.5 fb = 0.4	π = −0.5 fb = 0.4	π = −0.5 κ = −1.5

Note. Logically, the values should be chosen with equal intervals, e.g., chosen as 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 and 0.9 for f_b, but this is time-consuming. Therefore, we selected a few values within the range but with a larger interval. On the other hand, the step size for f_b is 0.04, and we selected the values of two-step sizes away from the boundary for the values closest to the boundary, i.e., we chose 0.28 and 0.92 (the same reason for the other two parameters).

Download table as:  ASCIITypeset image

The top two rows show the constrained results for the groups in T1 and T2 from the synthetic CDFs distribution, while the bottom two rows are for the groups in T3 and T4. We found that the results of stars in groups T3 and T4 (12 and 14 points) are better constrained than those from T1 and T2 (7 and 6 points), in which the point represents the constrained result. The results indicate that the method is applicable for sample stars in groups T3 and T4, but not for stars in groups T1 and T2, especially for κ. In order to explore the reasons for the unconstrained κ value, we repeat the analysis above using the observed CDFs of ΔRV and ΔMJD for stars in group T1. The results are displayed in the top panel of Figure 11. We note that the κ value is still not well constrained. We thus further explore the possible κ values by extending the range to −4.5 ∼ 0.5. The results are displayed in the bottom panel of Figure 11, where we see that κ is still unconstrained. The results for T2 are similar to those of T1. This indicates that the initially selected range of κ values is independent of whether it can be constrained or not. This might be explained by the small value of sample size (N) for groups T1 and T2, which is an initial input of the binomial distribution of the GMF.

Zoom In Zoom Out Reset image size

Considering that a large number of unconstrained κ values for stars in the groups T1 and T2 may affect the results of other parameters, we fixed the κ value to guarantee the reliability of other parameters (Sana et al. 2013). According to the initial mass function, the population of late-type stars shows an increasing trend toward small stellar mass. Therefore, the sample size of early-B type stars is larger than that of the O-type stars for the stars in groups T1 and T2. We thus adopt the value of κ from Dunstall et al. (2015), i.e., κ = −2.8 for early B-type stars and repeat the analysis above again to verify the applicability of the method. Based upon the validation tests, we then obtained the final results of π and f_b from the GMFs by repeating the procedures as mentioned in Section 5.1, except that we now fixed the value of κ to be −2.8 for sample stars in groups T1 and T2.

5.3.1. Results

The GMF projections of π, κ and f_b for stars in groups T1 to T4 are depicted in Figure 12. The red, blue and black contours indicate the loci of 10%, 50% and 95% confidence levels of GMF. The intrinsic binary fractions of these four groups are 68% ± 8%, 52% ± 3%, 44% ± 6% and 44% ± 6%, respectively. The π values of these groups are −1 ± 0.1, −1.1 ± 0.05, −1.1 ± 0.1 and −0.6 ± 0.05, respectively. The κ values of T3 and T4 are −2.4 ± 0.3 and −1.6 ± 0.3, respectively.

Zoom In Zoom Out Reset image size

Note that, although the threshold C-value does affect the observed binary fraction ${f}_{{\rm{b}}}^{{\rm{o}}}$, the bias has been corrected here for the intrinsic fraction since the same C-value is used in our Monte-Carlo simulations. Mahy et al. (2021) analyzed different threshold C-values using the Monte-Carlo method from Sana et al. (2013) to see whether different choices of C-value affect the final estimation for the intrinsic binary fraction or not. They found that the adopted threshold C does not affect the final results unless significant contamination exists by false-positive RV measurements.

5.3.2. Discussion

Our intrinsic binary fractions (red squares) for stars in groups T1 to T4 are shown in Figure 13, and the star mark is the result of Galactic O-type stars from Sana et al. (2012). The final result indicates that the intrinsic binary fractions decrease toward late-type stars.

Zoom In Zoom Out Reset image size

We are also interested in investigating the distribution of orbital period as a consequence of spectral type (on T_eff). We then plot the orbital period distribution in the left panel of Figure 14 for sample stars in groups T1 to T4 (shown as solid black line, green dashed line, black dashed line and black dotted line, respectively) as well as similar over-plotted distribution by adopting the results repeated from Sana et al. (2013) in blue and Dunstall et al. (2015) in red. We plotted each work's orbital period distribution with the upper and lower uncertainty ranges by adding and subtracting the uncertainty from the best-estimated values, respectively, in the right panel of Figure 14. It seems that the distribution of group T4 is more flat.

Zoom In Zoom Out Reset image size

Our κ (red squares) distributions are shown in Figure 15. The star mark is the result of Sana et al. (2013), and the triangle represents the result of Dunstall et al. (2015). All of them have large error bars, and no significant correlation is found.

Zoom In Zoom Out Reset image size