1 Introduction↩︎

Theoretical models suggest the Milky Way harbors billions of black holes (BHs) and neutron stars (NSs) [1], [2], driving decades-long searches for these objects. Historically, accreting BH X-ray binaries have been primarily discovered since the 1960s through their characteristic X-ray bursts, leading to the identification of \(\sim\)​30 accreting BH X-ray binaries [3]. NS searches have employed both X-ray bursts and radio pulse detections, yielding hundreds of detections via X-ray bursts [4] and thousands via radio pulses [5].

However, these classical methods sampled only a small fraction of the total BH and NS population. X-ray burst detections prefer accreting BHs/NSs, while radio pulse detections can only detect rapidly rotating NSs which emit electromagnetic pulse signals beaming towards us. Recent technological advances and data accumulation have spurred the development of novel search techniques. Microlensing remains the only currently feasible method for detecting isolated stellar-mass BHs, yielding several candidates [6] and one confirmed detection in the event OGLE-2011-BLG-0462 [7]–[9]. Future Galactic Bulge Time Domain Survey [10] using the Nancy Grace Roman Space Telescope, formerly called Wide Field Infrared Survey Telescope [11], is expected to discover many more such isolated BHs/NSs. With the detection of GW150914 [12], gravitational wave observations are dedicatedly revealing merging BH-BH and BH-NS binaries, and more than 300 events have been detected [13]. Large-scale spectroscopy surveys have amassed vast stellar spectral datasets, enabling the identification of non-accreting BHs and NSs in medium-to-short orbital period binaries through radial velocity (RV) monitoring [14]–[18]. Furthermore, photometry has also been adopted to search for dark massive companions (BHs/NSs) in ellipsoidal variables [19]–[22].

While discoveries from these novel methods have expanded the known populations of BHs and NSs, they have yet to significantly increase their total known numbers. Importantly, the third data release of the Gaia mission [23] offers a promising new avenue to overcome this limitation. Its initial data release included orbital solutions for 1.7\(\times\)​10\(^{5}\) binaries with full astrometric parameters. Utilizing these full astrometric solutions along with spectroscopic observations, BHs and NSs in AU-scale-orbit binaries are being revealed [24], [25]. In addition, several BH binaries associated with Gaia astrometric measurements have also been reported [26], [27]. However, only a small fraction of AU-scale-orbit binaries have such published data due to the stringent quality cuts, inevitably missing detections for many BHs and NSs in AU-scale-orbit systems. To address this, [28] enlarged their search to the extended Gaia DR3 dataset. They employed a forward-modeling framework, simulating Gaia observables for 21,028 red-giant branch (RGB) stars to target massive dark companions, yielding 556 RGB + BH candidates.

Assessing the risk of hidden hot companions is necessary for these candidates, given previous false alarms like J0521 and V723 Mon. Both systems were initially reported as RGB + BH binaries [29], [30]. However, later analyses revealed that they host stellar companions rather than BHs: ultraviolet (UV) features indicated a hot companion in J0521 [31], while for V723 Mon, spectral disentangling already demonstrated the presence of a stellar companion [32], with subsequent UV observations providing additional confirmation [33]. While UV excess provides a powerful diagnostic for identifying hot luminous stellar contaminants, kinematic signatures offer an independent probe of compact-object formation. Supernova explosions associated with the birth of BHs or NSs can impart substantial natal kicks to binary systems through asymmetric mass loss and explosion kinematics, leading to significantly heated binary orbits and high peculiar velocities (\(V_{\rm pec}\)) [34]–[37]. Kinematic information therefore provides a complementary avenue to distinguish genuine compact-object binaries from ordinary stellar systems.

In this work, we perform spectral energy distribution (SED) analysis for 1,328 sources from [28], along with Galactic kinematic analysis, as a reference for evaluating their reliability and priority selection of follow-up observation sources. The structure of this paper is organized as follows. In Section 2, we introduce the source selection for analysis in this work and corresponding data collection. In Section 3, we introduce our SED fitting method and use UV diagnostics to select sources. In Section 4, we analyze the Galactic kinematics of our sample. We summarize our results and make a discussion in Section 5.

2 source selection and data collection↩︎

Our initial sample comprises 3,773 sources with fit_companion_mass \(\geq\) 1.4 \(M_\odot\) and flag_quality == True from [28]. Due to the necessity of UV diagnostics to evaluate the risk of hidden hot companions, we refine sources by requiring available UV photometry. The Galaxy Evolution Explorer (GALEX) [38], [39] is uniquely useful for this purpose, being the only modern mission that performed a large-area, UV photometric imaging survey of the sky with high spatial resolution. Equipped with a 50-cm Ritchey–Chrétien telescope, GALEX can simultaneously observe in two broadband channels—the far-ultraviolet (FUV; 1,350–1,750 \(\mathrm{\AA}\)) and near-ultraviolet (NUV; 1,750–2,750 \(\mathrm{\AA}\)). GALEX General Releases 6 and 7 (GR6/GR7), which constitute the final public data products of the mission and are based on the GALEX Merged Catalog of Sources [40], contain 82,992,086 sources. We cross-match our sample with GALEX GR6/GR7 to acquire UV photometry, reducing the sample to 903 sources. We further retrieve the MCAT (AIS_\(\ast\)_mcat.fits) by python codes from astroquery.mast import Catalogs and Catalogs.query_region() [41], which provides UV photometry for another 645 sources.

To construct the broadband SED, we supplement photometry from existing Gaia G, BP and RP bands by retrieving the archive of The AAVSO Photometric All-Sky Survey (APASS) [42], Two Micron All-Sky Survey (2MASS) [43] and Wide-field Infrared Survey Explorer [44], to collect optical and infrared photometry. Photometric measurements with undefined or zero uncertainties are regarded as invalid measurements. Furthermore, for photometric measurements from APASS, we also regard photometric measurements with Uncertainty flag = 1 as invalid measurements. Sources retrieved in APASS may only have valid photometry for partial bands, and we only retain sources with at least three bands of valid photometry to ensure enough photometric data points. Finally, we only consider sources with parallax_over_err\(>\)5, leaving 1,349 sources.

Given the inherent degeneracies in SED fitting, incorporating priors on stellar parameters is essential. We adopt stellar parameters published by [45], who have developed a new pipeline, Gaia Net, for reprocessing Gaia XP spectra, to predict stellar parameters for 220 million stars released in Gaia DR3 and published their results on Vizier¹ [46], [47], due to the complete sample coverage and accessibility. 21 sources without uncertainties for all three stellar parameters have been ruled out in this process, leaving 1,328 sources for the following analysis.

3 SED fitting↩︎

We set an additional systematic photometric uncertainty of 0.03 mag as an uncertainty floor for each band, then we fit the SEDs with a single-star model. Since our intention is to conduct UV excess diagnostics, we do not include the GALEX photometry in SED fitting. The parameters for SED fitting include effective temperature \(T_\mathrm{eff}\), surface gravity \(\log g\), metallicity [Fe/H], stellar radius \(R\), distance \(D\), and \(V\)-band extinction \(A_{\rm V}\). We perform the SED fitting using an affine-invariant Markov Chain Monte Carlo sampler implemented in emcee [48] to explore the posterior distribution of the model parameters. We run 24 walkers for 5,000 steps, regarding the initial 2,500 steps of each walker as the burn-in and discarding them. The posterior probability of the fitting parameters \(P_\mathrm{rob}\) is denoted as: \[\ln P_{\rm rob} = \ln \mathcal{L} + \ln \mathcal{P},\] where \(\mathcal{L}\) is the likelihood and \(\mathcal{P}\) is prior. We adopt separable priors for all model parameters. Gaussian priors are applied to (\(T_{\rm eff}\), \(\log g\), [Fe/H]) and \(D\), based on Gaia Net estimates and the Gaia parallax, respectively, with widths set by the quoted uncertainties of these quantities. \(A_\mathrm{V}\) and \(R\) are assigned non-informative flat priors over their physically allowed ranges. The likelihood \(\mathcal{L}\) is denoted as: \[\ln L = -\frac{1}{2} \sum_\mathrm{band}^{} \left(\frac{F_{\rm obs,band}-F_{\rm mod,band}(\emph{T}_\mathrm{eff}, \log g, {\rm [Fe/H]}, \emph{R}, \emph{D}, \emph{A}_{\rm V})}{\sigma_{\rm obs}}\right)^{2}\,\] where \(F_{\rm obs,band}\) is the observed flux, \(F_{\rm mod,band}\) is the model flux, and \(\sigma_{\rm obs}\) is the error of the observed flux. For each set of fitting parameters, the model flux of each band is calculated by convolving the synthetic spectrum with the corresponding transmission curve provided by the virtual observatory SED analyzer² (VOSA) [49]. The synthetic spectrum is generated by pystellibs³ for a given set of (\(T_\mathrm{eff}\), \(\log g\), [Fe/H], \(R\)), using BTSettl library [50]–[52], then reddened according to the Fitzpatrick extinction law [53]. The total-to-selective extinction ratio \(R_\mathrm{V}\) is set to be 3.1 [54] and the reddening coefficients at different wavelengths for a given \(A_\mathrm{V}\) are calculated by extinction [55]. Finally, the synthetic spectrum is diluted by 4\(\pi D^2\). The zero-point used for converting magnitude to flux of each band is also provided by VOSA.

For each source, we adopt the median of the posterior distribution for each parameter as the best-fitting value and generate the corresponding synthetic spectrum. We then compute the \(\chi^2\) per data point (\(\chi^2\)/\(N\)) of the fit using all non-UV photometric bands. To ensure that the observed broadband photometry is adequately described by a single-star SED model, we retain only sources whose \(\chi^2\)/\(N\) \(<\) 2. As a result, 137 sources are excluded, yielding a cleaner sample containing 1,191 sources for subsequent investigations.

For the retained sources, we compute the model-predicted NUV fluxes and compare them with the observed values (see Figure 1 (a)) and show SED examples for a candidate with a UV excess and one without a UV excess, respectively, in Figure 2. We define the NUV flux ratio as \(R_{\rm NUV}=F_\mathrm{NUV,obs}/F_\mathrm{NUV,mod}\). 298 sources with \(R_{\rm NUV}\leq1.2\) are considered to show no significant NUV excess and are directly retained. We note that some sources have \(R_{\rm NUV}<1\), which is not physically expected and likely reflects measurement uncertainties or limitations of the model; no additional filtering is applied to these objects. For 585 sources with moderate excess, \(1.2<R_{\rm NUV}~\leq~3\), we perform further analysis to test whether the observed NUV emission can be explained by a non-degenerate companion or not.

To test this hypothesis, we model the companion in the moderate-excess systems as a main-sequence (MS) star. Stellar parameters for the companion are inferred from Isochrones⁴ [56], adopting the same [Fe/H] as the RGB in each system. The predicted parameters are used to generate synthetic spectra, which are then reddened using the SED-derived extinction and scaled by the distance to predict the observed fluxes and compute the corresponding NUV fluxes. We then compare the total model-predicted NUV fluxes of the RGB and the companion with the observed values (see Figure 1 (b)). We find a strong positive correlation between the predicted-to-observed NUV flux ratio and the fit_companion_mass, as expected if the companion were a normal luminous star. For a substantial fraction of sources, the total model flux exceeds the observed NUV flux by more than an order. Such large discrepancies are difficult to reconcile with measurement uncertainties or model limitations, indicating that a non-degenerate luminous companion is less possible. We therefore interpret 512 systems with \(F_\mathrm{NUV,tot\_mod}/F_\mathrm{NUV,obs}\geq10\) as more consistent with hosting dark companions and retain them as candidate BH/NS binaries. After applying both the SED quality selection and the UV consistency tests, 810 sources satisfy the SED diagnostics.

4 Galactic kinematic diagnostics↩︎

We characterize the kinematic properties of 1,191 sources, whose optical-to-infrared SEDs can be well fitted by single-star model, by deriving their \(V_{\rm pec}\) and maximum vertical heights above/below the Galactic plane (\(|Z|_\mathrm{max}\)) in the Galactic potential. We evaluate the \(V_{\rm pec}\) at the Galactic-plane crossing phase (\(Z=0\)). This choice provides a uniform kinematic reference across the sample, where the local circular velocity is well defined in the adopted Galactic potential model, and simultaneously reduces sensitivity to orbital-phase–dependent exchange between kinetic and gravitational potential energies. Systemic velocities needed for \(V_{\rm pec}\) calculation are taken from radial_velocity published by Gaia DR3, which are derived from multi-epoch RV observations and provided as a robust average value for each source. Therefore, the use of Gaia DR3 radial_velocity as systemic velocities is reasonable. We integrate the orbits backward for 1 Gyr in the Galactic potential under the Milky Way potential model McMillan17 [57], the most used Milky Way potential model, implemented by python package galpy [58], to determine the Galactic plane crossing phase and \(|Z|_\mathrm{max}\). Then we calculate the 3D velocities in Galactocentric Cartesian frame (\(V_{\rm x}\), \(V_{\rm y}\), and \(V_{\rm z}\)) at the Galactic plane crossing phase. In the calculation process, we place the Sun at (X, Z) = (-8.178, 0.025) kpc [59], [60]. At such place, the circular velocity is set to be 233.2 km s\(^{-1}\) according to McMillan17. \(V_{\rm pec}\) of the Sun relative to the Local Standard of Rest is set to be (\(U_\odot\), \(V_\odot\), \(W_\odot\)) = (7.01, 10.13, 4.95) km s\(^{-1}\) [61]. We convert \(V_{\rm x}\), \(V_{\rm y}\), and \(V_{\rm z}\) into 3D Galactic space velocity (\(U,~V,~W\)) by the following matrix transformation:

\[\begin{bmatrix} -\cos\alpha & -\sin\alpha & 0 \\ -\sin\alpha & \cos\alpha & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} V_{\rm x} \\ V_{\rm y} \\ V_{\rm z} \end{bmatrix} = \begin{bmatrix} U \\ V \\ W \end{bmatrix} ,\] where \(\alpha\) is the angle between the Galactic Center-to-source vector and the positive X-axis in Galactocentric Cartesian frame. From the resulting Galactic velocity components, we define the \(V_{\rm pec}\) as \(V_\mathrm{pec}=\sqrt{U^{2} + (V-V_{\rm c})^{2} + W^{2}}\), where \(V_{\rm c}\) is the local circular velocity at the disk-crossing radius.

[62] investigated \(V_{\rm pec}\) of compact-object binaries and major contaminating populations. By comparing their empirical cumulative distributions, they showed that the \(V_{\rm pec}\) of more than 90% of contaminating sources lie below 100 km s\(^{-1}\) (see their Figure 3). We therefore adopt \(V_{\rm pec}\geq100\) km s\(^{-1}\) to isolate systems probably accelerated by supernova natal kicks, yielding 227 sources. Among them, 182 sources also satisfy the SED-based selection criteria and 19 are BH candidates with fit_companion_mass \(\geq3\) \(M_\odot\). On the other hand, theoretical studies suggest that BHs may experience substantially reduced natal kicks due to strong fallback [63] or even form via direct collapse with negligible mass ejection [64], naturally leading to lower \(V_{\rm pec}\). Observations also support that some BHs must form with very weak kicks [65]. Therefore, 77 sources that satisfy the SED criterion alone and have fit_companion_mass \(\geq3\) \(M_\odot\) remain plausible BH binary candidates.

In Figure 3, we respectively compare \(V_{\rm pec}\) with fit_companion_mass, \(R_{\rm NUV}\), and \(|Z|_{\rm max}\), for sources with \(\chi^2\)/\(N\) \(<\) 2. According to the fit_companion_mass, we divide the systems into NS-like zone, transition zone, and BH-like zone, respectively corresponding to fit_companion_mass of \(\leq2\) \(M_\odot\), \(<2\) \(M_\odot\) and \(<3\) \(M_\odot\), and \(\geq3\) \(M_\odot\). These three zones contain 677, 322, and 134 sources, respectively. The NS-like zone extends to higher \(V_\mathrm{pec}\) tails than the BH-like zone; however, given the much larger sample size in the low-mass bin and possible selection effects, our data do not permit a quantitative comparison of natal kick distributions between the two populations. We further note a mild tendency for systems with higher \(V_{\rm pec}\) to exhibit lower NUV flux ratios, with the median \(R_{\rm NUV}\) shifting from \(\sim\)​2 at \(V_{\rm pec}< 100\) km s\(^{-1}\) to \(\sim\)​1.25 at higher \(V_{\rm pec}\). This behavior is consistent with kinematically heated systems preferentially hosting compact companions, while more significant NUV excess is more common among kinematically cold populations. Interestingly, the median \(R_{\rm NUV}\) shifting in the subsample with \(V_{\rm pec}\geq100\) km s\(^{-1}\), \(\sim\)​1.25, matches the threshold adopted in Section 3 for identifying sources showing no significant NUV excess, providing an independent consistency check between the UV and kinematic diagnostics. Furthermore, a clear positive correlation is observed between \(V_{\rm pec}\) and \(|Z|_{\rm max}\), indicating progressively heated orbits from disk-like to halo-like kinematics.

5 Summary and discussion↩︎

[28] employed a forward-modeling framework, simulating Gaia observables for 21,028 RGB stars to target compact object candidates. According to their results, we perform SED and Galactic kinematic diagnostics to provide additional vetting for 1,328 dormant BH/NS binary candidates and to prioritize targets for follow-up observations. A total of 810 sources satisfy the SED-based selection criteria. Their optical and infrared photometry can be well fitted by a single-star model and exhibit no UV excess or moderate UV excess that is completely insufficient to accommodate the existence of non-degenerate companions. Independent of the SED-based diagnostics, 227 sources meet the kinematic criterion with \(V_\mathrm{pec}\geq100\) km s\(^{-1}\). We identify the intersection of these two samples as the highest-priority candidates for follow-up observations. The highest-priority candidates contain 182 sources, of which 19 are BH candidates with fit_companion_mass \(\geq\) 3 \(M_\odot\). Besides, we designate 77 sources that satisfy the SED criterion alone and have fit_companion_mass \(\geq\) 3 \(M_\odot\) as a secondary-priority sample which only contains BH candidates. Both the highest-priority and secondary-priority candidates are listed in Table [tab:candidates].

Hierarchical systems may contaminate our sample, as their AU-scale orbits can accommodate inner binaries. In addition, for systems with moderate NUV excess but unlikely to host MS companions, contamination may arise from evolved non-degenerate companions in evolutionary phases where the NUV flux is significantly reduced. Based on Isochrones, we estimate that this phase occupies a non-negligible but limited fraction (typically \(\sim\)​14% to \(\sim\)​32%). Notably, such scenarios are only relevant to the subset of sources with moderate NUV excess, which accounts for \(\sim\)​63% of the systems that satisfy the SED diagnostics, suggesting that the overall impact of such contamination is further limited. For highest-priority candidates, the high \(V_{\rm pec}\) indicates a comparatively higher likelihood of having experienced supernova explosions. Consequently, the probability of such systems being contaminated by evolved non-degenerate companions becomes lower. Furthermore, even if such systems host inner binaries, they are expected to contain at least one compact object (BH or NS) rather than being composed purely of non-degenerate stars. In contrast, the secondary-priority sample, with lower \(V_{\rm pec}\), is more likely to be contaminated by hierarchical systems consisting of pure non-degenerate binaries and systems hosting evolved non-degenerate companions that have not undergone supernova events. This distinction reflects a fundamental physical difference between the two samples. Within the secondary-priority candidates, systems with higher fit_companion_mass or comparatively higher \(V_\mathrm{pec}\) are therefore comparatively more promising sources for follow-up observations. On the other hand, UV excess could also arise from the chromospheric activity of the RGB. Consequently, such systems may be removed by our SED diagnostic criterion. This is consistent with our goal of prioritizing sample purity over completeness.

References↩︎