Open Access
Issue
A&A
Volume 670, February 2023
Article Number A79
Number of page(s) 20
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/202244818
Published online 07 February 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1. Introduction

Double massive main-sequence (MS) binaries (OB+OB) are expected to play an important role in the formation of compact object mergers of neutron stars (NSs) and black holes (BHs, e.g. Mandel & Farmer 2022). A particularly important intermediate evolutionary stage before a merger is a single-degenerate phase (either OB+NS or OB+BH), which might be identified through the emission of X-rays (see Walter et al. 2015; Corral-Santana et al. 2016, for identified BHs), called a high-mass X-ray binary (HMXB). This occurs when the OB star has evolved, either almost filling its Roche lobe, leading to a wind-fed X-ray binary, or filling its Roche lobe after which copious mass is being accreted onto the BH. This phase in the evolution is, however, very short-lived compared to the OB+BH phase, which is determined by the MS lifetime of the OB star (Sen et al. 2021; Hirai & Mandel 2021). On top of that, OB+BH systems will not always evolve into HMXBs. For example, an OB+BH with P >  10 d would never become a HMXB.

The search for X-ray-quiescent OB+BH systems is proving to be highly challenging. In the past two years, reported BHs in spectroscopic studies have been contested by follow-up studies, often leading to explanations that do not involve BHs (see, for example, some discussions on LB-1: Liu et al. 2019, 2020; Abdul-Masih et al. 2020; El-Badry & Quataert 2020; Shenar et al. 2020; Irrgang et al. 2020; Lennon et al. 2021; on HR6819: Bodensteiner et al. 2020; Mazeh & Faigler 2020; Rivinius et al. 2020; Romagnolo et al. 2022; Frost et al. 2022; on NGC 1850: El-Badry & Burdge 2022; Saracino et al. 2022; Stevance et al. 2022 to name a few). Nonetheless, very recent spectroscopic investigations have revealed a handful of such OB+BH systems (e.g. Mahy et al. 2022; Shenar et al. 2022). However, their numbers remain small since these dedicated spectroscopic studies require long time series and high signal-to-noise spectra.

The Gaia space mission (Gaia Collaboration 2016) has great potential for the detection of single-degenerate binaries (e.g. Breivik et al. 2017, 2019; Mashian & Loeb 2017; Yalinewich et al. 2018; Yamaguchi et al. 2018; Andrews et al. 2019; Wiktorowicz et al. 2019; Shikauchi et al. 2020, 2022). A method for astrometrically identifying such single-degenerate binaries was developed by Shahaf et al. (2019, hereafter SMF19).

Using the results of detailed binary evolution simulations of Langer et al. (2020), Janssens et al. (2022, hereafter Paper I) showed that the identification method presented in SMF19 could be adapted to identify a predicted 100–200 OB+BH systems among the sources in the second Alma Luminous Star (ALS II) catalogue (Pantaleoni González et al. 2021). The authors further show that the distributions of the parameters of detected OB+BHs are expected to retain information on the BH-formation physics (the kick and the fallback mass).

In Paper I, the authors only considered the impact of unevolved non-degenerate OB+MS binaries on the number of false-positive detections – systems that do not contain a BH that are misidentified as having a BH. However, many OB stars are found in triples (Sana et al. 2014; Moe & Di Stefano 2017), which may also contribute to the amount of false positives. On top of that, the ALS II catalogue is magnitude limited and hence evolved massive stars such as supergiants are also abundant. Supergiants are much brighter than MS dwarfs, and their mass-magnitude relations are very different, impacting the used OB+BH-identification method. Hence, supergiants too are a potential cause of many false positives that has not yet been quantified. Finally, in Paper I the authors used much less stringent constraints than those that were implemented when constructing the recently published Gaia Data Release 3 (DR3) catalogue.

In this work we investigate the impact of a different mass-magnitude relation, from supergiant primaries, and additional companions on the fractions of identifiable OB+BH systems and false positives. We also explore how accurately the primary masses need to be known to limit the amount of false positives. Moreover, we investigate the effect of using the adopted Gaia selection criteria to publish astrometric orbits in DR3 and compare to this the results predicted in Paper I.

Section 2 gives an overview of the method used in Paper I for the identification of OB+BH systems, together with a summary of the possible biases and uncertainties. The impact of the adopted mass-magnitude relation on the identification fractions, higher-order multiple systems, and the accuracy of the primary mass measurements are discussed in Sects. 3 to 5. The results of OB+BH detections with DR3 in comparison to new predictions using the detection criteria made in DR3 are discussed in Sect. 6. We end with a summary in Sect. 7.

2. The identification method, possible biases, and uncertainties

2.1. The astrometric mass-ratio function

The astrometric mass-ratio function (AMRF; SMF19) was adapted to massive stars in Paper I. The AMRF is a unit-less parameter defined as

(1)

where, α is the observed astrometric signal (the observed semi-major axis traced by the photo-centre), ϖ is the parallax of the system, M1 is the mass of the primary or most luminous star, and P is the orbital period of the system in years. In the right-most term, q = M2/M1 is the current mass ratio of the secondary (least luminous) to the primary star and S = I2/I1 their intensity ratio in the pass band of the instrument measuring the astrometry, in this case the GaiaG-band. From Eq. (1), it is clear that the AMRF is sensitive to the adopted mass-luminosity relation through S. The right-most term in Eq. (1) will be referred to as the theoretical AMRF and the middle term as the observational AMRF or 𝒜.

Both SMF19 and Paper I considered three different scenarios for the AMRF method: single-degenerate binaries (MS+BH/NS), MS binaries (MS+MS), and triple systems where the outer star is more luminous than the inner binary (triple) since these systems result in the highest AMRF curves (see Sect. 4.1). Hence, the inner binary in the triple systems is considered the secondary. Since we are here focusing on massive primaries, the abbreviations above change to OB+BH and OB+MS.

The theoretical AMRF curves for systems with massive MS dwarf primaries are determined in Paper I for the three abovementioned kinds of systems (see Fig. 1) as a function of q. However, without a mass estimate for the luminous companion, the mass ratio of the system is unconstrained. We therefore work under the assumption that q is not a known quantity.

thumbnail Fig. 1.

Different theoretical AMRFs for three different primary masses: 10 M (blue), 20 M (orange), and 30 M (grey). The solid black line shows the theoretical AMRF for MS+BH systems, the dotted curves are for OB+MS binaries, and dashed curves are for triple MS dwarf systems. Different classes are also indicated for the 30 M primary by the horizontal dotted lines at the maxima of the triple AMRF and the binary AMRF (Class III: identifiable OB+BH systems, Class II: triple MS or OB+BH systems, Class I: OB+MS binaries, triple MS systems, or OB+BH binaries).

SMF19 introduced three classes, based on the principle that, without measurement uncertainties, the observational AMRF of any kind of system cannot exceed the theoretical curve for that kind of system (e.g. OB+MS binaries cannot exceed the theoretical curve for OB+MS binaries). Class III systems are systems that are found with an observational AMRF above the maximum of the triple curve. These are hence classified as OB+BH systems. Systems found with an observational AMRF below the maximum of the triple curve but above the maximum of the binary curve are classified as Class II. Here, we find a mix of triples and OB+BH systems. Class I systems are systems with an observational AMRF below the maximum of the binary curve. Class I contains OB+MS systems, but also triples and OB+BH systems. The three classes are marked in Fig. 1 only for a 30 M primary.

2.2. Biases previously considered in Paper I

The identification criterion for Class III imposed in Paper I is such that a system requires an observational AMRF above the maximum of the triple curve within one sigma or 𝒜 − σ𝒜 >  𝒜triple. Although the three classes are in principle well defined, measurement uncertainties might cause observational data points in the diagram to shift into a different class, which impact the identification of single-degenerate binaries through the AMRF method. A proper determination of the true numbers of the expected OB+BHs is important, as it can contain information on the survivability of single-degenerate systems through the collapse of their non-degenerate progenitors (Paper I).

On the one hand, OB+MS or triple systems might get falsely identified as OB+BHs. These false identifications are called false positives. On the other hand, measurement uncertainties can also cause some OB+BH systems to appear as Class II or I systems. Both false positives and non-detections can affect the number of identifications. As shown in Paper I, the expected contribution of OB+MS binaries to the amount of false positives is negligible.

2.3. Biases and uncertainties considered in this work

Stars spend most of their life as MS luminosity-class V (i.e. dwarfs). Main-sequence dwarf stars are hence expected to be the most abundant stars in the Universe. However, when considering a magnitude limited sample, bright stars are typically over-represented (unless the magnitude cut encompasses all stars), especially at larger distances. Non-dwarf stars have a considerable different mass-magnitude relation than dwarf stars and will thus change the AMRF curves. Paper I only considered the mass-magnitude relation for dwarfs. A consideration of massive binaries with non-dwarf primaries might not only impact the determined AMRF curves, but might also have a large impact on true negatives and false positives. In Sect. 3, we investigate the impact of non-dwarf primaries. Here, we also briefly touch upon the impact of systems with a stripped star, which resemble systems such as LB-1 and HR6819.

In the investigation of false positives, Paper I only investigated the contribution of OB+MS binaries. Although triple systems with outer periods shorter than the Gaia observation time – 3 yr for DR3 or 5 yr for the fourth data release (DR4) – are not expected to be abundant (Sana et al. 2014; Tramper et al., in prep.), their possible contribution to the false-positive identification should be investigated, since we are searching for systems that are intrinsically rare. This is done in Sect. 4, together with a discussion on higher-order hierarchical multiple systems.

In the ideal scenario, all four parameters in the middle term of Eq. (1) are known with good accuracy and precision. The astrometric signal α, the parallax ϖ, and the period P will be published with the astrometric Gaia orbits. The mass of the primary M1 is for some sources also provided in DR3. However, it is either lacking for most massive stars or inaccurate for most massive stars that have a published M1. Therefore, it must be determined using other data. We discuss the effect of the accuracy of the mass measurement in Sect. 5.

To determine the fraction of detectable OB+BH systems, Paper I required that the astrometric signal produced by a system is larger than three times the Gaia along-scan variance (see also Sect. 3.3). In Sect. 6, we present new predictions on the detectability of OB+BH systems using the basic selection criterion implemented by the Gaia Collaboration in DR3. We also compare these predictions with the ALS II sources having an astrometric DR3 orbit.

3. The impact of an OB supergiant mass-magnitude relation on the AMRF

The theoretical AMRF curves are sensitive to the adopted mass-magnitude relation through the band-specific intensity ratio S. Evolved non-dwarf primaries or supergiants are substantially brighter than MS dwarf stars. Therefore, non-dwarf primaries can significantly alter the theoretical AMRF curves presented in Fig. 1, potentially altering both the identification fractions as well as the amount of false positives.

In order to fully understand which systems might end up as false positives, systems with evolved supergiant primaries should also be included in the mock population. Since systems with cool red supergiant (RSG) primaries have very long periods (> 5 yrs) due to their large radii, they are not expected to have a binary solution in DR3 or DR4. Potentially, the large convection zones on the surface of RSGs can cause their photo-centre to change, creating a chaotic photo-centre motion, which is also not in favour of RSGs having a binary solution (Chiavassa et al. 2022). Therefore, we only considered blue supergiant (BSG) primaries.

In Sect. 3.1, we derive a mass-magnitude relation for BSGs appropriate for Gaia. We obtain the BSG theoretical AMRF curves in Sect. 3.2. We simulate a population of binaries with OB and BSG primaries in Sect. 3.3, where we also explore the effect of the newly obtained BSG curves. In Sect. 3.4, we briefly discuss the impact of BH imposters (i.e. stripped He stars) such as the systems LB-1 and HR6819.

3.1. The G-band magnitudes of BSGs

In order to create a mock population of BSGs and their theoretical Gaia AMRF curves, the magnitudes of BSGs need to be obtained. As a starting point, we used the parameters derived for stars with a luminosity-class I from Martins et al. (2005), corresponding to spectral types between O3.5 I and O9.5 I and masses 30 M ≲ M ≲ 66 M. For lower-mass BSGs (9 M ≲ M ≲ 30 M), we obtained parameters from the evolutionary tracks published by Brott et al. (2011) at the time the core-H fraction reaches 10%. We used the galactic tracks for stars with initial masses in between 9–60 M and initial rotational velocities in between 250–300 km s−1, corresponding to one track per mass in their model grid. The positions of the BSGs in the Hertzsprung-Russell diagram (HRD) are shown in Fig. 2.

thumbnail Fig. 2.

Positions in the HRD for BSGs. Cyan dots are luminosity-class I stars from Martins et al. (2005) and blue stars are the positions on the tracks corresponding to a 10% core-H fraction. The dashed line represents a second-degree interpolation between the data from Martins et al. (2005) and the four lowest-mass points with a 10% core-H fraction (see text for more details).

The data from Martins et al. (2005) are for O supergiants and do not match the indicated positions for massive stars with a core-H fraction of 10% on the tracks, where the most massive stars would appear as cool B supergiants. In order to connect both sets of data, we used a second degree spline to interpolate the luminosities between the data from Martins et al. (2005) and the four lowest-mass points with a core-H fraction of 10% as a function of effective temperature. This result is shown in Fig. 2 along with the corresponding positions on the evolutionary tracks.

The same fit was also performed using the first five low-mass points, slightly altering the fit parameters. Hence, the uncertainties in Table 1 in the low-mass range are dominated by this effect. The parameters for the BSGs are given in Table A.1.

Table 1.

Fit parameters for the mass-magnitude relation (G = alog(M/M)+b) of BSGs in different mass regimes.

The absolute G-band magnitudes of the BSGs were determined in the same way as described in Sect. 3.2 of Paper I. We derived a log-linear mass-magnitude relation in the same way as described in Sect. 4.1 of Paper I. The fitted relation is shown in Fig. 3 and the fit parameters are given in Table 1.

thumbnail Fig. 3.

Absolute G-band magnitudes for BSGs in the mass range 8.9–66.89 M obtained from Martins et al. (2005) and Brott et al. (2011) (see Table A.1). The solid red line shows the fit and the vertical grey lines indicate the different mass regimes. The masses are shown on a logarithmic scale. As a reference, the obtained magnitudes for MS dwarf stars from Paper I are also shown.

We also explored an additional scenario, using the positions on the HRD tracks where the star has a core-H fraction of 1%. Equivalents of Figs. 2 and 3 are shown in Appendix B when using a core-H fraction of 1% on the tracks (Figs. B.1 and B.2). The main difference between the two is for stars with M <  30 M.

3.2. The theoretical AMRF curves for BSGs

Combining the mass-magnitude relation obtained in Sect. 3.1 and the right-most term in Eq. (1), the theoretical AMRF curves for BSGs can be determined. Unless stated otherwise, the used mass-magnitude relation is that obtained from the interpolation between the data from Martins et al. (2005) and the four lowest-mass points with a core-H fraction of 10% on the tracks from Brott et al. (2011).

The BSG curves are shown in Fig. 4 for two different primary masses together with the curves for MS dwarf stars or MS curves. The companion for the BSG-binary curves is assumed to be an unevolved MS dwarf star and for the triple curve a double MS binary companion is assumed. As expected, the BSG curves are reaching higher AMRFs compared to the MS equivalent because BSGs are more luminous and hence yield larger astrometric signals (Eq. 1). This can also be seen from the horizontal lines, indicating the maximum of the curves and thus the separation between the Classes. For a 10, 20 (not shown), and 30 M primary, the Class III thresholds increase from 0.36, 0.31, and 0.26 to 0.51, 0.45, and 0.42, respectively. The Class II thresholds change from 0.23, 0.19, and 0.16 to 0.35, 0.30, and 0.34 for a 10, 20 (not shown), and 30 M primary, respectively.

thumbnail Fig. 4.

Different theoretical AMRF curves for two different primary masses: 10 M (left) and 30 M (right). In addition to the MS curves, AMRF curves for BSGs are also shown for the same primary masses. The different classes are now also indicated for both the OBs and the BSGs.

To investigate the impact of the adopted core-H fraction, we obtained theoretical AMRF curves using the mass-magnitude relation obtained using a 1% core-H fraction. A comparison between the 10% curves and the 1% curves is shown in Fig. B.3. The 1% curves can reach values that are even higher than the 10% curves.

3.3. A mock population that includes BSG primaries

We performed a cross-match between the ALS II and the Astronomical database SIMBAD1 to obtain spectral types for the sources in the ALS II. Of those that have a spectral type, about 50% had at least one companion that was not a dwarf (luminosity-class V) and hence more evolved. With this information, we simulated a population of binaries where 50% of the systems have a BSG primary.

We performed a Monte Carlo (MC) simulation to immediately include the effect of measurement uncertainties. Following Paper I, we take the uncertainties on α and ϖ equal to the derived along-scan variance data provided to us by A. Everall (priv. comm., based on a similar method described in Everall et al. 2021). We use a 10% uncertainty on P and a 50% uncertainty on M1. The published DR3 periods have on average a 3% precision, smaller than what is used here. However, in order to directly compare to the predictions of Paper I, a 10% uncertainty on the period is used. If the precision on P would be lowered, this could result in a slightly lower fraction of false positives, as it decreases the mobility of the data points in the AMRF diagram. For the systems, Ω (longitude of the ascending node) and ω (argument of periastron) were uniformly randomly distributed between 0 and 2π, and cosi (where i is the inclination) between −1 and 1. Distances were drawn from the distribution of distances in the ALS II.

In total we simulated 1000 samples of 150 OB+BH, 150 BSG+BH, 3500 OB+MS, and 3500 BSG+MS binaries. The numbers originate from the estimations done in Sect. 6.3 of Paper I. From the > 13 000 sources in the ALS II, 50% would be post-interaction binaries (de Mink et al. 2014), amounting to ∼7000 sources. Furthermore, 3% of OB stars in binaries are expected to have a BH companion (Langer et al. 2020), resulting in almost 300 OB+BH systems. Using now that 50% of sources in the ALS II are non-dwarf stars, we obtain the above-stated population.

To determine which systems are detectable by Gaia, we adopted the same assumptions as Paper I used for DR3: P <  3 yr and α >  3σG, DR3. We note that, due to α always being positive, the criterion of α >  3σG, DR3 may result in a larger fraction of single stars falsely identified as binaries in comparison to that expected in the case of a normal distribution (0.3%). This criterion is hence an optimistic one and, as discussed in Sect. 6, it is currently not used for DR3. However, we do adopt it here to facilitate a comparative study to Paper I.

Also similar to Paper I, identifiable systems need to be detectable and have 𝒜 − σ𝒜 >  𝒜triple. Here, 𝒜triple can be either from the MS or the BSG curves. Unless stated otherwise, these assumptions will hold for the remainder of the manuscript.

The average detection and identification fractions are given in the first four columns of Table 2 and the estimated numbers of sources in the ALS II are given in Table 3. In Table 2, fdet is the fraction of systems detectable by Gaia, fid is the fraction of the systems detected by Gaia that would be identified with the AMRF method, and ftot is the fraction of systems identified in the full simulated population, defined as

(2)

Table 2.

Detection, identification, and total identification fractions for different kinds of systems.

Table 3.

Estimated number of sources detected and identified.

First, we can see that the usage of the BSG curves greatly reduces the fraction of identifiable OB/BSG+BH systems – true positives – by a little more than half. For the OB+BH and BSG+BH systems, the identification fractions now decrease to 29% compared to 68% and 65%, respectively, when using the MS curves. Assuming around 300 OB/BSG+BH systems, about 85 BHs can be detected compared to almost 200 when using the MS curves. Hence, a significant fraction of the population is lost when using the BSG curves.

Second, using the BSG curves also reduces the amount of false positives. When using the MS curves, the fraction of false-positive OB+MS binaries is around 0.1% and it is well below 0.1% when using the BSG curves. With an estimated 3500 of OB+MS systems in the ALS II, the number of false-positive systems would be between ∼1–6 depending on the used curves. This is indeed a negligible amount just as shown in Paper I.

The false-positive fraction of BSG+MS systems is on average about 1.3% when using the MS curves. Given 3500 systems, this would translate to ∼50 systems in the ALS II. With an estimated 200 OB/BSG+BH systems, this would mean 20–25% of the identified sources are false positive. When using the BSG curves, this fraction decreases to 0.2%, which would result in ∼7 systems compared to ∼85 OB/BSG+BH binaries.

We thus see that the inclusion of BSGs changes the identification fractions. Ideally, the evolutionary stage of the primary is known and each system can then be compared to the curves corresponding to the evolutionary stage of the primary to benefit both from the high identification fraction of OB+BHs and the low false-positive fraction of BSG+MS systems. However, if the evolutionary stage of the primary is not known one might consider using the BSG curves to obtain the least amount of false positives despite the significant decrease in true positives.

3.4. BH imposters such as LB-1 and HR 6819

The systems LB-1 (Liu et al. 2019) and HR 6819 (Rivinius et al. 2020) are two of the most recent and well known BH imposters. Both these systems host a Be star bound to a bloated stripped star that is 5–10 times less massive than its Be companion, resulting in an extreme mass ratio. Such binaries are in an extremely rare evolutionary phase, caught as the stripped star contracts towards the helium MS, while it is still extended and visually bright. The fact that the stripped star is overluminous compared to its mass can make these systems appear as OB+BHs in the AMRF diagram. However, according to Bodensteiner et al. (2020), Shenar et al. (2020), and El-Badry & Quataert (2021), LB-1 and HR 6819 consist of two components with very similar visual brightness. This hence means that the visual photo-centre of these systems is located very close to the centre of mass and thus the photo-centre motion is very small. Therefore, we would not expect these systems to be astrometrically detected as binaries in DR3.

Nevertheless, it is possible that other ‘LB-1-like’ systems exist in which the stripped star is substantially brighter than its more massive companion. In such cases, the system could be erroneously identified as an OB+BH binary using our method. Given the overall rarity of such binaries and the associated uncertainties of binary evolution, it is difficult to assess the contamination of such systems among the OB+BH candidates. However, a rough estimate can be made using the prototypical example of HR 6819. This is the only binary within a ∼350 pc volume that is known to host a bloated stripped star. There are roughly 300 OB stars in the ALS II catalogue within this volume, and, according to Langer et al. (2020), ∼2% should host BHs, amounting to roughly six targets. Hence, the population of OB+BH might be expected to be larger by a factor of a few. This factor increases if one considers that the bloated stripped stars need to be significantly more luminous than their companions to be confused as OB+BH binaries astrometrically. An important caveat here, however, is that the estimated number of OB+BH binaries is a prediction, not an observation, while the estimated number of binaries hosting bloated stripped stars hinges on the assumption that HR 6819 is unique within this volume.

El-Badry & Rix (2022) analysed a sub-sample of single-lined spectroscopic binaries (SB1s) detected by Gaia in DR3 and reported they all contain a bloated stripped star. While no Gaia SB1 compact-object systems have been reported so far, the Gaia astrometric binary solutions revealed one BH (Gaia DR3 4373465352415301632, El-Badry et al. 2023) and one NS or BH candidate (Gaia DR3 6328149636482597888, Andrews et al. 2022; El-Badry et al. 2023), both orbiting a low-mass star. This could indicate that bloated stripped stars are more common than BHs. However, the stripped stars around massive stars most likely originate from another massive star, whereas for the low-mass binaries detected in DR3, the stripped stars most likely originate from other low-mass stars due to the selection biases of DR3. Thus, the stripped-star and BH populations in these DR3 low-mass binaries originate from different progenitors. In contrast, for massive stars, both stripped stars and BHs have the same progenitors. Hence, the apparent overabundance of stripped stars in low-mass binaries in DR3 compared to BHs cannot be directly extended to the massive star domain.

To conclude, the numbers currently indicate that, among a population of OB+BH candidates, not more than 10–20% should be LB-1-like systems erroneously classified as OB+BH. While relatively small, this number is not negligible. Moreover, it is possible that the true proportion of such imposters is larger (e.g. if these systems are more common, or if OB+BH binaries are rarer than predicted). Dedicated spectroscopic and interferometric surveys of the candidates would be needed to fully remove this uncertainty (e.g. El-Badry & Rix 2022; Frost et al. 2022; Shenar et al. 2020).

4. Higher-order multiple systems

The theoretical AMRF curves for OB+MS binaries cover much lower AMRF values than those for triples. Indeed, triple systems can produce larger astrometric signals, as a binary companion is less bright than a MS star of equal mass. Observational errors thus can more easily cause a triple system to become a false positive than for an OB+MS binary. Here, we investigate the impact of triple systems and higher-order multiple systems on the amount of false positives.

The number of hierarchical triples or higher-order multiple systems that have a period shorter than ∼3–5 years are limited (∼10%; Sana et al. 2014). On top of that, the binary orbits published in DR3 satisfy hard quality cuts (see Sect. 6.1 for some details). Chaotic photo-centre motions caused by higher-order multiple systems where all three components significantly contribute to the light will thus most likely not have a clear astrometric solution even though for DR4 the quality cuts will probably be slightly less conservative than those for DR3. Therefore, based on the Gaia detection alone, we do not expect higher-order multiple systems to greatly impact the amount of false positives. However, there might be systems with outer periods P ≲ 3 yr where the photo-centre motion is not extremely chaotic that will have an astrometric binary solution and the impact of these needs to be investigated.

First, we explored the effect of different triple configurations on the theoretical AMRF curves (Sect. 4.1). Then, in Sect. 4.2, we simulated a grid of triples to investigate which triple configuration has the largest impact on the amount of false positives. Next, in Sect. 4.3, we used the determined distributions for the triples in the Southern MAssive Stars at High angular resolution survey (SMASH+; Sana et al. 2014) to simulate a population of triples. Last, we discuss the impact of quadruples and higher-order multiple systems in Sect. 4.4.

4.1. Investigation of distinct triple configurations

The theoretical AMRF curves for triples in Fig. 1 represent the case where the outer star is more luminous than the inner binary. On top of that, the curves also assume that the inner binary has a mass ratio of unity (q1 = 1, not to be confused with the q on the x-axis, which is in this case the mass of the inner binary divided by that of the outer star). Only the curves for q1 = 1 are shown because they are the ‘worst-case’ scenario (i.e. they give the highest AMRF curves). Starting at q1 = 1, which coincides with the previously presented curves (Fig. 1), the AMRF curves keep on moving downwards for decreasing q1, with the limit being a binary system at q1 = 0 (Fig. 5). The inner binary appears under-luminous compared to a single star with a mass equal to the total mass of the inner binary. Hence, the photo-centre will be located closer to the outer star and produce a larger astrometric signal for any 0 <  q1 ≤ 1 than compared to a binary system.

thumbnail Fig. 5.

Different theoretical AMRF curves for a 20 M primary. Triple curves for different mass ratios of the inner binary q1 are shown for the case where the inner binary is the most luminous (‘Inner’) or the outer binary is more luminous (‘Outer’).

What is thought to be more common are triples where the most luminous component is located in the inner binary. For these systems, q is the mass of the outer star divided by that of the inner binary and M1 in Eq. (1) is the mass of the most luminous star in the binary. We can distinguish the two extremes: q1 = 1 and q1 = 0. The latter case again reduces the system to a binary system. In contrast to the other kinds of triples, where q1 = 1 is the worst-case scenario, here q1 = 1 is the ‘best-case’ scenario, giving the lowest AMRF curves (Fig. 5). Following the same reasoning as before, the inner binary looks under-luminous compared to a single star with a mass equal to the total mass of the inner binary. Hence, the photo-centre will be located farther away from the inner binary and, in this case, produce a smaller astrometric signal for any 0 <  q1 ≤ 1 than compared to a binary system.

From the theoretical AMRF curves presented in Fig. 5, we can argue that triples where the most luminous component is the inner binary will not contribute significantly to the amount of false positives. The reasons for this are that these kind of triples all fall below the OB+MS curves and that OB+MS binaries already give very few false positives (Paper I and Sect. 3.3). Thus, triple systems where the outer star is the most luminous are the only triples with the potential to cause many false positives, as these have AMRF curves above those of OB+MS binaries.

4.2. A grid approach

As previously discussed, the worst-case scenario is that in which the binary companion has a mass ratio equal to unity. Therefore, we only simulated triples with q1 = 1. We varied the outer eccentricity between 0 and 0.75 with a step size of 0.25. The outer period Pout ranged from 10 to 1000 days with a logarithmic step size of 0.25. For the masses, we varied the mass of each component in the inner binary between 4 to 40 M. We varied the mass ratio q of the triple system from 0.1 to 2.1 with a step size of 0.2, such that the mass of the outer star is equal to Mout = 2M1q. We put a lower and upper limit on Mout, which are 8 M and 60 M, respectively.

Again, we performed a MC simulation with the same uncertainties and randomly drawn orientations and distances as described in Sect. 3.3. For each combination of the above-stated parameters, 6000 systems were simulated with random orientations and distances.

The total identification fractions (i.e. the fraction of systems detected by Gaia and identified with the AMRF) are shown in Fig. 6 when using the MS curves, for the mass of one of the components of the inner binary M1a (M1a = M1b) or the outer period Pout as a function of the mass of the outer star Mout. More parameters are shown in Fig. C.1. The fraction of false-positive triples for OB stars is fairly low (< 3%), even at the grid points giving the highest false-positive fraction. However, for triple systems where the outer star is a BSG, the false-positive identification fraction can reach values of 50%, which is extremely high. This, however, does not result from the fact that the system is triple but because the AMRF curve for dwarfs are inappropriate for BSGs (Sect. 3).

thumbnail Fig. 6.

Total identification fraction for triple systems where the outer star is the most luminous component and the inner binary has q1 = 1. In the left panels, M1a is the mass of one of the components in the inner binary. The top panels are for OB systems, and the bottom panels are for systems with a BSG outer star.

In Fig. 7, the same false-positive fractions are shown but now when using the BSG curves (and more parameters in Fig. C.2). This decreases the false-positive fraction among systems with BSG primaries to < 3%. Hence, the usage of the BSG curves can significantly reduce the amount of false positives in case the luminosity class of the primary is not known. However, again, the usage of the BSG curves also significantly reduces the amount of true positives. Therefore, it is highly recommended and beneficial to establish the nature of the most luminous component to still obtain the high amount of true positives among the OB+BH systems.

thumbnail Fig. 7.

Same as Fig. 6, but when using the BSG curves.

4.3. A mock population of triples

The SMASH+ has observed a significant amount of triple systems with massive star components. The SMASH+ combines high angular resolution data of massive stars in the southern hemisphere with (most often spectroscopic) information found in the literature to investigate the multiplicity of massive stars over the full range of separation up to > 10 000 au. Tramper et al. (in prep.) have derived masses and physical separations for the companions detected with interferometry. Using this information, we can determine distributions for the periods and mass ratios of the triple systems.

Practically all of the triples in the SMASH+ database have an inner binary that is more luminous compared to an outer third component. As explained in Sect. 4.1, triples where the inner binary is more luminous theoretically have a smaller 𝒜 than binaries. Hence, we do not expect those to significantly contribute to the amount of false positives. To verify this, we simulated a population of triples using the distributions derived from SMASH+ (Tramper et al., in prep.). While doing this, we assumed that the photo-centre motion of the inner binary combined with the photo-centre motion of the triple is not chaotic, but resembles a binary (i.e. that the point-spread function of the triple would not be impacted by the inner-binary orbit).

For triples where the outer star is the most luminous, it is effectively more difficult to identify those as triples instead of binaries. Therefore, the SMASH+ database might be biased against such systems, though a few are observed. Nonetheless, such triples are not expected to be very abundant.

We simulated both triples with all MS dwarf components and triples where the most luminous star is a BSG. A similar approach was used as described in Sect. 3.3, where we used a MC simulation with similar uncertainties and random orientations and distances. In total, 1000 samples of 7000 systems were simulated both for OB and BSG primaries.

The detection and identification fractions can be found in the two last columns of Table 2 and the number of estimated sources in Table 3. The fraction of triples estimated to be detected by Gaia DR3 is only ≲4%. This is mostly because either the outer period of the system is much larger than the DR3 baseline or because the mass ratio of the outer star to the inner binary system is very low, resulting in a very small photo-centre motion.

As expected, the false-positive identification fraction of the SMASH+ triples detected in DR3 is extremely low (around 0.01%), resulting in potentially one false detection. Since most triples have fairly large mass ratios for the inner binary, these systems will produce much smaller photo-centre motions as binaries and end up at lower values for 𝒜, making it more difficult for these systems to be falsely identified. If instead we would use the BSG curves, the identification fractions become practically zero. Unlike with binaries, for the triple population as observed by the SMASH+, there is no direct need to use BSG curves instead of MS curves.

4.4. Quadruples and higher-order multiple systems

Besides triples, higher-order multiple systems exist and contribute to ∼5–10% of the sample studied by SMASH+. There are two kinds of hierarchical quadruple systems: (i) two binaries orbiting around each other (e.g. QZ Car Rainot et al. 2020) and (ii) hierarchical systems formed by a binary orbited by two other stars at increasing separation (e.g. Plaskett’s star Linder et al. 2008; Sana et al. 2014). In both cases, the binaries themselves most likely will be unresolved. In the hierarchical case, the photo-centre motion will be very chaotic, unless one or both of the outer stars are very low mass and thus have a relatively low luminosity compared to the inner binary. If only one star is very low mass, the photo-centre motion approximately reduces to that of a triple system, a scenario that is already discussed. In case both outer stars are low-mass stars, the photo-centre motion will be very small, because the luminous binary in the quadruple will not show large motion. Hence, it is likely that these sources will not have an astrometric binary solution in DR3 or DR4. Alternatively, the published astrometric solution is that of the inner binary, in which case the quadruples should not greatly impact the amount of false positives.

For the quadruples in case (i), there are three (extreme) configurations where the photo-centre motion will resemble closely that of a binary system : (a) two binaries each with a mass ratio equal to unity (e.g. (15+15) M + (10+10) M), (b) one binary has an extreme mass ratio (e.g. (29+1) M + (10+10) M), and (c) both binaries have extreme mass ratios (e.g. (29+1) M + (19+1) M).

An (a)-type quadruple gives the same AMRF as a binary and is most easily explained by an example. In the example of a (15+15) M + (10+10) M, this quadruple has the same mass ratio and intensity ratio as a (15+10) M binary. Hence, the theoretical AMRF reduces to that of a binary. As shown in Paper I, binaries are not expected to contribute to the false positives, and thus nor will these kinds of quadruples.

The case (b) quadruple resembles very closely a triple system. As shown in Sects. 4.3 and 4.2, the number of false positives expected from triples is very low and will even be lower for these kinds of quadruples as they are even less abundant.

The last case, case (c), closely resembles the case of a binary system. Due to the extreme mass ratios of both binary components, the AMRF of these kinds of systems will approach that of binaries. Hence, we also do not expect case (c) quadruples to affect the false positives.

For even higher-order multiple systems, like quintuples, the reasoning is very similar to that of quadruples. Either the photo-centre motion is very chaotic, and no solution for this system will be given in DR3, or the systems photo-centre motion can be approximated by any of the previous discussed cases (binaries, triples, and quadruples). Hence, higher-order multiple systems should also not impact the false-positive OB+BH identifications.

5. The non-necessity of accurate mass measurements

Three of the four parameters necessary to determine the observational AMRF, the middle term of Eq. (1), can be obtained from the Gaia astrometric binary solutions (the astrometric signal, the parallax, and the period). The masses of the primary stars M1 are also provided for a sub-sample of stars in DR3. However, the DR3 astrophysical parameters of massive stars are unreliable due to a lack of quality diagnostic lines in the Gaia spectral range for objects hotter than 10 to 15 kK. Hence, reliable M1 values still need to be determined. There are several ways these masses can be estimated.

Direct mass measurements can be obtained from eclipsing double-lined spectroscopic binaries. However, these are often short period and most will be too tight to be identified in DR3. For non-eclipsing spectroscopic binaries, only the mass function or minimal mass can be determined instead of the absolute masses (see e.g. Mahy et al. 2020). Mass estimates can also be obtained from atmospheric analysis (effective temperature and surface gravity) or by comparison to evolutionary tracks. Even the knowledge of spectral type and luminosity class can provide a first estimate by comparison with calibration tables (e.g. Martins et al. 2005). However, spectroscopic data are not available for the vast majority of the sources in the ALS II.

In principle, effective temperatures and luminosities can also be determined through fitting of a spectral energy distribution. However, this fitting is often ambiguous for massive stars as the available photometric data in the visual and ultraviolet is in the Rayleigh Jeans domain.

Since there is no spectrum available at hand for most stars in the ALS II, we explored the effect of inaccurate and imprecise mass measurements of the primary on the identification fractions of OB/BSG+BHs and false positives (i.e. OB/BSG+MS and triple systems). Through a MC simulation, the observational parameters of the simulated systems were varied with the same uncertainties as used in Sect. 3.3. We distinguish between several scenarios. Three scenarios assume that the mass of the system is determined, but with different mass errors of 50%, 100%, and 200% (i.e. 0.5 times the mass, one time the mass, and 2 times the mass, respectively). In figures, these scenarios are marked as ‘0.5Merr’, ‘1Merr’, and ‘2Merr’, respectively. The lower and upper limit on the simulated masses was 7 and 60 M, respectively. The other scenarios investigate what would happen to the fractions if we assumed a fixed mass for all systems, for example, fixing all the masses of the primaries to 10 M, labelled with ‘10 M’.

We discuss the results of the simulations using the MS and BSG curves in Sects. 5.1 and 5.2, respectively. While these two sections argue that M1 does not need to be known accurately, the derived mass of the compact object will depend on M1 and might even change the nature of the compact object. We discuss this in Sect. 5.3.

5.1. Total identification fractions using the MS curves

Figure 8 shows the distribution of total identification fractions (see Eq. 2) of 100 samples of 150 OB+BH, 150 BSG+BH, 3500 OB+MS, 3500 BSG+MS, and also two sets of 7000 triples from the SMASH+ with OB and BSG primaries using the MS curves. The total identification fractions show a trend for single-degenerate systems with a BH (OB/BSG+BH systems, top row of Fig. 8). For the three scenarios assuming a determined mass with an error (‘0.5Merr’, ‘1Merr’, and ‘2Merr’) there is a clear difference. The larger the error on the primary mass, the larger the error on 𝒜, hence the smaller the identification fraction (the identification criterion of Paper I was 𝒜 − σ𝒜 >  𝒜triple). The total identification fractions are even slightly higher when using a fixed mass for all primaries, since no error on the mass is assumed in this case.

thumbnail Fig. 8.

Total fraction of identified systems when using the MS curves for single-degenerate binaries with a BH (top row), non-degenerate binaries (middle row), and triples from the SMASH+ (bottom row). Left and right panels are for OB and BSG primaries, respectively.

The value of the fixed mass does not significantly affect the total identification fractions of the single-degenerate systems. When assuming a different primary mass, the curves in the AMRF diagram will change, going downwards for higher primary masses (see Fig. 1). A higher primary mass also decreases 𝒜 (see Eq. (1)). When altering the primary mass, the rate with which both the MS curves and 𝒜 move downwards in the AMRF diagram is comparable over most of the massive-star dwarf mass range that was derived here (8 <  M/M ≲ 60). Hence, the value of the fixed mass does not play a large role for the identification of OB/BSG+BHs.

The above described effect of larger uncertainties on the mass or a fixed mass is almost not present for the false-positive identification of OB+MS binaries (middle-left panel of Fig. 8). The reason for this is that OB+MS systems are located very low in the AMRF diagram. Therefore, it is already very difficult for these systems to be falsely identified as a single-degenerate binary. Hence, the effect of larger mass uncertainties or zero mass uncertainties becomes barely prominent.

For the false-positive BSG+MS binaries, we see that a larger uncertainty on the mass does have the same effect as for the single-degenerate binaries. However, we now also see a change in fractions for different values of the fixed mass. The mass-magnitude relation for BSGs (see Fig. 3) is rather steep for 8 <  M/M <  30 and flattens for M ≳ 30 M, but the mass-magnitude relation for MS stars (see the grey stars in Fig. 3 in this paper and Fig. 3 in Paper I) does not. Over almost the whole massive mass regime, the mass-magnitude relation for MS stars retains almost the same slope, except for 8 <  M/M <  15, where the slope changes more significantly. Because of the combination of the two mass-magnitude relations, BSG+MS binaries with different primary masses take up significantly different positions in the AMRF diagram compared to their OB+MS counterparts. The effect of this is non-trivial and is most likely the cause of the change in total identification fractions for different values of the fixed mass. However, each value of the fixed mass does still result in total identification fractions with the same order of magnitude.

Finally, we can see that a fixed value for the mass slightly increases the total identification fraction on average for SMASH+ triples with a BSG primary because of the same effect as with the BSG+MS binaries. Overall, the total identification fraction of the SMASH+ triples is extremely low, both for OB primaries as well as BSG primaries, as expected from Sect. 4.3.

When using the MS curves, an accurate mass estimate is not necessary to identify OB/BSG+BH binaries. In fact, a fixed value of the mass will lead to slightly higher total identification fractions because in that case no error on the mass is assumed. The order of magnitude of the false-positive identification among the BSG+MS binaries remains the same and below 5% for any explored scenario, but can vary depending on the value of the fixed mass.

5.2. Total identification fractions using the BSG curves

The total identification fractions for OB/BSG+BH binaries become much lower when using the BSG curves (Fig. 9, top panels), as already shown in Sect. 3.3. As for Sect. 5.1, a larger uncertainty on the mass results in lower fractions. This time, however, we can also see here that different values for the fixed mass somewhat change the fractions, even for the OB+MS systems. The decrease of the maximum of the theoretical BSG curves and 𝒜 is no longer comparable and hence causes a change in the fractions when using different values for the fixed mass.

thumbnail Fig. 9.

Same as Fig. 8, but when using the BSG curves.

The total false-positive identification fractions of BSG+MS binaries using the BSG curves (Fig. 9, middle-right panel) now behaves the same as those for OB+MS curves when using the MS curves (Fig. 9, middle-left panel). They are extremely small, reaching maximum values of ∼0.4%. The total identification fractions for OB+MS binaries are even lower when using BSG curves, with a maximum below 0.2%.

The triples in the SMASH+ are not resulting in any false positives (Fig. 9, bottom panels). Sometimes however, one system (∼0.014%) is falsely identified when using the BSG curves.

When using the BSG curves, the choice of the value of a fixed mass can slightly affect the total identification fractions. However, in all panels, we can see that the order of magnitude of the fractions remains the same. Especially for the false-positive identifications, the most important result is that the value of the fixed mass will not significantly increase nor decrease the fraction of false positives.

5.3. The effect on the mass of the compact object

In Sects. 5.1 and 5.2, we show that the value of the primary mass does not need to be accurately determined in order to identify a massive binary as hosting a compact object. However, an inaccurate mass measurement will alter the derived mass of the compact object, as also recently shown by El-Badry et al. (2023). By measuring a different metallicity for the luminous component of Gaia DR3 5136025521527939072, El-Badry et al. (2023) concluded that the primary was almost 50% less massive than was reported in DR3, leaving open the possibility that the compact companion may also be a 1.2 M white dwarf (WD) instead of a 1.5 M NS (Gaia Collaboration 2023a).

The mass ranges of WDs and NSs are close to one another, such that a small decrease in the mass of a low-mass NS can easily make it fall into the WD mass range. In contrast, for typical BH masses, a large change in the primary mass is necessary to alter the nature of the compact object. For example, in one of the worst-case scenarios, an estimated primary mass of 16 M will have a minimum observational AMRF value 𝒜 = 0.37, corresponding to a minimum q of 0.47. This results in a minimum BH mass of 7.5 M. If the primary mass is actually 8 M, the minimum mass of the compact object at the same AMRF value is 5 M. The object would thus still be identified as a BH (see Shahaf et al. 2023, for the determination of the minimum compact-object mass from 𝒜).

An important caveat, however, is the following. While we have just shown that the identification of the binary as hosting a BH is robust against inaccurate mass estimates of the visible primary star, the accuracy of the estimate of the BH mass itself depends on the accuracy of the primary mass. Upon identification, accurate masses for these objects are then best determined by follow-up spectroscopy (as discussed at the beginning of Sect. 5). The spectroscopic data may further allow for a determination of the metallicity, which was the cause for the mass discrepancy in the WD or NS candidate.

6. Results from Gaia DR3

On June 13, 2022, the third Gaia data release DR3 was published (Gaia Collaboration 2023b), including the first Gaia astrometric binary solutions (Halbwachs et al. 2023). In total, the DR3 catalogue contains astrometric binary solutions for 169 227 sources. The assumptions adopted in Paper I led to an estimated 2000 sources in the ALS II having an astrometric binary solution in DR3 (see Table 3). Performing a cross-match between the ALS II and the DR3 astrometric binary sources table resulted in 11 sources, none of which satisfy our identification criterion of 𝒜 − σ𝒜 >  𝒜triple.

Although we find no strong OB+BH candidates, the non-detection reported here cannot be interpreted as a result of the BH-formation scenario, even though for some BH-formation scenarios we expected barely any detections. Instead, the small number of massive binaries with an astrometric binary solution in DR3 is in line with new predictions using the selection criterion imposed by the Gaia Collaboration, which we elaborate on in Sect. 6.1. We investigate how much less conservative the selection criteria need to be in order to obtain information on the BH-formation scenario in Sect. 6.2. In Sect. 6.3, we discuss briefly the work done for low-mass single-degenerate systems.

6.1. The impact of the DR3 selection criteria on the number of published astrometric binaries

There are several selection criteria for publishing astrometric binary solutions listed in Halbwachs et al. (2023). Here, we investigate the impact of two of those.

The required precision on ϖ

One of the criteria for an astrometric binary solution to be published in DR3 is related to the relative parallax precision:

(3)

where ϖ and σϖ are the parallax and its error and Pday is the derived period in days. For orbital periods of 100, 500, and 1000 days, this criterion implies a relative precision on ϖ of 200, 40, and 20, respectively. The majority of the OB+BH systems are expected to have P = 100 − 500 days.

In order to investigate how many OB+BH sources we would expect in the DR3 astrometric binary catalogue assuming only the criterion in Eq. (3), we performed the same simulations presented in Paper I now with a different detectability criterion. We obtained parallaxes and parallax uncertainties from the Gaia catalogue for all sources in the ALS II. Instead of only drawing distances for the simulated OB+BH systems, we randomly drew a parallax and its uncertainty for each system in order to be able to determine ϖ/σϖ (this also equals the parameter parallax_over_error in the Gaia catalogue). The distribution of ϖ/σϖ of the sources in the ALS II is given in Fig. 10. There are no systems with ϖ/σϖ >  200.

thumbnail Fig. 10.

Histogram of ϖ/σϖ of the massive sources in the ALS II. The vertical lines indicate a relative precision on ϖ of 200, 40, and 20.

Using Eq. (3) as the detectability criterion for the simulated OB+BH systems, we find that only 0.14% are detected. Assuming 300 OB+BH systems in the ALS II, this would result in 0 or 1 detection. Similarly, for non-degenerate binaries, only 0.3% are detectable. Should there be 7000 binaries in the ALS II, this would lead to the detection of ∼20 astrometric binaries. This predicted amount is very close to the 11 ALS II sources with astrometric binary solutions in DR3. Hence, the non-detection of OB+BH candidates is expected from the adopted cut and does not teach us anything yet on the underlying BH-formation scenario.

The required precision on the derived astrometric signal

Another criterion is limiting the precision on the derived astrometric signal α0 or the significance s:

(4)

where σα0 is the uncertainty on α0. When using a significance of 10 (i.e. assuming a 10% uncertainty on α0), our simulations show that, after applying the criterion in Eq. (3), no additional non-degenerate massive binaries are filtered out. Hence, the detected non-degenerate massive-binary fraction in our simulations remains 0.3% or ∼20 astrometric binaries. On average, the significance of the astrometric DR3 binaries in the ALS II is ∼3 (i.e. a 3% uncertainty on α0). Hence, in our simulations, the criterion in Eq. (4) did not affect the number of published astrometric massive binaries in DR3.

6.2. Towards less conservative constraints

As shown in Sect. 6.1, the criterion in Eq. (3) is the most stringent one. This criterion has been imposed on the DR3 data to avoid many spurious solutions, which otherwise could result in the false-positive detection of many massive BH binaries. While the restriction also filtered out a lot of probably non-spurious solutions (see Figs. 3a and b in Halbwachs et al. 2023), a significant amount of false-positive sources could have led to wrong conclusions and many false reports. If this issue is resolved in DR4, Eq. (3) can be relaxed, opening up the possibility to discover a large population of BHs. Below we elaborate on how much less conservative the constraints need to be to obtain information on the BH-formation scenario.

We performed simulations with several constraints on the parallax precision while keeping the constraint on the significance as presented in Eq. (4). Here, ϖ and σϖ are randomly drawn from the ALS II catalogue and the relative uncertainty on α0 is assumed to be 3%.

As expected, we find that the constraints need to be significantly lowered to detect a significant amount of BHs. In order to be able to detect ∼20 OB+BH binaries or ∼7% of the simulated population, the constraint in Eq. (3) needs to be decreased by 75% (i.e. requiring that ϖ/σϖ >  5000/Pday). Here it is assumed that BHs form with no natal kick and no mass loss. However, there are several scenarios that lead to fewer OB+BH binaries when, for example, the BH receives a strong kick (see Paper I for the impact of different BH-formation scenarios on the predicted amount of OB+BHs).

In order to detect a sensible amount of OB+BH systems in other extreme scenarios involving strong kicks, we would argue for even lower constraints. We suggest using ϖ/σϖ >  1000/Pday or ϖ/σϖ >  5, resulting in three identifications in the ‘worst-case’ scenario (strong kick), while > 30 identifications in other scenarios.

When lowering the required relative precision on ϖ, and keeping the relative precision on α0 assumed to be 3%, the significance criterion (Eq. (4)) still does not impact the number of detectable OB+BHs. This shows that the precision on the parallax remains the most critical cut and the significance criterion does not severely impact the number of detectable OB+BHs.

We did not, however, take into account that the relative precision on ϖ improves in the astrometric binary solution and might still improve in DR4 when more data are available. In that case, our imposed constraints can easily be scaled.

To conclude, our simulations showed that the constraint on the relative parallax precision needs to be significantly lowered. If the relative parallax precision of the DR3 single star solution catalogue is used, a decrease by 95% is needed. Otherwise, if the relative precision improves due to more data or the binary solution, the maximum conservativeness the constraint can have is

(5)

with n the improvement on the relative parallax precision compared to that of the DR3 single star solutions. The constraint on the significance does not necessarily have to be lowered.

6.3. Lower-mass single-degenerate binaries with Gaia DR3

While the manuscript focuses on massive single-degenerate binaries with a BH component, it is worth mentioning that some candidates of low-mass single-degenerate binaries with a NS or BH component (i.e. NS or BH systems that have a luminous companion with an initial mass < 8 M) are reported in several other manuscripts. Already in Gaia Collaboration (2023a), one of the Gaia verification papers, the authors use the AMRF on a sample of DR3 astrometric binaries with primary masses derived by Gaia in the range of 0.2 − 2 M. While most of their candidates are likely WDs, they also report on a NS candidate, Gaia DR3 5136025521527939072. However, as mentioned in Sect. 5.3, El-Badry et al. (2023) argue that this candidate may alternatively be a WD due to an ambiguity in the mass of the primary star.

Currently, there has been one report on a candidate BH system orbiting around a solar-type star (Gaia DR3 4373465352415301632, El-Badry et al. 2023), which was almost filtered out by Eq. (3). While the solution showed signs of potentially being spurious, follow-up spectroscopic data showed that the astrometric orbit derived by Gaia is indeed correct. This candidate was also picked up by Shahaf et al. (2023) as a candidate system containing a BH with the AMRF method. Shahaf et al. (2023) also report on some NS and WD candidates. Furthermore, there is also a report of an astrometric NS or BH candidate (Gaia DR3 6328149636482597888, Andrews et al. 2022; El-Badry et al. 2023).

Different research teams also used the other available Gaia binary solutions (i.e. photometric and spectroscopic), in the search for compact objects. While some photometrically detected single-degenerate binary candidates are reported (Gomel et al. 2023), several studies conclude there are no spectroscopic candidates (El-Badry & Rix 2022; Fu et al. 2022; Jayasinghe et al. 2022).

7. Summary

We have investigated how different biases and uncertainties may affect the number of expected OB+BH identifications predicted in Paper I. We have expanded that work by including the effect of a different mass-magnitude relation on the total identification fractions (i.e. the fraction of systems detectable by Gaia and identifiable with the AMRF). For this, we derived a mass-magnitude relation for BSG stars, resulting in the BSG curves.

We found that the contribution of OB+MS binaries to the amount of false positives is negligible, and around 0.1% when using the mass-magnitude relation for OBs derived in Paper I, that resulted in the MS curves. This fraction decreases when using the BSG curves. With an estimated 3500 OB+MS binaries in the ALS II, the number of false-positive systems would be between ∼2 and 6 depending on the curves used. The fraction of false-positive BSG+MS systems is ∼1.3% when using the MS curves and drops to less than 1% when using the BSG curves. For an estimated 3500 BSG+MS binaries in the ALS II, this would lead to ∼50 or ∼2 false-positive identifications. Comparing this to the expected ∼200 or ∼85 OB/BSG+BH identifications using MS or BSG curves, we can see that when including BSG primaries, it is better to use the BSG curves to increase the ratio of true identifications to false positives.

We further widened the contribution of false positives to triple systems and higher-order multiple systems. First, we simulated a grid of triples where the outer star is more luminous than the inner binary. There is a large difference between using the MS curves and using the BSG curves, especially for systems with a BSG primary. When using the MS curves, the false-positive fraction in some parts of the explored grid can go up to 50%. Although these kinds of triples are not expected to be very abundant, 50% is a highly non-negligible false-positive fraction. Instead, when using the BSG curves, the maximum false-positive fraction drops to around ∼2.5%. Here, the usage of the BSG curves can thus highly reduce the amount of false positives.

We also simulated a mock population using the triples in the SMASH+ both for OB and BSG primaries. Most of the triples in the SMASH+ have inner binaries that are more luminous than the outer star. Hence, the false-positive fraction of triples in the SMASH+ is extremely low (≲0.01%) when using the MS or BSG curves, potentially leading to only one false-positive identification. This is a highly negligible amount.

The contribution of quadruple or higher-order multiple systems should be extremely low. This is mostly due to the chaotic motion of the photo-centre and the generally long outer periods (P >  5 yr).

Finally, we have shown that an accurate mass measurement is not very important for the identification of OB/BSG+BH binaries, nor for the amount of false positives. We experimented both with different errors on the mass and with fixing the mass of the primary to the same value for the whole population. As expected, larger errors result in lower identification fractions. A fixed mass, however, results in even slightly higher identification fractions for the OB/BSG+BH binaries and, most importantly, not in significantly different false-positive fractions. This result is not any different when using the MS curves or the BSG curves.

Ideally, the evolutionary stage and mass of the primary star are known. This way, the corresponding mass-magnitude relation can be used to establish the AMRF curve that corresponds to the correct primary mass. If the evolutionary stage is not known, which is the case for most of the stars in the ALS II, then the BSG curves can be used instead of the MS curves, especially to lower the false-positive identification fractions of triples where the outer star is a luminous BSG.

Eleven ALS II sources have an astrometric binary solution in DR3, of which none are good OB+BH candidates. Although this null-result was predicted for some BH-formation scenarios, the low number of ALS II sources with an astrometric binary solution is likely caused by the stringent selection criterion imposed on ϖ/σϖ by the Gaia Collaboration.

Few extra compact-object candidates are found among the large pool of data provided by DR3. It is clear that the search for single-degenerate binaries in any mass regime remains a hot topic among several research groups. While no detection criteria are known yet for DR4, the relaxation of the Gaia detection criteria could tremendously improve the Gaia discovery capability for single-degenerate binaries, especially among the massive stars. Our simulations show that the constraint on the relative parallax precision needs to be lowered by 95% to obtain information on the BH-formation scenario. We therefore suggest that, the most conservative the constraint can be in DR4 is n × (ϖ/σϖ)DR3,single >  n × 1000/Pday, with (ϖ/σϖ)DR3,single the relative parallax precision for the single source solution in DR3 and n the DR4 improvement of the relative precision compared to the DR3 single star solutions.


Acknowledgments

SJ acknowledges support from the FWO PhD fellowship under project 11E1721N. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement ndeg 772225/MULTIPLES). TS acknowledges support from the European Union’s Horizon 2020 under the Marie Skłodowska-Curie grant agreement No 101024605. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

References

  1. Abdul-Masih, M., Banyard, G., Bodensteiner, J., et al. 2020, Nature, 580, E11 [Google Scholar]
  2. Andrews, J. J., Breivik, K., & Chatterjee, S. 2019, ApJ, 886, 68 [NASA ADS] [CrossRef] [Google Scholar]
  3. Andrews, J.~J., Taggart, K., & Foley, R. 2022, AAS J., submitted [arXiv:2207.00680] [Google Scholar]
  4. Bodensteiner, J., Shenar, T., Mahy, L., et al. 2020, A&A, 641, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Breivik, K., Chatterjee, S., & Larson, S. L. 2017, ApJ, 850, L13 [CrossRef] [Google Scholar]
  6. Breivik, K., Chatterjee, S., & Andrews, J. J. 2019, ApJ, 878, L4 [NASA ADS] [CrossRef] [Google Scholar]
  7. Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Chiavassa, A., Kudritzki, R., Davies, B., Freytag, B., & de Mink, S. E. 2022, A&A, 661, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Corral-Santana, J. M., Casares, J., Muñoz-Darias, T., et al. 2016, A&A, 587, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. de Mink, S. E., Sana, H., Langer, N., Izzard, R. G., & Schneider, F. R. N. 2014, ApJ, 782, 7 [Google Scholar]
  11. El-Badry, K., & Burdge, K. B. 2022, MNRAS, 511, 24 [NASA ADS] [CrossRef] [Google Scholar]
  12. El-Badry, K., & Quataert, E. 2020, MNRAS, 493, L22 [NASA ADS] [CrossRef] [Google Scholar]
  13. El-Badry, K., & Quataert, E. 2021, MNRAS, 502, 3436 [NASA ADS] [CrossRef] [Google Scholar]
  14. El-Badry, K., & Rix, H.-W. 2022, MNRAS, 515, 1266 [NASA ADS] [CrossRef] [Google Scholar]
  15. El-Badry, K., Rix, H.-W., Quataert, E., et al. 2023, MNRAS, 518, 1057 [Google Scholar]
  16. Everall, A., Boubert, D., Koposov, S. E., Smith, L., & Holl, B. 2021, MNRAS, 502, 1908 [CrossRef] [Google Scholar]
  17. Frost, A. J., Bodensteiner, J., Rivinius, T., et al. 2022, A&A, 659, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Fu, J.-B., Gu, W.-M., Zhang, Z.-X., et al. 2022, ApJ, 940, 126 [NASA ADS] [CrossRef] [Google Scholar]
  19. Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Gaia Collaboration (Arenou, F., et al.) 2023a, A&A, in press, https://doi.org/10.1051/0004-6361/202243782 [Google Scholar]
  21. Gaia Collaboration (Vallenari, A., et al.) 2023b, A&A, in press, https://doi.org/10.1051/0004-6361/202243940 [Google Scholar]
  22. Gomel, R., Mazeh, T., Faigler, S., et al. 2023, A&A, in press, https://doi.org/10.1051/0004-6361/202243626 [Google Scholar]
  23. Halbwachs, J.-L., Pourbaix, D., Arenou, F., et al. 2023, A&A, in press, https://doi.org/10.1051/0004-6361/202243969 [Google Scholar]
  24. Hirai, R., & Mandel, I. 2021, PASA, 38, e056 [NASA ADS] [CrossRef] [Google Scholar]
  25. Irrgang, A., Geier, S., Kreuzer, S., Pelisoli, I., & Heber, U. 2020, A&A, 633, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Janssens, S., Shenar, T., Sana, H., et al. 2022, A&A, 658, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Jayasinghe, T., Rowan, D.~M., Thompson, T.~A., Kochanek, C.~S., & Stanek, K. Z. 2022, MNRAS, submitted [arXiv:2207.05086] [Google Scholar]
  28. Langer, N., Schürmann, C., Stoll, K., et al. 2020, A&A, 638, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Lennon, D. J., Maíz Apellániz, J., Irrgang, A., et al. 2021, A&A, 649, A167 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Linder, N., Rauw, G., Martins, F., et al. 2008, A&A, 489, 713 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Liu, J., Zhang, H., Howard, A. W., et al. 2019, Nature, 575, 618 [Google Scholar]
  32. Liu, J., Zheng, Z., Soria, R., et al. 2020, ApJ, 900, 42 [Google Scholar]
  33. Mahy, L., Almeida, L. A., Sana, H., et al. 2020, A&A, 634, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Mahy, L., Sana, H., Shenar, T., et al. 2022, A&A, 664, A159 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Mandel, I., & Farmer, A. 2022, Phys. Rep., 955, 1 [NASA ADS] [CrossRef] [Google Scholar]
  36. Martins, F., Schaerer, D., & Hillier, D. J. 2005, A&A, 436, 1049 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Mashian, N., & Loeb, A. 2017, MNRAS, 470, 2611 [NASA ADS] [CrossRef] [Google Scholar]
  38. Mazeh, T., & Faigler, S. 2020, MNRAS, 498, L58 [NASA ADS] [CrossRef] [Google Scholar]
  39. Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15 [Google Scholar]
  40. Pantaleoni González, M., Maíz Apellániz, J., Barbá, R. H., & Reed, B. C. 2021, MNRAS, 504, 2968 [CrossRef] [Google Scholar]
  41. Rainot, A., Reggiani, M., Sana, H., et al. 2020, A&A, 640, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Rivinius, T., Baade, D., Hadrava, P., Heida, M., & Klement, R. 2020, A&A, 637, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Romagnolo, A., Olejak, A., Hypki, A., Wiktorowicz, G., & Belczynski, K. 2022, A&A, 667, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Sana, H., Le Bouquin, J. B., Lacour, S., et al. 2014, ApJS, 215, 15 [Google Scholar]
  45. Saracino, S., Kamann, S., Guarcello, M. G., et al. 2022, MNRAS, 511, 2914 [NASA ADS] [CrossRef] [Google Scholar]
  46. Sen, K., Xu, X. T., Langer, N., et al. 2021, A&A, 652, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Shahaf, S., Mazeh, T., Faigler, S., & Holl, B. 2019, MNRAS, 487, 5610 [NASA ADS] [CrossRef] [Google Scholar]
  48. Shahaf, S., Bashi, D., Mazeh, T., et al. 2023, MNRAS, 518, 2991 [Google Scholar]
  49. Shenar, T., Bodensteiner, J., Abdul-Masih, M., et al. 2020, A&A, 639, L6 [EDP Sciences] [Google Scholar]
  50. Shenar, T., Sana, H., Mahy, L., et al. 2022, Nat. Astron., 6, 1085 [NASA ADS] [CrossRef] [Google Scholar]
  51. Shikauchi, M., Kumamoto, J., Tanikawa, A., & Fujii, M. S. 2020, PASJ, 72, 45 [NASA ADS] [CrossRef] [Google Scholar]
  52. Shikauchi, M., Tanikawa, A., & Kawanaka, N. 2022, ApJ, 928, 13 [NASA ADS] [CrossRef] [Google Scholar]
  53. Stevance, H. F., Parsons, S. G., & Eldridge, J. J. 2022, MNRAS, 511, L77 [NASA ADS] [CrossRef] [Google Scholar]
  54. Walter, R., Lutovinov, A. A., Bozzo, E., & Tsygankov, S. S. 2015, A&ARv, 23, 2 [Google Scholar]
  55. Wiktorowicz, G., Wyrzykowski, Ł., Chruslinska, M., et al. 2019, ApJ, 885, 1 [NASA ADS] [CrossRef] [Google Scholar]
  56. Yalinewich, A., Beniamini, P., Hotokezaka, K., & Zhu, W. 2018, MNRAS, 481, 930 [NASA ADS] [CrossRef] [Google Scholar]
  57. Yamaguchi, M. S., Kawanaka, N., Bulik, T., & Piran, T. 2018, ApJ, 861, 21 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Masses, effective temperatures, and radii

Table A.1 shows the literature data obtained for BSGs.

Table A.1.

Masses, effective temperatures, and radii of BSGs.

Appendix B: Using a 10% versus a 1% core-H fraction for the BSGs

Here, we show the mass-magnitude relation for BSG stars derived as described in Sect. 3.1 using a 1% core-H fraction instead of 10% (Figs. B.1 and B.2). The fit parameters are listed in Table B.1. We also show an example of the difference in theoretical AMRF curves between these two cases (Fig. B.3).

Table B.1.

Fit parameters for the mass-magnitude relation of BSGs in different mass regimes using a 1% core-H fraction.

thumbnail Fig. B.1.

Same as Fig. 2 but using a 1% core-H fraction.

thumbnail Fig. B.2.

Same as Fig. 3 but using a 1% core-H fraction. The absolute G magnitudes using a 10% core-H fraction are also shown as a reference.

thumbnail Fig. B.3.

Theoretical AMRF curves for a 10 M (left) and a 20 M (right) primary BSG star obtained from the mass-magnitude relation using a 10% (blue) or a 1% (orange) core-H fraction. The different classes in both cases are also marked.

Appendix C: Identification fractions for a grid of triples

Figures C.1 and C.2 show the total identification fractions for the triples investigated in Sect. 4.2, but for more parameters, such as the eccentricity, e.

thumbnail Fig. C.1.

Same as Fig. 6, but for more parameters.

thumbnail Fig. C.2.

Same as Fig. 7, but for more parameters.

All Tables

Table 1.

Fit parameters for the mass-magnitude relation (G = alog(M/M)+b) of BSGs in different mass regimes.

Table 2.

Detection, identification, and total identification fractions for different kinds of systems.

Table 3.

Estimated number of sources detected and identified.

Table A.1.

Masses, effective temperatures, and radii of BSGs.

Table B.1.

Fit parameters for the mass-magnitude relation of BSGs in different mass regimes using a 1% core-H fraction.

All Figures

thumbnail Fig. 1.

Different theoretical AMRFs for three different primary masses: 10 M (blue), 20 M (orange), and 30 M (grey). The solid black line shows the theoretical AMRF for MS+BH systems, the dotted curves are for OB+MS binaries, and dashed curves are for triple MS dwarf systems. Different classes are also indicated for the 30 M primary by the horizontal dotted lines at the maxima of the triple AMRF and the binary AMRF (Class III: identifiable OB+BH systems, Class II: triple MS or OB+BH systems, Class I: OB+MS binaries, triple MS systems, or OB+BH binaries).

In the text
thumbnail Fig. 2.

Positions in the HRD for BSGs. Cyan dots are luminosity-class I stars from Martins et al. (2005) and blue stars are the positions on the tracks corresponding to a 10% core-H fraction. The dashed line represents a second-degree interpolation between the data from Martins et al. (2005) and the four lowest-mass points with a 10% core-H fraction (see text for more details).

In the text
thumbnail Fig. 3.

Absolute G-band magnitudes for BSGs in the mass range 8.9–66.89 M obtained from Martins et al. (2005) and Brott et al. (2011) (see Table A.1). The solid red line shows the fit and the vertical grey lines indicate the different mass regimes. The masses are shown on a logarithmic scale. As a reference, the obtained magnitudes for MS dwarf stars from Paper I are also shown.

In the text
thumbnail Fig. 4.

Different theoretical AMRF curves for two different primary masses: 10 M (left) and 30 M (right). In addition to the MS curves, AMRF curves for BSGs are also shown for the same primary masses. The different classes are now also indicated for both the OBs and the BSGs.

In the text
thumbnail Fig. 5.

Different theoretical AMRF curves for a 20 M primary. Triple curves for different mass ratios of the inner binary q1 are shown for the case where the inner binary is the most luminous (‘Inner’) or the outer binary is more luminous (‘Outer’).

In the text
thumbnail Fig. 6.

Total identification fraction for triple systems where the outer star is the most luminous component and the inner binary has q1 = 1. In the left panels, M1a is the mass of one of the components in the inner binary. The top panels are for OB systems, and the bottom panels are for systems with a BSG outer star.

In the text
thumbnail Fig. 7.

Same as Fig. 6, but when using the BSG curves.

In the text
thumbnail Fig. 8.

Total fraction of identified systems when using the MS curves for single-degenerate binaries with a BH (top row), non-degenerate binaries (middle row), and triples from the SMASH+ (bottom row). Left and right panels are for OB and BSG primaries, respectively.

In the text
thumbnail Fig. 9.

Same as Fig. 8, but when using the BSG curves.

In the text
thumbnail Fig. 10.

Histogram of ϖ/σϖ of the massive sources in the ALS II. The vertical lines indicate a relative precision on ϖ of 200, 40, and 20.

In the text
thumbnail Fig. B.1.

Same as Fig. 2 but using a 1% core-H fraction.

In the text
thumbnail Fig. B.2.

Same as Fig. 3 but using a 1% core-H fraction. The absolute G magnitudes using a 10% core-H fraction are also shown as a reference.

In the text
thumbnail Fig. B.3.

Theoretical AMRF curves for a 10 M (left) and a 20 M (right) primary BSG star obtained from the mass-magnitude relation using a 10% (blue) or a 1% (orange) core-H fraction. The different classes in both cases are also marked.

In the text
thumbnail Fig. C.1.

Same as Fig. 6, but for more parameters.

In the text
thumbnail Fig. C.2.

Same as Fig. 7, but for more parameters.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.