On the formation of a $33\,M_{\odot}$ black hole in a low-metallicity binary

A $33\,M_\odot$ black hole (BH) was recently discovered in an 11.6-year binary only 590 pc from the Sun. The system, Gaia BH3, contains a $0.8\,M_\odot$ low-metallicity giant ($\rm [M/H]=-2.2$) that is a member of the ED-2 stellar stream. This paper investigates whether the system could have formed via isolated binary evolution. I construct evolutionary models for metal-poor massive stars with initial masses ranging from $35-55\,M_{\odot}$, which reach maximum radii of $1150-1800\,R_{\odot}$ as red supergiants. I then explore what combinations of initial orbit, mass loss, and natal kick can produce the period and eccentricity of Gaia BH3. Initial orbits wide enough to accommodate the BH progenitor as a red supergiant can match the observed period and eccentricity, but only if the BH formed with a significant natal kick ($v_{\rm kick}\gtrsim 20\, {\rm km\,s^{-1}}$). These models are disfavored because such a kick would have ejected the binary from the ED-2 progenitor cluster. I conclude that Gaia BH3 likely formed through dynamical interactions, unless the BH progenitor did not expand to become a red supergiant. Only about 1 in 10,000 stars in the solar neighborhood have metallicities as low as Gaia BH3. This suggests that BH companions are dramatically over-represented at low-metallicity, though caveats related to small number statistics apply. The fact that the luminous star in Gaia BH3 has been a giant -- greatly boosting its detectability -- only for $\sim$1% of the time since the system's formation implies that additional massive BHs remain to be discovered with only moderately fainter companions. Both isolated and dynamically-formed BH binaries with orbits similar to Gaia BH3 are likely to be discovered in Gaia DR4.

1. INTRODUCTION Gaia Collaboration et al. (2024) recently reported discovery of an astrometric binary containing a 0.8 M ⊙ star on the lower giant branch and a 33 M ⊙ dark companion that is presumably a black hole (BH).The system, which they named Gaia BH3, was discovered with pre-release DR4 astrometry, and its orbital motion was confirmed with (nearly) independent radial velocity measurements from the Gaia RVS spectrometer (Cropper et al. 2018).
The discovery is exciting for several reasons.First, Gaia BH3 contains the most massive robustly measured stellar-mass BH in the local Universe.Second, the star has the lowest metallicity of any known star orbiting a BH, [M/H] = −2.2.Third, the system provides the first empirical evidence that low-metallicity massive stars leave behind massive BHs, making them prime candidates for being progenitors of gravitational wave sources.And fourth, the system is bright and nearby (G = 11.2, d = 590 pc), suggesting that many similar systems that are fainter and/or more distant remain to be discovered.
Gaia BH3 is the widest BH binary discovered to date, with an orbital period of ∼ 4250 days and a semimajor axis of ∼ 16.5 au.The eccentricity is relatively high, e ≈ 0.73, such that the star and BH pass within 4.5 au of each other at periastron.How the system formed is uncertain.Gaia Collaboration et al. (2024) note that a typical red supergiant would not fit within the current orbit at periastron.They thus consider formation of the system from an isolated binary unlikely, and in-Corresponding author: kelbadry@caltech.eduvoke dynamical exchange as a possible alternative.Balbinot et al. (2024) report that Gaia BH3 is a member of the ED-2 stellar stream (Dodd et al. 2023), lending credence to the dynamical formation hypothesis.Because ED-2 stars do not display light element abundance variations characteristic of globular clusters, Balbinot et al. (2024) propose that ED-2 formed from a low-mass cluster (M cluster ≲ 4 × 10 4 M ⊙ ).
Here, I consider the possibility of formation from an isolated binary in more detail.While the association of Gaia BH3 with ED-2 appears robust, the progenitor cluster's initial density and mass are uncertain, leaving open the possibility that Gaia BH3 is a primordial binary as opposed to one formed by dynamical exchange.I explore (a) the expected mass and maximum radius of the BH progenitor, and (b) the range of possible orbits that could have produced the current orbit in the presence of possible impulsive mass loss and kicks to the BH during its formation.I ultimately conclude that, while orbits similar to Gaia BH3's can naturally be produced from isolated binaries if BHs are born with modest natal kicks, this channel is disfavored for Gaia BH3 in particular because such natal kicks would have ejected the binary from the ED-2 progenitor cluster.
The remainder of this paper is organized as follows.In Section 2, I describe evolutionary models for lowmetallicity massive stars representative of the BH progenitor.Section 3 uses Monte Carlo simulations of BH formation to assess the likely initial separation if Gaia BH3 formed from an isolated binary.Section 4 argues that BH companions appear to be strongly overrepre-sented at low metallicity, and that additional nearby BHs orbited by low-metallicity stars likely remain to be discovered.Future prospects and broader implications are briefly discussed in Section 5.

THE BH PROGENITOR
To estimate the maximum radius and expected mass loss history of the BH progenitor, I calculated a small suite of evolutionary models using MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019;;Jermyn et al. 2023).The calculations closely follow the setup described by Klencki et al. (2020), with MESA updated to version r23.05.1.The most important ingredients are summarized below, and I refer to Klencki et al. (2020Klencki et al. ( , 2021) ) for more details.
I evolve single stars with initial masses ranging from 30 to 55 M ⊙ and metallicities ranging from Convection is modeled following Henyey et al. (1965) (mlt option = 'Henyey') with a default mixing length parameter α = 1.5.The Ledoux criterion is used and semiconvection is modeled following Langer et al. (1983), with efficiency parameter α SC = 100 motivated by Schootemeijer et al. (2019).MLT++ is not used.
Step overshooting is assumed above the H and He burning cores, with overshoot f = 0.33 and overshoot f0 = 0.05 by default.I experiment with varying the mixing length and overshooting parameter, as described below.
Stellar winds are modeled following Brott et al. (2011): for hot, hydrogen-rich stars (T eff > 25 kK; X S > 0.7), the models take the wind mass loss rates from Vink et al. (2000Vink et al. ( , 2001, here X S is the surface hydrogen mass fraction).For hydrogen-poor stars (X S < 0.4), they use the Wolf-Rayet wind model from Hamann et al. (1995) and reduce the mass loss rate by a factor of 10 to account for wind clumping (Yoon et al. 2006).For stars with intermediate surface hydrogen abundances (0.4 < X S < 0.7), the mass loss rate is interpolated between the two prescriptions above.For stars cooler than T eff = 25 kK, the mass loss rate is set to the larger of the rates predicted by Nieuwenhuijzen & de Jager (1990) and Vink et al. (2001).These mass loss rates all include a Ṁ ∝ Z 0.85 scaling with metallicity.As a result, the predicted wind mass-loss rates are all rather low at the metallicities of interest for Gaia BH3.For consistency with Klencki et al. (2020), the models assume Z ⊙ = 0.017 when calculating mass loss rates, but I label the metallicities by [M/H] = log(Z/0.014)following Asplund et al. (2009).
For rotating models, I include the effects of Eddington-Sweet circulation (D ES factor = 1.0), secular shear instabilities (D SSI factor = 1.0), and the Goldreich-Schubert-Fricke instability (D GSF factor = 1.0), with an efficiency factor am D mix factor = 1/30 following Heger et al. (2000).I initialize models using the MESA create pre main sequence model routine and run until the end of carbon burning, by which time the model is within days of core collapse and there is no time for further radius evolution.I limit the timestep with varcontrol target = 0.0001 and delta HR limit = 0.0005 and the mesh resolution with max dq = 0.001.Sensitivity to these choices was explored by Klencki et al. (2021).
To verify that the problem setup is consistent with the one used by Klencki et al. (2020), I began by reproduc-ing their 35 M ⊙ model at metallicity Z = 0.00017 and comparing to their published tracks.The model reaches a maximum radius of 1331 R ⊙ shortly before core collapse.Their corresponding track has a maximum radius 1340 R ⊙ at the same point, and its evolution in the HR diagram is nearly indistinguishable from that of the model calculated here.
Figure 1 shows the evolution of several models in the HR diagram.Core hydrogen exhaustion, the onset of carbon burning, and the end of carbon burning are marked in each panel with a diamond, circle, and star symbol, respectively.As a fiducial model -shown in black in each panel -I take the nonrotating model with initial mass of 35 M ⊙ and metallicity [M/H] = −2.2.The mass of this model at the end of carbon burning is 34.7 M ⊙ , reflecting the fact that winds are predicted to be very weak at these metallicities.This model could be expected to produce a ∼ 33 M ⊙ BH if the SN shock fails and it undergoes core collapse with little mass loss, since loss of up to a few M ⊙ is still expected to occur as a result of neutrino mass loss and the subsequent shock (Fernández et al. 2018).The fiducial model reaches a maximum radius of 1150 R ⊙ shortly before core collapse.
The four panels in Figure 1 show the effects of varying metallicity, mass, rotation, and convective mixing length and overshooting.The maximum radius reached by our models varies relatively weakly with metallicity over −2.2 < [M/H] < −1.6.The model with [M/H] = −2.5 (slightly lower than Gaia BH3) is significantly smaller than other models, reaching a maximum radius of ∼ 550 R ⊙ .The maximum radius increases with initial mass, ranging from ≈ 850 R ⊙ for the 30 M ⊙ model to ≈ 1800 R ⊙ for the 50 M ⊙ model at [M/H] = −2.2.The maximum radius is sensitive to overshooting, as explored by Schootemeijer et al. (2019): smaller overshooting parameters lead to more compact supergiants.For the choice of parameters explored here, sensitivity to the mixing length is relatively weak: increasing α MLT from 1.5 to 3.0 reduces the maximum radius by ∼ 10%.
The shaded regions in Figure 1 show R > 800 R ⊙ .In Section 3, I argue that this is the largest progenitor that could plausibly have fit inside the orbit of Gaia BH3 before the BH's formation.Most of the models exceed this radius shortly before core collapse.A notable exception is provided by the two rapidly-rotating models shown in the lower left panel of Figure 1, which go through chemically homogeneous evolution (CHE).At high rotation rates, large-scale meridional circulations can efficiently mix helium into stellar envelopes, preventing the buildup of a chemical gradient and core/envelope dichotomy that leads to expansion in nonrotating models (Maeder 1987;Yoon & Langer 2005;Woosley & Heger 2006).Under CHE, stars evolve from the main sequence toward the helium main sequence, where they have effective temperatures above 10 5 K and radii smaller than 1 R ⊙ .These models eventually ignite carbon burning and terminate their evolution without ever reaching R > 10 R ⊙ .At the low metallicities considered here, CHE is predicted to occur with initial rotation rates of Ω/Ω crit ≥ 0.45.
One of the main differences between the models shown here and higher-metallicity massive star models is that low-metallicity models remain compact during core helium burning: although they are red supergiants just before core collapse, they spend the large majority of co m m on en ve lo pe lik el y -MESA evolutionary models for low-metallicity massive stars.The model shown in black is the same in all panels: a 35 M ⊙ star with metallicity [M/H] = −2.2 and no rotation.Colored lines show models varying metallicity (upper left), initial mass (upper right), initial rotation rate (lower left), and overshooting or mixing length (lower right).Diamonds, circles, and stars mark core hydrogen exhaustion, the onset of carbon burning, and core collapse.Dashed diagonal lines show constant radius.Shaded gray region marks R > 800 R ⊙ , beyond which the star would have been unlikely to fit inside the pre-SN orbit of Gaia BH3 (Section 3).Most of the models exceed this limit at their maximum radius and thus would have overflowed their Roche lobes.The only exceptions are rapidly-rotating models that undergo chemically homogeneous evolution (CHE; lower left) and models with reduced overshooting (lower right).
their post-main sequence evolution as blue or yellow supergiants with R ≲ 100 R ⊙ , only expanding to red supergiant dimensions during their final ∼ 1000 years (see Klencki et al. 2020, for further discussion).This is illustrated in Figure 2, which shows the late-stage radius evolution for models with a range of masses (top panel) and metallicities (bottom panel).The fiducial model increases its radius by a factor of 2 in its final 800 years.This behavior is not found in higher metallicity models, even at "low" metallicites comparable to those found in the SMC and LMC.The upshot is that if Gaia BH3 formed as a primordial binary, there may not have been sufficient time for tides to circularize the orbit when the BH progenitor was a red supergiant.
-Radius evolution of MESA models in the final years before core collapse.Mass is varied in the top panel, and metallicity in the bottom panel.At low metallicity, the models expand significantly in the final ∼ 1000 yr of their evolution, only appearing as red supergiants shortly before core collapse.Expansion occurs earlier at higher metallicity.This means that wide low-metallicity binaries are less likely to be tidally circularized before the formation of a BH.
BH3: P orb,final = 4000 − 4500 d, and e final = 0.68 − 0.78.Narrowing these ranges reduces the number of surviving Monte Carlo samples but has no significant effect on their distribution.
3.1.Circular orbits I first consider pre-SN orbits that have zero eccentricity, as would be expected if tides circularized the orbit when the BH progenitor expanded.In this case, the system's current eccentricity must be a result of natal kicks and/or instantaneous mass loss during the SN.I use the formalism from Brandt & Podsiadlowski (1995) to predict post-SN parameters.
I begin with a uniform distribution of pre-SN orbital periods between 0 and 8500 d (twice the currently observed orbital period of Gaia BH3).I simulate kicks assuming orientations distributed uniformly on the sky and velocities distributed uniformly between 0 and 100 km s −1 .I take a uniform distribution of ∆m, the mass loss of the BH progenitor during the SN, between 0 and 67 M ⊙ , implying a progenitor mass of 33-100 M ⊙ .This is not a population synthesis simulation -pre-SN orbits are unlikely to be uniformly distributed in all parameters -but should be viewed as an exploration of what combinations of orbits, kicks, and mass loss could have produced the observed orbit.
The results of this experiment are shown in Figure 3. Black histograms show all simulations that produce post-SN periods and eccentricities similar to Gaia BH3 today, while red histograms show the subset of those that predict a post-SN systemic velocity v sys < 10 km s −1 .This is a conservative estimate of the largest post-SN systemic velocity for which Gaia BH3 could plausibly have remained bound to the low-mass cluster that formed the ED-2 stream.
When there is no restriction on v sys (black histograms), initial separations ranging from ∼ 800 to ∼ 6000 R ⊙ can produce the observed orbit.Initial separation closer than ∼ 800 R ⊙ cannot, because widening these to the separation observed today via a kick or mass loss would result in a higher final eccentricity than is observed.Wider initial separations could match the orbit for suitably chosen kicks, but these are excluded by the prior of P orb, init < 8500 d.
The upper right panel of Figure 3 shows the pre-SN separations and post-SN systemic velocities of all simulations that match the observed period and eccentricity.The combination of low v sys (≲ 20 km s −1 ) and large initial separation (≳ 1200 R ⊙ ) is ruled out because such models produce too low eccentricities.The red histograms show that among simulations with v sys < 10 km s −1 , only relatively tight pre-SN orbits with a init < 1200 R ⊙ can match the observed period and eccentricity.These calculations all have weak kicks (v kick ≲ 10 km s −1 ) and substantial mass loss, such that the final eccentricity is mainly a result of the Blaauw (1961) kick.Initial masses M 2,init ≳ 70 M ⊙ are ruled out because they become unbound due to mass loss.
The largest star that could fit inside a pre-SN orbit of semimajor axis a init has radius R max = a init f q , where f q (M 1 /M 2 ) is given by Eggleton (1983) and varies from 0.68 to 0.71 over M 1 = 35 − 70 M ⊙ , assuming M 2 = 0.8 M ⊙ .This means that for the pre-SN orbits that can produce the observed period and eccentricity of Gaia BH3 without escaping the ED-2 progenitor, the ±1σ range of maximum radii is R max ≈ 680 +73 −87 R ⊙ .As shown in Section 2, most plausible models for the progenitor of the BH reach maximum radii larger than this value, so it seems unlikely that the present orbit of Gaia BH3 could have resulted from isolated binary evolution if the pre-SN orbit was circular.
3.2.Eccentric orbits I next consider pre-SN orbits that are eccentric.These may be expected if the BH progenitor only expands late in its evolution, leaving insufficient time for tidal circularization (e.g. Figure 2).I again simulate N = 10 7 orbits, now with a uniform distribution of eccentricities between 0 and 1.I assume that the SN occurs at a random phase, corresponding to a uniform distribution of the pre-SN mean anomaly, and use the formalism from Hurley et al. (2002) to predict the post-SN orbital eccentricity, period, and center-of-mass velocity.
The results are shown in Figure 4.In this case, I show the distributions of pre-SN periastron distance, d peri, init , rather than semimajor axis, since this is the key parameter for determining whether the BH progenitor would have overflowed its Roche lobe.As in the circular case, a broad range of initial orbits are possible for arbitrary kicks, but a relatively narrow range of d peri, init can produce the observed period and eccentricity subject to the restriction of v sys < 10 km s −1 .Most of these orbits have weak kicks and little mass loss, and thus are similar to the orbit of Gaia BH3 today, which has d peri ≈ 800 R ⊙ .-Constraints on kicks and the pre-SN orbit of Gaia BH3-like binaries.I assume that the binary had a circular orbit before the SN and that the BH received a kick with velocity v kick when it formed.I simulate 10 7 pre-SN orbits with a uniform period distribution P orb ∼ U (0, 2P orb,BH3 ), a pre-SN mass distribution M 2, init /M ⊙ ∼ U (33, 100), and a kick velocity distribution v kick /(km s −1 ) ∼ U (0, 100).Black histograms show the properies of all the systems whose post-SN orbits have periods and eccentricities similar to Gaia BH3 (upper left).Red histograms show the subset of these that result in a post-SN systemic velocity vsys < 10 km s −1 , as required to remain bound to the ED-2 stream progenitor.Orbits with a broad range of pre-SN separations can produce the observed period and eccentricity.However, most of these require significant kicks and result in vsys > 10 km s −1 .Simulations with vsys < 10 km s −1 can only produce the observed period and eccentricity for a init ≲ 1200 R ⊙ .The BH progenitor would overflow its Roche lobe in such orbits for R ≳ 800 R ⊙ .
Orbits with significantly larger d peri, init and weak kicks end up with longer periods and/or lower eccentricities than the values observed for Gaia BH3 today.
Among orbits that reproduce the observed eccentricity and period of Gaia BH3 with v sys < 10 km s −1 , 90% have d peri, init < 1120 R ⊙ , such that they would overflow their Roche lobes at periastron for radii R > 775 R ⊙ .99% would overflow their Roche lobes at periastron for R > 1125 R ⊙ .Given that the MESA models shown in Figure 1 predict maximum radii of ≳ 1000 R ⊙ , I conclude that the BH progenitor is unlikely to have avoided Roche lobe overflow in any orbit that could have produced the  3, but for eccentric pre-SN orbits.I simulate 10 7 orbits with initial eccentricity e init ∼ U (0, 1) and the SN occurring at a random phase.Other parameters are distributed as in Figure 3. d peri, init is the separation between the BH progenitor and companion star at periastron in the pre-SN orbit.Most orbits with vsys < 10 km s −1 (red) have d peri, init ≲ 1100 R ⊙ and thus would overflow their Roche lobes if the BH progenitor reached R ≳ 800 R ⊙ .
current orbit of Gaia BH3, unless it formed through chemically homogeneous evolution or had a maximum radius significantly smaller than predicted by the evolutionary models shown in Figure 1.Gaia BH3 has a much lower metallicity than typical stars in the solar neighborhood.It is expected that lowmetallicity massive stars may form higher-mass BHs due to their weaker winds (Woosley et al. 2002), and Gaia BH3 was published in advance of DR4 primarily because of the BH's high mass.However, I argue here thateven ignoring the BH's mass -the fact that a BH of any mass was found with such a low-metallicity companion suggests that wide BH companions are overrepresented -Metallicity distribution of red giants within 1 kpc of the Sun, calculated from the Gaia XP metallicity catalog of Andrae et al. (2023, black) and APOGEE DR17 (Abdurro'uf et al. 2022, red).The metallicity of Gaia BH3 is marked with a shaded region: [M/H] = −1.82according to Andrae et al. (2023), and [M/H] = −2.2according to Gaia Collaboration et al. (2024).The fraction of giants with metallicities as low as Gaia BH3 is only ∼ 0.0001.at low metallicity.
Figure 5 shows differential and cumulative metallicity distributions of giants within 1 kpc of the Sun.In black, I show giants with high-quality metallicities from Gaia XP spectra, which I take from Table 2 of Andrae et al. (2023).In red, I show giants with high-SNR spectra from DR 17 of the APOGEE survey (Majewski et al. 2017;Abdurro'uf et al. 2022).Here I exclude stars targeted non-randomly, since several of the APOGEE targeting cartons intentionally selected low-metallicity stars.The shaded gray region shows the plausible metallicity uncertainty of Gaia BH3, which goes from [M/H] = −2.2 as estimated by Gaia Collaboration et al. (2024) to [M/H] = −1.82 as estimated for the source by Andrae et al. (2023).I assume [M/H] = −2.0 in the discussion below.
Both the XP and APOGEE samples suggest that only 1 in ∼ 10 4 giants within 1 kpc have [M/H] < −2.0.This very small fraction strongly implies that wide BH companions are overrepresented at low metallicity.If they were not, we could expect to find ∼ 10 4 BH companions -just to red giants -within 1 kpc.The true number of BH + giant binaries within 1 kpc is unknown, but 0 have been discovered so far: the nearest such system is Gaia BH2 at d = 1.16 kpc (El-Badry et al. 2023b).It seems quite unlikely that 10 4 remain to be discovered (although published Gaia data are not yet sensitive to orbital periods as long as that of Gaia BH3), and so a more plausible explanation is that wide BH companions are more common at low metallicity.Early results from searches for neutron star companions in the Gaia data similarly imply that they are overrepresented at low metallicity (El-Badry et al. 2024a).More work is needed to understand why.

More low-metallicity BHs are expected in DR4
There is also reason to believe that more BHs orbited by metal-poor stars exist in the solar neighborhood.Such stars are uniformly ≳ 12 Gyr old.This means that giants sample a narrow mass range of roughly (0.78 − 0.80) M ⊙ .For a Kroupa IMF, there are ∼165 stars with M = (0.1 − 0.78) M ⊙ for every star with M = (0.78 − 0.80) M ⊙ , suggesting that even within 1 kpc, other massive BHs are likely orbited by lower-mass metal-poor stars that are still on the main sequence.The mass distribution of companions to BHs of course may be different from the IMF, but it seems unlikely that 0.8 M ⊙ companions would be significantly more common than companions with mass of 0.7 M ⊙ or 0.6 M ⊙ .At the distance and extinction of Gaia BH3, such companions would have G < 14 for M > 0.75 M ⊙ , or G < 16 for M > 0.61 M ⊙ .These magnitudes are too faint for Gaia to have measured multi-epoch RVs from RVS spectra, but bright enough that the astrometric orbit would be resolved with SNR ≳ 100.With careful processing of the astrometric data, a source like Gaia BH3 could likely be detected even at G = 18−19 -where individual astrometric measurements have uncertainties of 1-3 mas in the along-scan direction, a factor of 10-30 smaller than Gaia BH3's photocenter ellipse (Lindegren et al. 2018(Lindegren et al. , 2021;;Holl et al. 2023) -although spurious solutions will be more abundant at low SNR.Gaia BH3-like systems can thus quite plausibly be detected with significantly lowermass secondaries, and to larger distances, than Gaia BH3 itself.
Unless Gaia BH3 is a major statistical fluke, these considerations suggest that (a) more massive BHs remain to be discovered around nearby metal-poor stars, and (b) the occurrence rate of BH companions is considerably higher at low metallicity than at solar values.
5. DISCUSSION 5.1.Formation from an isolated binary I have shown that wide BH + stellar binaries with orbits similar to Gaia BH3 can in principle be produced through isolated binary evolution, without the BH progenitor ever interacting with its companion.However, the orbit of Gaia BH3 in particular appears unlikely to have formed from this channel if we accept that the binary was born in the progenitor of the ED-2 stream and remained bound to it when the BH formed: matching the system's observed period and eccentricity while accommodating the BH progenitor as a red supergiant requires a natal kick of velocity ≳ 20 km s −1 , and such a kick would have ejected the binary from a low-mass cluster.Given the old age of the systems and the short lifetime of the BH progenitor, it seems unlikely that the binary was ejected from the stream with a significant kick and still remains closely associated with the stream in phase space The main uncertainty in these calculations is in the maximum radii predicted by the evolutionary models.For example, Figure 1 shows that rapidly-rotating BH progenitors could have undergone CHE and never overflowed their Roche lobes in any plausible pre-SN orbit of Gaia BH3.The radii of red supergiants are also sensitive to a variety of uncertain inputs, such as the convective mixing length and the treatment of overshooting and semiconvection (Chun et al. 2018;Schootemeijer et al. 2019;Goldberg & Bildsten 2020;Klencki et al. 2021;Romagnolo et al. 2023).Indeed, Iorio et al. (2024) perform similar calculations to those presented here but employ evolutionary tracks from the MIST project (Choi et al. 2016), which predict maximum radii of only ∼ 200 R ⊙ .Based on these models, they find that Gaia BH3's properties are consistent with formation from an isolated binary.The modeling choices adopted here produce red supergiant radii in reasonably good agreement with observations of massive stars in the LMC and SMC (e.g.Schootemeijer et al. 2019), which are the lowestmetallicity galaxies currently observable with a large population of evolved massive stars.The maximum radii they predict at [M/H] = −2.2 are not significantly different at fixed mass from those predicated at [M/H] = −0.7, the metallicity of the SMC.No one has ever observed a massive star with a metallicity as low as Gaia BH3, so the evolutionary models are untested in the lowestmetallicity regime.

Dynamical formation
With formation from a primordial binary unlikely, the alternative is dynamical assembly in the cluster that formed ED-2.Cluster formation models have also previously been found promising for forming the other two BH binaries discovered so far from Gaia data (Rastello et al. 2023;Di Carlo et al. 2023;Tanikawa et al. 2024b), and for BH binaries found in globular clusters with spectroscopic surveys, some of which have similar parameters to the Gaia BHs (Kremer et al. 2018;Gieles et al. 2021).
Recent modeling by Marín Pina et al. (2024) predicts a population of binaries with properties similar to Gaia BH3 that escape from Monte Carlo cluster models.They predict a separation distribution for dynamically-formed binaries that peaks at close separations but is relatively flat in log period over P orb = 10 0−4 d.All the clusters in the library analyzed by Marín Pina et al. (2024) have total masses significantly larger than Balbinot et al. (2024) estimate for the progenitor of the ED-2 stream, so it would be useful to model the formation of Gaia BH3-like binaries in lower-mass clusters.The most extensive effort in this direction thus far was undertaken by Tanikawa et al. (2024b,a), who also predict a period distribution that is relatively flat in log space for clusters with M cluster = 10 3 M ⊙ .A full exploration of the low-mass cluster regime is still warranted.
Marín Pina et al. (2024) predict that the formation efficiency of wide BH binaries per unit mass increases as ∼ M −1.2 cluster , such that low-mass clusters will contribute the large majority of BH+star binaries and the predicted number of detectable binaries depends sensitively on the behavior of the adopted cluster mass function at the lowmass end.Integrating down to M cluster = 10 2 M ⊙ (which may be optimistic, since the lowest-mass clusters may not form any BHs or be dense enough for dynamical interactions to be important), they predict ≲ 1 low-metallicity BH binary within 1 kpc.The empirical rate considerations in Section 4.2 suggests that the true number of such binaries may be larger than this, but the Poisson uncertainty associated with N = 1 detection still allows for the possibility that Gaia BH3 was a lucky discovery.

The possibility of pollution
It may be possible to discriminate between dynamical and isolated binary evolution channels for future BH discoveries on the basis of the abundances of the companion stars.Gaia Collaboration et al. (2024) and Balbinot et al. (2024) show that Gaia BH3's companion does not show any significant abundance anomalies, and they suggest that this favors a dynamical formation channel.In support of this possibility, some of the neutron star binary candidates identified by El-Badry et al. (2024a) do show abundance anomalies, though their origin is still uncertain.
For Gaia BH3, I find that pollution from SN ejecta would likely not be detectable, even if the system had formed from isolated binary evolution.For a plausible separation of 10 au and companion star radius of 1 R ⊙ at the time of the SN, the companion would have subtended a fraction f ∼ 5 × 10 −8 of the sky as viewed from the BH progenitor.Assuming 20 M ⊙ of isotropic ejecta, only ∼ 10 −6 M ⊙ of ejecta -likely consisting mostly of oxygen (Thielemann et al. 1996) -would have been deposited on the star.Since the star is now a giant, its convective zone extends deep into the interior: from a MESA model, I find that the outer ∼ 0.35 M ⊙ is likely convective.If ∼ 7 × 10 −7 M ⊙ of oxygen had been accreted, the ejecta would have been mixed with ∼ 3 × 10 −5 M ⊙ of oxygen already in the envelope, changing the star's surface abundances by at most a few percent, which is unlikely to be detectable.Chances of detecting pollution would be much higher if the companion star were still on the main sequence, where the outer convective zone is predicted to contain < 10 −3 M ⊙ .

Future prospects
Gaia BH3 is just at the separation beyond which BH+star binaries can avoid a common envelope.Irrespective of whether Gaia BH3 in particular formed from an isolated binary, the models presented here imply that isolated binary evolution produces binaries with orbits similar to Gaia BH3, and future Gaia data releases are likely to find these systems.This makes the system different from Gaia BH1 and BH2, which most studies have found difficult to produce through isolated binary evolution under standard assumptions (El-Badry et al. 2023a;Rastello et al. 2023;Tanikawa et al. 2024b, though see Kotko et al. 2024 for a different view).
Constraints on the separation distribution of BH binaries from DR4 will provide further insights into their dominant formation channel.If these systems form mainly by isolated binary evolution, models predict the BH + low mass star period distribution to rise steeply at P orb ≳ 10 yr (e.g.Breivik et al. 2017;Chawla et al. 2022), the period beyond which a significant fraction of binaries will have avoided a common envelope inspiral.For dynamical formation channels, the period distribution will depend on the mass of the clusters in which most systems are formed, but there is little reason to expect a steep increase at P orb ≳ 10 yr, and models published so far predict a period distribution that is flat or declines at longer periods (Tanikawa et al. 2024b;Marín Pina et al. 2024).Fortunately, Gaia will obtain an 11-year observational baseline (Gaia Collaboration et al. 2016), enabling robust orbital constraints at periods up to ∼ 20 years, and perhaps even significantly longer (e.g.Andrews et al. 2023).
The metallicity distribution of BH binaries may also help constrain how they formed.The vast majority of clusters that have contributed stars to the solar neighborhood had metallicities near solar.If dynamical interactions form most wide BH binaries, most systems should likely have metallicities near solar, similar to Gaia BH1 and BH2 (El-Badry et al. 2023a,b).For example, Tanikawa et al. (2024a) predict that the formation efficiency per unit mass of dynamically formed BH + star binaries increases by less than a factor of ∼ 5 between [M/H] = 0 and [M/H] = −2.Curiously, the empirical data so far suggests a significant increase in the prevalence of both BH and neutron star companions at low metallicity (El-Badry et al. 2024b,a).
We do not know how many other BHs in the pre-release Gaia data were not published.Because Gaia BH3 has too long of an orbital period to have been accessible in DR3 -and the BH binary period distribution may rise steeply at long periods, where binaries can avoid interaction and a common envelope -population inference is still uncertain.But prospects for constraining this population in the coming years are bright.
on en ve lo pe lik el y 10 Fig.1.-MESAevolutionary models for low-metallicity massive stars.The model shown in black is the same in all panels: a 35 M ⊙ star with metallicity [M/H] = −2.2 and no rotation.Colored lines show models varying metallicity (upper left), initial mass (upper right), initial rotation rate (lower left), and overshooting or mixing length (lower right).Diamonds, circles, and stars mark core hydrogen exhaustion, the onset of carbon burning, and core collapse.Dashed diagonal lines show constant radius.Shaded gray region marks R > 800 R ⊙ , beyond which the star would have been unlikely to fit inside the pre-SN orbit of Gaia BH3 (Section 3).Most of the models exceed this limit at their maximum radius and thus would have overflowed their Roche lobes.The only exceptions are rapidly-rotating models that undergo chemically homogeneous evolution (CHE; lower left) and models with reduced overshooting (lower right).
3. CONSTRAINTS ON THE PROGENITOR ORBITI use Monte Carlo simulations to study the combinations of pre-supernova 1 (SN) orbits and natal kicks that could have produced an orbit similar to that of Gaia BH3.This approach followsTauris et al. (2017) andEl-Badry et al. (2024b); seeTauris & van den Heuvel (2023) for additional details about SN kicks.In brief, I simulate a large number (N = 10 7 ) of orbits and kicks, assume that the SN occurs at a random phase with a kick oriented in a random direction, and predict post-SN parameters for all orbits.Finally, I analyze the initial parameters of the simulations for which the final period and eccentricity are close to the observed values for Gaia Fig.3.-Constraints on kicks and the pre-SN orbit of Gaia BH3-like binaries.I assume that the binary had a circular orbit before the SN and that the BH received a kick with velocity v kick when it formed.I simulate 10 7 pre-SN orbits with a uniform period distribution P orb ∼ U (0, 2P orb,BH3 ), a pre-SN mass distribution M 2, init /M ⊙ ∼ U (33, 100), and a kick velocity distribution v kick /(km s −1 ) ∼ U (0, 100).Black histograms show the properies of all the systems whose post-SN orbits have periods and eccentricities similar to Gaia BH3 (upper left).Red histograms show the subset of these that result in a post-SN systemic velocity vsys < 10 km s −1 , as required to remain bound to the ED-2 stream progenitor.Orbits with a broad range of pre-SN separations can produce the observed period and eccentricity.However, most of these require significant kicks and result in vsys > 10 km s −1 .Simulations with vsys < 10 km s −1 can only produce the observed period and eccentricity for a init ≲ 1200 R ⊙ .The BH progenitor would overflow its Roche lobe in such orbits for R ≳ 800 R ⊙ .
Fig.4.-Similar to Figure3, but for eccentric pre-SN orbits.I simulate 10 7 orbits with initial eccentricity e init ∼ U (0, 1) and the SN occurring at a random phase.Other parameters are distributed as in Figure3.d peri, init is the separation between the BH progenitor and companion star at periastron in the pre-SN orbit.Most orbits with vsys < 10 km s −1 (red) have d peri, init ≲ 1100 R ⊙ and thus would overflow their Roche lobes if the BH progenitor reached R ≳ 800 R ⊙ .