© The Authors 2024

1. Introduction

Traditionally, it is difficult to observe non-accreting black holes (BHs) because they do not emit in the electromagnetic spectrum. The first indirect measurements of stellar mass BHs are from X-ray binaries, where hot material is accreted onto an unseen but very compact object (e.g. Webster & Murdin 1972). These binaries are well observed due to their high-energy emission (e.g. Remillard & McClintock 2006, and references therein). The compact object itself is usually hidden behind the strong emission of the accreting material, which can outshine the donating companion as well. Additionally, the best-studied stellar X-ray sources are traditionally found to be Galactic sources (Casares & Jonker 2014).

In recent years, observations of gravitational waves showed a broader existence of BHs in binaries throughout the Universe. These observations question our understanding of BHs because of the large variety of BH masses (Abbott et al. 2019, 2021, 2023a,b, 2024), while we infer only moderate BH masses for X-ray sources (Özel et al. 2010; Farr et al. 2011).

Most recently, Gaia DR3 data has been searched for unseen massive companions to stars in the Milky Way (MW). These investigations revealed some candidates, named Gaia BH1 (Chakrabarti et al. 2023; El-Badry et al. 2023a), Gaia BH2 (El-Badry et al. 2023b; Tanikawa et al. 2023), and Gaia BH3 (Gaia Collaboration 2024, using pre-release DR4 data), as well as other possible candidates of unseen compact objects (e.g. Andrews et al. 2022; Wang et al. 2024). Follow-up radial velocity measurements have verified three individual Gaia BHs. Gaia BH3 (Gaia Collaboration 2024) is rather peculiar among BHs in the MW in terms of its mass, but within expectations from gravitational waves. It has a mass of 32.7 ± 0.82 M_⊙ and is thought to be formed within a tidally disrupted globular cluster at lower metallicity (Balbinot et al. 2024; El-Badry 2024a; Iorio et al. 2024; Marín Pina et al. 2024). In contrast, Gaia BH1 and BH2 have masses of 9.62 ± 0.18 M_⊙ and 8.94 ± 0.34 M_⊙, respectively, and they both lie in the MW disk with companion stars of roughly solar metallicity. Some other formation scenario is required.

From classical binary evolution, binaries consisting of a BH with a stellar companion in orbits up to a few thousand days are expected as shown in previous studies (Breivik et al. 2017; Yalinewich et al. 2018; Chawla et al. 2022). The binaries in these studies favoured companion masses above 1 M_⊙ (Mashian & Loeb 2017; Yamaguchi et al. 2018). Depending on their mass, OB or Wolf-Rayet (WR) star companions can stably accrete mass upon Roche-lobe overflow from the BH progenitor (Shenar et al. 2019; Langer et al. 2020). Even when this mass transfer becomes unstable, these binaries are more likely to survive. Low-mass companions, on the other hand, such as those in Gaia BH1 and BH2, will most certainly experience unstable mass transfer and either merge within a common envelope, because they lack the orbital energy to successfully eject a common envelope, or result in tight orbits, < 100 days, like the progenitors of low-mass X-ray binaries (e.g. Podsiadlowski et al. 2003; Fragos & McClintock 2015). Previous studies have suggested that Gaia BH1 and BH2 were formed dynamically in an open cluster (Rastello et al. 2023; Di Carlo et al. 2024; Tanikawa et al. 2024). However, this is not the only possible formation scenario. In their discovery paper for Gaia BH1, El-Badry et al. (2023a) proposed four other possible scenarios, including classical binary evolution involving at least one mass-transfer phase, triple system dynamics (e.g. Generozov & Perets 2024), and formation in a hierarchical triple (e.g. Li et al. 2024). We explore one of these possible scenarios in depth: The idea that the progenitors of Gaia BH1 and BH2 were initially very massive stars (M > 50 M_⊙), as recently studied by Gilkis & Mazeh (2024). The wind mass-loss rates of these massive stars are so extreme on the main sequence that the stars lose their entire envelopes prior to core hydrogen exhaustion, and they therefore never expand onto the giant branch. Hence, these systems avoid mass transfer. Subsequent wind mass loss during their evolution as WR stars further reduces their initially high masses (see, e.g., Bavera et al. 2023; Vink et al. 2024). The resulting BHs are of sufficiently low mass that they could potentially explain the BHs in Gaia BH1 and BH2.

Our approach is to use the POSYDON binary population synthesis framework (Fragos et al. 2023; Andrews et al. 2024) to evolve a population of stellar binaries with extreme mass ratios. The primary star is sufficiently massive that it never expands into a giant star, while the companion is ≲1 M_⊙ to explain the luminous companions to Gaia BH1 and BH2. In Sect. 2 we describe our adaptations to POSYDON and associated single-star models calculated with MESA, made specifically to model Gaia BHs. In Sect. 3 we provide our results from our single-star model grids and binary population synthesis models. Finally, in Sect. 4 we discuss our results and provide some conclusions.

2. Methods

All results presented in this paper were derived using the binary population synthesis framework POSYDON (Fragos et al. 2023; Andrews et al. 2024). POSYDON uses extensive detailed single- and binary-star model grids, calculated with the MESA stellar structure and binary evolution code (Paxton et al. 2011, 2013, 2015, 2018, 2019), to self-consistently follow the entire evolution of both stars in a stellar binary. Specifically, we employed the POSYDON infrastructure to calculate additional single-star grids with MESA. We analysed them with POSYDON tools and used POSYDON to perform population synthesis calculations.

2.1. Single star models

In addition to the default single-stars grids in POSYDON v2 (Andrews et al. 2024), which cover different metallicities, we generated several additional models at solar metallicity (we use Z_⊙ = 0.0142 following Asplund et al. 2009) with variations in our default wind prescriptions. We present them in Sect. 3.3. POSYDON modifies the default wind scheme in MESA and is described in Appendix A.

While multiple wind prescriptions can affect stars throughout their lifetimes depending on their evolutionary state, we found that two winds dominate for massive stars. We call them Vink and WR winds throughout the paper. The Vink wind (Vink et al. 2001) dominates the early evolution of the massive stars during the core-H burning, while the WR wind (Nugis & Lamers 2000) is usually optically thick and kicks in after the star loses most of its H-rich envelope. Therefore, we tested models in which both of these winds were altered with multiplicative factors. Separately, winds are well known to be affected by metallicity. Therefore, although we focused on models with Z_⊙ (observations suggest that the companions to Gaia BH1 and BH2 have ≲Z_⊙), we additionally show models with 0.45 Z_⊙ and 2 Z_⊙ without additional modifications in Sect. 3.3.

2.2. Binary population synthesis models

To compare our single-star models to an observed population, we used POSYDON to simulate the intrinsic population of non-interacting binaries with a low-mass star (0.5–2.0 M_⊙) orbiting a BH. We modified the default setup of POSYDON by adding a user module¹ to allow for the modelling of binaries with an extreme initial mass ratio (mass ratio < 0.05). Specifically, we modified the ‘flow chart’ to start our simulations in the so-called detached step (solving a simplified set of ordinary differential equations to model the effects of wind mass loss, tides, magnetic braking, and gravitational radiation on the orbital and rotation parameters; for more details see Sect. 8.1 of Fragos et al. 2023) to make use of our various single-star grids and allow initial non-zero eccentricity at the zero-age main sequence (ZAMS).

Our population synthesis models contained one million binaries each. These binaries were generated in a starburst and weighted in a post-processing step to resemble the MW, assuming a constant star formation history. For the results presented in Sect. 3.4, we used a thermal eccentricity distribution at the ZAMS². To focus on potential progenitors of Gaia BH1- and BH2-like systems, we only sampled primary masses ∈[20, 150] M_⊙ and secondary masses ∈[0.5, 2.0] M_⊙. Since we lack robust observational constraints for these extreme mass ratios (for example, see Moe & Di Stefano 2017), we adopted a canonical flat mass ratio distribution. We stopped the evolution of all binaries that initiated Roche-lobe overflow mass transfer. For our extreme mass ratios chosen at the ZAMS, it is expected that mass transfer always becomes dynamically unstable, leading to a merger during the common-envelope phase. Compared to the POSYDON default², we used zero kicks for simplicity and the delayed prescription from (Fryer et al. 2012) to form BHs at all masses. In this prescription, a proto-neutron star is first created, and the fallback of envelope material onto it determines the mass of the remaining compact remnant (see Sect. 8.3 of Fragos et al. 2023, for details). Finally, we allowed for a wide range of initial orbital periods ∈[0.75, 6000] days to cover binaries with and without interaction at any of our initial masses.

3. Results

First, we highlight some relevant properties of our single-star models calculated in MESA, including those where we varied our default wind prescriptions. We then present the binary population synthesis results based on these single-star models.

3.1. Maximum radius of single massive stars

Massive stars have larger radii than their less massive counterparts on the ZAMS. While very low-mass stars are often referred to as dwarfs, very massive stars are giants early on in their evolution. During their evolution, stars usually expand, especially past the main sequence. Most of these expansion phases are a consequence of a mirroring effect that burning shells have through their stabilising effect on the burning itself. The contraction of the core material leads to an expansion of the outer envelope layers above the burning region (for a review see e.g. Lamers & Levesque 2017).

For very massive stars, there is another effect which drives the expansion. The high luminosity produced by vigorous nuclear burning pushes the outer layers of a star to larger radii. At the surface of these massive stars, this luminosity launches enhanced winds. These winds remove the outer stellar material, leading to a movement of the stellar photosphere to a lower radial coordinate. The combination of these two effects settles to an equilibrium radial coordinate for the stellar surface. In the case of very massive stars, especially at higher metallicity, the winds are so strong that material that would have led to an expansion of the stellar surface instead escapes from the star. This wind causes these stars to remain at a radius not far above the initial radius at the ZAMS (for a review see e.g. Lamers & Cassinelli 1999).

Figure 1 shows the maximum radial extension of a star of a given ZAMS mass during its evolution (limited to the age of the Universe ∼13.8 Gyr, which causes a gap at very low masses, where stars would need more time to become giants) as a function of its initial mass. Depending on the mass at the ZAMS, different physical effects are at work to determine the maximum radius. First, the more massive a star, the larger its radius. Second, the stars that evolve through shell burning will significantly expand during contraction phases of the core. These phases lead to the tip of the red giant branch. Subsequent burning phases can cause even larger radii on the asymptotic giant branch. This difference in the giant branches is the first and small step on the very left of Fig. 1 at a few M_⊙ and differentiates between a few hundred R_⊙ and several hundred or thousand R_⊙. This feature exists at all metallicities, while a higher abundance of heavier elements leads to a larger radius because the resulting opacity is higher. Depending on the metallicity, for stellar masses in the range from 30 to several hundreds of M_⊙, the winds become so strong that the maximum radius is reached before the shell burning can lead to a maximum radius > 1000 R_⊙. The most massive stars reach their maximum extension already before they complete core-H burning. Transitions between wind prescriptions lead to several smaller features in Fig. 1. For the Gaia BHs to form without mass transfer, the BH progenitor needs a maximum radius below a few hundred R_⊙.

Fig. 1. Maximum radius of single stars as a function of their initial mass. The colour and symbol show different metallicities.

3.2. Final mass of single massive stars

Fig. 2 shows the final compact object mass using a simple direct collapse model to form the BH (alternative supernova prescriptions are shown in Appendix B). When we compare the ZAMS mass to the final BH mass, Fig. 2 exhibits a peak in the remnant mass that varies as a function of metallicity. Variations in the location of this peak are due to the amount of mass that is lost prior to collapse, which varies with metallicity. At high metallicity, the stellar winds are so strong that large parts, if not all, of the stellar envelope are lost during the stellar evolution, leading to relatively lower-mass BHs (Bavera et al. 2023; Vink et al. 2024). However, with their weaker winds, lower-metallicity stars remain more massive and can enter the regime of pair instability, which ejects the envelope or even disrupts the whole star so that no remnant is left (zero mass in Fig. 2). At metallicities below one-fifth of the solar value, the early winds are too weak to strip massive stars and enable the WR winds. Here, strong winds connected to luminous blue variables have a similar effect, allowing the existence of a plateau in the BH mass above the pair-instability regime (see Fig. 2).

Fig. 2. Remnant mass function of single stars depending on their initial mass. The colour and symbol show different metallicities. We used the direct collapse prescription in POSYDON. At Z ≤ 0.2 Z⊙ stars with initial masses above 100 M⊙ might not form a BH because of pair instability and are shown here with a zero mass. Stars with initial masses above the pair-instability gap form BHs outside the plotted range, with masses above 100 M⊙ at the lowest shown metallicity.

In the remainder of the paper, we concentrate on the stellar winds. Winds produce a plateau in the BH mass over a wide range of high ZAMS masses. The typical BH mass at this plateau depends on the strength of the wind and therefore indirectly on metallicity, as shown in Fig. 2. Although the plateau is above 20 M_⊙ for less than half the solar metallicity, plateau BHs would have typical masses below 10 M_⊙ for twice the solar metallicity for the default wind prescription in the MESA models from POSYDON v2. Since Gaia BH1 and BH2 are slightly sub-solar, as are most systems in the MW, we focus on solar-metallicity models below.

3.3. Variations in the wind for massive stars

During the evolution of a massive star at high metallicity, there are two main wind regimes we wish to investigate in more detail. These are the Vink and the WR winds (see Sect. 2.1 and Appendix A for more details).

The top panel of Fig. 3 shows that the Vink wind strength determines the position of the maximum remnant mass peak. The stronger the Vink wind, the lower the maximum mass and the lower the ZAMS mass leading to the peak value. In addition to the mass, the Vink wind determines when the star turns to become bluer, and it therefore defines its maximum radius. Sufficiently strong winds remove the step at about 500 R_⊙, and this characterises the predominant difference between hot and cold winds in the bottom panel of Fig. 3. This step feature is also modified by other transitions of wind prescriptions, for instance, stronger winds when the star tries to enter the regime of luminous blue variables. Enhancements of the Vink winds have a minor effect on the most massive stars that characterise the plateau region because these stars never evolve to be redder than the main sequence. Their maximum radius is therefore very close to their ZAMS radius. A decrease in the Vink winds (like suggested by Smith 2014; Gormaz-Matamala et al. 2022) would have effects similar to approaching lower metallicity (see Sects. 3.1 and 3.2).

Fig. 3. Same as Fig. 2 (top) and Fig. 1 (bottom), but here we show different variations of the wind at solar metallicity (see legend). The dotted grey line shows the ZAMS radius, where the discontinuity at about 130 M⊙ is a numerical artefact in the ZAMS models due to the detection criterion of ZAMS. The WR wind variations all overlap in the bottom panel.

The WR wind does not affect the ZAMS mass that leads to the peak in BH mass, but it can lower the resulting BH mass at the peak slightly. None of the stars that keep their H-rich envelope throughout their evolution experiences this type of wind. Additionally, the WR wind causes the formation of the plateau due to its luminosity scaling. Therefore, the typical BH mass in the plateau regime strongly depends on the WR wind prescription adopted. We point out that the maximum radius of a star is typically reached before this wind kicks in, and it therefore has no effect on whether binary interactions are triggered (see the bottom panel of Fig. 3).

In Fig. 4 we show the variation in the plateau mass with different wind prescriptions and metallicities. We defined the plateau as the mass range between the minimum compact-object mass at higher ZAMS masses than the peak and the maximum at even higher masses than the aforementioned minimum, which usually coincides with the MESA model with the highest ZAMS mass. The plateau mass depends only weakly on the Vink wind that is applied during the main sequence of massive stars. The optically thick WR wind that is applied after the H-rich envelope is lost has a clear impact on the typical BH mass in the plateau regime (see Fig. 4). Additionally, its strength causes the plateau mass range to become narrower. For a plateau in the mass range of the Gaia BHs, our simulations of solar metallicity would require a boost of the WR winds by a factor in the range [1.15, 1.35] depending on the supernova prescription. Some trends are slightly modified for different prescriptions for the BH formation (see Appendix B).

Fig. 4. Range of BH masses in the plateau regime depending only on the metallicity in Z⊙ (orange, 0.45 and 2 Z⊙), Vink wind factor (green), or WR wind factor (blue), while the others are at default. The default at solar metallicity and without any wind modification is shown in black.

3.4. Population models of non-interacting binaries

Using our adapted version of POSYDON, we performed several population synthesis simulations to inspect the distribution of detached BH binary properties with a low-mass companion in wide-period orbits, such as the Gaia BH systems. Figure 5 shows an example intrinsic population, without accounting for any Gaia selection effects, using the setup outlined in Sect. 2. The observed Gaia BHs are over-plotted in colour. We did not try to reproduce Gaia BH3 here because it clearly formed at lower metallicity, which is required to form BHs of this mass without binary interactions. Nevertheless, we included Gaia BH3 in Fig. 5 to show that systems with similar orbital parameters and only lower BH mass are formed. Hence, it is expected that Gaia BH3 might have formed at lower metallicity without binary interactions as well. On the other hand, very low-mass BHs, neutron stars, or compact objects in the mass gap (like the one recently proposed by Wang et al. 2024) would require a very high metallicity or extremely strong winds to form without binary interactions.

Fig. 5. Corner plot showing the BH mass, the companion mass, the orbital period, and the eccentricity of binaries in a population synthesis run with an enhanced WR wind (factor = 1.2). While the grey histograms show the full population, the green histograms are limited to binaries with orbital periods up to 5000 days. The observationally determined values for the Gaia BHs are shown in colour (El-Badry et al. 2023a,b; Gaia Collaboration 2024). The histograms are weighted to represent the number of these binaries in the MW, which has large uncertainties as given by the unit $[1_{-1.0}^{+3.3}]$ (for details, see Appendix D).

A typical progenitor of Gaia BH2 started in our simulation with stars of 92.07 and 1.092 M_⊙ in an orbit of 15.13 days and an eccentricity of 0.4796 at the ZAMS. Prior to collapse, the primary loses most of its mass in winds that reduce its mass to 11.12 M_⊙, which forms a 9.446 M_⊙ BH. The mass loss due to winds and the BH formation causes the system to widen to 878.9 days prior and 1 275 days post collapse. While the eccentricity stays nearly unchanged during the progenitor evolution, it becomes 0.5168 after the BH formation. While the companion evolves, the system is observable for > 8 Gyr, where tides start to circularize the binary after the companion left its main sequence. Finally, the system becomes a wide low-mass X-ray binary.

For the population, a peak is caused in the BH mass function that originated from the plateau region (see BH masses ≲10 M_⊙ in the top left panel of Figs. 5 and C.1). Depending on the strength assumed for the kicks at BH formation, there can be a second, usually smaller, peak at the BH mass that corresponds to the lowest-mass stars (at initial mass of ≲50 M_⊙, leading to a ≈15 M_⊙ BH), which reach a maximum radius of about 500 R_⊙. Hence, this peak requires a longer orbital period than the one of the plateau. In the two-dimensional heat map (first column of the third row in Figs. 5 and C.1) between BH mass and orbital period, the two sub-populations are easily distinguished. None of the Gaia BHs can belong to the second sub-population, which requires a period of about 20 000 days.

The companion mass shows a clear preference of lower masses. This preference is simply caused by their longer lifetimes, which make lower-mass systems statistically more likely. Slightly below 1 M_⊙, this trend levels off because the stellar lifetime reaches the age of the Universe, which is the maximum allowed in our population synthesis. In observations, companions slightly below 1 M_⊙ are most likely to be found first because they are the brightest of the most common companions. Companions with masses of 2 M_⊙ are less common by more than an order of magnitude. The power-law behaviour suppresses companions with higher masses even more strongly, which results in a very small sample that will be very rare to observe in the MW.

The only other two-dimensional representation that shows a correlation is the representation between the orbital period and the eccentricity. The structure is a result of the mass loss at BH formation. The reported BHs observed by Gaia are at the lower end of the period distribution.

3.5. Detectability by Gaia

Although a detailed modelling and analyses of the observability of our sample by Gaia is beyond the scope of this paper, it is instructive to examine the Gaia sensitivity function in the visible star mass-range we simulated. The vast majority of orbital companions with orbital periods > 10 days will be identified in the astrometric data of Gaia (see e.g. El-Badry 2024b). We illustrate the expected Gaia astrometric sensitivity function in Fig. 6 for 5 yr (∼DR4) and 10 yr (∼DR5), showing the sky-averaged detectability horizon for unseen objects on a circular orbit around companion stars with masses of 0.65, 1.02, and 1.59 M_⊙ (following Appendix A of Holl et al. 2022). As expected, the sensitivity is higher (i.e. a more distant detection horizon) for longer periods, but limited to roughly the time-span of observations of 1000, 2000, and 4000 d for DR3, DR4, and DR5, respectively. Comparing this to the orbital period distribution of Fig. 5, we see that with each succeeding data release, Gaia will cover increasingly more of the expected orbital period distribution. Thus, each release is expected to increase the number of Gaia BHs. Additionally, each successive data release will double the number of data points, which increases the overall signal-to-noise ratio, as shown by the offset between the 5 and 10 yr lines of Fig. 6.

Fig. 6. Gaia detection limits for a non-luminous companions of 10 M⊙ on a circular orbit around a F0V, G2V, and K2V star. The lines represent the (sky-uniform) median distance limit at which the Δχ2 = 100 between a single star and a circular Keplerian orbital model fit, corresponding to (typically) well-constrained orbital fits. Estimates were provided for the 5 yr nominal mission (∼DR4) and 10 yr extended mission (∼DR5) as a function of period. Successful detections drop quickly for orbital periods longer than the data time-span, and the graphs are therefore truncated at 5 and 10 yr. At the top, the semi-major axis (a) is specified for a G2V star, although the difference for the other masses would be hardly noticeable. On the right side, we provide a rough estimate of the MW coverage for a specific distance limit, assuming a disk of height 1.15 kpc and radius 15 kpc, and we therefore obtain a full coverage at 23 kpc. All detection values assume zero extinction.

We note that the eccentricity distribution of the systems shown in Fig. 5 is far from zero and has a typical value of 0.4 − 0.8, like the first Gaia BHs. As illustrated in Figure A.1 of Holl et al. (2022), a higher eccentricity reduces the detectability slowly: with a distance limit drop of about 10% at e = 0.5 and only exceeding a drop of more than 20% beyond e = 0.7. This and Fig. 6 describe the all-sky median behaviour, while in reality, the detectability as function of all astrometric orbital parameters can differ strongly as a function of the sky position due to the Gaia scanning law.

4. Discussion and conclusions

We showed that BHs found in wide binaries, such as Gaia BH1 and Gaia BH2, can form via a non-interacting isolated binary evolution channel at higher solar-like metallicities. The clearest signature of the formation channel we presented is a preferred BH mass in wide binaries with low-mass companions at a given metallicity. Additionally, if the eccentricity distribution of systems like the Gaia BHs is similar to samples of young systems, a formation channel with weak interaction (e.g. tides) and low kicks might be favoured. The verification that this channel is dominant requires a significantly larger observational sample.

If it can be demonstrated that the Gaia BHs formed without binary interaction, the masses of the BH and the metallicity of the companion can be used to constrain the winds of massive stars, especially the WR wind. We find a best match of Gaia BHs for a boosting factor of about 1.2 on the WR winds. Nevertheless, there are uncertainties that are not covered by our simple models, for instance, a different scaling of the winds, or the impact of the BH formation prescription. The latter is less important for BHs plateaus above 10 M_⊙ (see Appendix B). These degeneracies can only be overcome with a large sample of wide binaries containing a BH and a low-mass companion. Finally, we only measured the location at which the final mass settles. We therefore cannot distinguish between stronger mass loss during a shorter period of time versus moderate mass loss over a longer time period, that is, different prescriptions for BH formation favour different wind strengths to accumulate at the same mass range.

Our simulations predict a second smaller peak in the BH mass distribution of wide non-interacting BH binaries (see ≈15 M_⊙ in Fig. 5). This local overabundance should always coincide with periods above 1000 days. There, the progenitors are stars with ZAMS masses below the mass of the progenitors that form the most massive BHs at a given metallicity. These stars experienced a very short phase at most with WR winds in their evolution. Thus, the final BH mass from this sub-population that causes the second peak probes the early winds and radial expansions of stars. However, the weak contribution of these secondary sub-populations and the long periods of their binaries make it unlikely that this final BH mass can be directly constrained by Gaia observations.

Other modifications such as enhanced overshooting (Gilkis & Mazeh 2024) have effects very similar to our changes in the Vink wind. Changes in the overshooting or other mixing processes can enhance the luminosity on the main sequence and therefore cause stroger winds as well. As shown in Fig. 4 and Gilkis & Mazeh (2024), rather extreme changes are required to obtain BHs in the reported mass range of the first two Gaia BHs. No such strong push of the main sequence towards higher luminosities has been observed so far.

It is very difficult to compare this channel to other binary formation channels, which probably coincide in nature, because of the extreme mass ratio it requires already early on the ZAMS. In this regime, observational samples are notoriously incomplete, which precludes reliable estimates of the prevalence of initial binaries at the ZAMS with the extreme mass ratios required to form Gaia BH-like binaries. Nevertheless, adopting a flat mass-ratio distribution, which is a pure extrapolation, we estimate that hundreds of systems like Gaia BH1 and BH2 probably exist in the MW (see the bottom rows of Table D.1). Furthermore, our models predict many more systems at longer orbital periods, which are difficult to find given the limited lifetime of Gaia. Nevertheless, we expect many more systems to be found in future Gaia data releases, which will cover a longer time baseline.

Acknowledgments

This work was supported by the Swiss National Science Foundation (project number PP00P2_211006). The POSYDON project is supported primarily by two sources: the Swiss National Science Foundation (PI Fragos, project numbers PP00P2_211006 and CRSII5_213497) and the Gordon and Betty Moore Foundation (PI Kalogera, grant award GBMF8477). J.J.A. acknowledges support for Program number (JWST-AR-04369.001-A) provided through a grant from the STScI under NASA contract NAS5-03127. M.B. acknowledges support from the Boninchi Foundation. K.K. acknowledges support from the Spanish State Research Agency, through the María de Maeztu Program for Centers and Units of Excellence in R&D, No. CEX2020-001058-M. K.A.R. is also supported by the Riedel Family Fellowship and thanks the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant No. 1829740, the Brinson Foundation, and the Moore Foundation; their participation in the program has benefited this work. Z.X. acknowledges support from the China Scholarship Council (CSC). E.Z. acknowledges funding support from the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “3rd Call for H.F.R.I. Research Projects to support Post-Doctoral Researchers” (Project No: 7933).

References

Appendix A: Wind prescriptions

Here we summarise the wind prescriptions used in POSYDON’s default adoptions to winds in MESA.

We calculated a ‘cool’ wind for low photospheric temperatures, T_photo ≤ 12 000 K. If the central He abundance is smaller than 10⁻⁶ and the difference between He and C core is small, M_He, core − M_C, core < 0.1 M_⊙; hence for cool asymptotic giant branch stars, wed use the ‘Blocker’ wind scheme. Otherwise, for cool red giant branch stars, we use the ‘Dutch’ wind scheme. For T_photo ≥ 8 000 K we calculate a ‘hot’ wind using the ‘Dutch’ wind scheme (initial mass M_ini > 8 M_⊙) or ‘Reimers’ wind (M_ini ≤ 8 M_⊙). For T_photo < 8 000 K we use the cold wind only and for T_photo > 12 000 K we use the hot wind only. In the intermediate range we interpolate between the two winds linearly with T_photo.

In the ‘Blocker’ wind scheme, we use the maximum of 0.2 times the value of Eq. (2) in Bloecker (1995) and our ‘Reimers’ wind – including the scaling factor. In the ‘Dutch’ wind scheme (Glebbeek et al. 2009), for T_photo ≤ 10 000 K we use ‘de Jager’ wind, for T_photo ≥ 11 000 K we apply ‘Nugis & Lamers’ wind if the surface H abundance, X < 0.4 and ‘Vink’ wind scheme for larger H abundances, otherwise we linearly interpolate with T_photo.

In the ‘Reimers’ wind prescription, we adopt a wind mass loss rate according to Reimers (1975) with an additional scaling factor of 0.1. In the ‘de Jager’ wind prescription, we use the linear approximation of the empirical interpolation given in Eq. (1) of de Jager et al. (1988). In the ‘Vink’ wind prescription, we use Eqs. (24) and (25) of Vink et al. (2001) with Eqs. (14) and (15) to make the transition between the high- and low-temperature regime. Throughout this work, we modify the ‘Vink’ wind prescription by including a multiplicative factor to the final wind mass loss rate, see Sect. 3.3. In the ‘Nugis & Lamers’ wind prescription, we use Eq. (1) of Nugis & Lamers (2000), where we use the surface values of the He abundance and Z = 1 − X − Y. In our wind variations we refer to this wind as WR wind and use a simple multiplicative factor in front of the usual wind value, see Sect. 3.3.

If a star crosses the Humphreys-Davidson limit (L > 6 × 10⁵ L_⊙ and R R_⊙⁻¹ > 10⁵ L^−0.5 L_⊙^0.5, Humphreys & Davidson 1979, 1994), is not on the thermally pulsating asymptotic giant branch, and has a surface H abundance X > 0.1, then we apply luminous blue variable winds which is taken to be the maximum of 10⁻⁴ M_⊙yr⁻¹ and the otherwise calculated wind.

Appendix B: Other prescriptions for the compact object formation

In Figs. B.1 to B.3 we show the impact of alternative prescriptions for the BH formation. Following the delayed prescription of Fryer et al. (2012) looks very similar to the simple direct collapse model (Fig. 2 and top panel of Fig. 3), where only the lower compact object masses are shifted due to differences in the amount of fallback. On the other hand the prescription based on explodability (Patton & Sukhbold 2020) may form neutron stars instead of BHs for some models in the plateau range. It should be noted that it is unclear whether this prescription is valid for very massive stars which got highly stripped.

Fig. B.1. Same as Fig. 2 but showing the Fryer+12-delayed (top; Fryer et al. 2012) and Patton+Sukhbold20-engine (bottom; Patton & Sukhbold 2020) prescription for the compact object formation.

Fig. B.2. Same as the top panel of Fig. 3 but showing the Fryer+12-delayed (top; Fryer et al. 2012) and Patton+Sukhbold20-engine (bottom; Patton & Sukhbold 2020) prescription for the compact object formation.

Fig. B.3. Same as Fig. 4 but showing the Fryer+12-delayed (top; Fryer et al. 2012) and Patton+Sukhbold20-engine (bottom; Patton & Sukhbold 2020) prescription for the compact object formation. It should be noted that here the plateau spans over the mass gap between neutron stars and BHs for the strongest wind variations.

B.1. The maximum black hole mass

Figure B.4 shows how the maximum BH mass scales with the metallicity and additional factors on the Vink or WR wind. This is in agreement with the findings of wind variations on the maximum mass of a BH formed at solar metallicity by Bavera et al. (2023). Because the peak masses are always well above 10 M_⊙ the different prescriptions of BH formation in POSYDON lead to the same results.

Fig. B.4. Maximum BH mass depending on the metallicity in Z⊙ (orange), Vink-wind factor (green), or WR-wind factor (blue). The default at solar metallicity and without any wind modification is shown in black. The different symbols indicate the prescription for the BH formation.

B.2. The black hole progenitor

In Fig. B.5 we show the final profile at core carbon depletion similar to the example star in Sect. 3.4 as a typical case from the BH mass plateau. This profiles is used to determine the BH after collapse. There is no hydrogen left in the star and the winds have even stripped off most of the helium rich material. The final BH consists of all the carbon-oxygen core and a small amount of helium will fall back into the final remnant of about 9.4 M_⊙.

Fig. B.5. Final profile of an initially 92.42 M⊙ star with WR winds enhanced by a factor 1.2. The dashed black line is the radius, and the coloured lines are the mass fractions of hydrogen (X), helium (Y), and metals (Z).

Appendix C: Another population run

Figure C.1 show the results from a population synthesis run with several changes compared to Fig. 5. Those are: (1) using stronger WR winds (factor = 1.3), (2) using Patton+Sukhbold20-engine (Patton & Sukhbold 2020) for BH formation, (3) using kicks scaled by one over the BH mass, and (4) using zero eccentricity distribution at ZAMS.

Fig. C.1. Corner plot showing the BH mass, the companion mass, the orbital period, and the eccentricity of binaries in a default POSYDON population synthesis run with enhanced WR wind (factor = 1.3); see the text for more details. The grey histograms are the full population, and the green ones are limited to binaries with orbital periods of up to 5 000 days. The observationally determined values for the Gaia BHs are shown in colour (El-Badry et al. 2023a,b; Gaia Collaboration 2024). It should be noted that the histograms are weighted to represent the number of such binaries in the MW, which has large uncertainties as given by the unity unit $[1_{-1.0}^{+3.3}]$; for details see Appendix D.

With a different engine to calculate the final BH mass, a different wind strength is needed to reproduce the BH masses of Gaia BH1 and BH2. Still, Gaia BH3 is above the maximum mass of BHs formed at solar metallicity. Using kicks at BH formation leads to significant orbital changes when the BH is born. Hence, the choice of initial eccentricity distribution is less important and higher eccentricities are more common. Additionally, the kick leaves an imprint in the period-eccentricity plane. The other key findings of Sect. 3.4 hold here as well.

Appendix D: Normalising the populations

To get numbers of binaries in the MW from our simulations, we have to normalise our counts:

$$\begin{array}{c}\hfill N_{\mathrm{MW}}=f_{\mathrm{bin}}{f}_{\mathrm{mass}}{f}_{\mathrm{MW}}{N}_{\mathrm{sim}},\end{array}$$(D.1)

where f_bin = 0.7 is the assumed initial overall binary fraction³. For the initial primary mass distribution we follow Kroupa (2001), ξ(m)∝m^−α_i with

$$\begin{array}{ccc}& {\alpha }_0=0.3\pm 0.7\hfill & \hfill \text{for}0.01M_{\odot }\le m<0.08M_{\odot }\\ \hfill & {\alpha }_1=1.3\pm 0.5\hfill & \hfill \text{for}0.08M_{\odot }\le m<0.50M_{\odot }\\ \hfill & {\alpha }_2=2.3\pm 0.7\hfill & \hfill \text{for}0.50M_{\odot }\le m<1.00M_{\odot }\\ \hfill & {\alpha }_3=2.3\pm 0.7\hfill & \hfill \text{for}1.00M_{\odot }\le m.\end{array}$$(D.2)

The secondary mass we infer from a mass ratio distribution, p_q ∝ q^γ. This we assume to be flat, hence γ = 0 for all mass ratios. Thus we get a correction factor³

$$\begin{array}{c}\hfill f_{\mathrm{mass}}=\frac{{\int }_{20}^{150}{\int }_{0.5/m}^{2.0/m}\xi (m)p_q(q)\mathrm{d}q\mathrm{d}m}{{\int }_{0.01}^{200}{\int }_0^1\xi (m)p_q(q)\mathrm{d}q\mathrm{d}m}.\end{array}$$(D.3)

We obtain our uncertainty for the dominating factor by having different assumptions for the distribution of q < 0.1. First, there are no binaries with a low mass ratio, p_q = 0 ⇒ N_MW, min = 0. Second, all single stars are actually binaries with a low mass ratio, ∫₀^0.1p_q dq = 0.3 and f_bin = 1. Thus we have N_MW, max ≈ 4.3 N_MW.

We normalised to a MW like galaxy by comparing the simulated mass to the stellar mass of the MW, M_MW = 5 × 10¹⁰ M_⊙ (Hammer et al. 2007; Yin et al. 2009, where we optimistically assume that all binaries have formed at solar metallicity):

$$\begin{array}{c}\hfill f_{\mathrm{MW}}=\frac{M_{\mathrm{MW}}}{M_{\mathrm{sim}}}.\end{array}$$(D.4)

Finally, we corrected the simulated number of binaries already to account for the lifetime of systems like the Gaia BHs:

$$\begin{array}{c}\hfill N_{\mathrm{sim}}=\sum _{\mathrm{sim}}\frac{t_{\mathrm{selected}}}{t_{\max }}.\end{array}$$(D.5)

We selected Gaia BH-like systems as binaries containing a BH with a detached companion. This provides us with the time those binaries spend in this phase, t_selected. Each binary is in maximum evolved for the age of the Universe³, t_max = 13.8 Gyr, hence t_selected < t_max.

Table D.1 summarises counts of systems in the MW for different selections on our simulates samples. The first one includes all simulated binaries resulting in a detached binary with a BH and a low-mass companion, which did not interact during the evolution. The second line tries to mimic the selection done in El-Badry et al. (2023a). The third selection unifies the parameter ranges to just include Gaia BH1 and BH2 at the same time. The fourth line represents the green histograms in Figs 5 and C.1. The last two join the mass estimates of Gaia BH1 and BH2, while being more generous on the orbital parameters.

Appendix E: Milky Way disk coverage

For the volume coverage, we used a simple geometric calculation to convert the distance to a coverage of the MW disk with a thickness of 2.3 kpc and a radius of 15 kpc (Hammer et al. 2007; Yin et al. 2009). It should be noted that Gaia is not in the centre of the disk and we used a solar system distance of about 8 kpc to the disk centre. Thus, we can use a spherical volume observed by Gaia for distances up to half the thickness. At larger distances we had to remove the two spherical sections above and below the disk, which hold up to about 7 kpc. Finally, we know that we would have the full disk covered within about 23 kpc. All this still assumes no extinction, which does not hold towards the galactic centre. To combine the population synthesis results (normalised by mass) with the volume coverage, one would need to assume a constant average density.

All Tables

All Figures

	Fig. 1. Maximum radius of single stars as a function of their initial mass. The colour and symbol show different metallicities.
In the text

Fig. 2. Remnant mass function of single stars depending on their initial mass. The colour and symbol show different metallicities. We used the direct collapse prescription in POSYDON. At Z ≤ 0.2 Z⊙ stars with initial masses above 100 M⊙ might not form a BH because of pair instability and are shown here with a zero mass. Stars with initial masses above the pair-instability gap form BHs outside the plotted range, with masses above 100 M⊙ at the lowest shown metallicity.

In the text

	Fig. 3. Same as Fig. 2 (top) and Fig. 1 (bottom), but here we show different variations of the wind at solar metallicity (see legend). The dotted grey line shows the ZAMS radius, where the discontinuity at about 130 M⊙ is a numerical artefact in the ZAMS models due to the detection criterion of ZAMS. The WR wind variations all overlap in the bottom panel.
In the text

	Fig. 4. Range of BH masses in the plateau regime depending only on the metallicity in Z⊙ (orange, 0.45 and 2 Z⊙), Vink wind factor (green), or WR wind factor (blue), while the others are at default. The default at solar metallicity and without any wind modification is shown in black.
In the text

Fig. 5. Corner plot showing the BH mass, the companion mass, the orbital period, and the eccentricity of binaries in a population synthesis run with an enhanced WR wind (factor = 1.2). While the grey histograms show the full population, the green histograms are limited to binaries with orbital periods up to 5000 days. The observationally determined values for the Gaia BHs are shown in colour (El-Badry et al. 2023a,b; Gaia Collaboration 2024). The histograms are weighted to represent the number of these binaries in the MW, which has large uncertainties as given by the unit $[1_{-1.0}^{+3.3}]$ (for details, see Appendix D).

In the text

Fig. 6. Gaia detection limits for a non-luminous companions of 10 M⊙ on a circular orbit around a F0V, G2V, and K2V star. The lines represent the (sky-uniform) median distance limit at which the Δχ2 = 100 between a single star and a circular Keplerian orbital model fit, corresponding to (typically) well-constrained orbital fits. Estimates were provided for the 5 yr nominal mission (∼DR4) and 10 yr extended mission (∼DR5) as a function of period. Successful detections drop quickly for orbital periods longer than the data time-span, and the graphs are therefore truncated at 5 and 10 yr. At the top, the semi-major axis (a) is specified for a G2V star, although the difference for the other masses would be hardly noticeable. On the right side, we provide a rough estimate of the MW coverage for a specific distance limit, assuming a disk of height 1.15 kpc and radius 15 kpc, and we therefore obtain a full coverage at 23 kpc. All detection values assume zero extinction.

In the text

	Fig. B.1. Same as Fig. 2 but showing the Fryer+12-delayed (top; Fryer et al. 2012) and Patton+Sukhbold20-engine (bottom; Patton & Sukhbold 2020) prescription for the compact object formation.
In the text

	Fig. B.2. Same as the top panel of Fig. 3 but showing the Fryer+12-delayed (top; Fryer et al. 2012) and Patton+Sukhbold20-engine (bottom; Patton & Sukhbold 2020) prescription for the compact object formation.
In the text

	Fig. B.3. Same as Fig. 4 but showing the Fryer+12-delayed (top; Fryer et al. 2012) and Patton+Sukhbold20-engine (bottom; Patton & Sukhbold 2020) prescription for the compact object formation. It should be noted that here the plateau spans over the mass gap between neutron stars and BHs for the strongest wind variations.
In the text

	Fig. B.4. Maximum BH mass depending on the metallicity in Z⊙ (orange), Vink-wind factor (green), or WR-wind factor (blue). The default at solar metallicity and without any wind modification is shown in black. The different symbols indicate the prescription for the BH formation.
In the text

	Fig. B.5. Final profile of an initially 92.42 M⊙ star with WR winds enhanced by a factor 1.2. The dashed black line is the radius, and the coloured lines are the mass fractions of hydrogen (X), helium (Y), and metals (Z).
In the text

Fig. C.1. Corner plot showing the BH mass, the companion mass, the orbital period, and the eccentricity of binaries in a default POSYDON population synthesis run with enhanced WR wind (factor = 1.3); see the text for more details. The grey histograms are the full population, and the green ones are limited to binaries with orbital periods of up to 5 000 days. The observationally determined values for the Gaia BHs are shown in colour (El-Badry et al. 2023a,b; Gaia Collaboration 2024). It should be noted that the histograms are weighted to represent the number of such binaries in the MW, which has large uncertainties as given by the unity unit $[1_{-1.0}^{+3.3}]$; for details see Appendix D.

In the text