Infrared Excess of a Large OB Star Sample

The infrared excess from OB stars are commonly considered as contributions from ionized stellar wind or circumstellar dust. With the newly published LAMOST-OB catalog and GOSSS data, this work steps further on understanding the infrared excess of OB stars. Based on a forward modeling approach comparing the spectral slope of observational Spectral Energy Distributions (SED) and photospheric models, 1147 stars are found to have infrared excess from 7818 stars with good-quality photometric data. After removing the objects in the sightline of dark clouds, 532 ($\sim7\%$) B-type stars and 118 ($\sim23\%$) O-type stars are identified to be true OB stars with circumstellar infrared excess emission. The ionized stellar wind model and the circumstellar dust model are adopted to explain the infrared excess, and Bayes Factors are computed to quantitatively compare the two. It is shown that the infrared excess can be accounted for by the stellar wind for about 65\% cases in which 33\% by free-free emission and 32\% by synchrotron radiation. Other 30\% sources could have and 4\% should have a dust component or other mechanisms to explain the sharply increase flux at $\lambda>10\mu$m. The parameters of dust model indicate a large-scale circumstellar halo structure which implies the origin of the dust from the birthplace of the OB stars. A statistical study suggests that the proportion with infrared excess in OB stars increases with stellar effective temperature and luminosity, and that there is no systematic change of the mechanism for infrared emission with stellar parameters.


INTRODUCTION
The infrared (IR) excess in early-type stars was firstly detected by Geisel (1970) who found the observed color index K(2.2 µm)-N (10.2 µm) being redder than that corresponding to the star's spectral type. Geisel (1970) explained the infrared excess in the early-type stars by dust grains formed in the mass loss process like in the evolved low-mass stars. Then Allen & Swings (1972) detected several forbidden emission lines from some of these stars and argued that the K − N excess comes from free-free emission by hot and ionized circumstellar gas instead of dust thermal emission. Hovhannessian & Hovhannessian (2001) studied 58 O, B, A, and F-type stars (including 45 OB stars) observed by the Infrared Astronomical Satellite (IRAS ). They explained the IR excess of 34 stars as the contribution by both blackbody emission from dust and free-free emission from ionized gas, and described the structure as 'gas-dust shell' or 'gas-dust disk'. Then, Siebenmorgen et al. (2018) identified twelve stars with IR excess from a sample of 22 OB stars with the Spitzer /IRS spectrum available (Houck et al. 2004) alongside the 2MASS (Skrutskie et al. 2006) and WISE (Wright et al. 2010;Cutri et al. 2013) photometric data. Similar to the arguments of Hovhannessian & Hovhannessian (2001), the observational results can be successfully explained either by free-free emission from ionized gas or dust thermal emission.
The debate on the mechanism of infrared excess in early type stars has been continuing since its detection. Theoretically, early-type (OB) stars have strong stellar winds which produce hot and dense ionized gas to bring about IR excess by free-free emission (Hartmann & Cassinelli 1977), which is a power-law continuum emission from IR to radio bands. Meanwhile, the IR excess can also be explained by optically thin dust emission. To discriminate the two mechanisms, a widewavelength-range spectral energy distribution to well define the profile can be helpful. Nevertheless, previous works agreed that both free-free emission from ionized gas and thermal blackbody emission from dust can explain the infrared excess. But we are not clear how much proportion of the IR excess can be explained by free-free emission or dust emission and how the mechanism depends on stellar parameters.
By its unique design of large field-of-view with four thousand of fibers, the Large Sky Area Multi-Object Fiber Spectroscopy Telescope (LAMOST; Cui et al. 2012) has acquired over 10 million stellar spectra in the Galaxy. This huge database brings the possibility to significantly expand the scale of OB star sample. From the LAMOST/LRS (Low Resolution Spectra), about 16,000 OB stars were identified by Liu et al. (2019) and from the LAMOST/MRS (Mid Resolution Spectra), other ∼9,000 OB stars were then identified by Guo et al. (2021a). For these LAMOST-OB stars, the basic stellar parameters including effective temperature (T eff ) and surface gravity (log g) were determined by Guo et al. (2021b) by the data-driven technique Stellar LAbel Machine (SLAM) with the non-LTE TLUSTY synthetic spectra as the training dataset. Moreover, the 1st catalog based on LRS is used by Deng et al. (2020) to determine the intrinsic color indexes of these OB stars and consequently an accurate measurement of the interstellar extinction to them, therefore the infrared excess can be calculated with high precision.
Based on this biggest sample of OB stars ever, this work tries to solve the following problems. How many of them have IR excess? What are their typical emission characteristics? How much do stellar wind and dust contribute? The data will be introduced in Section 2 followed by the method to detect IR excess in the objects in Section 3. To explain this phenomenon, the models of ionized stellar wind and circumstellar dust will be presented in Section 4. The results are shown in Section 5 and more details will be discussed in Section 6.

The OB star Sample
The preliminary catalog contains more than 20,000 OB stars with stellar parameters by Guo et al. (2021b) from the LAMOST survey (Cui et al. 2012) that is a reflective Schmidt telescope located at the Xinglong Station of the National Astronomical Observatory of China. In order to obtain more uniform properties from the statistical sample, and most of the parent samples are dwarfs, so only the dwarfs with log g ≥ 3.5 is included in our sample. Because this catalog lacks Otype stars, the Galactic O-Star Spectroscopic Survey (GOSSS, a project dedicated to O-stars;Maíz Apellániz et al. 2016) which identifies more than 1000 O-stars is supplemented.

Photometric Data
The photometric data to define the spectral energy distribution (SED) covers the optical-to-infrared wavebands. The optical photometry from the Gaia EDR3 (Gaia Collaboration et al. 2016; Gaia Collaboration 2020), Pan-STARRS1 DR2 (Chambers et al. 2016) and APASS DR9 (Henden et al. 2016) surveys were adopted. Gaia takes high-quality photometry in three passbands: G, G BP and G RP . Pan-STARRS1 uses a 1.8 m telescope located in Hawaii to observe in five bands: g, r, i, z and y, and APASS works in traditional B and V bands. In infrared, the data from 2MASS, Spitzer /SEIP (Spitzer Enhanced Imaging Products), Spitzer /GLIMPSE (Churchwell et al. 2009) and WISE surveys were adopted. 2MASS is an infrared full-sky survey in the J, H and K s bands. Both Spitzer /SEIP and Spitzer /GLIMPSE perform photometry in Spitzer /IRAC (Infrared Array Camera; Fazio et al. 2004) bands at 3.6, 4.5, 5.8, 8 µm respectively, and the 24 µm data from Spitzer /SEIP and Spitzer /MIPS (Multiband Imaging Photometer; Rieke et al. 2004) are searched as well. WISE, though not as sensitive as Spitzer /GLIMPSE, surveyed all the sky in the W 1, W 2, W 3 and W 4 bands.
In addition, the photometry at mid-and far-infrared wavelengths is searched from the observations by the MSX satellite (Price et al. 2001), which surveyed the entire Galactic plane within |b| ≤ 5 • in four mid-infrared spectral bands between 6 and 25 µm, by the IRAS satellite (Neugebauer et al. 1984) which surveyed more than 96% of the sky at 12, 25, 60 and 100 µm as the first infrared space telescope, by the AKARI satellite (Murakami et al. 2007) which covered more than 90% of the sky at 9, 18, 65, 90, 140 and 160 µm, and by the Herschel Space Observatory (Pilbratt et al. 2010) which mapped nearly 8% of the far-infrared up to sub-millimeter sky as the latest infrared space facility. The OB stars are cross-identified in the IRAS (Helou & Walker 1988), AKARI (Ishihara et al. 2010) and Herschel /PACS point source catalogs (Herschel Point Source Catalogue Working Group et al. 2020). This results in 24 objects with IR excess which is detected in at least one band among these long wavelengths, which will be exclusively discussed later. Consequently, the SED of most of the sample stars extends from optical to about 22-24 µm.
The cross-identification between catalogs are performed within a radius of 3 (extending to 5 for MSX, IRAS, AKARI and Herschel ), which is about three times the positional uncertainties of LAMOST. In cases where there is more than one object within the 3 radius, the nearest one is selected. The object is required to have the Gaia, 2MASS and WISE data available for a wide coverage of the SED. With these highsensitivity all-sky surveys, the final LAMOST sample contains 20,551 stars, i.e. about 95% of the preliminary sample, while the GOSSS sample is left with 589 sources, i.e. about 60%.

Extinction Correction
Prior to searching for infrared excess, the interstellar extinction and reddening are corrected to get stellar intrinsic SED. Since Gaia provides the best photometry quality, the intrinsic color index Gaia/G BP − G RP is calculated by its relation with T eff for early-type stars derived by Deng et al. (2020), which brings about an accurate determination of the color excess E(G BP − G RP ) that is proportional to absolute interstellar extinction. The interstellar extinction in each photometric band is then calculated by the extinction law. Specifically, this color excess is first converted to the classical reddening (Wang et al. 2018), and then to extinction in other bands according to the extinction law by Fitzpatrick & Massa (2007) in optical and by Xue et al. (2016) in infrared.
The interstellar extinction is checked by comparing E(B − V ) with the widely used Bayestar 3D extinction map by Green et al. (2019) at the Gaia distance of the object. For most of the objects, the two extinctions agree with each other very well. But they differ for about 10% stars by E(B − V ) > 0.3 mag, in particular for the stars with heavy reddening. This is because Green et al. (2019) take the average extinction in a range of area while our estimation dealt with the photospheric properties for individual star. Thus, the interstellar extinction is corrected according to the intrinsic color index by Deng et al. (2020). One may question that the circumstellar dust can also redden the star and be taken into the interstellar extinction. This may be true. But no correction is made to this point, because (1) we have no idea how much this might be, and (2) as will be shown later, this is a very small value on the order of 0.001 and negligible.
The photospheric emission should be calculated in order to extract the infrared excess. The stellar atmospheric model from the Kurucz ATLAS9 (Kurucz 1979;Castelli et al. 1997) and Tlusty (Hubeny & Lanz 1995;Lanz & Hubeny 2003, 2007 grids are both examined. Though the Tlusty model is more suitable for the massive stars such as OB stars, its grid begins from T eff = 15, 000 K. Instead, the ATLAS9 model includes the range of T eff ∈ [10000, 15000] K. In addition, no apparent difference appears in these two models. So, the ATLAS9 model is adopted for all the sources.
The output flux from the stellar atmospheric model is converted to the observed flux by a normalization factor C: where R refers to the radius of the stellar photosphere and D is the distance. For this purpose, only the optical bands of the Gaia, APASS and PS1 missions are adopted because the infrared emission may come from circumstellar matter in addition to the photosphere. In practice, the coefficient C is calculated for each band and the mean value of C is finally adopted. The dispersion of C from about 10 bands is typically about 10%. Together with the distance D measured by Gaia (Bailer-Jones et al. 2021), the stellar radius is derived by this mean C and will be used later. Besides of C, the typical fractional uncertainty of D is about 10%. Therefore, based on the error propagation, the fractional uncertainty of stellar radius is also about 10%.
Since T eff and log g are already given in the LAMOST-OB catalog, the closest model was simply matched to each star. As for stars from GOSSS, the best model was selected by the least-χ 2 method from all ATLAS9 grid for O-type stars (T eff > 32000 K) and a pair of T eff and log g is given. Figure 1 shows two examples, including the photometric brightness before (open circle) and after (filled circle) correcting for interstellar extinction, and the model spectrum with the stellar parameters from the LAMOST spectroscopy. Both the ATLAS9 and Tlusty models are displayed with no visible difference.

Identifying Infrared Excess by Forward Modeling
The aforementioned early detection of infrared excess in early type stars by Geisel (1970)    β β obs, ± σ β obs, ± 2σ β obs, ± 3σ β ATLAS9 Figure 1. Two typical spectral energy distributions for a star of T eff ∼ 20, 000 K without (left) or with IR excess (right). The photometric data from a few surveys are presented, where the open and filled circles denote the flux before and after correction for interstellar extinction respectively. The photometry error is plotted in grey shades behind the data points. The orange and brown lines show the ATLAS9 and Tlusty models respectively for the given effective temperature and surface gravity from the LAMOST survey. The dashed black and blue lines show the best-fit power-law (F obs − F model ) ∝ ν β for observational SED and ATLAS9 model, respectively. The likelihood distribution regarding to spectral index β is presented in the lower panels. The dashed lines demonstrate 68% (1σ), 95% (2σ) and 99.7% (3σ) confidence intervals and the solid blue lines represent the spectral slope index for the ATLAS9 photospheric data (2018) detected the infrared excess by the ratio between observed flux and the blackbody photosphere flux as FIRS−FBB FBB > 0.1 to find out the IR excess among their samples. These two methods are simple to use, but strongly affected by the uncertainties of specific photometry data adopted for IR excess identification.
To avoid such issue, instead of doing a simple comparison on flux, a power-law fit on the long wavelengths SED is carried out with a spectral slope index β (also see Figure 1): The IR excess could be detected by comparing the spectral slope index β obs measured from the observational SED with the β ATLAS9 from ATLAS9 photospheric model. Because the T eff of our sources ranges from 10,000 K to nearly 50,000 K, the spectral slope of their theoretical photosphere radiation varies from approximately β ATLAS9 ∼ 1.8 to ∼ 2.0, a universal threshold simply comparing these two β is not appropriate. Hence, a forward modeling method based on Bayesian statistic framework is adopted. The goodness of model fitting is calculated by the following likelihood function L = − exp ( χ 2 IR 2 ) which takes only the measurements at λ > 2.15 µm (2MASS/K s band) into account, since the emission at shorter wavelength comes from stellar photosphere. The χ 2 IR is taken as: Quality control of the observational data is firstly conducted for accurate photometry data. The photometric error is limited to be smaller than 0.03 mag, 0.05 mag, 0.05 mag and 0.1 mag in the 2MASS/K s , WISE /W 1, W 2, and WISE /W 3 band respectively, which keeps 7671 sources in the LAMOST sample and 147 sources in the GOSSS sample.
Then, the infrared excess detection is carried out in following steps: 1. For each star, the spectral index of its photospheric model β ATLAS9 is computed by leastsquare method from the lg F ATLAS9 ∝ β ATLAS9 × lg ν. Synthetic photometry on the ATLAS9 model is performed based on the filters of each passband to get F ATLAS9 , and no uncertainties are assumed.
2. A grid of β obs points ranging from −5.0 to 2.0 with a step size of 0.01 is created for finding the best one for observational data. For each of the β obs,t value, an additional simulated flux F add = C add ν β obs,t is added to the photospheric model to get the full modeling flux F model = F ATLAS9 + F add . To best demonstrate any potential infrared excess at long wavelengths, the constant C add is determined by forcing the F model to match the WISE /W 3 observation.
3. A distribution of likelihood function L(β obs ) is compiled by going through the β obs,t in the grid. The β obs,t that gives the highest likelihood value is the best-fit point.
5. If the spectral slope from ATLAS9 photospheric model is outside of the 99.7% confidence interval, i.e., β ATLAS9 > β obs,+3σ , this star is recognised as having infrared excess.
Top panels of Figure 1 show two stars as examples of the typical SEDs with or without IR excess. The best-fit spectral index on the observations and ATLAS9 are shown as the dashed black and blue lines respectively. The β obs and its likelihood distribution for these two examples are presented in the lower panels. The dashed lines demonstrate the boundaries of different confidence intervals and the solid blue lines show the spectral slope index for the ATLAS9 photospheric data. Smaller percentage for confidence intervals such as 68% (1σ) or 95% (2σ) could result in more detections of IR excess, but most of those additional samples are highly likely to be mis-identifications as the β ATLAS9 goes into the uncertainty range. Besides, the stellar parameters of these massive stars are with high uncertainty (err(T eff )/T eff ∼ 10%), resulting in some uncertainties of err(β ATLAS9 ) ∼ 0.05, but this uncertainty is not considered in this work. Therefore, as a safe choice, the upper bound of 99.7% confidence interval is adopted aiming for clear IR excess identification.
For the entire sample of 7818 stars, 1147 stars are identified with IR excess, giving a percentage of ∼ 15%.

Spectral Lines
LAMOST provides low-resolution (R1800) and medium-resolution (R7500) spectra, from which the spectral lines could be measured and analyzed. The H α line index was measured by integrating a continuumsubtracted flux within 12Å for the 1071 of 1088 stars with clear mid-IR excess and available spectrum from LAMOST sample. The continuum was subtracted by a second-order polynomial fit using 30Å of data on either side of the H α line. We found 238 of them (∼ 22%) have strong emission (H α line index > 10). Three kinds of profiles appear as single-peak emission, self-absorption in the center and center emission with wing-absorption, which are expected from a circumstellar disk. For stars that have multiple-epoch observations, though the time variation of line profile is obvious, no periodicity is visible, which may indicate the change is not due to periodic phenomena such as binary, rotation or pulsation. But the observations are not numerous enough, we will not investigate further based on the presently available data.

Association with Dark Clouds
Born in a dusty environment, many OB stars at their youth are still immersed in their birthplaces with significant amount of dust. It is possible that background sources bring about the IR excess. Though a singletemperature modified blackbody of f ν ∝ ν 2 B ν (T ) radiation cannot fit the SED, the contribution by some surrounding clouds causes additional infrared radiation. Dark clouds (DCs) are nearby members of the densest and coldest phase in the Galactic interstellar medium, and represent the sites where stars are currently being born. Early-type massive stars are young and likely to be associated with those dark clouds with gas temperatures > 10 K, or the so-called IR dark clouds (Bergin & Tafalla 2007). Both the spectrum and photometry of these stars in the dark clouds are very likely to be affected, and the observed infrared excess may then come from the cloud instead of the circumstellar matter.
In order to exclude the objects associated with dark clouds, the above sources that show infrared excess are cross-matched with the catalog of dark clouds from the Atlas and Catalog of Dark Clouds 1 (Dobashi et al. 2005;Dobashi 2011) based on the optical Digitized Sky Survey (DSS) and infrared 2MASS images. It is found that 463 stars from LAMOST and 34 stars from GOSSS are in the sightlines of dark clouds whose infrared excess is very likely to be caused by the radiation of dust in the cloud. They are excluded in further analysis of the mechanism for infrared emission. After removing these stars from the preliminary sample, 625 stars from LAMOST and 25 GOSSS stars (650 stars in total) are left for further study.

Stellar Wind Model
As massive stars, OB stars normally blow out strong hot ionized stellar wind. Seaquist & Gregory (1973) pointed out that for the isotropic spherical condition, the flux distribution of the free-free emission of electrons in an ionized gas which is generally thinner in the outer region following the power law F ν ∝ ν α , where spectral index α ∈ [−0.1, 2]. Theoretically, α = 2 corresponds to the optically thick case that simulates the blackbody radiation, while α = −0.1 refers to extremely optically thin case. Wright & Barlow (1975) derived that ionized stellar wind under spherical isotropic isothermal expansion, will have a flux of F ν ∝ ν 0.6 . This conclusion was then adopted to describe the radiation property in infrared to radio bands by many works such as Crowther (2007) and Fogerty et al. (2016). Furthermore, Barlow (1979) summarized the observational results of free-free emission from early-type stars as F ν ∝ ν 0.7 , which was adopted by Siebenmorgen et al. (2018) to explain the IR excess of OB stars. Apparently, the results of Wright & Barlow (1975) and Barlow (1979) are highly consistent, with the index α ∼ 0.6 − 0.7, implying a rather optically thin free-free emission.
In our model, the power law F ν = C F ν α is adopted, where the spectral index α and scaling constant C F are varied to fit the observational IR excess of each star. The upper limit of α is set to 2.0, coincide with both the Rayleigh-Jeans approximation of the high-temperature photospheric radiation and the free-free emission of the electrons. On the other hand, the lower limit of α is free. Though the lower limit of α is −0.1 for free-free emission as mentioned above, synchrotron radiation, which also follows a power law F ν ∝ ν α but with a much larger negative index, i.e. α < −0.1, is a possible source of infrared excess. Shchekinov & Sobolev (2004) argue that the interaction of stellar wind with the surface of a circumstellar accretion (or protoplanetary) disk around massive stars can result in the acceleration of relativistic electrons in an external layer of the disk and produce synchrotron radiation. Leaving the lower limit of α free opens the possibility of identifying synchrotron radiation as the source of infrared excess. While the modelling simply takes the power law, the free-free and synchrotron radiation will be discriminated by the power law index yielded from the modelling.
Similar to the approach in Section 3.2, a grid of α and C F is created to compile the likelihood distribution of L(α, C F ). The spectral index α is ranging from -10.0 to 2.0 with a step of 0.1, and 10 C F grid points with uniform interval sampled between the value matching the lowestflux photometry and the highest-flux at the wavelengths > 2.15 µm. Then, the best-fit parameters are identified with the highest likelihood value and the uncertainties for both are determined from 68% confidence intervals.

Circumstellar Dust Model
Other than the stellar wind, thermal radiation from circumstellar dust could also be responsible for the IR excess around massive stars.
The code DUSTY (Ivezic & Elitzur 1997;Ivezic et al. 1999) is adopted to analyze the properties of circumstellar dust. DUSTY solves the radiative transfer problem in a dusty environment and offers many options for input radiation, dust types and density distributions. Because for most of the stars in our sample, they are lack of long wavelengths observations, thus hard to constrain their dust parameters. A simple Bayesian approach considering a prior distribution is adopted.
There are 27 stars that contain long wavelengths observations at λ ≥ 60 µm from IRAS, AKARI or Herschel, standing out from most of our samples that the observation only reaches to the wavelengths of W 3 band of 12 µm. These stars are chosen as the training sample for dust model fitting to obtain the priors of dust parameters. Using the training sample, the parameters of the DUSTY models library are chosen to cover the reasonable range for OB stars as following. First, the isotropic central radiation of a single heat source is adopted, specifically the input stellar spectrum is taken from the ATLAS9 model at given stellar parameters. Some dust-related parameters are fixed as well: the size distribution follows the MRN power law i.e. n(a) = a −q for a ∈ [0.005, 0.25] µm and q = 3.5; the upper limit of the dust temperature is set to be 1500 K; the outer radius of this dust shell is set to be 10 3 times of the inner radius. The grids of models are built within a range of parameters. For the temperature of inner dust shell (T d,inner , hereafter T d ), 29 equally spaced points with an interval of 50 K are sampled from 100 K to 1,500 K. The optical depth at 550 nm (τ dust V , hereafter τ ) is explored at 33 equal logarithmic interval from 10 −7 to 10. Various combinations of chemical composition (e.g. silicate, amorphous carbon and graphite) and dust density distributions (e.g. inverse square attenuation with ra-dius and the AGB stellar wind density model) are also tested. The experimental running of DUSTY found that the optical depth of dust τ is very low that the discrimination of dust species is meaningless. Thus, the dust composition is fixed as a mixture of 53% silicate and 47% graphite, i.e. the average interstellar dust composition from Draine & Lee (1984). This option presumes that the dust around these OB stars comes from their birthplace -molecular clouds, which is evidenced by the sub-parsec-scale dust structure (∼ 0.1 pc) since stellar wind can hardly reach such a distance. Correspondingly, the dust density distribution is set to be constant. All physical parameters adopted are shown in Table 1. Also from this training sample, a simple prior probability distribution of τ and T d is set up, as shown in Figure 2.
From this training sample, a wind-dust model is established by combining the DUSTY model and the stellar wind model of F ν = C F ν α in following steps: 1. Compute a grid of dust models from DUSTY. For each combination of T d and τ from dust parameter grid described in Table 1, dust thermal radiation is first calculated for each star as F dust .
2. For each source, compute the residual for each combination of dust parameters (T d , τ ) and its photosphere, then fit the residual with a powerlaw by a simple least-squared method: 3. Compile a distribution of likelihoods L test (T d , τ ) by going through all the T d and τ in the grid.
4. Multiply the likelihoods with the prior found by training samples to get the posterior: L post (T d , τ ) = L test (T d , τ )P prior . From this posterior distribution, the values (T d , τ ) given highest probability is chosen as the best-fit parameters, and the uncertainties are also presented accordingly as 68% confidence intervals.
As the same as the above method, only photometry data with wavelength > 2.15 µm are considered in computing L test (T d , τ ) to exaggerate the differences among the grids. Differently, the spectral index in this winddust model is limited within α ∈ [−0.1, 2] to represent free-free emission from the ionized stellar wind, while synchrotron radiation is not taken into account. Though the stellar wind component is added to all stars, it is unnecessary for some of them for which C F is very small. The post likelihood distribution and the chosen fit parameters are highly sensitive to the prior distribution based on the training sample of 27 stars. But due to the lack of long wavelengths observations on those OB stars in the total sample, this is currently the best estimation that could be chosen on wind-dust model.

Model Comparison
To quantitatively understand which model explains observations better, Bayes Factors (BF, Jeffreys 1961;Kass & Raftery 1995) are computed between two models for each star: where M i refers to two different models and P is the marginalized probability. There are two grids as described above for each star: (a) stellar wind model with likelihood distribution L(α, C F ), (b) wind-dust model with L(T d , τ ). From these grids, the marginalized posterior probability P (M i | data) is computed by: where λ 1 and λ 2 refer to the two free parameters in two models, and L post is the posterior probability. The computed BF value can tell which model is better quantitatively between the two models.

RESULTS
Previous works usually found that nearly 50% of the sample OB stars show IR excess. Excluding the stars in the dark clouds sightlines, there are 532/5634 (∼ 9.4%) B-type and 118/301 (∼ 39.2%) O-type stars with true infrared excess. This apparently lower percentage of OB stars with IR excess may be attributed to the fact that our OB star sample is not biased to any objects potentially having IR excess, instead it mainly comes from the LAMOST survey in optical. Moreover, the infrared excess that may come from surrounding medium is exclusively removed. If taking the ratio of observed flux in the infrared band to the photosphere flux as Siebenmorgen et al. (2018) did, a larger proportion of IR excess would be obtained. However, the additional stars usually have marginal IR excess. Our results may represent the portion with IR excess among normal OB stars.
The Bayes Factor (BF, described in Section 4.3) is computed for each star. In this work, the stellar wind model is called M 1 while the wind-dust model called M 2 . When the BF > 100 (or < 0.01), which means M 1 model is 100 times more (or less) likely to explain the observations than the M 2 (Kass & Raftery 1995), the M 1 (M 2 ) would be assigned as the model that best explains the data, and when 0.01 ≤ BF ≤ 100, both are appropriate models. The BF distribution for all the stars are presented in Figure 3. The prior probability distribution of dust parameters τ and T d based on the 27 training samples with long wavelengths observations, so that the dust parameters can be better constrained. The columns are the parameters fitted best for individual star, while the red columns refer to the one off the sightlines of dark clouds. The black lines represent the selected prior probability distribution later used for analyzing the whole sample. Since the stellar wind component is also in the winddust model, ideally the likelihood distribution should be complied from L(T d , τ, α, C F ) with four variables together. However, it is not practical as the computation time for each star with only two variables (T d , τ ) is already ∼ 1 minute. If a typical grid with ∼ 100 points for (α, C F ) were added to the whole grids, the whole process would be too computational expensive. A Markov chain Monte Carlo (MCMC) simulation would be helpful, but it is not necessary here. Because the best stellar wind component (α best , C F,best ) is already included in the wind-dust model for each combination of the dust parameters (T d , τ ), the marginalized probability computed for wind-dust model here is larger than the marginalized probability from the comprehensive one: With an overestimated marginalized probability for wind-dust model, the comparison would be inclined to it. Even though, as it shows below, the overall comparison tells that only dozens of stars are clearly better explained by wind-dust model, while stellar wind model is still a better choice or at least as good as the wind-dust model for most of the stars in our sample. Thus, the adopted simpler grid with only two variables, L(T d , τ ), would not give a different result comparing with the comprehensive one.
For 532 B-type and 118 O-type stars off the dark clouds sightlines and with IR excess, the stellar wind model with theoretical predicted free-free emission (α ∈ [−0.1, 2]) can satisfactorily explain the infrared excess of 178 (∼ 19%) and 39 (∼ 19%) stars, respectively. In addition, there are 169 (∼ 18%) B-type and 40 (∼ 19%) O-type stars, respectively for which steep increase at λ > 10 µm can be fitted by a single power-law radiation, but requiring a larger negative spectral index ure 4. Such large negative index is generally an indicator of synchrotron radiation. Figure 5 shows the distribution of the spectral index α fitted in the stellar wind model for all the stars with IR excess and off dark clouds. There are clearly two groups on either side of α = −0.1, suggesting there are two separate mechanisms. In addition, there are 10 stars fitted by α = −10 which is on the lower edge of the grid and the upper bound of 68% confidence interval could even reach to α +σ ∼ 1.3, giving a huge uncertainty of ∆α > 10. It is because there is only one photometry data point (WISE /W 3) showing obvious IR excess in their SED. Hence, the constraint on the spectral index α is so weak that it is not even possible to pin down the best value.
It should be mentioned that when a star can be explained well by the power-law model no matter the value of the spectral index α, it is possible that the IR excess comes from an exceptionally optically thin dust shell. There are a clear group suggesting such trend in the second row of Figure 5, where both models fit the observational SED as good as each other. Moreover, a dust component of τ < 10 −7 , which is even smaller than the lower limit in our grid could even fit the stars that are currently fitted better by stellar wind model. A practical problem is that the observations available to constrain the dust model are scarce and lacking at long wavelengths, and an obligatory dust model would be unreliable. In addition, for those with α ∈ [−0.1, 2], the free-free emission from stellar wind model already works quite well, and it is theoretically plausible. Hence, the free-free emission from ionized stellar wind is set to be the best model for those stars with α ∈ [−0.1, 2] and BF > 100. From them, the mean spectral index in F ν ∝ ν α can be concluded as α = 1.603 with a standard deviation of σ α = 0.492 for B-type stars, and α = 1.857 with a standard deviation of σ α = 0.313 for O-type. This mean index for B-type stars of 1.603 is higher to the expected value of ∼ 0.7 from Wright & Barlow (1975) and Barlow (1979), because the IR excess identification methodology adopted here is so sensitive that even stars with very weak IR excess are identified.
There are 216 stars that could be explained by the wind-dust model, in which 23 of them is fitted better by it. Figure 6 displays the results of best-fit T d,inner and τ for those stars. The distributions of T d,inner and τ from both samples are similar. As mentioned in Section 4.2 about the training samples, most of OB star sample are fitted by very low T d,inner (∼ 500 K) with small τ (∼ 10 −5 ). Most of the stars in our sample are lack of long wavelengths observations, and there is a very high degree of degeneracy among all the dust parameters, and it is impossible to choose which parameter is better by individual star. For example, the second star shown in Figure 4 only contains data points that reach to W 3 and W 4 bands and its IR excess only can be seen in these two bands. Before applying the prior probability distribution, there is a very high degree of degeneracy in the 'test' runs, i.e., there are multiple maxima in the likelihood function in the parameter space (left panel of Figure 7). The overall likelihood function tends to bias to the parameters on the edge of grids for many of stars under this circumstance. With this simple Bayesian method, the parameter degeneracy is eliminated and the optimal parameters can be found in the posterior within the parameter space we set up. Though the final fitted dust parameters are highly sensitive to the prior probability distribution based on the training sample, these are the best parameters we could estimate for wind-dust model due to the lack of the observational data at long wavelengths. This incompleteness of wind-dust model grids for some special stars can also be noticed by the 8 stars that couldn't be fitted by either models. For those targets, a single power-law radiation is not appropriate for explaining their IR excess, and the dust component needed is outside of the current grids. More variables, such as another set of dust compositions, are needed, but it is beyond the scope of this study.
The results of fitting as well as stellar parameters are presented in Tables 2 and 3   GOSSS sample respectively. Because of the asymmetry nature of fitted parameters, the lower and upper bounds of the 68% (1 σ) confidence intervals, λ −σ and λ +σ , are adopted to show their uncertainties. Only the stars off the sightlines of dark clouds are presented. The 'Best' column shows which model is the best: 'S' or 'F' is assigned when the stellar wind model with synchrotron radiation or free-free emission is better fitted as BF > 100; 'B' is for the case both models work and couldn't be distinguished for which is better as 0.01 ≤ BF ≤ 100; 'D' is used when more likely a wind-dust model is needed as BF < 0.01; 'U' is used for the stars that none of the model could explain their SED; and 'X' is marked when the SED is better fitted by synchrotron radiation in stellar wind but its spectral index is fitted at the lower edge of the grid α = −10.0. There is also an online version of these tables including all of the fitting results. Meanwhile, Table 4 shows the samples of various classes.

Be-stars
Be-stars are non-supergiant B-type stars with at least one Balmer line emission. For classical Be-stars, which are fast-rotating main-sequence B-type stars, the outflowing material forms a gaseous dust-free Keplerian circumstellar disk (Rivinius et al. 2013). The IR excess is often observed on them and it is mostly due to free-free emission from that disk. For this Be phenomenon, its emission is most likely to be a power-law with a spectral index α ∈ [0.6, 2] (Klement et al. 2017). Carciofi & Bjorkman (2006) predicted that for fully ionized gaseous disks near Be-stars, the IR excess should show up at wavelengths ≤ 1 µm for disks inclinations ≤ 60 • and at ∼ 10 µm for edge-on condition. Rivinius et al. (2013) summarizes the H α profile under the present classical Be-star model with a gaseous disk. The double peaks around the core correspond to the edge-on case. The single peak represents the pole-on condition, in which the wing-absorption might exist due to the strong stellar photospheric absorption. All these line profiles are observed in our LAMOST spectra. Similar to samples from Siebenmorgen et al. (2018), there are many stars in our sample that contain far-IR emission which only appears at wavelengths > 10 µm. Most of them also have double-peak H α emission line, which well support the edge-on condition of the model for classical Be-star with gas disk. Among the 532 stars with IR excess from the LAMOST-OB catalog off dark clouds sightlines and which are identified as B-type stars, 153 sources exhibit a clear H α emission line (∼ 29%) indicating that they are Be-stars. For 109 sources of them, the SED could be best fitted by a stellar wind model (BF > 100) and free-free emission works for 57 of them. It leads to a relatively higher proportion (109/153 ∼ 71%) of successful wind model than for the entire sample of B-type stars ( 178+169 532 ∼ 65%), which is consistent with the scenario that a Be-star usually has an ionized gaseous disk.

Possibility of Synchrotron Radiation
For the stars whose SED can be fitted by a powerlaw but with α < −0.1, the infrared excess can hardly be attributed to free-free emission which in general decreases with wavelength. On the other hand, such an SED that increases with wavelength resembles the synchrotron radiation. Although some early-type in particular B-type stars are found to have strong magnetic field, there must be some mechanism to obtain relativistic electrons if synchrotron radiation is to occur. White (1985) proved that electrons can be accelerated to relativistic energies by chaotic stellar winds in hot stars, and Shchekinov & Sobolev (2004) argued that the interaction of stellar wind with the surface of a circumstellar disk can result in the acceleration of relativistic radiation.
27 objects (the training sample in dust model discussed in Section 4.2) are detected at relatively longer wavelengths by either IRAS, AKARI or Herschel with observations at λ ≥ 60 µm, and 10 of them are out of the dark clouds sightlines. For all these 27 stars, their infrared excess can be fitted by a steep power law (α < −0.1) except one with α ∼ 0.2, and several with a plausible turn-over point at ∼ 100 µm (see Figure 8). Besides, the search for a counterpart in the NVSS catalog (Condon et al. 1998) within 5 for each star results in only 3 stars in spite that the extrapolation of the power law predicts an intensity at 21 cm greatly higher than the sensitivity of the NVSS survey. In addition, for those 3 stars, the NVSS observed flux is too low comparing with the predictions from synchrotron radiation starting from IR wavelengths (more than 3 orders of magnitude lower in the unit of Jy), so that their radio emission is more likely to be originated from circumstel-  Note-'in DC': inside the sightline of Dark Clouds. 'off DC': off the sightline of Dark Clouds. 'Both': both stellar wind and wind-dust models could explain the observations. 'Unknown': none of the two could explain the observations. lar dust emission or other mechanisms. In the total sample of 650 stars, the IR excess of 275 (24%) sources can be fitted by the synchrotron radiation and 209 (18%) of then even favor this model. But it is difficult to explain such large proportion deserving synchrotron radiation. Since the dust in dark cloud is usually cold and emits in far-infrared, the infrared excess originates very possibly from the dark cloud, and the variation of the SED within the studied wavelength range is caused by the difference in temperature. Alternatively, the wind-dust model is a better explanation for many of them. Thus, the final solution lies in more brightness measurements at longer wavelengths which can help distinguish these models.

Possibility of Debris Disk
Disk is one of the ubiquitous dynamic structures in the universe. As pre-main-sequence stars, Herbig Be stars with protoplanetary disks have now been studied extensively. Dusty debris disks might also exist around main-sequence massive stars as they evolved. Roberge & Weinberger (2008) presented SED for 16 nearby main-sequence massive stars including one Bestar (HD142926), one early F-type and 14 A-type stars. Both of the Spitzer /MIPS/24 and 70 µm photometry are included in their observations. It is found that the mid-IR excess in that Be star does not like the classical Be stars, whose fluxes should be a power-law decline with increasing wavelength. A debris disk model with blackbody dust grains or 1 µm silicate grains can both fit the observational data but with different parameters such as dust temperature, indicating a possible existence of the dusty debris disk around it. For a cold dust disk model, its SED also shows a continuous increase in fluxes from 10 to 100 µm, which is very similar to our current IR observation.
It is impossible to distinguish the disk structure from the spherical circumstellar dust described in this work (Section 4.2) until a directly resolved image is acquired. Siebenmorgen et al. (2018) performed near-IR highcontrast imaging of three O-type stars with far-IR excess in their sample, and they didn't find any significant disk structures except the stellar halo, which might be the scatter-light from disks. Also, it is basically impossible to reach a reasonable constraint on the disk model due to the lack of observation data at longer wavelength, and there is no direct evidence that debris disks exist around main-sequence OB stars. Therefore, further research on the disk model is not conducted in this work. Normally, the strong stellar wind of the star is accompanied by violent material ejection, which then brings considerable circumstellar matter. However, the harsh environment near OB stars brought challenges to the survival of circumstellar dust. The strong stellar wind  blows away dust, and the Poynting-Robertson drag also causes dust near inner radius to lose angular momentum and fall into the stellar atmosphere (Draine 2011).

Dust Source
From the grids generated by DUSTY (Section 4.2), a parsec-scale dusty sphere with a very small optical depth (τ ≤ 10 −5 ) and rather low dust temperature (T d,inner ∼ 500 K) can best explain the observational data. A possible condition should be: as an OB star born in its molecular cloud, the strong stellar wind quickly blows away the surrounding cloud and a huge structure was constructed, making a parsec-scale dusty envelope. That is to say, this envelope is the molecular cloud blown larger by the stellar wind. Thus, it has the same chemical composition and dust density distribution as this cloud, which are both represented in our wind-dust model (see Section 4.2). Far from the central star, the temperature is low, and the optical depth is small, which looks like a dusty halo as suggested by Siebenmorgen et al. (2018). This dusty circumstellar halo, together with the photosphere and stellar wind radiation inside, can interpret the observational SED.
This circumstellar halo strucutre is similar to the scenario of a Young Stellar Objects (YSO) growing in a dark molecular cloud. Molinari et al. (2008) described how the SEDs of massive YSOs evolve in the star forming region. They adopted both DUSTY code and 3dimentional model by Whitney et al. (2003) to estimate the SED and fit the observational data points focusing on wavelengths from ∼ 10 µm to > 1000 µm, c.f., Figures in Appendix A in Molinari et al. (2008). Although their objects are massive YSO, their observational SEDs are very similar to our samples of main-sequence OBtype stars requiring dust components: there is a sharp increase of IR excess to 100 µm. Their fitted DUSTY parameters are also similar to our results in the winddust model. Besides, shown in Figure 8, there are turnover points at 100 µm for these two stars staying in dark clouds sightlines, which makes it even closer to the YSO in dark cloud scenario. Therefore, the dust components in the IR excess actually reflects the conditions of the birthplaces of these massive stars: as they grow extremely fast, the dark clouds as the birthplaces still surround them when they are already in the main-sequence phase. For those stay out of the currently identified dark clouds sightlines, they might stay in a very small or faint dark cloud that has not been identified.

Dependence of infrared excess on stellar parameters
The LAMOST sample is taken to investigate the dependence of IR excess on stellar parameters, for which the stellar parameters are available. The GOSSS sample provides no stellar parameters and is dropped off. The stellar luminosity is derived with the effective temperature from the LAMOST spectroscopy and the radius calculated by Equation 1. Given that the errors of the effective temperature and the radius are both 10%, the luminosity is with ∼ 45% uncertainty. Figure 9 shows the distribution of stars with IR excess and best models in the Hertzsprung-Russell diagram (HR diagram). The percentage with IR excess increases apparently with stellar effective temperature (T eff ), and the luminosity since the objects are all dwarfs. For Btype stars , when T eff > 20, 000 K, more than 25% of them presents IR excess, while this percentage reduces to < 5% at T eff ∼ 10, 000 K. This trend agrees with the origin of the IR excess that is mainly stellar wind since luminous OB stars generally have strong stellar wind. For the mechanisms explaining the IR excess (right panels of Figure 9), all three cases generally stay at the same level at different T eff , except for slightly changing trends T eff (K) Figure 9.
Dependence of the IR excess and its model on stellar effective temperature and luminosity. The distribution of the stars with (pink dots) and without IR excess (grey dots) in the HR diagram are shown in the upper-left panel, and the stars with IR excess off the dark clouds are shown in the red dots. The change of the proportion with infrared excess is displayed in the lower-left panel. The right column is for the SED that is best fitted by free-free emission (light blue dots), synchrotron radiation (dark blue dots) of stellar wind model, or could be fitted by wind-dust model (orange dots). The T eff = 32, 000 K is shown in black lines to divide the O-type and B-type stars in our sample.
in lower temperature end. Free-free emission from ionized stellar wind and synchrotron radiation are more common for stars with higher T eff , rising to ∼ 40% at T eff ∼ 20, 000 K. Meanwhile, the dusty structures are more likely to survive around stars with lower temperatures, thus the proportion of sources requiring dust component reaches to ∼ 40% on the low temperature end with a small increase. However, since the variations in both trends are just above the Poisson uncertainties ( √ N bin N total ∼ 5% in ratio) in each bin, those two trends are questionable. Thus, we conclude that there is no obvious systematic change of the mechanism for infrared emission with stellar parameters.

SUMMARY
The infrared excess of OB stars is systematically studied based on the largest OB star catalog with stellar parameters from the LAMOST survey and the GOSSS O-type star sample. After a precise extinction-correction with the intrinsic color indexes from our previous work (Deng et al. 2020), IR excess are identified by comparing their spectral index in the infrared SED with photospheric model in a forward modeling approach. It is found that 939 and 208 stars show infrared excess among the 7303 B-type and 515 O-type stars, respectively. To better analyze the circumstellar condition of them, 407 B-type and 90 O-type stars in the dark clouds sightlines are eliminated, which leaves 532 (∼ 7%) and 118 (∼ 23%) stars respectively with true circumstellar infrared excess. Afterwards, the observational SED from optical bands (Gaia, PS1 and APASS) to infrared (2MASS, MSX, Spitzer and WISE ) is interpreted by synchrotron radiation or free-free emission in stellar wind or together with dust thermal radiation. Bayes Factors (BF) are computed for both models to quantitatively compare which one fitting observations better. The IR excess in one third of the OB-type stars (∼ 33%) can be better explained by free-free emission in ionized stellar wind, another one third (∼ 32%) of them is explained better by synchrotron radiation, both models work well for other 193 stars (∼ 30%), while winddust model works better for only 23 sources (∼ 4%) and 8 stars (∼ 1%) couldn't find a proper model. For those objects that wind-dust model could fit the observations (∼ 34% of the total sample), a parsec-scale dusty envelope with a low dust temperature and exceedingly small optical depth is identified, which implies a largescale circumstellar dust halo possibly originated from the birthplace cloud.