Young Star Clusters Dominate the Production of Detached Black Hole-Star Binaries

The recent discovery of two detached black hole-star (BH-star) binaries from Gaia's third data release has sparkled interest in understanding the formation mechanisms of these systems. We investigate the formation of these systems by dynamical processes in young open star clusters (SCs) and via isolated binary (IB) evolution, using a combination of direct $N$-body models and population synthesis simulations. By comparing dynamical and isolated systems created using the same model of binary stellar evolution, we find that dynamical formation in SCs is nearly 40 times as efficient per unit of star formation at producing BH-star binaries compared to IB evolution. We expand this analysis to the full Milky Way (MW) using a FIRE-2 hydrodynamical simulation of a MW-mass galaxy. Even assuming that only $10\%$ of star formation produces SCs with masses $>1000\,\mathrm{M_{\odot}}$, we find that the MW contains $\sim 2 \times 10^5$ BH-star systems, with approximately 4 out of every 5 systems being formed dynamically. Many of these dynamically-formed systems have larger orbital periods, eccentricities, and black hole masses than their isolated counterparts. For binaries older than 100 Myr, we show that any detectable system with $e\gtrsim0.5$ or $M_{\rm BH}\gtrsim 10\,\mathrm{M_{\odot}}$ can only be formed through dynamical processes. Our MW model predicts between 61 and 210 such detections from the complete DR4 Gaia catalog, with the majority of systems being dynamically formed in massive and metal-rich SCs. Finally, we compare our populations to the recently discovered Gaia BH1 and Gaia BH2, and conclude that the dynamical scenario is the most favorable formation pathway for both systems.


INTRODUCTION
Current estimates anticipate that the Milky Way (MW) harbors a population of ∼ 10 7 − 10 9 stellar black holes (BHs) and ∼ 10 4 massive stars which are likely BH progenitors (Garmany et al. 1982;Reed 2003).While the exact binary fraction of stellar BHs is unknown, the high binary fraction of BH-progenitor stars (e.g., Sana et al. 2012) suggests that a substantial quantity of their BH descendants likely exist in binaries as well.Most of the current observational evidence for BH binary systems in the MW comes from X-ray binary systems, consisting of 20 dynamically confirmed BHs in X-ray binaries, with an additional ∼ 50 X-ray sources identified as strong candidates for containing BHs (McClintock et al. 2006;Remillard & McClintock 2006;Corral & Santana 2016).However, population synthesis models indicate udicarlo@sissa.itthat the majority of BHs in binary systems within the MW are likely to be dormant BHs, with large orbital periods that preclude their involvement in X-ray binaries (Portegies Zwart et al. 1997;Corral & Santana 2016;Breivik et al. 2017;Chawla et al. 2022), making them harder to be detected.Over the past few years, it has become possible to identify dormant BHs through the motion of their luminous companions, with the first of these binaries being identified through the radial velocity of their luminous components in the globular cluster (GC) NGC 3201 (Giesers et al. 2018(Giesers et al. , 2019)).More recently, proper motion data from ESO Gaia's Data Release 3 (DR3, Gaia Collaboration et al. 2022) has enabled the discovery of two dormant BHs in BH-star binary systems in the galactic field (Chakrabarti et al. 2022;El-Badry et al. 2023a,b).Combined with followup radial velocity measurements to fully characterize the orbit (El-Badry et al. 2023a,b), these systems offer new insights into the formation, evolution, and characteristics of BHs in our galaxy.The first system, Gaia BH1, consists of a Sun-like main sequence (MS) star orbiting a black hole with mass M BH = 9.62 M ⊙ , eccentricity e = 0.45 and orbital period P = 185.6 days (El-Badry et al. 2023b).Gaia BH2 is a ∼ 1 M ⊙ red giant orbiting a a black hole with mass M BH = 8.9 M ⊙ , eccentricity e = 0.52 and orbital period P = 1277 days (El-Badry et al. 2023a).According to El-Badry et al. (2023b) and El-Badry et al. (2023a), the formation histories of both Gaia BH1 and Gaia BH2 are extremely difficult to explain by isolated binary (IB) evolution, due to the incompatibility between their orbital properties and the outcomes predicted for common envelope evolution.Star clusters (SCs), on the other hand, are highly favored environments for the formation of such systems; not only can dynamical interactions significantly alter the orbital properties of primordial binaries allowing them to follow other evolutionary pathways, but BHstar binaries and their progenitors may also dynamically assemble through binary exchanges and/or multi-body encounters.Despite this, the thin-disk-like Galactic orbits and near-solar metallicities of the Gaia BHs suggest that they most likely did not form dynamically in an old, metal-poor GC.Young/open SCs, on the other hand, are dynamical systems with Galactic orbits and metallicities perfectly compatible with those of Gaia BHs.Additionally, they are the birthplace of the vast majority of massive stars, which are the progenitors of compact objects (e.g.Lada & Lada 2003a;Portegies Zwart et al. 2010).Hence, the majority of BHs in the MW have probably spent the first part of their life in young/open SCs undergoing dynamical interactions.
In this paper, we investigate the dynamical formation of binary systems consisting of a star and a BH in young/open SCs in the MW and via IB evolution.We use a large set of 3 × 10 3 N -body simulations of young/open SCs (from Di Carlo et al. 2020b), along with population synthesis simulations of IB evolution.We describe the population of BH-star binaries in the MW, predicting the number and properties of the intrinsic population, the expected detections in Gaia DR3, and in future Gaia data releases.Finally, we compare our systems with Gaia BH1 and Gaia BH2 to understand their possible formation pathway, and conclude that both systems (and especially BH2) are more compatible with a dynamical formation scenario.Throughout this paper, we use "star" to refer to every luminous stellar type except for white dwarfs and neutron stars (further discriminating between MS stars and giants when appropriate).
In this work, we analyze 3 × 10 3 simulations of young/open SCs, which are extensively discussed in Di Carlo et al. (2020aCarlo et al. ( ,b, 2021)).The initial distribution of stars in space is modeled with fractal initial conditions (Küpper et al. 2011;Ballone et al. 2020Ballone et al. , 2021;;Torniamenti et al. 2021), to mimic the asymmetry and clumpiness of star forming regions (Goodwin & Whitworth 2004).The initial stellar mass M SC of every SC ranges from 10 3 M ⊙ to 3 × 10 4 M ⊙ .Each cluster mass is drawn from a dN/dM SC ∝ M −2 SC distribution, i.e. the initial mass function of YSCs in the MW described in Lada & Lada (2003b).We calculate the initial half-mass radius r h using the relationship from Marks & Kroupa (2012).
The SC initial conditions are generated using mcluster (Küpper et al. 2011).We adopt the initial mass function from Kroupa (2001), with a minimum stellar mass of 0.1 M ⊙ and a maximum stellar mass of 150 M ⊙ .The initial total binary fraction is f bin = 0.4.Mass ratios are drawn from a distribution P(q) ∝ q −0.1 (where q = m 2 /m 1 ∈ [0.1, 1], Sana et al. 2012).All the stars with m ≥ 5 M ⊙ are in binary systems, while the stars with m < 5 M ⊙ are paired stochastically until the initial binary fraction is reached.This is consistent with the multiplicity properties of O/B-type stars (Sana et al. 2012;Moe & Di Stefano 2017).Eccentricities and orbital periods are drawn from the distributions from Sana et al. (2012).
We set the efficiency of common envelope ejection to α = 5.We adopt the rapid core-collapse supernovae model, described in Fryer et al. (2012).Natal kicks are randomly drawn from a Maxwellian velocity distribution.A one-dimensional root mean square velocity σ = 15 km/s is adopted for both core-collapse supernovae and electron-capture supernovae.We have simulated SCs with three different metallicities: Z = 0.02, 0.002 and 0.0002.Each SC is simulated for 100 Myr in a static solar neighbourhood-like tidal field (Wang et al. 2016).From each simulation, we extract all the binaries composed of a BH and a luminous companion which are still in the SC at the end of the simulations, as well as all the binaries of the same type which escape the SC.

Synthetic MW of Binaries and Clusters
Following Chawla et al. (2022), we seed both IBs and SCs following the star-formation history of a cosmological zoom-in simulation of a MW-mass galaxy (specifically the m12i galaxy from the Latte simulations, Wetzel et al. 2016) from the Feedback In Realistic Environments (FIRE-2 ) suite of galaxies (Hopkins 2015;Hopkins et al. 2018).We extract from the publicly-available data every star particle formed in the simulation (which have masses at formation of 7070M ⊙ ), including the particle's metallicity and formation time.These are then binned into low ([Fe/H] < −1.5), intermediate (−1.5 < [Fe/H] < −0.5) and high (−0.5 < [Fe/H]) metallicity bins, designed to match the metallicities of the original nbody6++gpu cluster simulations.Each metallicity bin is then further divided into 100 equal intervals of cosmic time between 0 and 13.78 Gyr, providing us with the metallicity-dependent star formation rate of our MWanalogue galaxy.
To convert this star formation per unit time into stars and clusters, we begin by assuming that all star formation occurs in clusters (Lada & Lada 2003b).Since the lowest-mass clusters from Di Carlo et al. (2021) begin at 1000M ⊙ , we assume that all star formation clumps between 100 − 1000M ⊙ rapidly dissolve into isolated stars and binaries, while clusters above 1000M ⊙ remain bound and undergo dynamical encounters.1 Along with our cluster initial mass function (dN/dM SC ∝ M −2 SC ), this suggests that 10% of star formation should occur in open clusters.
We seed our clusters in the m12i galaxy in the following way: for each bin of star formation, we draw randomly with replacement from the cluster sample, subtracting that cluster's mass from the total mass of star particles formed in that bin, until we have completely turned that 10% of the star formation into clusters.Each time we draw a cluster, we randomly pick an m12i star particle from the same time/metallicity bin, and assign that particle's exact formation time and present-day position in the m12i galaxy to the BH-star binaries produced by that cluster.Note that this procedure generates more clusters that present in the original catalog, meaning that many of our dynamical sources are sampled repeatedly.
The remaining 90% of star formation is turned into massive binaries in a similar fashion: for each time/metallicity bin, the remaining star formation is converted directly into star particles, assuming the primary masses (that will become BHs) follow a 1/M 2.3 initial mass function (IMF) starting at 18M ⊙ (Kroupa 2001) and up to 150 M ⊙ , while the secondaries are drawn from a uniform mass ratio distribution (Mazeh et al. 1992;Goldberg & Mazeh 1994) from 0.1M ⊙ up to the primary masses.The binary orbital periods are drawn randomly from a flat-in-log distribution (Abt 1983), while the eccentricities are drawn from a thermal distribution.Note that these initial conditions (other than the IMF) are different from those assumed in §2.1: they are instead chosen to be consistent with Chawla et al. (2022), with a specific focus on forming BH-star binaries.These systems are then evolved to the present day as described in the next section.For the m12i galaxy, this procedure yields approximately 9.4 × 10 7 binaries with at least one BH progenitor.

Evolving Binaries to the Present Day
For consistency with El-Badry et al. (2023a,b), we evolve both IB and SC binaries to the present day using the binary population synthesis code COSMIC (Breivik et al. 2020).The binaries are evolved using the same stellar and binary evolution parameters and recipes used in Di Carlo et al. (2020b), described in §2.1, to ensure that binary evolution is treated identically in both populations.IBs are evolved starting from the Zero-Age Main Sequence (ZAMS) with absolute metallicities of 0.02, 0.002 and 0.0002 (to match the original cluster simulations), from the time when they were formed (T form ) in the galaxy until the current age of the Universe (13.78 Gyr).
The evolution of SC binaries is restarted from the last recorded state of the binary in the N -body simulationthat is, either the time of binary ejection or the end of the cluster simulation at 100 Myr.Restarting the evolution of cluster binaries yields additional complications, as the binary evolution in COSMIC for systems having undergone mass loss or accretion depend on the current mass of the star, its effective initial mass, and the change in its effective age due to the mass loss/transfer (referred to as M (t), M 0 , and the stellar "epoch" in Hurley et al. 2000, respectively).M 0 and epoch are not part of the default nbody6++gpu output, but since we are primarily interested in BHs and MS stars, this did not present a major difficulty: in COSMIC, M 0 is equal to M t for both BHs and MS stars, while the epoch is only relevant for the evolution calculations of MS stars.We assume epoch = 0 Myr for MS stars when restarting their evolu-tion.As a test, we ran the entire population with epoch values of 0 and 100 Myr (the minimum and maximum values it could have given the SC integration time) and found the difference in the population-level statistics to be negligible (only 1% systems show a relative difference greater than ∼1% in their final masses and orbital periods between epoch = 0 Myr and epoch = 100 Myr).While restarting post-MS stars in COSMIC is substantially more complicated, none of the luminous companion stars evolved beyond the MS before the end of the SC simulations (as opposed to evolving onto the giant branch later in the galactic field).In our analysis, we distinguish between BH-MS binaries and BH-Giant binaries.We follow the criterion from Drout et al. (2012) and classify post-MS stars with T eff ≤ 4800 K as giants.This includes stars with BSE stellar types associated with red giants, core helium burning stars and asymptotic giant branch stars (types 3, 4 and 5).
For our analysis, we only select systems with an age t age ≥ 100 Myr.This choice is justified for three reasons: first, our SCs are evolved for 100 Myr, and therefore we are unable to accurately track systems within SCs which formed in the past 100 Myr.Second, the star-formation history of the FIRE-2 simulations can be larger than that observed in MW-mass galaxies at late times by a factor of a few ( §4.5 of Hafen et al. 2022).This is true in particular for the m12i galaxy we consider, whose z = 0 star-formation rate (6M ⊙ /yr, Wetzel et al. 2016) is ≳ 3 − 4 times the best observational estimate of star formation in the MW (1.65 ± 0.19M ⊙ /yr, Licquia & Newman 2015).Finally, one of the main focuses of this paper is to understand the formation channels of Gaia BH1 and Gaia BH2, which are both older than ∼ 1 Gyr (and likely significantly older, especially in the case of Gaia BH2, El-Badry et al. 2023b, §3.8).We briefly discuss the implications of this cutoff (and what our results would look like including such systems) in §5.3.

Gaia Detectability
Following Chawla et al. (2022), we further determine whether Gaia can astrometrically resolve each orbit based on the motion of the luminous companion based on an optimistic and pessimistic criteria.In the optimistic case, we assume that Gaia can resolve any orbit for which the star's motion is at least as large as Gaia's single pointing precision, σ G , when it is projected on the 2d plane of the sky.In the pessimistic case, we require the projected star's orbit to be at least three times as large as σ G .
We also consider the cut in the parallax signal-tonoise ratio (SNR), ω/σ ω , which was applied in in Gaia DR3.This cut was applied as a function of orbital period as shown in Figure 1.The results are shown for the SC BH-star binaries in red and IB binaries in blue.Dynamically-formed BH-star binaries are preferentially located above the SNR cut while IB BH-star binaries which lie above the cut are mostly due to random-chance placements that are unrealistically close to Earth.
To increase our statistical sample of "detected" systems, we consider different positions for Earth in the m12i simulation and combine the results.We use 100 different equally-spaced starting positions along a ring in the disk with radius 8.5 Kpc.For what follows, the results for our synthetic Gaia population are divided by 100 when quoting the number of detectable binaries and all other relevant quantities.

Intrinisic Population
We first explore the intrinsic population of BH-star binary systems in the galaxy; that is, the full population before considering our Gaia detectability criterion.In all our BH-star systems, the luminous member is either a MS or giant star.The exact numbers of systems per formation channel and stellar type is reported in Table 1.Based on the m12i simulation, we predict that the MW hosts a total of ∼ 2×10 5 BH-star systems from the IB and the young/open SC channels, with ∼ 81% of the systems formed dynamically and ∼ 19% formed in isolation.This means that according to our models, SCs produce nearly 4 times more systems than IBs.However, we must emphasize that we have assumed that only 10% of the star formation occurs in young SCs.If we take this into account, it is obvious that SCs are dramatically more efficient than IBs at producing BHstar systems, with SCs producing ∼ 40 times as many systems overall, per unit mass, as IB evolution.This is driven largely by the production of BH-MS star systems; for evolved systems however, the situation is reversed, with IBs producing 2.4 times more BH-Giant binaries than SCs.
Figure 2 shows the distribution of the orbital parameters of the systems in the intrinsic population.SCs are much more efficient at producing systems with large orbital periods for both BH-MS and BH-Giant binaries.IB systems have eccentricities up to ∼ 0.84, while the most eccentric dynamical system has an eccentricity of ∼ 0.99.IBs produce fewer systems with more massive stars: the maximum value of M * for IB is ∼ 5 M ⊙ , while it is ∼ 7 M ⊙ for SCs (though this does depend on the 100 Myr cutoff we have employed, as described in §5.3).The maximum BH mass is ∼ 50 M ⊙ for IBs and ∼ 67 M ⊙ for the SC channel.The SC channel is thus able to form systems with eccentricities and masses not accessible from IB evolution.

Population Detectable by Gaia
The Gaia mission provides a unique opportunity to detect a subset of the intrinsic population of BH-star binary systems in the Milky Way.Table 1 presents the number of systems in our data-set according to their detectability by Gaia.The table displays two subsets: the first containing the number of systems detectable by Gaia DR3 in our model, while the second includes systems detectable by both Gaia DR3 and future Gaia data releases.It is evident how SCs are significantly more efficient at producing detectable binaries.In particular, we see that all the 7 expected Gaia DR3 detections are expected to come from the dynamical formation channel.If we also consider detections with future Gaia data releases, we predict that Gaia will detect between 8 and 32 binaries from IBs, and between 53 and 178 binaries from the dynamical channel, i.e. a factor of ∼ 6 more dynamical systems.In general, our models produce more detectable BH-MS binaries than BH-Giant binaries.
Figure 3 shows the distributions and correlations of the orbital parameters for the detectable population.It is evident that only SCs are able to produce detectable systems with large eccentricities and with large  BH masses.In particular, the detection of a system with e ≳ 0.5 or with M BH ≳ 10 M ⊙ would be a smoking gun of dynamical formation, as IBs are unable to form detectable systems with these characteristics in our models.The SC channel is more efficient at producing detectable binaries with large orbital periods as well.However, we reiterate that these results are for systems older than 100 Myr; we explore the implications of this cutoff in §5.3.

Dynamical Formation and Star Cluster Properties
The results from our analysis of the intrinsic and detectable populations strongly suggest that the majority of these systems are formed dynamically in SCs.Given that, it is informative to understand how many detectable systems produced by SC evolution come from primordial binaries (i.e.binaries present in the SC since the beginning of the simulation that may have been altered through dynamical processes) or from binaries that were assembled later dynamically.We find that approximately 94% of detections stem from these dynamically-assembled binaries, while the remaining 6% arise from primordial binaries; ∼ 20% of the dynamically assembled binaries initially assemble as a star-star system (that later evolves to a BH-star system), while the remaining ∼ 80% form dynamically from a BH and a star.
We can also ask which SCs produce the majority of the BH-star binaries.In Figure 4, we show the number of BH-star binaries produced by SCs of a given mass.Even though the cluster catalog is sampled following a ∝ M −2 SC cluster mass function (containing many more low-mass clusters), the intrinsic population exhibits a relatively uniform distribution across SC masses, while the detectable population predominantly arises from SCs with M SC ≳ 3 × 10 3 M ⊙ .For insights into the formation of BH-MS binaries in SCs with lower masses (M SC = 300 − 1000 M ⊙ ), Rastello et al. (2023) show a comprehensive analysis that compares results from lowand high-mass SCs.
The overwhelming majority of observable binary systems originate from SCs with high metallicity (Z=0.02), with ∼ 99.84% of them forming in such clusters.This is because of the metallicity-dependent star formation history of the MW.Massive metal-rich SCs are thus very efficient at producing detectable binaries, and we expect the majority of BH-star systems in the MW and current/future Gaia detections to have formed in such environments.In a recent study, Tanikawa et al. (2023) estimated that the MW hosts around 1.6 × 10 4 BH-star systems by simulating SCs with M SC = 1000 M ⊙ .In contrast, our analysis suggests a much larger number of ∼ 1.9 × 10 5 BH-star systems, highlighting the crucial role of more massive SCs to comprehend the formation and evolution of these systems.
We also highlight that ∼ 85% of the BH-star binaries formed in SCs are retained by their host SC at the end of the simulations.Despite the relatively low initial mass of the SCs considered in this study, most of them have not undergone complete disruption by the end of the simulations and would likely survive longer (see e.g.Torniamenti et al. 2022).While it is possible that some of the retained binaries may eventually disrupt due to dynamical encounters, it is also highly likely that new BH-star binaries will form as a result of ongoing dynamical interactions.

FORMING GAIA BH1 AND GAIA BH2
In order to understand how the BH-star binary systems detected by Gaia might have formed, we compare their properties to those of our simulated populations.For both the IB and the SC channel, we find the most similar system to Gaia BH1 and Gaia BH2 in the detectable population, and describe their formation history.The identification of the most similar system entails assessing the fractional differences in key parameters, namely the orbital period, primary mass, and secondary mass, between the simulated binaries and the Gaia BHs.The fractional difference is defined as |x − x GAIA |/x GAIA , where x is the parameter of the object in the simulated population and x GAIA is the corresponding parameter of the Gaia BH.The system that exhibits the smallest maximum fractional difference across these parameters (i.e. the one that deviates the least across all the parameters considered) is regarded as the most similar to the Gaia BH.Due to the relatively limited statistical data available for the SC channel, we make the assumption that eccentricities are subject to randomization through dynamical encounters in SCs, where the eccentricity distribution is expected to follow a thermal profile.Thus, we exclude the eccentricity in the evaluation of the most similar systems.The formation histories of these systems is summarized in Figure 5.

Gaia BH1
The parameters of the most similar IB system to Gaia BH1 fall within approximately 25% of those of Gaia BH1.The system has an orbital period P ∼ 232 days, M BH ≃ 8.4M ⊙ and M * ≃ 1.05M ⊙ .The initial conditions of the system are P ≃ 3486 days, M 1 ≃ 55.1M ⊙ , M 2 ≃ 1.04M ⊙ and a metallicity of Z = 0.02.The system undergoes a common envelope episode that significantly shrinks its orbit, and then turns into a BH-MS binary after 4.8 Myr.The binary then evolves unperturbed to the present day.The system has an age of ∼ 1.5 Gyr.
The closest system from the SC channel has the parameters within ∼ 32% of Gaia BH1.The system has an orbital period P ∼ 130 days, M BH ≃ 6.6M ⊙ and M * ≃ 0.6M ⊙ .The system assembles dynamically inside a SC with metallicity Z = 0.02 and initial mass M SC = 3204 M ⊙ .The binary forms ∼ 5 Myr after the beginning of the simulation.The primary turns into a giant, triggering a common envelope episode that shrinks the orbit.Afterwards, the envelope is ejected, and the primary turns into a BH.∼ 10 Myr after its formation, the system is dynamically ejected from the SC.The age of this system is ∼ 3 Gyr.
Both the IB and SC models struggle to reproduce a system that very closely resembles Gaia BH1.Rastello et al. (2023) suggests that the formation of Gaia BH1 might be explained by less massive SCs, with masses ranging from 300M ⊙ to 1000M ⊙ .These lower mass star clusters could provide a favorable environment for the formation of a system with properties similar to Gaia BH1.Based on our analysis of the age distribution and of our detection rates for BH-MS systems in Gaia DR3, we draw the conclusion that the dynamical formation channel is the most favorable for Gaia BH1, as we better explain in §5.2.

Gaia BH2
The closest IB system to Gaia BH2 is a BH-Giant binary with P ∼ 1044 days, M BH ≃ 8.0M ⊙ and M * ≃ 1.3M ⊙ .The values lie within ∼ 19% of the parameters of Gaia BH2.The initial conditions of the system are P ≃ 3967 days, M 1 ≃ 45.8M ⊙ , M 2 ≃ 1.36M ⊙ .The system has a metallicity of Z = 0.02, and undergoes a common envelope episode throughout its evolution.The system has an age of ∼ 4.1 Gyr, which is not compatible with the age of Gaia BH2.
The best candidate from the SC channel has the parameters within 13% of the Gaia BH2 parameters.The system has an orbital period P ∼ 1197 days, M BH ≃ 8.1M ⊙ and M * ≃ 1.2M ⊙ .The binary is dynamically assembled in a SC with metallicity Z = 0.02, initial mass M SC = 14284 M ⊙ and an age of ∼ 6 Gyr.The system forms as a BH-MS binary ∼ 39 Myr after the beginning of the simulation, and it later evolves to a BH-Giant.The system is retained by its host SC at the end of the simulation.
Both channels produce systems that exhibit close parameter matches to Gaia BH2.Among these, the dynamical scenario yields the closest resemblance.If we also take into account that Gaia BH2 is estimated to be older than 5 Gyr, and that our Gaia DR3 detections for BH-Giant systems come entirely from SCs, we conclude that Gaia BH2 has likely formed via the dynamical formation channel.We provide a more extensive analysis on this in §5.2.

Dynamical versus Isolated
We see from Table 1 that the dynamical channel produces ∼ 4.3 times more BH-star systems than the isolated channel.Since we consider that only 10% of the total star formation occurs into SCs, we find that the number of systems produced per unit stellar mass by SCs is ∼ 40 times larger than the number of systems produced by isolated binary evolution.This discrep- The time axis and the size of the objects and orbits are not to scale.Even though Gaia BH-like systems can be formed by both the formation channels, our results on rates and ages presented in §5.2 strongly suggest that the dynamical scenario is the primary formation pathway for the Gaia BHs.
ancy can be attributed to several mechanisms.Stellar dynamics in SCs significantly affects the formation and evolution of binaries.Multi-body encounters can dynamically assemble BH-star or star-star binaries which can later evolve into BH-star binaries, resulting in combinations of orbital parameters that are inaccessible in IB evolution.Furthermore, three-body encounters between a binary and a single object can modify the eccentricity and orbital period of binaries, triggering or preventing mass transfer episodes (e.g., Roche-lobe overflow and common envelope evolution) that may or may not occur if the binary evolved in isolation.Common envelope episodes significantly reduce the orbital period of binaries.From Figure 2 and Figure 3, we see how the dynamical channel is extremely more efficient at producing binaries with larger orbital periods.Indeed, we find that 99.4% of all the simulated IB systems which become BH-star binaries undergo at least one common envelope episode throughout their evolution.This can be easily avoided in SCs; as we said in §3.3, the majority of detectable SC binaries dynamically assembles as a BH-star binary already, which means that they do not undergo a common envelope episode2 .The dynamical channel is thus not only more efficient, but can also leave unique fingerprints on the properties of the binaries.

Gaia BHs formation channels
Our models provide several hints on how the dynamical scenario is a more favorable formation pathway for Gaia BH1 and Gaia BH2.According to our results in Table 1, all the expected Gaia DR3 detections come from SCs, suggesting that Gaia BHs have more likely formed dynamically.In order to form systems like Gaia BH1 or Gaia BH2 through IB evolution, El-Badry et al. (2023b) and El-Badry et al. (2023a) highlight the necessity of assuming an exceedingly large and potentially unrealistic value for the common envelope efficiency (α ≈ 12.8 for Gaia BH1).Lower values of α would have led to the merger of the binary components, preventing the formation of Gaia BHs through this channel.As discussed in §5.1, dynamics offers an alternative pathway for the formation of BH-star systems that completely bypasses the common envelope phase.
The age distribution of the detected systems is shown in Figure 6.The efficiency of both channels is comparable at 0.1 Gyr, but the dynamical channel is much more efficient at forming detectable systems with older ages.Gaia BH1 has an age ≳1 Gyr, while Gaia BH2 has an age ≳5 Gyr, indicating that both Gaia BHs, particularly Gaia BH2, have likely formed dynamically, according to our models.Rastello et al. (2023) show that both lower mass (300-1000M ⊙ ) and higher mass (1000-30000 M ⊙ ) young SCs are efficient in forming Gaia BH1-like systems.Gaia BHs may also have formed via other channels not taken into account in this study, such as formation in isolated hierarchical triples or hierarchical triples formed in SCs (see e.g.Trani et al. 2022).Another plausible channel could be the dynamical formation in a globular cluster (GC); while the Galactic orbit of both Gaia BHs is not aligned with any GC in the MW (El-Badry et al. 2023b,a), these systems might have dynamically formed in GCs which have already completely disrupted.

BH-star ages and the 100 Myr age cutoff
Throughout this paper, we have argued that SCs dominate the production of old BH-star systems in the MW.We show this explicitly in Figure 6, which displays the distribution of ages of the systems, and reveals that the dynamical formation channel appears to be more efficient at producing older systems, particularly those older than 1 Gyr.Of course, it is also obvious in Figure 6 that our 100 Myr cutoff has removed a large number of young systems, particularly from the IB channel, that contribute significantly to the detectable population.For comparison, we present the results obtained without applying the 100 Myr age cutoff (which was justified in §2.3), starting with Table 2 which reports the number of systems in the intrinsic and detectable populations.Although there is an increase in the number of systems in the isolated intrinsic population, the dynamical channel remains dominant, constituting approximately 80% of the total intrinsic population.
However, we observe a significant increase in the number of isolated detectable systems.The expected number of IB detections by Gaia DR3 increases from 0 to 87, while the total number of IB Gaia detections increases from 8-32 to 286-455.Conversely, the number of detectable systems from the dynamical channel barely increases.This disparity can be attributed to the fact that systems in the isolated population are predominantly very young (see Figure 6) and preferentially host bright massive stars that are significantly easier to detect.This trend is evident from Figure 7 and Figure 8, which illustrate the significant impact of the cutoff on the M * distribution.The difference comes almost entirely from a huge population of young massive O-stars.The IB detectable systems younger than 20 Myr have an average BH mass of 8 M ⊙ and an average star mass of 41 M ⊙ , and constitute ∼ 86% of the total number of detectable systems in the whole IB catalogue without the 100 Myr age cutoff.We emphasize that our predictions for the formation channels of both Gaia BH1 and Gaia BH2 should remain unaffected by this cutoff, as they both are older than 1 Gyr.

CONCLUSIONS
The MW is very likely populated by a large number of stellar BHs, many of which probably reside in binary systems.The recent discovery of two BH-star binaries in Gaia DR3 has sparkled new interest in understanding the formation channels and population characteristics of such systems.In this paper, we have investigated the formation of BH-star binaries in young SCs and via IB evolution.
According to our simulations, the MW harbors a total of ∼ 2×10 5 BH-star systems of which ∼ 81% formed dynamically and ∼ 19% formed in isolation.We find that dynamical formation in SCs is nearly 40 times more efficient per unit stellar mass at producing BH-star binaries compared to isolated binary evolution.Dynamical systems tend to have larger orbital periods and eccentricities than the isolated ones.We expect that a total of 7 BH-star systems are present in the Gaia DR3 data, all of which come from the dynamical channel.We also predict between 61 and 210 detections from the whole Gaia mission, ∼ 83% of which come from SCs.Overall, dynamics enhances dramatically the number of BH-star binaries, both in the intrinsic and the detectable populations.
We compare our detectable populations with Gaia BHs, and we conclude that our models support the dynamical scenario as the primary formation pathway for Gaia BH1 and Gaia BH2.Our results suggest that identifying systems with an an eccentricity greater than ∼ 0.5, or with black hole mass M BH exceeding 10 M ⊙ would provide compelling evidence for their dynamical formation.With new detections from future Gaia data releases expected over the next few years, these results will help disentangle the formation channels of BH-star systems and the characteristics of the BH population in the MW.
Finally, our analysis also revealed the presence of some BH-Giant systems with 3 M ⊙ < M BH < 5 M ⊙ , i.e. below the minimum BH mass allowed by the rapid supernova mechanism (Fryer et al. 2012).These peculiar systems form through accretion of matter by neutron stars from their binary companions, leading to the formation of BHs through accretion-induced collapse (AIC).None of these systems are detectable by Gaia according to our analysis, but we have not yet compared these systems  (identified only in the IB channel) to their dynamical counterparts.Efforts to better characterize these systems are currently underway.

Figure 1 .
Figure 1.The orbital period vs parallax SNR, ω/σω are shown for BH-star binaries formed in clusters (red) and in isolated binaries (blue) from all of our 100 MW-like simulation realizations.Dots show BH-MS binaries and crosses show BH-Giant binaries.The dashed black like shows the parallax SNR cut placed as a function of orbital period for Gaia DR3 such that only binaries above the line are included in the data set.

Figure 2 .
Figure 2. Orbital parameters of systems in the intrinsic population.From top to bottom we show the distributions of binary orbital periods, eccentricities, stellar masses M * , and BH masses M BH , with SC binaries shown in red lines and IB shown in blue lines.Solid lines represent BH-MS binaries, while dashed lines show BH-Giant binaries.The yellow and purple vertical lines represent the values of Gaia BH1 and Gaia BH2 respectively.

Figure 3 .
Figure 3. Correlations and distributions of orbital period, eccentricity, mass of the star and mass of the BH of detectable systems.Dynamical systems are shown in red, while isolated systems are shown in blue.Gaia BH1 is shown in yellow, while Gaia BH2 is shown in violet.BH-MS systems are represented by filled circles in the scatter panels and by solid lines in the histogram panels.BH-Giant systems are represented by crosses in the scatter panels and by dashed lines in the histogram panels.Scatter plots show systems from all the 100 observers around the Galaxy (see §2.4), while histograms show the number of detections N det = N/100.

Figure 4 .
Figure 4. Distributions of the initial SC masses M SC which produce BH-star binaries.The intrinsic population is shown by the solid red line, while the detectable population is shown by the dotted red line.The initial masses of the simulated SCs are drawn from a dN/dM SC ∝ M −2SC distribution(Lada & Lada 2003b), meaning that we simulated more low mass SCs and less high mass SCs.

Figure 5 .
Figure 5. Formation history of the most similar systems to Gaia BH1 (left panels) and Gaia BH2 (right panels) in our simulations.Top panels show systems from the isolated binary evolution channel; bottom panels show systems from the star clusters channel.Main sequence stars (with label MS) are represented as yellow stars; core helium burning stars (label cHeB) are visualized as red stars; naked helium stars (label Naked He) are represented as blue stars; black holes (label BH) are shown as black circles.The mass of each object is shown next to them.The orbital period in days of the final binaries is shown underneath.The time axis and the size of the objects and orbits are not to scale.Even though Gaia BH-like systems can be formed by both the formation channels, our results on rates and ages presented in §5.2 strongly suggest that the dynamical scenario is the primary formation pathway for the Gaia BHs.

Figure 6 .
Figure 6.Distributions of the ages of systems in our datasets.The intrinsic population is shown by solid lines, while Gaia detectable systems are shown by dotted lines.Dynamical systems are shown in red, while isolated systems are shown in blue.The two dashed vertical lines show the lower limits of the age of Gaia BH1 and Gaia BH2.The grey shaded area on the left represents the systems not taken into account due to the age cutoff described in §2.3.

Figure 7 .
Figure 7. Distributions of orbital parameters of systems in the intrinsic population, including systems with age < 100 Myr.From top to bottom: orbital period, orbital eccentricity e, star mass M * , BH mass M BH .SC binaries are shown by red lines, IB are shown by blue lines.Solid lines represent BH-MS binaries, while dashed lines show BH-Giant binaries.The yellow and purple vertical lines represent the values of Gaia BH1 and Gaia BH2 respectively.
Figure 8. Same as Figure 7, but for systems detectable in future Gaia data releases.