11institutetext:  Physics and Astronomy Department Galileo Galilei, University
of Padova, Vicolo dell’Osservatorio 3, I–35122, Padova, Italy 22institutetext:
INFN - Padova, Via Marzolo 8, I–35131 Padova, Italy 33institutetext: INAF –
Padova, Vicolo dell’Osservatorio 5, I-35122 Padova, I 44institutetext:
Institut für Theoretische Astrophysik, ZAH, Universität Heidelberg, Albert-
Ueberle-Straße 2, D-69120, Heidelberg, Germany 55institutetext: Niels Bohr
International Academy, Niels Bohr Institute, Blegdamsvej 17, 2100 Copenhagen,
Denmark 66institutetext: Departament de Física Quàntica i Astrofísica,
Institut de Ciències del Cosmos, Universitat de Barcelona, Martí i Franquès 1,
E-08028 Barcelona, Spain 77institutetext: SISSA, via Bonomea 365, I-34136
Trieste, Italy 88institutetext: Univ Lyon, Univ Lyon1, ENS de Lyon, CNRS,
Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-Laval,
France 99institutetext: INAF, Osservatorio Astronomico di Bologna, Via Gobetti
93/3, I-40129 Bologna, Italy 1010institutetext: Gran Sasso Science Institute,
Via F. Crispi 7, I-67100 L’Aquila, Italy 1111institutetext: INFN - Laboratori
Nazionali del Gran Sasso, I-67100 L’Aquila, Italy 1212institutetext: INAF -
Osservatorio Astronomico d’Abruzzo, Via Mentore Maggini, s.n.c., I-64100
Teramo, Italy

# The boring history of Gaia BH3 from isolated binary evolution

Giuliano Iorio<EMAIL_ADDRESS>Stefano Torniamenti
<EMAIL_ADDRESS>Michela Mapelli mapelli@uni-
heidelberg.de44112233 Marco Dall’Amico 1122 Alessandro A. Trani 55 Sara
Rastello 66 Cecilia Sgalletta 2277 Stefano Rinaldi 1144 Guglielmo Costa 8833
Bera A. Dahl-Lahtinen 11 Gastón J. Escobar 1122 Erika Korb 1122 M. Paola
Vaccaro 44 Elena Lacchin 11224499 Benedetta Mestichelli 4410101111 Ugo N. Di
Carlo 77 Mario Spera 77 Manuel Arca Sedda, 101011111212

(Received September 15, 1996; accepted March 16, 1997)

Gaia BH3 is the first observed dormant black hole (BH) with a mass of
$\approx{30}$ M⊙ and represents the first confirmation that such massive BHs
are associated with metal-poor stars. Here, we explore the isolated binary
formation channel for Gaia BH3 focusing on the old and metal-poor stellar
population of the Milky Way halo. We use the mist stellar models and our open-
source population synthesis code sevn to evolve $3.2\times 10^{8}$ binaries
exploring 16 sets of parameters, that encompass different natal kicks,
metallicities, common-envelope efficiencies, and models for the Roche-lobe
overflow. We find that systems like Gaia BH3 form preferentially from binaries
initially composed of a massive star ($40-60$ M⊙) and a low mass companion
($<1$ M⊙) in a wide ($P>10^{3}$ days) and eccentric orbit ($e>0.6$). Such
progenitor binary stars do not undergo any Roche-lobe overflow episode during
their entire evolution, so that the final orbital properties of the BH-star
system are essentially determined at the core collapse of the primary star.
Low natal kicks ($\approx$ 10 km/s) significantly favour the formation of Gaia
BH3-like systems, but high velocity kicks up to $\approx 220$ km/s are also
allowed. We estimate the formation efficiency for Gaia BH3-like systems in old
($t>10$ Gyr) and metal-poor ($Z<0.01$) populations to be $4\times 10^{-8}$
M${}_{\odot}^{-1}$ (for our fiducial model), representing $\approx 3\%$ of the
whole simulated BH-star population. We expect up to $\sim 3000$ BH-star
systems in the Galactic halo formed through isolated evolution, of which $\sim
100$ are compatible with Gaia BH3. Gaia BH3-like systems represent a common
product of isolated binary evolution at low metallicity ($Z<0.01$), but given
the steep density profile of the Galactic halo we do not expect more than one
at the observed distance of Gaia BH3. Considering the estimated formation
efficiency for both the isolated and dynamical formation channel, we conclude
that they are almost equally likely to explain the origin of Gaia BH3.

###### Key Words.:

black hole physics – methods: numerical – binaries: general – stars: black
holes – stars: massive – Galaxy: halo

##  1 Introduction

Gaia data (Gaia Collaboration et al. 2023) offer a unique opportunity to
search for dormant black holes (BHs), i.e. X-ray and radio-quiet BHs that are
members of a binary system with a star (Guseinov & Zel’dovich 1966; Trimble &
Thorne 1969). Gaia BH3, the third dormant BH retrieved from these data, has a
mass of $\sim{33}$ M⊙ and is member of a loose eccentric binary system
(orbital period $P\approx{4000}$ days, eccentricity $e\approx{0.72}$) together
with a $\approx{0.76}$ M⊙ very metal-poor giant star
($[\textrm{Fe}/\textrm{H}]=-2.56\pm{0.12}$ from the UVES spectrum, Gaia
Collaboration et al. 2024).

Gaia BH3 is the first ever discovered dormant BH with mass $\gtrsim{}30$ M⊙,
and the first BH with such a high mass unequivocally observed in the Milky
Way. The other two Gaia BHs (BH1 and BH2) both have a mass of $\sim{10}$ M⊙
(El-Badry et al. 2023b, a; Tanikawa et al. 2024). All of the other dormant BH
candidates (Casares et al. 2014; Ribó et al. 2017; Giesers et al. 2018, 2019;
Thompson et al. 2019; Shenar et al. 2022a, b; Mahy et al. 2022; Lennon et al.
2022; Saracino et al. 2022, 2023) and the BHs in X-ray binary systems with a
dynamical measurement have mass between a few and $\approx{20}$ M⊙ (Oëzel et
al. 2010; Miller-Jones et al. 2021).

BHs with mass $\gtrsim{30}$ M⊙ have been observed for the first time with
gravitational waves (Abbott et al. 2016a, b, 2021a, 2021c, 2021b). According
to the main formation scenario, such massive BHs are expected to form from
metal-poor stars, which do not lose much mass by stellar winds and collapse to
a BH directly (Heger & Woosley 2002; Mapelli et al. 2009, 2010, 2013; Zampieri
& Roberts 2009; Belczynski et al. 2010; Fryer et al. 2012; Ziosi et al. 2014;
Spera et al. 2015). Gravitational-wave data have proved the existence of such
”oversize” BHs, but they do not carry any information about the progenitor’s
metallicity, because we did not observe any electromagnetic counterpart.
Hence, Gaia BH3 is the first direct confirmation that BHs with mass in excess
of 30 M⊙ are associated with metal-poor stellar populations, supporting the
main theoretical scenario.

For both its chemical composition and Galactic orbit, Gaia BH3 is associated
with the ED–2 stream, composed of stars that are at least 12 Gyr old (Balbinot
et al. 2024). Thus, this is the first BH unambiguously associated with a
disrupted star cluster with a mass $2\times{}10^{3}-4.2\times{}10^{4}$ M⊙
(Balbinot et al. 2024). Preliminary numerical models show that Gaia BH3 might
have assembled dynamically from a low-mass star and a BH that have not evolved
in the same binary star system (Marín Pina et al. 2024, Arca Sedda et al.
2024, in prep.).

While the association of Gaia BH3 with the ED–2 stream is indisputable, here
we show that this system is compatible with the nearly unperturbed evolution
of a primordial binary star: our population-synthesis simulations show that we
do not need to invoke a dynamical exchange or a dynamical capture to explain
its orbital properties (see also El-Badry 2024).

Gaia BH3 is very different from both Gaia BH1 and Gaia BH2, which challenge
our models of binary evolution (El-Badry et al. 2023b, a; Tanikawa et al.
2023; Rastello et al. 2023; Di Carlo et al. 2024; Generozov & Perets 2024;
Kotko et al. 2024): its wide orbital separation does not require evolution via
Roche-lobe overflow or common envelope. Rather, it is consistent with a wide
system that never went through mass transfer and survived the formation of the
BH because this happened via direct collapse, with a mild natal kick. We show
that the measured orbital eccentricity helps to constrain the phase of the
binary at the time of the formation of the BH. If formed in a star cluster,
Gaia BH3 might have evolved nearly unperturbed inside it, until it evaporated
or was ejected. The paper is organized as follows. We describe our methodology
in Section 2, present the main results in Section 3, and discuss the
implications in Section 4. Section 5 summarizes our main takeaways.

##  2 Methods

In this work, we compare the observational features of the Gaia BH3 system
with a synthetic population of stellar binaries produced with the Stellar
EVolution for N-body (sevn) code (Spera & Mapelli 2017; Spera et al. 2019;
Mapelli et al. 2020; Iorio et al. 2023). In sevn, the evolution of single
stars is estimated by interpolating on the fly across pre-calculated stellar
evolution tracks stored in custom lookup tables. The models presented in this
study are computed using the mist stellar tracks by Choi et al. (2016). We
used trackcruncher (Iorio et al. 2023) to produce the sevn tables from the
mist stellar tracks. We produce eight sets of tables from metallicity
$Z=1.4\times 10^{-5}$ to $4.5\times 10^{-2}$; each set contains stellar tracks
from 0.71 M⊙ to 150 M⊙.

In sevn, binary interactions are described through analytical and semi-
analytical prescriptions, encompassing stable mass transfer by Roche-lobe
overflow and winds, common envelope evolution, angular momentum dissipation by
magnetic braking, tidal interactions, orbital decay by gravitational-wave
emission, dynamical hardening, chemically homogeneous evolution, and stellar
mergers. In our models, we treat tidal forces with the analytical formalism
developed by Hut (1981), as implemented by Hurley et al. (2002). For more
details, we refer to Iorio et al. (2023).

In this work, we use the sevn version V 2.8.0 (commit 3375978) publicly
available at the GitLab repository https://gitlab.com/sevncodes/sevn. The code
trackcruncher used to produce the mist table is available at
https://gitlab.com/sevncodes/trackcruncher (commit 0827337).

###  2.1 Initial conditions

We have generated the initial conditions for $2\times 10^{7}$ binary systems
by extracting the zero-age main-sequence mass of the primary (more massive)
star, $m_{1}$, from the Kroupa initial mass function (Kroupa 2001) from 5 to
150 M⊙. The mass of the secondary (less massive) star is set from the binary
mass ratio $q$ following the distribution PDF$(q)\propto q^{-0.1}$ with
$q\in\left[0.71\,{}\textrm{M}_{\odot}/m_{1},1\right]$. We have drawn the
initial orbital period ($P$) and eccentricity ($e$) from the following
distributions: PDF$\left(\log P/\mathrm{days}\right)\propto\left(\log
P\right)^{-0.55}$, PDF$(e)\propto e^{-0.42}$. These distributions of mass
ratios, orbital periods and eccentricities are fitting formulas to
observations of massive binary stars in Galactic young clusters (Sana et al.
2012). The periods range from 1.4 days to 866 years, while the eccentricities
range from 0 to $e_{\mathrm{max}}=1-\left(P/2\ \mathrm{days}\right)^{-2/3}$,
according to Moe & Di Stefano (2017). The total mass of the simulated
population is $M_{\rm{sim}}=4.2\times 10^{8}\ \mathrm{M}_{\odot}$
corresponding to a mass of $M_{\rm{tot}}=1.4\times 10^{9}\ \mathrm{M}_{\odot}$
when accounting for the fact that we do not simulate low-mass primary stars
and systems with secondary stars less massive than 0.71 M⊙111We assume a
parent population with primary masses $m_{1}$ between 0.08 M⊙ and 150 M⊙, and
mass-ratios in the range $q\in[0.08\,{}\textrm{M}_{\odot}/m_{1},1]$.. We use
sevn to evolve all the generated binaries from the zero age main sequence up
to to the moment in which the two stars are compact remnants (white dwars,
neutron stars, BHs).

###  2.2 Fiducial model

Figure 1:  Hertzsprung–Russell diagram of a sample of stars at metallicity
$Z=4.1\times 10^{-5}$ interpolated by sevn from the mist stellar tracks. The
blue and red lines show the likely progenitors of the BH and stellar companion
in Gaia BH3, respectively.

In our fiducial model, we use the default sevn options as described by Iorio
et al. (2023) (see also Sgalletta et al. 2023; Costa et al. 2023) except for
the circularisation option. Here, we use the new default circularisation
option _angmom_ in which at the onset of the Roche-lobe overflow the orbit is
circularised conserving the angular momentum, so the new semi-major axis is
$a_{\mathrm{new}}=a_{\mathrm{old}}(1-e^{2}_{\mathrm{old}})$.

The stability of the Roche-lobe mass transfer is checked by comparing the
current binary mass ratio to a critical mass ratio that depends on the stellar
evolutionary phase. In this work, we use the default option of sevn (Iorio et
al. 2023), in which Roche-lobe overflow episodes involving stars with
radiative envelope (main sequence stars and low-mass stars before climbing the
red giant branch) are always stable, while in all the other cases we use the
value of the critical mass ratio by Hurley et al. (2002). When the mass
transfer is unstable, a common envelope evolution phase begins (Ivanova et al.
2013). We use the standard $\alpha-\lambda$ prescription for the common
envelope (Webbink 1984; Hurley et al. 2002, and reference therein) with the
ejection efficiency parameter $\alpha=1$, whereas we estimate $\lambda$ using
the formalism by Claeys et al. (2014), as implemented in sevn (see Appendix
A1.4 in Iorio et al. 2023).

After a core-collapse supernova, we decide the mass of the BH based on the
delayed model by Fryer et al. (2012). We adopt the supernova kick model by
Giacobbo & Mapelli (2020), where the natal kick is drawn from a Maxwellian
distribution with one-dimensional root-mean-square $\sigma_{\mathrm{rms}}=265\
\mathrm{km/s}$ (Hobbs et al. 2005), re-scaled by the ejecta mass and the
compact-remnant mass resulting in a much lower ($\lessapprox$ 10 km/s)
effective kick for the massive BHs ($>30$ M⊙).

In the fiducial model, we set the metallicity to $Z=4.1\times 10^{-5}$,
corresponding to the observed value of $[\textrm{Fe}/\textrm{H}]=-2.56$
considering the mist stellar models. The mist stellar tracks with the fiducial
metallicity are shown in Figure 1.

###  2.3 Additional models: parameter-space exploration

To explore the parameter space for the formation of Gaia BH3-like systems, we
ran an additional set of 15 simulations varying one parameter of the fiducial
simulation per time. For the initial conditions, we used the same set of
$2\times{}10^{7}$ binary systems as for the fiducial model.

We have simulated eight additional metallicities: $Z=0.0001$ (model Z1E-4),
$Z=0.0004$ (model Z4E-4), $Z=0.0008$ (model Z8E-4), $Z=0.001$ (model Z1E-3),
$Z=0.004$ model (Z4E-3), $Z=0.006$ (model Z6E-3), $Z=0.008$ (model Z8E-3), and
$Z=0.01$ (model Z1E-2). We have also explored the effect of the supernova kick
by running four different simulations in which the kicks are drawn from a
Maxwellian distribution with $\sigma_{\mathrm{rms}}=1\ \mathrm{km/s}$
(essentially no kick, model VK1), $50\ \mathrm{km/s}$ (VK50), $100\
\mathrm{km/s}$ (VK100), and $265\ \mathrm{km/s}$ (VK265). In these models, the
kick is not re-scaled for the ejected mass and the mass of the compact
remnant. We also test the importance of the common-envelope phase by running
two simulations varying the common-envelope efficiency: $\alpha=3$ (model
CE$\alpha$3) and $\alpha=5$ (model CE$\alpha$5). Higher common-envelope
efficiency enhances the likelihood of system survival without merger, leading
to increased post-common envelope orbital separation.

Finally, we have run an additional model (RLOp) in which we allow the system
to start a Roche-lobe overflow in an eccentric orbit if the condition of
Roche-lobe overflow is satisfied at the periastron, i.e. if the stellar radius
is larger than $R_{\mathrm{L}}=r_{\mathrm{peri}}f(q)$, where
$r_{\mathrm{peri}}$ is the periastron and $f(q)$ is a function that depends on
the binary mass ratio (Eggleton 1983).

##  3 Results

Figure 2: Properties of non-interacting BH-star systems in the fiducial model.
Here, we only show stars in the giant phase, i.e. outside the main sequence
and that have not started the core helium burning phase yet. Also, we do not
show systems that undergo a Roche-lobe overflow episode. From top to bottom
and from left to right: orbital period versus companion star mass; BH mass
versus companion star mass; eccentricity versus orbital period; eccentricity
versus companion star mass. The color map shows the integrated time spent by
the systems in a given bin, thus it is a proxy of the likelihood to observe a
system with such properties. The map does not include information about the
stellar luminosity, which is also important for detectability. We assume a
metallicity $Z=4.1\times 10^{-5}$. The star indicates the properties of Gaia
BH3, while the light blue boxes show the criteria used to select Gaia BH3-like
systems (hereafter, Gaia BH3box). The color maps saturate at 98% of the
maximum value.

Figure 2 shows the properties of the BH-star systems in the fiducial model. We
consider all the simulation outputs in which a binary system is composed of a
BH and a giant star, i.e. a star that left the main sequence and has not
started the core helium burning phase yet. We impose the additional condition
that the system is not undergoing a common-envelope or Roche-lobe overflow
episode. The color map of Fig. 2 represents the time spent by all the systems
in such configuration, which is a proxy of the likelihood to observe a given
system.

Binaries hosting a BH and a massive companion are, in general, more common in
the simulations, but the short lifetime of the massive star makes them less
likely to be observed. In contrast, BH-low mass star binaries are rarer, but
once formed, they can remain in the same configuration for several Gyrs.

Gaia BH3 is located in a region of the parameter space well covered by the
fiducial simulation, especially if we consider binary systems older than 10
Gyr. In addition, the Gaia BH3 properties (mostly period and BH mass) cover a
range in which Gaia is particularly sensible to astrometric signals of
binarity (see e.g. Penoyre et al. 2022). From Figure 2, we can conclude that
if a BH-star system is observed in an old and metal-poor population, it is not
surprising to find it with the Gaia BH3 configuration. In this sense, Gaia BH3
is not an outlier, unlike Gaia BH1 and BH2 (El-Badry et al. 2023b; Rastello et
al. 2023).

It is convenient to define a region of the observed parameter space to
identify Gaia BH3-like systems. For the sake of simplicity, and to help the
comparison with previous work, we use the same selection cut used in Marín
Pina et al. (2024): $P\in[400,40000]\ \mathrm{days}$, $m_{s}\in[0.2,1]\
\mathrm{M}_{\odot}$, $m_{\mathrm{bh}}\in[25,45]\ \mathrm{M}_{\odot}$,
$e\in[0.63,0.83]$. The box represents about 1% of the volume of the parameter
space that includes BH-star systems at ages older than 10 Gyr.

###  3.1 Formation history

All the Gaia BH3-like systems in our models have the same simple formation
history. Their progenitors consist of a massive star ($>40$ M⊙) and a low mass
companion ($<1$ M⊙). The binary components never interact through Roche-lobe
overflow, so the orbital properties at BH formation are almost the same as
their initial ones. The final properties of the BH-star system are then
determined by the binary phase at the moment of the supernova explosion, the
amount of mass lost during the evolution, and the BH natal kick.

In our fiducial model, the initial BH progenitor mass is about 40–50 M⊙ and
the ejected mass at BH formation is 10–20 M⊙. The mild BH natal kick ($\sim
10$ km/s) cannot drastically change the binary orbital parameters. As a
consequence, the large majority of BH progenitors originated with very similar
orbital properties to Gaia BH3: a period between $10^{3}$ and $10^{4}$ days
and an eccentricity in the range $0.6-0.9$. The system more similar to Gaia
BH3 has an initial period of about 1500 days, an eccentricity of 0.65, stellar
mass of 0.83 M⊙ and initial BH progenitor mass of 41 M⊙. We find a sub-
population ($<10$% of the total) that produces Gaia BH3-like systems starting
from almost circular orbits ($e<0.2$) and lower initial periods ($<10^{3}$
days). For these systems, the final orbital properties are determined by the
high value of the BH natal kick. A few systems ($\approx{}3\%$) have very
massive BH progenitors of about 120 M⊙, which produce BHs with a mass of
$\sim{30}$ M⊙ because of mass loss during pulsational pair instability (Heger
& Woosley 2002; Belczynski et al. 2016; Woosley 2017, 2019; Farmer et al.
2019; Costa et al. 2021; Tanikawa et al. 2021).

Since Gaia BH3-like progenitors do not experience any binary interaction,
there are no differences for the model CE$\alpha 3$ and CE$\alpha 5$. Also, we
obtain the same results for the model RLOp, which indicates that including or
not the possibility to trigger a stable mass transfer at the periastron is not
important for the production of such binary systems. We find the same
progenitors in the case of the kick model VK1, as a confirmation that in the
fiducial model the implemented kick model produces very low kicks for the case
of Gaia BH3-like progenitors.

Models with higher kick velocities show significant differences in the initial
orbital properties. The BH natal kicks in models VK50 and VK100 are able to
wash out any correlation between the initial and final period and
eccentricity. The Gaia BH3-like progenitors have initial periods and
eccentricities almost uniformly distributed in the range 100–4000 days, and
0–0.8, respectively. At very high velocities (model VK265) the formation of
Gaia BH3-like binaries is strongly suppressed. Most systems break during BH
formation and just a few of them (4 over 20 million simulated binaries) had a
lucky configuration of binary phase and/or velocity kick that let them survive
the supernova event.

We find a trend with metallicity for the initial orbital period, eccentricity
and BH progenitor mass. From low metallicity (model Z1E-4) to intermediate-
high metallicity (model Z1E-2) the lower limit of the initial period moves
toward higher values ($P>1000$ days for model Z1E-2), while the sub-
populations with initial low eccentricity ($e<0.6$) tend to disappear. As for
the BH progenitor mass, the range of values widens with increasing
metallicity, up to initial masses of 100 M⊙. This is due to the non-monotonic
trend of the initial BH mass relation of the mist tracks used in this work
when paired with the rapid and delayed supernova models by Fryer et al.
(2012). When assuming initial masses lower than 60 M⊙, the combination of the
Fryer et al. (2012) supernova models and the mist tracks produces BHs more
massive than 25 and 30 M⊙ up to a metallicity $Z=0.02$ and $0.01$,
respectively. Relaxing the limit on the initial masses, BHs more massive than
30 M⊙ are excluded only for $Z>0.02$.

Figure 3: Distribution of the initial conditions of the binaries that
reproduce the Gaia BH3 properties after the BH formation. From top left to
bottom right, the panels show the initial binary period, initial eccentricity,
the module of the BH natal kick velocity, and initial mass of the BH
progenitor. We show the distributions obtained assuming different Maxwellian-
distribution priors for the kick velocity: $\sigma_{v}=1\ \mathrm{km/s}$
(blue, thick solid), $\sigma_{v}=10\ \mathrm{km/s}$ (orange,thick dashed),
$\sigma_{v}=50\ \mathrm{km/s}$ (green, dotted), $\sigma_{v}=100\
\mathrm{km/s}$ (red, thin solid), $\sigma_{v}=256\ \mathrm{km/s}$ (violet,
thin dashed).

###  3.2 Bayesian sampling

The sevn simulations indicate that the region of Gaia BH3 in the observed
parameter space ($P$, $e$, $m_{\mathrm{BH}}$, $m_{\mathrm{ms}}$) tends to be
populated by systems in which the binary components never interact. Therefore,
the final binary properties are determined by the initial binary parameters
and the details of the BH formation only. In order to focus on Gaia BH3
progenitors independently of the chosen selection cut, we implement a simple
Bayesian sampling technique in which we sample the posteriors of the binary
initial condition using as likelihood a multivariate normal distribution
centred on the Gaia BH3 value and with standard deviation equal to the
reported errors (assumed independent).

For the parameter prior we use an uniform distribution for the mean anomaly
(between $0$ and $2\pi$) at the moment of the BH formation, an uniform
distribution for the mass of the BH progenitor (between 20 and 100 M⊙) and
eccentricity (between 0 and 1), a log-uniform distribution for the period
(between 10 and $10^{6}$ days), and a normal distribution for the three
spherical components of the velocity kick, corresponding to a Maxwellian
distribution of the velocity module.

We test five different values of the one-dimensional root-mean-square velocity
for the Maxwellian distribution ($1,10,50,100,256$ km/s). We also exclude
interacting binary systems, by setting the prior to be equal to 0 when the
maximum radius of the BH progenitor is larger than either the separation at
periastron or the Roche-lobe radius at periastron (see Section 2). To estimate
the final properties of the binary we use the formalism in Appendix of Hurley
et al. (2002) and we estimate the maximum radius and the final mass of the BH
progenitor by using sevn. The posteriors are sampled using the Monte-Carlo
Markov Chain (MCMC) Ensemble Sampler available in the emcee python package
(Foreman-Mackey et al. 2013). This method is similar to the grid search
presented by El-Badry (2024), but our method is tailored to reproduce exactly
the observed properties of Gaia BH3, exploring widely the parameter space. In
particular, we do not assume initial circular orbits and this has a
significant impact on the final results. Since we include the reduction of the
parameter space due to stellar interactions directly in the priors, the
posteriors do not need to be corrected a posteriori.

We show the posterior samples in Figure 3. Overall, they are consistent with
what obtained from the sevn simulations for the wider space of the selection
cut assumed in Sect. 3.1. Assuming low kick velocities, the Gaia BH3
progenitor is very likely to originate from a wide and eccentric binary
($P\approx 3500$ days, $e\sim 0.75$). The smaller peak at around $10^{3}$ days
and $e\sim 0.4$ is due to systems that form the BH at periastron, where it is
more favourable to produce large final eccentricities due to the impact of
mass loss during the supernova explosion. Larger kicks widen the range of
initial periods and eccentricities, with most of the initial periods around
500–600 days. The period distribution is sharply cut below 300 days because
the objects are tight enough to interact. It is possible to form a Gaia BH3
system with initial kicks up to 200 km/s. For higher values, the kick is
larger than the escape velocity at periastron in the whole range of periods
consistent with a non-interacting binary.

Independently of the explored parameters, a clear picture emerges. In presence
of low supernova kicks ($<10$ km/s), Gaia BH3 like stystems form from non
interacting binary progenitors preferentially in a wide
($500<P/\mathrm{days}<2\times{}10^{4}$) and eccentric configuration.
Alternatively, on the off-chance that the binary survives to higher supernova
kicks, systems with Gaia BH3 properties can form preferentially from tighter
($P/\mathrm{days}\sim 500$) systems almost independently of their initial
eccentricity if $e<0.7$. For this reason, the initial properties of the binary
can play a key role in setting the formation efficiency of Gaia BH3-like
systems.

###  3.3 Halo population

Model | $\eta_{\mathrm{GaiaBH3-like}}$ |  | $\eta_{\mathrm{GaiaBH3-like}}$
---
(t>10 Gyr)
| $\eta_{\mathrm{BHStar}}$
---
(t>10 Gyr)
fiducial | $(1.1\pm 0.1)\times 10^{-7}$ | $(4.7\pm 0.6)\times 10^{-8}$ | $(1.5\pm 0.1)\times 10^{-6}$
RLOp | $(1.1\pm 0.1)\times 10^{-7}$ | $(4.0\pm 0.5)\times 10^{-8}$ | $(1.5\pm 0.1)\times 10^{-6}$
VK1 | $(1.1\pm 0.1)\times 10^{-7}$ | $(4.3\pm 0.6)\times 10^{-8}$ | $(2.2\pm 0.1)\times 10^{-6}$
VK50 | $(9.0\pm 0.8)\times 10^{-8}$ | $(3.9\pm 0.5)\times 10^{-8}$ | $(5.0\pm 0.2)\times 10^{-7}$
VK100 | $(4.0\pm 0.5)\times 10^{-8}$ | $(1.3\pm 0.3)\times 10^{-8}$ | $(5.0\pm 0.2)\times 10^{-7}$
VK265 | $(6.4\pm 2.1)\times 10^{-9}$ | $(2.1\pm 1.2)\times 10^{-9}$ | $(5.4\pm 0.2)\times 10^{-8}$
CE$\alpha$3 | $(1.1\pm 0.1)\times 10^{-7}$ | $(4.7\pm 0.6)\times 10^{-8}$ | $(2.1\pm 0.1)\times 10^{-6}$
CE$\alpha$5 | $(1.1\pm 0.1)\times 10^{-7}$ | $(4.7\pm 0.6)\times 10^{-8}$ | $(2.4\pm 0.1)\times 10^{-6}$
Z1E-4 | $(1.3\pm 0.1)\times 10^{-7}$ | $(4.8\pm 0.6)\times 10^{-8}$ | $(1.4\pm 0.1)\times 10^{-6}$
Z4E-4 | $(1.1\pm 0.1)\times 10^{-7}$ | $(4.8\pm 0.6)\times 10^{-8}$ | $(1.4\pm 0.1)\times 10^{-6}$
Z8E-4 | $(9.4\pm 0.8)\times 10^{-8}$ | $(3.9\pm 0.5)\times 10^{-8}$ | $(1.5\pm 0.1)\times 10^{-6}$
Z1E-3 | $(9.2\pm 0.8)\times 10^{-8}$ | $(4.0\pm 0.5)\times 10^{-8}$ | $(1.6\pm 0.1)\times 10^{-6}$
Z4E-3 | $(9.1\pm 0.8)\times 10^{-8}$ | $(5.3\pm 0.6)\times 10^{-8}$ | $(2.4\pm 0.1)\times 10^{-6}$
Z6E-3 | $(9.4\pm 0.8)\times 10^{-8}$ | $(5.8\pm 0.6)\times 10^{-8}$ | $(2.7\pm 0.1)\times 10^{-6}$
Z8E-3 | $(8.2\pm 0.8)\times 10^{-8}$ | $(6.2\pm 0.7)\times 10^{-8}$ | $(3.1\pm 0.1)\times 10^{-6}$
Z1E-2 | $(8.2\pm 0.8)\times 10^{-8}$ | $(7.1\pm 0.7)\times 10^{-8}$ | $(3.4\pm 0.1)\times 10^{-6}$

Table 1: Formation efficiency (M${}^{-1}_{\odot}$) obtained from different
sevn simulations. The column reports the formation efficiency for the the Gaia
BH3 selection cut (see text), the same selection cut but only for systems that
are still in the BH-star configuration after 10 Gyr, and for all the BH-star
systems with age larger than 10 Gyr. The formation efficiency assumes a 100%
binary fraction. The errors are estimated as the square root of the counts
over the total mass.

In the explored models, Gaia BH3-like systems represent one of the standard
outcomes of binary evolution. We now investigate how common Gaia BH3-like
systems and BH-star binaries are in the Galactic halo environment,
characterised by low metallicity and old ages. We also assume $t>10$ Gyr.
First, we estimate the formation efficiency, defined as the ratio of the
number of objects to the total simulated mass. We correct the total mass in
the initial conditions for the incomplete mass sampling of the primary, and
assume an initial binary fraction of 1 for the sake of simplicity. Table 1
reports the formation efficiency for all the systems in the Gaia BH3 selection
box (Figure 2), and for those that are in the selection box and in the BH-star
configuration after 10 Gyr.

In the fiducial model, the population of Gaia BH3-like objects accounts for
about 3% of the whole population at $t>10$ Gyr. This number is consistent with
the volume fraction sampled by the selection box. Models with higher kicks
significantly suppress the formation of both Gaia BH3-like and BH-star
systems, while a higher metallicity increases their formation efficiency up to
a factor of two. Given the formation history of Gaia BH3-like objects, the
efficiency is not affected by the choice of the common-envelope ejection
efficiency. However, larger values of the common-envelope efficiency parameter
$\alpha$ slightly increase the formation efficiency of the whole BH-star
systems resulting in a lower fraction of Gaia BH3-like systems.

Even if the Galactic stellar halo is the least massive component of our
Galaxy, the most recent estimates revise its mass from few $10^{8}$ M⊙ (see
e.g. Bell et al. 2008) up to $0.7-1.2\ \times 10^{9}$ M⊙ (Deason et al. 2019;
Lane et al. 2023). Combining the total halo mass with the formation
efficiencies in Table 1, we predict a number of BH-stars in the halo in the
range $\approx 1800-3000$. Among these systems, about $50-100$ fall in the
Gaia BH3 selection cut used in this work (Figure 2). The Galactic stellar halo
is a diffuse component extended up to 100 kpc with a steep density profile
($\rho\propto r^{-\gamma}$, with $2<\gamma<5$, see e.g. Medina et al. 2024;
Iorio et al. 2018, and references therein).

Deason et al. (2019) report a local density profile at Sun location of
$\rho_{\mathrm{h,\odot}}\approx 8\times 10^{4}\
\mathrm{M}_{\odot}/\mathrm{kpc}^{3}$. Multiplying this density by the
formation efficiencies of the old BH-star binaries we obtain a local density
of $n_{\mathrm{BHStar}}=\eta_{\mathrm{BHStar}}\,\rho_{\mathrm{h,\odot}}$
ranging from 0.1 $\mathrm{kpc}^{-3}$ to 0.3 $\mathrm{kpc}^{-3}$. The values
correspond to an expected number of objects between 0.12 and 0.27 at the
distance of Gaia BH3 ($\approx 0.6\ \mathrm{kpc}$). At face value, we do not
expect any BH-star binary system from the halo in the solar neighbourhood.
However, we can account for the statistical noise assuming that the number of
systems follows a Poissonian distribution with rate
$\lambda_{\mathrm{BHStar}}$ ranging between 0.12 and 0.27. In this case, the
probability to observe at least one object is between 11% and 25%. By simply
using the formation efficiency for the old systems in the Gaia BH3 selection
box (Figure 2) we obtain a Poisson rate between 0.003 and 0.006, corresponding
to the probability to observe at least one Gaia BH3-like system of 0.3–0.6%.
Although these numbers correspond to the very low probability tail of the
distribution, they cannot be rejected at 3$\sigma$ level. However, a single
future detection of another of Gaia BH3-like binary system in the Galactic
halo within the solar neighbourhood will be enough to significantly rule-out
the isolated formation channel, according to the models presented in this
work. However, this will represent a challenge also for the alternative
dynamical formation channel (see Section 4.3).

##  4 Discussion

###  4.1 Initial eccentricity and stellar tides

In all the simulations assuming our fiducial kick model (Giacobbo & Mapelli
2020), Gaia BH3-like systems form preferentially from binary systems with
initial eccentric orbit ($e>0.6$). Therefore, we expect that initial
distributions that favour eccentric orbits, like the thermal distribution,
boost the number of Gaia BH3-like systems. In contrast, circular orbits
significantly suppress the formation of such systems unless high kicks models
are used.

In some cases, initial circular binaries are justified by the effect of tides.
The sevn simulations used in this work implement the tides following the
formalism by Hurley et al. (2002) (see also Iorio et al. 2023) and show that
for Gaia BH3-like systems tides are extremely inefficient. The tidal
efficiency is at its maximum after the main sequence phase of the BH
progenitor, when the star expands and develops a convective envelope. Assuming
that the dominant tidal dissipation mechanism is due to convective motions in
the envelope (Zahn 1977; Hurley et al. 2002), and considering the mist models
at $Z=4.1\times 10^{-5}$ (Figure 1), we estimate a circularisation timescale
due to tides of $\tau_{\rm circ}\simeq 10^{25}$–$10^{38}\,\rm Gyr$. This
timescale is orders of magnitude longer than the lifetime of the binary system
and especially longer than the post-main-sequence phase, when tides are more
efficient.

###  4.2 Systematic uncertainties

Table 1 shows that the formation efficiency of BH-star systems drops
significantly for models with high supernova kicks. In the models VK100 and
VK265, we can exclude the hypothesis that the observation of Gaia BH3 is
consistent with the isolated channel at a high level of significance
($>3\sigma$). Estimates of BH kicks rely on a few tracers (e.g. X-ray
binaries) and are subjected to many uncertainties and systematics (see e.g.
Andrews & Kalogera 2022; Atri et al. 2019; Shenar et al. 2022a). The
association of Gaia BH3 with a stellar stream (Balbinot et al. 2024) indicates
that the binary system did not experience a kick much higher than the typical
escape velocity from a star cluster ($<10$ km/s) at the time of BH formation.

In our analysis, we are not considering binaries with companion stars less
massive than $0.71\ \mathrm{M}_{\odot}$. As discussed by El-Badry (2024), such
a population can represent a significant portion of the stellar companions of
Gaia BH3-like systems. To asses its relevance on the estimate of the formation
efficiency, we generate new initial conditions as in Section 2.1, but assuming
an extended mass ratio range from $0.08\,{}\textrm{M}_{\odot}/m_{1}$ to 1. We
find that the percentage of binaries hosting a BH progenitor and a secondary
star less massive than $0.71\ \mathrm{M}_{\odot}$ is 40% of similar binaries
with a secondary in the mass range $0.71-1.0\ \mathrm{M}_{\odot}$. Therefore,
we expect an overall increase of the formation efficiency of about 40% with
respect to the value reported in Table 1. This factor is balanced by the
commonly assumed binary fraction in stellar populations ($40-60\%$, see e.g.
Moe & Di Stefano 2017).

Although the formation of Gaia BH3 does not depend on binary processes such as
common envelope and mass transfer, the uncertainties in the stellar evolution
models can have an important impact on this study. For instance, differences
in the maximum radius of the BH progenitor are fundamental to select what
objects interact during binary evolution. Different wind prescriptions are
fundamental to set the final relation between the initial mass and the BH
mass. El-Badry (2024) explored the main uncertainties about stellar evolution
models showing that systems such as Gaia BH3 can be produced within a
reasonable range of stellar parameters including initial mass, metallicity,
stellar rotation, and convection models.

Many population synthesis studies on BH-star binaries are based on stellar
evolution models by Pols et al. (1998) (see e.g. Breivik et al. 2017; Shao &
Li 2019; Chawla et al. 2022; Shikauchi et al. 2023). In general, the stars in
these tracks tend to expand to larger radii compared to the mist models. Most
of these works focus on the Galactic population of which the stellar halo is a
negligible component, hence a direct comparison is not straightforward. We
notice that the properties of the few metal-poor stars reported by Chawla et
al. (2022) seem to be qualitatively consistent with the distribution of BH-
star systems in terms of period and BH mass.

###  4.3 Alternative formation channel

Gaia BH3 is associated with the ED-2 stream (Balbinot et al. 2024), therefore
it possibly spent a part of its life in a dynamical environment. The variety
of interactions in dynamical environments have two competitive effects. On the
one hand, they can pair BHs and stars or create binary systems with orbital
properties (e.g., high eccentricities) that boost the production of Gaia
BH3-like binaries. On the other hand, dynamical interactions could destroy
original and newly formed BH-star systems, especially if their orbit is wide
and the companion star is low mass.

The other two BH-star systems retrieved from Gaia data (BH1 El-Badry et al.
2023b, and BH2 El-Badry et al. 2023a) are difficult to reconcile with isolated
binary evolution models (El-Badry et al. 2023a), while they are more easily
formed in star clusters (Rastello et al. 2023; Di Carlo et al. 2024; Tanikawa
et al. 2024; Kotko et al. 2024).

Marín Pina et al. (2024) show that BH-star systems dynamically form in an
efficient way in globular clusters with mass $10^{5}-10^{6}\
\mathrm{M}_{\odot}$. They estimate a formation efficiency of $10^{-5}-3\times
10^{-4}$ M${}^{-1}_{\odot}$ by considering all the BH-star systems formed in
their simulations (both escaped or retained in the cluster). Their proposed
dynamical channel is one to two orders of magnitude times more efficient than
the isolated channel presented in this work. However, when considering the
halo (Deason et al. 2019) and the cluster radial profile (e.g., Gieles &
Gnedin 2023), the difference on the expected formation efficiency of BH-star
systems in the solar neighbourhood reduces to a factor of 3–10. If Gaia BH3
remains the only observed BH-star system in the halo, none of the two
formation channels can be excluded or significantly favoured.

When focusing specifically on Gaia BH3-like systems, Marín Pina et al. (2024)
find that only 0.6% of their BH-star systems fall within the same selection
box used in our work (Figure 2). This fraction is a factor of 3–5 smaller than
what we found in our isolated-binary simulations. Taken at face value, these
numbers seem to indicate that the dynamic environment does not considerably
boost the formation of Gaia BH3-like systems. Therefore, even if it was
clearly born in a cluster, Gaia BH3 is compatible with a primordial binary
star which escaped from the cluster without experiencing significant dynamical
interactions.

We point out that our models can easily reproduce the physical properties of
Gaia BH3 with low or negligible BH natal kicks. Thus, they also point towards
an efficient retention of the so-formed binary within the star cluster where
it was born. As a consequence, our isolated formation scenario can account for
the formation of Gaia BH3 in a low-density and dynamically-quite star cluster,
while being consistent with its membership in a stellar stream like ED-2
(Balbinot et al. 2024).

A more quantitative comparison and model selection requires a detailed
analysis of systematic, selection effects, channel definitions, halo
population models, and robust statistic tool (see e.g. Mould et al. 2023) that
are beyond the scope of this paper.

##  5 Summary

We have explored the formation of Gaia BH3 (Gaia Collaboration et al. 2024), a
non interacting binary system composed of a massive BH ($\approx 30\
\mathrm{M}_{\odot}$) and a low-mass giant star ($\approx 0.8\
\mathrm{M}_{\odot}$), by means of population synthesis simulations.

We have exploited the flexibility of our code sevn (Iorio et al. 2023) to
perform the first population study of BH-star systems using the mist stellar
models (Choi et al. 2016). We explored the parameter space by considering 9
different metallicities, 4 different supernova kick models, 2 models for the
common envelope efficiency parameter $\alpha$, and 2 models for the onset of
Roche-Lobe overflow. Also, we implemented a Bayesian sampling analysis to
study the distribution of initial conditions that reproduce exactly the
observed properties of Gaia BH3.

The main takeaways of our analysis are as follows.

  * •

Systems like Gaia BH3 represent a common product of isolated binary evolution
at low metallicity ($Z<0.01$). They form preferentially from a massive star
(40–60 M⊙) in an initially wide ($P>10^{3}$ days) and eccentric ($e>0.6$)
orbit with a low mass companion. The progenitor binary system never underwent
Roche-lobe overflow during its life.

  * •

Our fiducial supernova kick model (Giacobbo & Mapelli 2020) produces low natal
kicks, leaving the initial orbital properties of the binary systems almost
unchanged. The formation of Gaia BH3 is possible up to kicks as high as 220
km/s, but this requires fine tuning for the conditions at BH formation, and
the overall formation efficiency drops by about one order of magnitude.

  * •

Focusing on systems older than 10 Gyr, we find a formation efficiency
$\eta_{\rm BHStar}\sim 10^{-6}$ M${}_{\odot}^{-1}$ for BH-star binaries in
general, and $\eta_{\rm GaiaBH3-like}\sim 5\times 10^{-8}$
M${}_{\odot}^{-1}$for Gaia BH3-like systems.

  * •

Considering the Galactic stellar halo model by Deason et al. (2019), we expect
up to $\approx$ 3000 BH-star systems born in the stellar halo from isolated
evolution. About 50–100 of them are Gaia BH3-like systems in terms of mass and
orbital properties. Given these numbers, the hypothesis that Gaia BH3 formed
through isolated evolution cannot be ruled out at 3$\sigma$ significance
level.

  * •

Comparing our formation efficiencies with the ones derived by Marín Pina et
al. (2024) for the dynamical formation scenario, we conclude that none of the
two channels is strongly favoured for the formation of Gaia BH3.

Overall, Gaia BH3 is compatible with a primordial, metal-poor binary system
that never evolved via Roche-lobe overflow. If it formed in a star cluster, it
might have evolved nearly unperturbed until it quietly evaporated from it or
was dynamically ejected. Unlike Gaia BH1 and BH2, its orbital properties are
easy to explain with current stellar and binary evolution models.

###### Acknowledgements.

We thank Daniel Marín Pina for useful comments. Part of this work was
supported by the German _Deutsche Forschungsgemeinschaft, DFG_ project number
Ts 17/2–1. MM, GI, ST, MD, StR, GE, EK, EL and MPV acknowledge financial
support from the European Research Council for the ERC Consolidator grant
DEMOBLACK, under contract no. 770017. MM, ST, StR and MPV also acknowledge
financial support from the German Excellence Strategy via the Heidelberg
Cluster of Excellence (EXC 2181 - 390900948) STRUCTURES. MD also acknowledges
financial support from Cariparo foundation under grant 55440. AAT acknowledges
support from the European Union’s Horizon 2020 and Horizon Europe research and
innovation programs under the Marie Skłodowska-Curie grant agreements no.
847523 and 101103134. SaR acknowledges support from the Beatriu de Pinós
postdoctoral program under the Ministry of Research and Universities of the
Government of Catalonia (Grant Reference No. 2021 BP 00213). MAS acknowledges
funding from the European Union’s Horizon 2020 research and innovation
programme under the Marie Skłodowska-Curie grant agreement No. 101025436
(project GRACE-BH, PI: Manuel Arca Sedda). MAS acknowledge financial support
from the MERAC Foundation.

## References

  * Abbott et al. (2016a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, Phys. Rev. Lett., 116, 061102 [arXiv:1602.03837]
  * Abbott et al. (2016b) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016b, Phys. Rev. Lett., 116, 241102 [arXiv:1602.03840]
  * Abbott et al. (2021a) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021a, Physical Review X, 11, 021053 [arXiv:2010.14527]
  * Abbott et al. (2021b) Abbott, R., Abbott, T. D., Acernese, F., et al. 2021b, arXiv e-prints, arXiv:2111.03606 [arXiv:2111.03606]
  * Abbott et al. (2021c) Abbott, R., Abbott, T. D., Acernese, F., et al. 2021c, arXiv e-prints, arXiv:2108.01045 [arXiv:2108.01045]
  * Andrews & Kalogera (2022) Andrews, J. J. & Kalogera, V. 2022, ApJ, 930, 159 [arXiv:2203.15156]
  * Atri et al. (2019) Atri, P., Miller-Jones, J. C. A., Bahramian, A., et al. 2019, MNRAS, 489, 3116 [arXiv:1908.07199]
  * Balbinot et al. (2024) Balbinot, E., Dodd, E., Matsuno, T., et al. 2024, arXiv e-prints, arXiv:2404.11604 [arXiv:2404.11604]
  * Belczynski et al. (2010) Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010, ApJ, 714, 1217 [arXiv:0904.2784]
  * Belczynski et al. (2016) Belczynski, K., Heger, A., Gladysz, W., et al. 2016, A&A, 594, A97 [arXiv:1607.03116]
  * Bell et al. (2008) Bell, E. F., Zucker, D. B., Belokurov, V., et al. 2008, ApJ, 680, 295 [arXiv:0706.0004]
  * Breivik et al. (2017) Breivik, K., Chatterjee, S., & Larson, S. L. 2017, ApJ, 850, L13 [arXiv:1710.04657]
  * Casares et al. (2014) Casares, J., Negueruela, I., Ribó, M., et al. 2014, Nature, 505, 378 [arXiv:1401.3711]
  * Chawla et al. (2022) Chawla, C., Chatterjee, S., Breivik, K., et al. 2022, ApJ, 931, 107 [arXiv:2110.05979]
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102 [arXiv:1604.08592]
  * Claeys et al. (2014) Claeys, J. S. W., Pols, O. R., Izzard, R. G., Vink, J., & Verbunt, F. W. M. 2014, A&A, 563, A83 [arXiv:1401.2895]
  * Costa et al. (2021) Costa, G., Bressan, A., Mapelli, M., et al. 2021, MNRAS, 501, 4514 [arXiv:2010.02242]
  * Costa et al. (2023) Costa, G., Mapelli, M., Iorio, G., et al. 2023, arXiv e-prints, arXiv:2303.15511 [arXiv:2303.15511]
  * Deason et al. (2019) Deason, A. J., Belokurov, V., & Sanders, J. L. 2019, MNRAS, 490, 3426 [arXiv:1908.02763]
  * Di Carlo et al. (2024) Di Carlo, U. N., Agrawal, P., Rodriguez, C. L., & Breivik, K. 2024, ApJ, 965, 22 [arXiv:2306.13121]
  * Eggleton (1983) Eggleton, P. P. 1983, ApJ, 268, 368
  * El-Badry (2024) El-Badry, K. 2024, arXiv e-prints, arXiv:2404.13047 [arXiv:2404.13047]
  * El-Badry et al. (2023a) El-Badry, K., Rix, H.-W., Cendes, Y., et al. 2023a, MNRAS, 521, 4323
  * El-Badry et al. (2023b) El-Badry, K., Rix, H.-W., Quataert, E., et al. 2023b, MNRAS, 518, 1057 [arXiv:2209.06833]
  * Farmer et al. (2019) Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, ApJ, 887, 53 [arXiv:1910.12874]
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [arXiv:1202.3665]
  * Fryer et al. (2012) Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91 [arXiv:1110.1726]
  * Gaia Collaboration et al. (2024) Gaia Collaboration, Panuzzo, P., Mazeh, T., et al. 2024, arXiv e-prints, arXiv:2404.10486 [arXiv:2404.10486]
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1 [arXiv:2208.00211]
  * Generozov & Perets (2024) Generozov, A. & Perets, H. B. 2024, ApJ, 964, 83 [arXiv:2312.03066]
  * Giacobbo & Mapelli (2020) Giacobbo, N. & Mapelli, M. 2020, ApJ, 891, 141 [arXiv:1909.06385]
  * Gieles & Gnedin (2023) Gieles, M. & Gnedin, O. Y. 2023, MNRAS, 522, 5340 [arXiv:2303.03791]
  * Giesers et al. (2018) Giesers, B., Dreizler, S., Husser, T.-O., et al. 2018, MNRAS, 475, L15 [arXiv:1801.05642]
  * Giesers et al. (2019) Giesers, B., Kamann, S., Dreizler, S., et al. 2019, A&A, 632, A3 [arXiv:1909.04050]
  * Guseinov & Zel’dovich (1966) Guseinov, O. K. & Zel’dovich, Y. B. 1966Sov. Ast., 10, 251
  * Heger & Woosley (2002) Heger, A. & Woosley, S. E. 2002, ApJ, 567, 532 [arXiv:astro-ph/0107037]
  * Hobbs et al. (2005) Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974 [arXiv:astro-ph/0504584]
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [arXiv:astro-ph/0201220]
  * Hut (1981) Hut, P. 1981A&A, 99, 126
  * Iorio et al. (2018) Iorio, G., Belokurov, V., Erkal, D., et al. 2018, MNRAS, 474, 2142 [arXiv:1707.03833]
  * Iorio et al. (2023) Iorio, G., Mapelli, M., Costa, G., et al. 2023, MNRAS, 524, 426 [arXiv:2211.11774]
  * Ivanova et al. (2013) Ivanova, N., Justham, S., Chen, X., et al. 2013, A&A Rev., 21, 59 [arXiv:1209.4302]
  * Kotko et al. (2024) Kotko, I., Banerjee, S., & Belczynski, K. 2024, arXiv e-prints, arXiv:2403.13579 [arXiv:2403.13579]
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231 [arXiv:astro-ph/0009005]
  * Lane et al. (2023) Lane, J. M. M., Bovy, J., & Mackereth, J. T. 2023, MNRAS, 526, 1209 [arXiv:2306.03084]
  * Lennon et al. (2022) Lennon, D. J., Dufton, P. L., Villaseñor, J. I., et al. 2022, A&A, 665, A180 [arXiv:2111.12173]
  * Mahy et al. (2022) Mahy, L., Sana, H., Shenar, T., et al. 2022, A&A, 664, A159 [arXiv:2207.07752]
  * Mapelli et al. (2009) Mapelli, M., Colpi, M., & Zampieri, L. 2009, MNRAS, 395, L71 [arXiv:0902.3540]
  * Mapelli et al. (2010) Mapelli, M., Ripamonti, E., Zampieri, L., Colpi, M., & Bressan, A. 2010, MNRAS, 408, 234 [arXiv:1005.3548]
  * Mapelli et al. (2020) Mapelli, M., Spera, M., Montanari, E., et al. 2020, ApJ, 888, 76 [arXiv:1909.01371]
  * Mapelli et al. (2013) Mapelli, M., Zampieri, L., Ripamonti, E., & Bressan, A. 2013, MNRAS, 429, 2298 [arXiv:1211.6441]
  * Marín Pina et al. (2024) Marín Pina, D., Rastello, S., Gieles, M., et al. 2024arXiv e-prints, arXiv:2404.13036 [arXiv:2404.13036]
  * Medina et al. (2024) Medina, G. E., Muñoz, R. R., Carlin, J. L., et al. 2024, arXiv e-prints, arXiv:2402.14055 [arXiv:2402.14055]
  * Miller-Jones et al. (2021) Miller-Jones, J. C. A., Bahramian, A., Orosz, J. A., et al. 2021, Science, 371, 1046 [arXiv:2102.09091]
  * Moe & Di Stefano (2017) Moe, M. & Di Stefano, R. 2017, ApJS, 230, 15 [arXiv:1606.05347]
  * Mould et al. (2023) Mould, M., Gerosa, D., Dall’Amico, M., & Mapelli, M. 2023, MNRAS, 525, 3986 [arXiv:2305.18539]
  * Oëzel et al. (2010) Oëzel, F., Psaltis, D., Narayan, R., & McClintock, J. E. 2010, ApJ, 725, 1918 [arXiv:1006.2834]
  * Penoyre et al. (2022) Penoyre, Z., Belokurov, V., & Evans, N. W. 2022, MNRAS, 513, 2437 [arXiv:2111.10380]
  * Pols et al. (1998) Pols, O. R., Schröder, K.-P., Hurley, J. R., Tout, C. A., & Eggleton, P. P. 1998, MNRAS, 298, 525
  * Rastello et al. (2023) Rastello, S., Iorio, G., Mapelli, M., et al. 2023, MNRAS, 526, 740 [arXiv:2306.14679]
  * Ribó et al. (2017) Ribó, M., Munar-Adrover, P., Paredes, J. M., et al. 2017, ApJ, 835, L33 [arXiv:1701.07265]
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [arXiv:1207.6397]
  * Saracino et al. (2022) Saracino, S., Kamann, S., Guarcello, M. G., et al. 2022, MNRAS, 511, 2914 [arXiv:2111.06506]
  * Saracino et al. (2023) Saracino, S., Shenar, T., Kamann, S., et al. 2023, MNRAS, 521, 3162 [arXiv:2303.07369]
  * Sgalletta et al. (2023) Sgalletta, C., Iorio, G., Mapelli, M., et al. 2023, arXiv e-prints, arXiv:2305.04955 [arXiv:2305.04955]
  * Shao & Li (2019) Shao, Y. & Li, X.-D. 2019, ApJ, 885, 151 [arXiv:1909.11328]
  * Shenar et al. (2022a) Shenar, T., Sana, H., Mahy, L., et al. 2022a, Nature Astronomy, 6, 1085 [arXiv:2207.07675]
  * Shenar et al. (2022b) Shenar, T., Sana, H., Mahy, L., et al. 2022b, A&A, 665, A148 [arXiv:2207.07674]
  * Shikauchi et al. (2023) Shikauchi, M., Tsuna, D., Tanikawa, A., & Kawanaka, N. 2023, ApJ, 953, 52 [arXiv:2301.07207]
  * Spera & Mapelli (2017) Spera, M. & Mapelli, M. 2017, MNRAS, 470, 4739 [arXiv:1706.06109]
  * Spera et al. (2015) Spera, M., Mapelli, M., & Bressan, A. 2015, MNRAS, 451, 4086 [arXiv:1505.05201]
  * Spera et al. (2019) Spera, M., Mapelli, M., Giacobbo, N., et al. 2019, MNRAS, 485, 889 [arXiv:1809.04605]
  * Tanikawa et al. (2024) Tanikawa, A., Cary, S., Shikauchi, M., Wang, L., & Fujii, M. S. 2024, MNRAS, 527, 4031 [arXiv:2303.05743]
  * Tanikawa et al. (2023) Tanikawa, A., Hattori, K., Kawanaka, N., et al. 2023, ApJ, 946, 79 [arXiv:2209.05632]
  * Tanikawa et al. (2021) Tanikawa, A., Susa, H., Yoshida, T., Trani, A. A., & Kinugawa, T. 2021, ApJ, 910, 30 [arXiv:2008.01890]
  * Thompson et al. (2019) Thompson, T. A., Kochanek, C. S., Stanek, K. Z., et al. 2019, Science, 366, 637 [arXiv:1806.02751]
  * Trimble & Thorne (1969) Trimble, V. L. & Thorne, K. S. 1969, ApJ, 156, 1013
  * Webbink (1984) Webbink, R. F. 1984, ApJ, 277, 355
  * Woosley (2017) Woosley, S. E. 2017, ApJ, 836, 244 [arXiv:1608.08939]
  * Woosley (2019) Woosley, S. E. 2019, ApJ, 878, 49 [arXiv:1901.00215]
  * Zahn (1977) Zahn, J. P. 1977A&A, 57, 383
  * Zampieri & Roberts (2009) Zampieri, L. & Roberts, T. P. 2009, MNRAS, 400, 677 [arXiv:0909.1017]
  * Ziosi et al. (2014) Ziosi, B. M., Mapelli, M., Branchesi, M., & Tormen, G. 2014, MNRAS, 441, 3703 [arXiv:1404.7147]

  *[BH]: black hole