# The enigmatic origin of two dormant BH binaries: Gaia BH1 and Gaia BH2.

Kotko,I.,1 Banerjee,S.,2,3 Belczynski,K.,1
1Nicolaus Copernicus Astronomical Centre of the Polish Academy of Sciences,
ul. Bartycka 18, 00-716 Warszawa, Poland
2Helmholtz-Instituts für Strahlen- und Kernphysik (HISKP),
Nussallee 14-16, D-53115 Bonn, Germany
3Argelander-Institut für Astronomie (AIfA), Auf dem Hügel 71, D-53121, Bonn,
Germany
E-mail<EMAIL_ADDRESS>
###### Abstract

The two systems that have been confirmed as dormant (no X-ray emission
detected) black hole (BH) - low-mass star binaries (Gaia BH1 and Gaia BH2) in
the latest Gaia mission data release (DR3) are very intriguing in the context
of their evolution. Both systems consist of $\sim 9\,\mathrm{M_{\odot}}$ BH
and $\sim 1\,\mathrm{M_{\odot}}$ star orbiting each other on a wide, eccentric
orbit ($e\sim 0.5$). We argue that the formation of Gaia BH-like systems
through the isolated binary evolution (IBE) channel under standard common
envelope assumptions and from dynamical interactions in young massive and open
clusters is equally probable, and the formation rate of those binaries is on
order of $10^{-7}\,\mathrm{M_{\odot}}^{-1}$ in both channels.
What plays an important role in formation of Gaia BH-like systems in the case
of IBE channel is the mutual position of the natal kick velocity vector and
the binary angular momentum vector. We find that the natal kicks of median
magnitude $\sim 40$ km/s are preferred for the formation of Gaia BH1-like
binaries. $\sim 94\%$ of those binaries is formed with BH spin misaligned to
orbital axis by less than $40^{\circ}$. Gaia BH2-like binaries form if the low
velocity natal kick (of median magnitude $\sim 20$ km/s) is directed within
$15^{\circ}$ about the orbital plane.
We follow the subsequent evolution of the binaries once formed as Gaia BH1 and
Gaia BH2 systems to investigate their connection with the low-mass X-ray
binary population.

###### keywords:

binaries: general – stars: black holes – stars: evolution – methods: numerical
– open clusters and associations: general – stars: kinematics and dynamics

††pagerange: The enigmatic origin of two dormant BH binaries: Gaia BH1 and
Gaia BH2.–B

##  1 Introduction

The dormant black hole (BH) binaries are difficult to observe as they do not
emit X-ray radiation and can be discovered only by astrometric and/or
spectroscopic analysis. Although such systems were predicted to exist already
in the 1960s (Guseinov & Zel’dovich, 1966) only recently the new generation of
the astronomical instruments and new data analysis methods have allowed to
identify the first candidates for such objects. To our knowledge the first
strong candidate for a binary system with a heavy companion unseen in
electromagnetic spectrum was found in the globular cluster NGC 3201 Giesers et
al. (2018). Within last two years another four systems have been discovered:
(i) binary HD 130298 where the minimum mass of an unseen companion has been
estimated to be $7.7\pm 1.5$ $\mathrm{M_{\odot}}$ (Mahy et al., 2022), (ii)
VFTS 243 found in Large Magellanic Cloud which is a binary of massive star and
most probably black hole companion with mass $M=10.1\pm 2.0$
$\mathrm{M_{\odot}}$ (Shenar et al., 2022) and (iii) two systems which were
identified in the most recent Gaia catalogue DR3 for which there is the strong
evidence that they host a non-accreting BH and a low-mass companion star (El-
Badry et al., 2023a, b; Tanikawa et al., 2023b). In our studies we focus only
on the last two systems and following El-Badry et al. (2023b) we refer to
those binaries as Gaia BH1 and Gaia BH2.
Once Gaia BH1 and Gaia BH2 had been detected, the question arose about the
origin of those systems. Already several studies have been made which point
towards evolution through the dynamical interactions in young stellar clusters
as the most probable evolution scenario to create BH-low-mass star system on a
moderately wide and eccentric orbit (El-Badry et al., 2023a; Tanikawa et al.,
2023a; Rastello et al., 2023; Niccolò Di Carlo et al., 2023). However, as we
want to point out in this work, the isolated binary evolution (IBE) channel
may contribute to the total population of dormant BH binaries at a rate
comparable to the dynamical interactions channel. The important factor that
impacts the number of the binaries of interest that may be born through IBE
channel is the magnitude and geometry of the natal kick imparted to the newly
formed compact object after supernova (SN) explosion. One can infer the
information about the natal kick from the peculiar velocity111The peculiar
velocity corresponds to systemic velocity in the model, as is explained in
Sec. 2.1. of the system and the system orbital parameters.
There are two components that change the systemic velocity of a binary after
SN event in the model: (i) the Blaauw kick that impacts the velocity of the
binary centre of gravity due to symmetric mass ejection during SN (Blaauw,
1961), (ii) the natal kick that is received by the compact object due to
asymmetries of mass ejecta and/or neutrino emission during SN (Chugai, 1984;
Janka, 2013, 2017, e.g.). The observations of proper motions of pulsars and
X-ray binaries give constrains on the natal kicks received by neutron stars
(NSs) or BHs upon their formation. The Maxwellian distribution with dispersion
velocity $\sigma=265$ km/s derived by Hobbs et al. (2005) from observations of
pulsars is widely used in the binary population synthesis simulations. This
distribution is applied also to BHs but with additional assumption that the
kick velocity is reduced by the velocity of the matter that falls back pulled
by the gravitation field of BH forming during SN. However, recent observations
indicate that this approach may underestimate the number of lower velocity
systems and overestimate the number of the systems in higher velocity range
and that it should be revisited (Atri et al., 2019; O’Doherty et al., 2023;
Zhao et al., 2023). The magnitude of the natal kicks is only one part of a
puzzle, the other is the kick velocity vector orientation. In the codes based
on BSE code from Hurley et al. (2000) the direction of natal kicks is drawn
from isotropic spherical distribution. The question is whether the specific
direction of natal kicks that BHs receive in SN explosion could lead to the
formation of Gaia BH1 and Gaia BH2 (under the standard assumptions of isolated
binary evolution model). The reason we may attempt to answer this question in
the case of the two systems under considerations is that Gaia BH1 and Gaia BH2
are not impacted by the orbit changes due to the mass transfer events after
SN. Gaia BH1 has the orbital separation and eccentricity as it was shaped by
SN and Gaia BH2 orbit is only modulated by weak tidal forces.
Our work consists of two main parts: in Sec. 2 we consider IBE channel and
follow the evolution of the binaries which we classify as Gaia BH1-like and
Gaia BH2-like to find the probability of formation of such binaries from IBE,
we then run the evolution of a two well-matching Gaia BH1 and Gaia BH2
binaries till the Hubble time is reached and finally we investigate the
information that Gaia BH1 and Gaia BH2 may bring us about the natal kicks; in
Sec. 3.1 we present and discuss the formation rates of Gaia BH-like binaries
from young massive and open star cluster.

Table 1: The parameters of Gaia BH1 and Gaia BH2: the orbital period
$P_{\mathrm{orb}}$ in days, the eccentricity $e$, BH mass $M_{\mathrm{BH}}$
and companion star mass $M_{\star}$ in solar units, companion evolutionary
stage and age in Gyrs, adopted from El-Badry et al. (2023a, b). The peculiar
velocity at birth $v_{\mathrm{pec}}$ is taken from Table 3 in Zhao et al.
(2023).

| Gaia BH1 | Gaia BH2
---|---|---
$P_{\mathrm{orb}}$ [d] | $185.59\pm 0.05$ | $1276.7\pm 0.6$
$e$ | $0.451\pm 0.005$ | $0.5176\pm 0.0009$
$M_{\mathrm{BH}}$ [$\mathrm{M_{\odot}}$] | $9.62\pm 0.18$ | $8.94\pm 0.34$
$M_{\star}$ [$\mathrm{M_{\odot}}$] | $0.93\pm 0.05$ | $1.07\pm 0.19$
Star type | MS | RG
Age [Gyr] | $>4.0$ | $8.0-12.0$
$v_{\mathrm{pec}}\,[\mathrm{km/s}]$ | $71.3^{+9.2}_{-10.8}$ | $34.1^{+5.0}_{-5.1}$

##  2 Gaia BH binaries from isolated binary evolution

###  2.1 Method

We performed the population synthesis calculations using StarTrack. The
detailed description of the code can be found in Belczynski et al. (2002);
Belczynski et al. (2008, 2016). Here we highlight the assumptions that are
most important in the context of the present paper.
The adopted initial mass distribution is a power law function with three
exponents Kroupa et al. (1993); Kroupa (2001): $\alpha_{1}=-1.3$ for stars of
masses $0.08<M<0.5\,\mathrm{M_{\odot}}$, $\alpha_{2}=-2.2$ for $0.5\leq
M<1.0\,\mathrm{M_{\odot}}$ and $\alpha_{3}=-2.3$ for stars more massive than
$1.0\,\mathrm{M_{\odot}}$. To maximise the effectiveness of calculations we
focused on the systems that may be the progenitors of Gaia BH binaries: the
mass of BH-progenitor star was drawn from the mass range
$18.0-200.0\,\mathrm{M_{\odot}}$ and the mass of its companion from the range
$0.2-1.5\,\mathrm{M_{\odot}}$. The orbital period ($P/\mathrm{orb}$) (in units
of days) and eccentricity $e$ were drawn from the distributions:
$f(\log{P/\mathrm{orb}})=(\log{P/\mathrm{orb}})^{-0.55}$ with
$\log{P/\mathrm{orb}}$ in the range $[0.15,4.0]$ and $f(e)=e^{-0.42}$ with $e$
in the range $[0.0,0.9]$ respectively (Sana et al., 2012, 2013).
In our calculations the common envelope phase is parametrized by CE parameter
$\alpha_{\mathrm{CE}}=\alpha_{\mathrm{ce}}\times\lambda$ where $\lambda$ is a
parameter describing the donor central concentration (de Kool, 1989) and
$\alpha_{\mathrm{ce}}$ is parametrizing the efficiency at which the orbital
energy is transferred into the envelope (Webbink, 1984). We assume standard
value $\alpha_{\mathrm{ce}}=1.0$ and for $\lambda$ we use the fits from Xu &
Li (2010a); Xu & Li (2010b) and Wang et al. (2016).
The massive stars (of type O/B) stellar winds were described by the equations
for radiation driven mass loss from (Vink et al., 2001) which take into
account the Luminous Blue Variable mass loss (Belczynski et al., 2010). For
Wolf-Rayet stars we took into account the wind clumping (Hamann & Koesterke,
1998) and the metallicity dependence of the winds (Vink & de Koter, 2005). The
winds of less massive stars were calculated according to the formulae of
Hurley et al. (2000).
We assumed the delayed supernova engine (Fryer et al., 2012) in our
calculations. There are two components that make up the final kick received by
the system in SN explosion: the kick from asymmetric neutrino and/or ejecta
emission (natal kick) and the Blaauw kick resulting from the symmetric
envelope ejection from the exploding star. We tested two distributions of the
natal kick velocities: i) the widely used Maxwellian distribution with
$\sigma=265\,\mathrm{km\,s^{-1}}$ (Hobbs et al., 2005) which for BH is
decreased by the fallback of the fraction of the stellar envelope that is
initially ejected in SN explosion, this fraction is defined by the fallback
parameter $f_{\mathrm{fb}}$ which takes the values between $0$ and $1$ (model
V1), ii) the two-component Maxwellian with
$\sigma_{1}=21.3\,\mathrm{km\,s^{-1}}$ for low-velocity component and
$\sigma_{2}=106.7\,\mathrm{km\,s^{-1}}$ for high-velocity component derived
from the sample of 85 BH and NS X-ray binaries by Zhao et al. (2023) also
decreased by the fallback for BHs (model V2). The directions of the kicks are
drawn from isotropic distribution. For each system the drawing of natal kick
magnitude and orientation was repeated $200$ times.
The systemic velocity $v_{\mathrm{sys}}$ calculated in StarTrack is the
velocity of the binary centre of mass resulting from Blaauw kick and natal
kick. It is assumed that the center of mass is initially at rest.
$v_{\mathrm{sys}}$ is consistent with the potential peculiar velocity at birth
$v_{\mathrm{p,z=0}}$ (Atri et al., 2019), that is the velocity that the binary
gained from SN which caused its motion out of the Galactic plane. The
potential peculiar velocity is inferred from the observational present-day
peculiar velocity (for details see Atri et al. (2019); Zhao et al. (2023)).
We position the orbital plane of a binary before the first SN on x-y plane.
All three vectors, the orbital angular momentum and the spins of a star and
BH, are aligned and pointing in $z$ direction. Setting binary in such position
does not influence our results while it makes the analysis of the post-SN
relative orientation of those vectors much easier. The position of the
exploding star on the orbit is found by means of mean anomaly $M$, eccentric
anomaly $E$, eccentricity $e$ and the time of periapsis passage $\tau$.
Taking into account the age constrains derived from the observations of the
luminous component in Gaia BH1 and Gaia BH2 we assumed that both binaries are
$\sim 4$ Gyr old (El-Badry et al., 2023a, b) and, therefore, were formed in
the environment of metallicity $Z=0.7Z_{\odot}$ (see Fig.3 in Olejak et al.,
2020) where$Z_{\odot}=0.014$ (Asplund et al., 2009).
It is hard to expect that one may get many systems of the parameters exactly
corresponding to those inferred from observations of Gaia BH1 and Gaia BH2 in
the simulations. Therefore, by putting constrains on the orbital separation
$a$, the eccentricity $e$, the BH mass $M_{\mathrm{BH}}$, the companion star
mass $M_{2}$ and the companion star evolutionary stage, we defined three
groups of binaries which are of our interest: Gaia BH1-like, Gaia BH2-like and
Gaia BH-like. Gaia BH1- and Gaia BH2-like are the systems assumed to
correspond to observed Gaia BH1 and Gaia BH2 while Gaia BH-like is the general
group of all systems which may be classified as dormant BH binaries. To
investigate the interconnection between the natal kicks and the post-SN
quantities: $a$, $e$, $v_{\mathrm{sys}}$, orbital angular momentum and BH spin
vectors orientation we defined the fourth group - the "optimistic" Gaia BH-
like binaries - for which we allow the wider range of $a$, $e$ and both
components masses. The classification criteria for each group are given in
Table 2.
We define two angles which help to complement our analysis: (i) $\gamma$ \- a
misalignment between the binary angular momentum vector and the natal kick
velocity vector, and (ii) $\theta$ \- a misalignment between the binary
angular momentum vector and post-SN BH spin vector.

Table 2: The parameters of the binaries that we classify as Gaia BH1-, Gaia
BH2-, and Gaia BH-like in our standard models V1 and V2. We take into account
both components masses ($M_{\mathrm{BH}}$ and $M_{2}$), the evolutionary stage
of the companion (main sequence (MS) /red giant (RG)), the orbital separation
$a$, the eccentricity $e$ and the systemic velocity $v_{\mathrm{sys}}$ (km/s).

|  |  | standard | optimistic
---|---|---|---|---
| Gaia BH1-like | Gaia BH2-like | Gaia BH-like | Gaia BH-like
$M_{\mathrm{BH}}$ [$\mathrm{M_{\odot}}$] | 8.0-12.0 | 7.0-11.0 | <15.0 | < 15.0
$M_{2}$ [$\mathrm{M_{\odot}}$] | 0.8-1.2 | 0.8-1.2 | 0.8-1.2 | < 1.5
Companion star | MS | RG | MS/RG | MS/RG
$a$ [$\mathrm{R_{\odot}}$] | 200-400 | 950-1150 | 200-1200 | 100-1500
$e$ | 0.4-0.6 | 0.4-0.6 | 0.3-0.7 | 0.1-0.9
$v_{\mathrm{sys}}\,[\mathrm{km/s}]$ | no limit | no limit | no limit | no limit

###  2.2 Results

####  2.2.1 The formation rate of Gaia BH-like systems from IBE

The formation rate of Gaia BH-like binaries depends on the criteria according
to which one classifies the binary as Gaia BH-like and on the adopted natal
kick velocity distribution and it may differ as much as three orders of
magnitude as can be seen from Table 3. The formation rate for model V1 is an
order of magnitude smaller than formation rates for model V2. There is also
one order of magnitude difference between the standard and optimistic groups
in both models: for optimistic and standard Gaia BH-like groups in model V2
the formation rate is $\sim 1.4\times 10^{-7}\,\mathrm{M_{\odot}}^{-1}$ and
$\sim 1.0\times 10^{-8}\,\mathrm{M_{\odot}}^{-1}$ respectively and in model V1
it is $\sim 1.6\times 10^{-9}\,\mathrm{M_{\odot}}^{-1}$ and $\sim 1.1\times
10^{-8}\,\mathrm{M_{\odot}}^{-1}$ respectively.
The formation rates of Gaia BH1- and Gaia BH2-like systems, which are the two
binary categories with the most stringent classification criteria among the
binary types defined in Table 2, are on order of
$10^{-10}\,\mathrm{M_{\odot}}^{-1}$. The formation rate of Gaia BH2-like
systems is almost an order of magnitude higher in model V2 than in model V1,
while the formation rate of Gaia BH1-like binaries does not depend on the
natal kick distribution adopted (model V1 vs V2).
Additionally, to make the comparison with the results of dynamical channel
described in Sec. 3.2, we calculated the general population of BH-star
binaries with the only constrain imposed on the companion star mass ($M_{\rm
comp}<3.0\,\mathrm{M_{\odot}}$) and we found that the formation rate of such
binaries from IBE channel is $2.2\times 10^{-6}\,\mathrm{M_{\odot}}^{-1}$.

Table 3: The formation rates of Gaia BH-like binaries for different classification criteria and natal kick velocity distributions V1 and V2. The criteria adopted to calculated the formation rate for case "V2 low-mass comp" correspond to the population of BH-star binaries with low-mass companion ($M_{\rm comp}<3\,\mathrm{M_{\odot}}$) from young massive and open clusters described in Sec. 3.2. We present the formation rates of Gaia BH1- and Gaia BH2-like binaries in model V1 and V2 in the last three rows. Case | Rate [$\mathrm{M_{\odot}}^{-1}$]
---|---
V1 standard | $\sim 1.6\times 10^{-9}$
V1 optimistic | $\sim 1.1\times 10^{-8}$
V2 standard | $\sim 1.0\times 10^{-8}$
V2 optimistic | $\sim 1.4\times 10^{-7}$
V2 low-mass comp | $\sim 2.2\times 10^{-6}$
V1 standard | $\sim 1.6\times 10^{-9}$
V1 $\&$ V2 Gaia BH1 | $\sim 3.2\times 10^{-10}$
V1 Gaia BH2 | $\sim 1.4\times 10^{-10}$
V2 Gaia BH2 | $\sim 9.6\times 10^{-10}$

Among Gaia BH1-like an Gaia BH2-like binaries we present two which are our
closets match to the observational systems in Table 4.

Table 4: The parameters of the two Gaia BH1-like and Gaia BH2-like binaries
that we find closely matching the observed systems: the orbital period
$P_{\mathrm{orb}}$ in days, the eccentricity $e$, BH mass $M_{\mathrm{BH}}$
and companion star mass $M_{\star}$ in solar units and the systemic velocity
$v_{\mathrm{sys}}$.

| Gaia BH1-like | Gaia BH2-like
---|---|---
$P_{\mathrm{orb}}$ [d] | $181.07$ | $1591.63$
$e$ | $0.48$ | $0.56$
$M_{\mathrm{BH}}$ [$\mathrm{M_{\odot}}$] | $8.22$ | $8.55$
$M_{\star}$ [$\mathrm{M_{\odot}}$] | $1.38$ | $1.14$
$v_{\mathrm{pec}}\,[\mathrm{km/s}]$ | $57.73$ | $35.57$

####  2.2.2 The evolution of Gaia BH1- and Gaia BH2-like systems

The conditions that we set on BH masses for Gaia BH1- and Gaia BH2-like
binaries limit the ZAMS masses of BH progenitors to the range
$M_{\mathrm{ZAMS}}\sim 34.0-40.0\,\mathrm{M_{\odot}}$ (for $Z=0.01$). We find
that the initial orbital separations may reach up to $\sim
36\,000\,\mathrm{R_{\odot}}$ for very high initial eccentricities ($e\sim
0.9$). For the moderate initial eccentricities $e<0.4$ the initial orbital
separations are in the range $\sim 6000-9000\,\mathrm{R}_{\odot}$.
Whatever their starting parameters are, all Gaia BH1-like binaries follow the
same evolutionary path which we describe on one example also shown on Fig. 1.
The evolution starts with two main sequence stars with masses
$39.7\,\mathrm{M_{\odot}}$ and $1.2\,\mathrm{M_{\odot}}$ on a wide
($a=7\,942.3\,\mathrm{R_{\odot}}$), moderately eccentric ($e=0.31$) orbit. In
the course of the nuclear evolution of the more massive star (the primary) the
orbit widens due to the strong stellar winds and it reaches $a=8\,542.1$
$\mathrm{R_{\odot}}$ at the time when the star leaves the Hertzsprung gap and
starts to burn helium in its core. The still expanding radius of the core-
Helium burning giant (CHeB giant) is large enough ($\sim
1550\,\mathrm{R_{\odot}}$)that the tides starts to act between the two stars
causing the orbit to gradually tighten and circularize. By the time the orbit
becomes circular the orbital separation shrinks to
$a=5\,840.9\,\mathrm{R_{\odot}}$. But the tidal interactions are still at work
trying to speed up the slowly rotating massive star to synchronise it with its
companion which is spinning 3 orders of magnitude faster than the giant. The
synchronization cannot be reached and Darwin instability sets in leading to
further orbital decay until the giant star overflows its Roche lobe. Up to
this moment the evolved star has lost $\sim 57\%$ of its initial mass (its
mass is $17.2\,\mathrm{M_{\odot}}$), the mass of MS star has remained almost
unchanged and the orbital separation between the stars decreased to
$a=3\,256.2\,\mathrm{R_{\odot}}$. The mass transfer that follows is
dynamically unstable even if one considers the more restrictive conditions for
unstable mass transfer than in our standard model (see Pavlovskii et al.
(2017); Olejak & Belczynski (2021) for details). The evolved donor which is
near the end of the core He burning has the developed outer convective
envelope and the low effective temperature $T_{\mathrm{eff}}\sim 3.55\times
10^{3}$ K. The sufficiently low binding energy of such an envelope leads to
its successful ejection during the common envelope (CE) phase.
The giant star loses its envelope becoming a Wolf-Rayet star on a close orbit
($a=140.3\,\mathrm{R_{\odot}}$) with MS star. In a subsequent evolution
lasting $\sim 0.22$ Myr strong winds drive further mass loss from the primary
star what widens the orbit to $a=164.4\,\mathrm{R_{\odot}}$ just prior the
supernova (SN) explosion. The SN explodes $\sim 5.7$ Myr from the onset of the
binary evolution. The primary star becomes $9.7\,\mathrm{M_{\odot}}$ BH and
the effective kick (the composition of Blaauw and natal kicks) imparted to the
system changes its orbital parameters to $a=318.8\,\mathrm{R_{\odot}}$ and
$e=0.49$. Together with the MS primary of mass $1.2\,\mathrm{M_{\odot}}$ the
system becomes Gaia BH1-like binary.
Our synthetic Gaia BH1 binary evolution proceeds on the nuclear timescale of
the secondary which slowly evolves towards the end of hydrogen burning in its
core. It takes $\sim 5$ Gyr for the star to leave the main sequence and become
a red giant what marks the end of Gaia BH1 stage of its evolution. As the
radius of red giant expands the tides act to circularize the orbit. Around
$\sim 600$ Myr after the star became a red giant its radius reaches the radius
of the star Roche lobe and the mass transfer onto BH begins. The system turns
into low-mass X-ray binary (LMXB) of orbital period $P_{\mathrm{orb}}\sim
139.8$ days ($a=251.2\,\mathrm{R_{\odot}}$) and donor mass
$M_{2}=1.2\,\mathrm{M_{\odot}}$. The mass transfer proceeds from less to more
massive binary component, therefore, the orbit widens and so does the radius
of the secondary star Roche lobe. When the Roche lobe extends to the point
that the star fits inside it again the mass transfer ceases. It happens when
the orbital separation is $a=891.7$ $\mathrm{R_{\odot}}$, $M_{\mathrm{BH}}\sim
10.2$ $\mathrm{M_{\odot}}$ and $M_{2}\sim 0.6$ $\mathrm{M_{\odot}}$. The LMXB
phase lasts for $\sim 4$ Myr. From there the secondary star continues its
nuclear evolution. Because the mass of the star is below the limit for helium
ignition the red giant transforms into the helium white dwarf (WD)
($M_{\mathrm{WD}}=0.5$ $\mathrm{M_{\odot}}$ ) after shedding its diluted
hydrogen envelope in stellar wind. At the Hubble time $t_{\mathrm{Hub}}=13.7$
Gyr the distance between the two compact objects in the system is $a=1244.3$
$\mathrm{R_{\odot}}$. For the rest of its life BH-helium WD binary slowly
spirals in due to the gravitational radiation. The system will not merge
within Hubble time.

Figure 1: The evolution of Gaia BH1-like binary from the isolated binary
evolution channel. The chosen example has the parameters close to the
observational ones.

To illustrate the evolution of Gaia BH2 binary we choose the binary which
starts as the system consisting of massive MS star of mass
$22.6\,\mathrm{M_{\odot}}$ and MS star of mass $1.1\,\mathrm{M_{\odot}}$
orbiting each other on the moderately eccentric ($e=0.66$) wide orbit
($a=7\,474.6\,\mathrm{R_{\odot}}$). As in the case of Gaia BH1, the system
first evolves on the more massive (primary) star nuclear evolution time. At
first the orbital separation increases due to the stellar wind from the
primary but when the primary becomes a core-helium burning giant and its
radius reaches $R_{2}=520.4\,\mathrm{R_{\odot}}$ the tidal force starts to
control the orbital evolution. At this point the orbital separation is
$a=7\,794.7\,\mathrm{R_{\odot}}$ and it stars shrinking. It takes $0.323$ Myr
for the tidal force to circularize the orbit. While the orbit keeps shrinking
the primary evolves onto the early asymptotic giant branch (AGB). The mass of
the primary just before it explodes as SN has decreased to
$M_{1}=18.5\,\mathrm{M_{\odot}}$ and the orbital separation has contracted by
over a half of its initial size, to $a=3043.1\,\mathrm{R_{\odot}}$. The
massive star becomes a $9.2\,\mathrm{M_{\odot}}$ BH after $\sim 9.53$ Myr from
the start of the evolution. The large amount of mass ($\sim 50\%$) that it
loses in the explosion results in high eccentricity ($e=0.77$) that the orbit
gains mainly due to the Blaauw kick but at the same time the systemic velocity
remains low $v_{\mathrm{sys}}\sim 32.8\mathrm{km\,s^{-1}}$ due to low natal
kick velocity ($v_{\mathrm{kick}}\sim 34.8\mathrm{km\,s^{-1}}$).
The system consisting of BH-RG forms $\sim 6.2$ Gyr after the onset of the
binary evolution. While RG expands the tides act to tighten and circularize
the orbit. When the orbital separation and the eccentricity decrease to
$a=1175.3\,\mathrm{R_{\odot}}$ and $e=0.59$ respectively, they fall into the
ranges which define the binary as Gaia BH2-like (see Table 2). This phase of
the binary evolution lasts only $\sim 2$ Myr. In next $\sim 4$ Myr the
eccentricity drops to zero and the orbital separation reaches its minimum,
i.e. $a=757.7\,\mathrm{R_{\odot}}$. At this point the wind from the expanding
RG comes into play: the orbit widens again, and BH accrets part of the matter
lost by the evolved companion. By the time He is ignited in the core ($\sim
6.7$ Myr after orbit circularization) the star mass drops to
$M_{2}=0.75\,\mathrm{M_{\odot}}$, BH mass slightly increases to
$M_{\mathrm{BH}}=9.4\,\mathrm{M_{\odot}}$ and orbit widens to
$a=1165.1\,\mathrm{R_{\odot}}$. The exhaustion of He in the core occurs after
$\sim 112$ Myr and the companion becomes an early AGB star (E-AGB) and after
next $\sim 5.7$ Myr it evolves to become the thermally-pulsing AGB star (TP-
AGB). The AGB stage, when the star loses roughly $25\%$ of its mass and the
orbit widens by almost a factor of two reaching
$a=2216.9\,\mathrm{R_{\odot}}$, is followed by to the formation of carbon
oxygen WD (CO WD) of mass $M_{\mathrm{WD}}=0.57\,\mathrm{M_{\odot}}$ orbiting
BH of mass $M_{\mathrm{BH}}=9.5\,\mathrm{M_{\odot}}$ at a distance
$a=2246.6\,\mathrm{R_{\odot}}$.

Figure 2: The main stages of Gaia BH2-like binary evolution from the isolated
binary evolution channel. The details are described in Sec. 2.2.2.

####  2.2.3 Systemic velocities and natal kicks

We present $v_{\mathrm{sys}}$ distributions of standard Gaia BH-like binaries
in models V1 (Hobbs et al. (2005) natal kick velocity distribution) and V2
(Zhao et al. (2023) natal kick velocity distribution) in green in Fig. 3. On
top of them we show the distribution of $v_{\mathrm{sys}}$ for subsets of
binaries which we identify as Gaia BH1-like (magenta) and Gaia BH2-like
(blue), the systems which have $v_{\mathrm{sys}}$ exactly in the observational
range for those systems are marked in pattern. The observed $v_{\mathrm{sys}}$
of Gaia BH1 and Gaia BH2 are within the ranges predicted by isolated binary
evolution for Gaia BH1-like and Gaia BH2-like binaries.
For all standard Gaia BH-like binaries the range of achievable
$v_{\mathrm{sys}}$ is wide (up to $\sim 300$ km/s in model V1 and $\sim 200$
km/s in model V2) but the median values are well below $100$ km/s: it is
$v_{\mathrm{sys,med,V1}}=72.8$ km/s for all Gaia BH-like binaries in model V1
and $v_{\mathrm{sys,med,V2}}=45.3$ km/s in model V2. For Gaia BH1-like
binaries $v_{\mathrm{sys,med,V1}}\sim 54.8$ km/s in model V1 and
$v_{\mathrm{sys,med,V2}}\sim 35.4$ km/s in model V2, for Gaia BH2-like
binaries $v_{\mathrm{sys,med,V1}}\sim 19.2$ km/s in model V1 and
$v_{\mathrm{med,sys,V2}}\sim 16.7$ km/s in model V2. For Gaia BH2-like
binaries the median systemic velocity is comparable between the two natal kick
velocity distributions while Gaia BH1-like binaries move on average by $\sim
55\%$ faster in model V1 than in model V2.
Overall there is $25\%$ more Gaia BH-like binaries in model V2 than in model
V1. Considering the two specific subsets of binaries it appears that the
number of Gaia BH1-like systems is almost the same in both models ($\pm 3$
systems) but Gaia BH2 are almost seven times more abundant in model V2. It is
a straightforward consequence of the lower mean velocity of the natal kick
velocities distribution of Zhao et al. (2023) and low velocities of natal
kicks allow the wide and eccentric binaries like Gaia BH2 to survive SN.
The corresponding median magnitude of the natal kick is not very different
from $v_{\mathrm{sys}}$ due to the low mass of the companion star ($\sim
1\,\mathrm{M_{\odot}}$) and for standard Gaia BH-like binaries it is
$v_{\mathrm{kick,med,V1}}=82.9$ km/s in model V1 and
$v_{\mathrm{kick,med,V2}}=51.0$ km/s in model V2.
The distributions of the natal kick velocities received by Gaia BH1-like and
Gaia BH2-like binaries for model V2 are shown on top and bottom panels of Fig.
4 correspondingly. The median natal kick magnitude is $v_{\mathrm{kick}}\sim
39.3$ km/s for Gaia BH1-like systems and $v_{\mathrm{kick}}\sim 18.7$ km/s for
Gaia BH2-like systems. We mark the distributions of natal kick velocities for
Gaia BH1- and Gaia BH2-like binaries which have the systemic velocities
exactly in the observational limits (Zhao et al., 2023) in dark-blue.
In our further analysis in Sec. 2.3.2 we focus on model V2 as it is well
motivated by observations, it gives us better statistics of Gaia BH2 binaries
and the conclusions of Sec. 2.3.2 do not depend on the choice of
$v_{\mathrm{kick}}$ distribution.

Figure 3: The distribution of the systemic velocity of all standard Gaia BH-
like binaries (green) in model V1 (top) and model V2 (bottom). The magenta
distribution stands for the subset of Gaia BH1-like binaries and the
distribution in blue is for Gaia BH2-like binaries. The binaries in both
groups which have $v_{\mathrm{sys}}$ exactly in the observational ranges are
shown as boxes with pattern. The criteria that we adopt for each Gaia BH-like
group are listed in Table 2.

Figure 4: The distribution of the natal kick velocity received by the
synthetic Gaia BH1-like (top) and Gaia BH2-like (bottom) in model V2. The
dark-blue distribution stands for the subset of the binaries that have
$v_{\mathrm{sys}}$ in the range corresponding to the observational
estimations: $v_{\mathrm{sys}}\in[60.0-80.0]$ km/s for Gaia BH1-like and
$v_{\mathrm{sys}}\in[30.0-40.0]$ km/s for Gaia BH2-like.

###  2.3 Discussion

####  2.3.1 The formation rate of Gaia BH-like binaries

The total stellar mass of one IBE model is $M_{\mathrm{sim}}=3.917\times
10^{11}\,\mathrm{M_{\odot}}$. It comes from our choice to evolve only the
binaries which may be the progenitors of Gaia BH-like binaries, that is we
take primary ZAMS mass to be $M_{\mathrm{1,ZAMS}}>18\mathrm{M_{\odot}}$ (so
that it becomes BH) and the second star ZAMS mass to be
$M_{\mathrm{2,ZAMS}}<1.5\mathrm{M_{\odot}}$. We then evolve $N=10^{6}$ such
binaries. To calculate $M_{\mathrm{sim}}$ we draw the binaries from the whole
IMF and sum up their masses and we count the number of the subset of drawn
binaries that have $M_{\mathrm{1,ZAMS}}>18\mathrm{M_{\odot}}$ and
$M_{\mathrm{2,ZAMS}}<1.5\mathrm{M_{\odot}}$. We draw the binaries (summing up
their masses) until we have a subset of N ($10^{6}$) binaries. We also take
into account the median number of natal kick draws that are needed to get a
Gaia BH-like binary. Therefore, the calculated formation rates are corrected
for using truncated IMF for both of the binary members and for the number of
kick draws.
If we scale down IBE simulations to match the total simulated initial cluster
mass of star clusters simulations $M_{\rm cl,tot}$$=3.6\times
10^{6}\,\mathrm{M_{\odot}}$ (see Sec. 3.2) we find one Gaia BH-like binary per
simulation. This corresponds to the result of the dynamical interaction
channel where also one Gaia BH0like binary ejected from the cluster was found
in a given model cluster set. In conclusion, both channels give the same
formation rate that is $\sim 2.7\times 10^{-7}\,\mathrm{M_{\odot}}^{-1}$.

####  2.3.2 The impact of natal kicks on Gaia BH-like binaries

The natal kick velocity $v_{\mathrm{kick}}$ is the vector that plays an
important if not the decisive role in the Gaia BH-like systems formation. What
counts is not only the magnitude of the kick but also its direction and
relative position of the stars on their mutual orbit when SN takes place.
The natal kick of a given magnitude determines $v_{\mathrm{sys}}$ but not the
eccentricity of the binary what is demonstrated in Fig. 5. However, there is a
correlation between the eccentricity, post- and pre-SN orbital separation (see
e.g. Kalogera, 1996) which is followed by our population of Gaia BH-like
binaries as is shown in Fig. 6. For Gaia BH1, where the companion star is
still on its main sequence, the evolution has not have enough time to change
the binary orbit from its immediate post-SN configuration. The system may have
widen or shrunken due to the natal kick but the total amount of the orbital
separation change is tightly connected to the system post-SN eccentricity. If
the natal kick causes binary orbit to widen the orbital separation may
increase no more than twice ($a_{\mathrm{post-SN}}/a_{\mathrm{pre-SN}}\sim 2$)
for binaries that have eccentricities $e\in[0.4;0.6]$. If the orbit shrinks
due to the natal kick then it can shrink by no more than $40\%$
($0.6<a_{\mathrm{post-SN}}/a_{\mathrm{pre-SN}}<1.0$) (see Fig. 6) for such
binaries.
There is no similar dependency for Gaia BH2-like binaries. They harbor an
already evolved star (RG) and the tidal forces between the binary components
influence the orbital separation and eccentricity during Gaia BH2-like phase.
In the simulations we catch most of Gaia BH2-like binaries when they reach the
upper limits for $a$ and $e$ that we set for Gaia BH2-like classification
during their post-SN evolution. The relation between $a_{\mathrm{post-SN}}$,
$a_{\mathrm{pre-SN}}$ and $e$ for those systems is not straightforward.

Figure 5: The relation between the binary eccentricity $e$, systemic velocity
$v_{\mathrm{sys}}$ (km/s) and natal kick $v_{\mathrm{kick}}$ (km/s) magnitudes
in model V2 for "optimistic" Gaia BH-like binaries. We note that a given
$v_{\mathrm{kick}}$ may result in a whole range of eccentricities.

Figure 6: The relation between the change of the binary orbital separation
after and just before SN ($a_{\mathrm{post-SN}}$/$a_{\mathrm{pre-SN}}$) and
the eccentricity ($e$). The blue circles are Gaia BH1-like binaries and the
green crosses are Gaia BH2-like binaries. For binaries which tend to form Gaia
BH1-like systems the fraction by which the orbit may change during first SN is
strictly limited.

Whether the binary shrinks or widens after SN depends on the components of the
vector position of the exploding star on the orbit around its companion at the
moment of SN explosion and on the components of natal kick vector, where the
components of both vectors are drawn from distributions described in Sec. 2.1.
Apart from the shrinkage/widening of the orbit the natal kick changes the
angle between the direction of BH spin and the direction of the binary angular
momentum ($\theta$). We show the relation between angles $\theta$, $\gamma$
and $v_{\mathrm{sys}}$ (top) and $a_{\mathrm{pre-SN}}$/$a_{\mathrm{post-SN}}$
(bottom) for Gaia BH1-like and Gaia BH2-like binaries in Fig. 7. On both plots
Gaia BH1-like binaries are marked with stars and Gaia BH2-like binaries are
marked with filled circles. The formation of Gaia BH1-like binaries depends
differently on $\gamma$ than the formation of Gaia BH2-like binaries and so
does $\theta$.
Most of Gaia BH1-like binaries form when the natal kick is directed out of the
binary orbital plane. If $\gamma\sim 90^{\circ}$ only a few ($<20$) Gaia
BH1-like binaries form and only with BH spin almost aligned with binary
angular momentum vector ($\theta<2^{\circ}$). The higher is the natal kick
velocity the further from the binary orbital plane the kick has to be oriented
and the more BH spin is misaligned from binary angular momentum in Gaia
BH1-like binary, specifically $\sim 94\%$ of Gaia BH1-like population forms
with $\theta<40^{\circ}$ and the median is $\theta_{\mathrm{med}}\sim
10.5^{\circ}$ which can be seen in Fig. 8.
In the case of Gaia BH2-like binaries $\gamma$ is a major factor determining
their population. Most of them ($\sim 95\%$) form if natal kick direction is
close to the orbital plane ($\gamma\sim 90^{\circ}\pm 15$) and
$v_{\mathrm{kick}}<41.6$ km/s. The binaries that survive SN and form Gaia
BH2-like systems have $\theta$ in the whole range of values. It can be clearly
seen from the $\gamma$ and $\theta$ distributions shown in Fig. 9: $\theta$
distribution is almost flat and $\gamma$ distribution is centred around
$90^{\circ}$.
There is a dependence between the natal kick direction and the
tightening/expanding of binary orbit after SN. All systems marked in black in
lower plot of Fig. 7 have expanded due to the natal kick and the majority of
them are Gaia BH1-like binaries. The only Gaia BH1-like binaries that have
formed through tightening of the orbit after SN are those that form a
pronounced narrow blue strip in bottom plot of Fig. 7. There is almost one-to-
one relation between $\gamma$ and $\theta$ which allows the binary to shrink
after SN and to form Gaia BH1-like system.
The natal kick has to tighten the orbit if the binary is to form Gaia BH2-like
system. The wider is the pre-SN binary the lower velocity has to have the kick
and the closer to the orbital plane it has to be oriented to become Gaia
BH2-like.

Figure 7: The relation between the misalignment angle between the binary
angular momentum vector and BH spin vector after SN ($\theta$), the angle
between the natal kick velocity vector and orbital angular momentum of the
system $\gamma$ and natal kick velocity magnitude (top) or the the ratio of
pre-SN to post-SN orbital separation (bottom). Both angles are in degrees
units and $v_{\mathrm{nk}}$ is in km/s. Gaia BH1-like systems are marked with
stars and Gaia BH2-like systems with filled circles.

The distribution of angles $\theta$ and $\gamma$ for Gaia BH1-like and Gaia
BH2-like binaries are presented on top and bottom plots in Fig. 8 and Fig. 9
respectively. The subsets of Gaia BH1-like and Gaia BH2-like binaries that
have $v_{\mathrm{sys}}$ matching the observations are marked in dark colors.
There is a specific range in $\gamma$ distribution which results in Gaia
BH1-like systems having observed $v_{\mathrm{sys}}$: most of them form if the
natal kick acts at the angle $30^{\circ}-45^{\circ}$ to the orbital plane. In
general the systems with observational $v_{\mathrm{sys}}$ fit into the
distributions of $\theta$ and $\gamma$ for Gaia BH1-like and Gaia BH2-like
populations.

Figure 8: The distributions of $theta$ (top) and $\gamma$ (bottom) for Gaia
BH1-like systems. In dark color is marked the part of the distribution
corresponding to Gaia BH1-like binaries with $v_{\mathrm{sys}}$ in
observational ranges.

Figure 9: The distributions of $theta$ (top) and $\gamma$ (bottom) for Gaia
BH2-like systems. In dark color is marked the part of the distribution
corresponding to Gaia BH2-like binaries with $v_{\mathrm{sys}}$ in
observational ranges.

####  2.3.3 Common envelope parameters $\lambda$ and $\alpha_{\mathrm{ce}}$

The standard version of StarTrack has been using the $\lambda$ formulae from
Xu & Li (2010a); Xu & Li (2010b). In course of work on this paper we found
that for the stars with ZAMS masses in the range $\sim
35-45\,\mathrm{M_{\odot}}$, which are the progenitors of $\sim
8-10\,\mathrm{M_{\odot}}$ BHs, Xu & Li (2010b) prescription underestimated the
binding energy stars at the core He-buring stage of their evolution, resulting
in the high $\lambda$ values ($\lambda>4.0$). The formulae were fitted to the
models calculated with the old Eddington’s code EV and we learnt (Xiang-Dong
Li private communication) that they were updated using Modules for Experiments
in Stellar Astrophysics code (MESA, Paxton et al. (2011)) in 2016 by Wang et
al. (2016) and should be used in their most recent form. In calculations done
for the purpose of this paper we incorporated the formulate from Wang et al.
(2016) for stars more massive than $18.0\,\mathrm{M_{\odot}}$ (BH
progenitors). We chose the formulae which were calculated using the
prescription for the wind mass loss rates from Hurley et al. (2000) and Vink
et al. (2001) (Wind1 in Wang et al. (2016)) and we take $\lambda$ as the
average of $\lambda_{\mathrm{g}}$ (calculated taking into account
gravitational binding energy only), $\lambda_{\mathrm{b}}$ (calculated taking
into account the total energy, that is both gravitational and internal energy
which consists of thermal energy, radiation energy, ionization and
dissociation energy) and $\lambda_{\mathrm{h}}$ (which in addition to the
total energy takes into account the enthalpy of the stellar envelope).
One should be aware that the results may vary depending which $\lambda$
prescription is used in the calculations. To see the evolution of stellar
envelope binding energy for different $\lambda$ prescriptions we refer to
Iorio et al. (2023) where envelope binding energy is plotted in Fig. A1 in
Appendix. Unfortunately, it does not take into account the updated Wang et al.
(2016) formulae.
We would like to point out that our results are obtained with standard value
of $\alpha_{\mathrm{ce}}=1.0$. There are analytical estimations done for the
final orbital separation after common envelope phase in pre-Gaia BH1 systems,
however they are based on the assumption that the star has to loose
$20\,\mathrm{M_{\odot}}$ of mass in the envelope during CE which high
overestimation when stellar winds are taken into account. The mass of BH-
progenitor prior to CE stage is much lower: a star that starts on ZAMS with
$\sim 40\,\mathrm{M_{\odot}}$ looses $\sim 50\%$ of its initial mass in the
stellar winds by the time CE phase starts. Consequently, the envelope that has
to be unbound during CE has mass of only $\sim 3\,\mathrm{M_{\odot}}$. If
$M_{\mathrm{env}}=3.0\,\mathrm{M_{\odot}}$ is taken for calculations the post-
CE orbital separation of the binary is $>100\,\mathrm{R_{\odot}}$ instead of a
few $\mathrm{R_{\odot}}$. Invoking $\alpha_{\mathrm{ce}}>1$ is not necessary
to form Gaia BH-like binaries in IBE channel.

####  2.3.4 Gaia BH1 as Low Mass X-ray Binary

Although there is no X-ray emission detected from Gaia BH1 and Gaia BH2
(Rodriguez et al., 2023), it was found that during the post-dormant BH binary
evolution of both of those binaries they may become the symbiotic X-ray
binaries with the subsequent phase of stable Roche-lobe overflow (RLOF)
Rodriguez et al. (2023) becoming a a low mass X-ray binary (LMXB). The
population of LMXBs that originated from Gaia BH1-like binaries would have
long orbital periods ($P_{\mathrm{orb}}>100$ days).
Of the two binaries that we consider as the best examples of Gaia BH1 and Gaia
BH2 in our IBE simulations and which evolution we describe in Sec. 2.2.2 only
Gaia BH1-like system undergoes RLOF becoming LMXB in its subsequent evolution
while Gaia BH2-like binary remains detached.
As of today, all known LMXBs with BH accretors are transient sources, that is
they show periodic rapid increase in their luminosity up to
$L_{\mathrm{X}}\sim 10^{39}$ erg/s lasting from days to years, followed by
long periods (tens of years) of quiescence when their luminosity drops below
$L_{\mathrm{X}}\sim 10^{33}$ erg/s. The outbursts are triggered by an
instability in accretion disk formed around BH. The instability arises from
the changes in opacities when ionisation/recombination sets in, what in turn
entails the changes in the disk viscosity (for details of disk instability
model (DIM) see e.g. Lasota (2001); Hameury (2020)). According to the model
the system is persistently bright (has high level of luminosity with no
outbursts) if the mass transfer rate $\dot{M}_{\mathrm{tr}}$ is higher than
critical value of mass accretion rate at the outer disk radius
($\dot{M}^{+}_{\mathrm{crit}}(R_{\mathrm{d}})$) and the system is transient
(undergoes outbursts) if
$\dot{M}_{\mathrm{tr}}<\dot{M}^{+}_{\mathrm{crit}}(R_{\mathrm{d}})$ and
$\dot{M}_{\mathrm{tr}}>\dot{M}^{-}_{\mathrm{crit}}(R_{\mathrm{in}})$ where
$\dot{M}^{-}_{\mathrm{crit}}(R_{\mathrm{in}})$ is the critical mass accretion
rate at the inner disk radius. The formulae for $\dot{M}_{\mathrm{crit}}(R)$
are taken from Lasota et al. (2008):

$\begin{split}&\dot{M}^{+}_{\mathrm{crit}}(R)=8.07\times
10^{15}\alpha_{0.1}^{-0.01}\,R_{10}^{2.64}\,M_{1}^{-0.89}\,&\mathrm{g\,s^{-1}}\\\
&\dot{M}^{-}_{\mathrm{crit}}(R)=2.64\times
10^{15}\alpha_{0.1}^{0.01}\,R_{10}^{2.58}\,M_{1}^{-0.85}\,&\mathrm{g\,s^{-1}}\end{split}$
(1)

where $\alpha_{0.1}=0.1\alpha$ is the viscosity parameter in units $0.1$,
$R_{10}$ is the radius in units of $10^{10}$ cm and $M_{1}$ is the primary
mass in solar masses.
In following calculations we assumed that the disk outer radius extends to
$90\%$ of BH Roche lobe radius in the outburst and that the inner disk radius
$R_{\mathrm{in}}$ does not coincide with the last stable orbit in the
quiescence but instead it is truncated due to the transition from optically
thick to optically thin radiatively inefficient flow (e.g. Lasota et al.,
1996; De Marco et al., 2021). We further assumed that the truncation happens
at $3.4\times 10^{4}R_{\mathrm{g}}$ where $R_{\mathrm{g}}$ is gravitational
radius $R_{\mathrm{g}}=GM_{\mathrm{BH}}/c^{2}$.
$\dot{M}^{-}_{\mathrm{crit}}(R)$ decreases with $R$ and the assumption about
$R_{\mathrm{in}}$ is the upper limit for the inner disk truncation that is
justified by the observations of V404 Cyg (a transient LMXB) (Bernardini et
al., 2016). The viscosity parameter in the hot disk is taken to be
$\alpha_{0.1}=1.0$ and in the cold disk $\alpha_{0.1}=0.2$ (Lasota, 2001). We
calculated the critical mass accretion rates at the outer and inner disk radii
for parameters that we found for exemplary Gaia BH1-like systems described in
Sec. 2.2.2 at the onset of LMXB phase:

$\begin{split}&\dot{M}^{+}_{\mathrm{crit}}(R_{\mathrm{d}})=7.95\times
10^{-6}\,&\mathrm{M_{\odot}\,yr^{-1}}\\\
&\dot{M}^{-}_{\mathrm{crit}}(R_{\mathrm{in}})=3.93\times
10^{-10}\,&\mathrm{M_{\odot}\,yr^{-1}}\end{split}$ (2)

for $R_{\mathrm{d}}=20.44$ $\mathrm{R_{\odot}}$ and $R_{\mathrm{in}}=0.71$
$\mathrm{R_{\odot}}$.
The mass transfer rate for the considered binary in its LMXB state is $\sim
3.53\times 10^{-8}\,\mathrm{M_{\odot}}\,\mathrm{yr}^{-1}$, therefore,
according to the disk instability model, the system should become a transient
LMXB. If LMXBs that formed from Gaia BH1-like binaries exist, they would be
observed through their outbursts. It is of interest to predict what would be
the outburst recurrence time of such systems. The majority of the disk mass is
accreted onto the compact object during the long outburst. The disk depleted
from matter needs to refill so that another outburst may develop. Following
Hameury & Lasota (2020) the refilling time corresponding to the outburst
recurrence time is:

$t_{\mathrm{recc}}=\frac{\Delta M}{\dot{M}_{\mathrm{tr}}},$ (3)

where $\Delta M$ is the amount of the disk mass accreted during outburst and
$\dot{M}_{\mathrm{tr}}$ is the mass transfer rate.
We assumed that the amount of mass retained in the disk at the end of outburst
is $10\%$ of the maximum disk mass (just before the outburst is triggered) and
for the purpose of this calculations we took $\Delta M=0.9M_{\mathrm{d}}$,
where the disk mass is:

$M_{\mathrm{d}}=3.3\times
10^{26}\psi\alpha_{0.2}^{-0.8}\dot{M}_{19}^{0.7}M_{1}^{0.25}R_{12}^{1.25}\,\,\mathrm{g}$
(4)

where $\psi$ is a factor connected to opacity (for details see Hameury &
Lasota (2020)) and we used $\psi=1.3$ following the authors, $\alpha$ is the
disk viscosity parameter in units 0.2 (Lasota, 2001), $\dot{M}$ is the
accretion rate in units of $10^{19}\,\mathrm{g/s}$, $M_{1}$ is the compact
object mass in solar units and $R_{12}$ is the outer disk radius in outburst
in units of $10^{12}$ cm.
The mass transfer rate (here equal to mass accretion rate) 222The mass
accretion rate on compact object is assumed to be equal to the mass transfer
rate during the stable Roche lobe overflow in StarTrack, but in general this
is not true for the discs of transient LMXBs. However, it is a sufficiently
good approximation for our purpose of the recurrence time estimation. changes
in time with the nuclear evolution of the donor and so does the Roche lobe
radius of the accretor. Therefore, we estimate two limiting values of
$t_{\mathrm{recc}}$: at the onset (minimum $t_{\mathrm{recc}}$) and at the end
(maximum $t_{\mathrm{recc}}$) of LMXB phase of the considered Gaia BH1-like
system. The parameters used for calculations are give in Table 6 in Appendix.
The obtained recurrence times are $\sim 47.3$ yr and $\sim 159.8$ yr
respectively. This suggests that there may be a population of long orbital
period transient LMXBs that originated from dormant BH binaries and that are
hard to observe because their luminosity rises above the detection threshold
only once every few tens of years.

##  3 Gaia BH binaries from young massive and open star clusters

###  3.1 Computed model clusters

It would be worth considering and comparing our isolated-binary formation rate
of Gaia BH1- and Gaia BH2-like binaries with those from star clusters.
Dynamical interactions among stars, stellar binaries, and BHs in star
clusters, particularly in low mass open clusters, is widely proposed to be a
channel for assembling detached BH-star binaries (Banerjee, 2018b; El-Badry et
al., 2023a; Tanikawa et al., 2023a; Rastello et al., 2023; Niccolò Di Carlo et
al., 2023). In the present work, we consider the evolutionary model star
cluster set of Banerjee (2021a, hereafter Ba21) to estimate the formation rate
of dynamical BH-star binaries in the galactic field.

We waive a detailed description and discussion of the computed models in Ba21,
which are available in the Ba21 paper itself and are discussed further in
several follow-up studies (e.g., Banerjee, 2020, 2021b; Belczynski et al.,
2022). To summarize, the Ba21 model set comprise 72 model star clusters that
evolve from a young massive phase to an old, open cluster-like phase. The
clusters’ initial mass cover the range of $10^{4}\mathrm{M_{\odot}}\leq M_{\rm
cl}\leq 10^{5}\mathrm{M_{\odot}}$ and initial size (half mass radius) range
between $1-3$ pc. About half of the clusters have an initial _overall_
primordial binary fraction between $5\%-10\%$; however, the clusters’ massive-
star members, above $m_{\ast}>16\mathrm{M_{\odot}}$, are paired separately
among themselves with an initial binary fraction of 100%. The clusters cover
the metallicity range $0.0001\leq Z\leq 0.02$. All models are subjected to a
solar-neighbourhood-like galactic tidal field.

The cluster models are computed with the direct, star-by-star N-body
integration code NBODY7 (Aarseth, 2012), with various updates as described in
Ba21. In particular, stellar and binary evolution and remnant formation take
place in the models via the coupling of NBODY7 with an updated version of the
fast binary evolution code BSE (Hurley et al., 2002; Banerjee et al., 2020).
In most models, stellar remnants (NSs and BHs) form via the “rapid” remnant
mass model of Fryer et al. (2012), as implemented in Banerjee et al. (2020).

As seen in the above mentioned followup works of Ba21, these cluster models
produce GW mergers with properties and rate consistent with those inferred
from GW observations. The main motivation for choosing the Ba21 models for
this study is that they are evolved either at least until all the in-cluster
BHs are practically depleted or until $\approx 11$ Gyr (see Appendix C of
Ba21). Therefore, the population of escaped BH-star binaries (and other
escaped BH-containing systems) from these clusters can be considered to be the
complete yield, allowing for robust rate estimates for such field objects from
clusters. This is unlike any other recent study addressing the dynamical
formation of detached BH-star binaries (see above). Therefore, in the
following, we shall focus on only those BH-star binaries that have genuinely
escaped from the clusters, i.e., those which have crossed their parent
cluster’s tidal radius and have become tidal-tail members. This would allow
for a fair comparison with the yield from isolated binaries (see above) and
with Gaia BH1 and Gaia BH2, both of which are field objects.

###  3.2 Results: formation rate of field BH-star binaries

Fig. 10 shows the orbital period, $P_{\rm orb}$, eccentricity, $e$, and
companion mass, $M_{\rm comp}$, against BH mass, $M_{\rm BH}$, of all the BH-
star binaries that have escaped into the galactic field from the model
clusters of Ba21. In all panels, the individual data points are colour coded
according to the velocity of ejection, $v_{\rm ej}$, which is the velocity of
the escaping binary’s center of mass, at ejection. Unless otherwise stated,
all quoted or plotted values of a BH-star binary are those when the escaping
binary crosses the instantaneous tidal radius of its parent cluster; i.e., at
the moment of its ‘formal’ escape form the cluster. These confirmed escaped
binaries and their parameters are recovered from the NBODY7 simulation logs,
with the help of the “BINARY ESCAPE” records. In Fig. 10, the four rows
distinguish between assembly mechanism (dynamically or primordially assembled)
and the stellar companion’s evolutionary status (main sequence or evolved) of
the binary.

The figure demonstrates that escaped BH-star binaries with a low mass
companion ($M_{\rm comp}\lesssim 3\mathrm{M_{\odot}}$) are most likely
dynamically assembled, as anticipated by El-Badry et al. (2023a). In the Ba21
clusters, such pairing between a low mass star and a BH can happen only via
dynamical channel since, initially, massive, BH-progenitor stars are always
primordially paired among themselves (and never with a low mass star). The
dynamically assembled BH-MS binaries generally cover the $P_{\rm orb}$,
eccentricity, $M_{\rm comp}$, and peculiar velocity range of Gaia BH1.
However, the escaped BH-evolved binaries (both dynamically and primordially
paired) struggle to meet the observed properties of Gaia BH2. As such, there
are only a handful of BH-evolved binaries that have escaped from the entire
Ba21 set.

If all orbital parameters are compared together, there is only one ejected BH-
star binary that is somewhat comparable to Gaia BH1: this binary has $M_{\rm
comp}=1.1\mathrm{M_{\odot}}$ (main sequence), $M_{\rm
BH}=17.2\mathrm{M_{\odot}}$, $P_{\rm orb}=290$ day, $e=0.31$. The total
computed initial star cluster mass of the Ba21 set being $M_{\rm
cl,tot}=3.6\times 10^{6}\mathrm{M_{\odot}}$, this gives an unweighted (plain)
formation rate of Gaia BH1-like system, per unit initial cluster mass, of
$2.7\times 10^{-7}\mathrm{M_{\odot}}^{-1}$. This particular binary is formed
via a dynamical exchange encounter in one of the $M_{\rm cl}=2\times
10^{4}\mathrm{M_{\odot}}$ cluster and is also ejected dynamically from the
cluster. In the Ba21 set, there is no ejected BH-star system that is or can
evolve (after the ejection) anywhere similar to Gaia BH2.

Table 5 (middle column) shows the plain formation rates for escaped BH-star
systems with a low mass companion ($M_{\rm comp}\leq 3\mathrm{M_{\odot}}$),
with a solar mass companion ($0.8\mathrm{M_{\odot}}\leq M_{\rm comp}\leq
1.1\mathrm{M_{\odot}}$), and of Gaia BH1-type, as obtained from the Ba21 set.
The table (right column) also states the corresponding weighted rates when the
initial masses of the clusters are assumed to follow a power-law distribution
of index $-2$, as motivated by observations (e.g., Portegies Zwart et al.,
2010; Bastian et al., 2012). The weighted rates are similar to the plain rates
since there are only a handful of field low-mass-companion systems (21, 5, and
1 for the three cases, respectively) from a number of model clusters spanning
over a wide initial mass range (see Sec. 3.1). Fig. 11 shows the mass, orbital
parameter, and $v_{\rm ej}$ distributions of the escaped BH-star binary
population in Ba21, with low mass (orange line and circles) and solar mass
(blue line and circles) companion. The distributions are constructed assuming
a $\propto M_{\rm cl}^{-2}$ initial cluster mass distribution. Note that there
is only one specimen that is comparable to a Gaia-detected BH-star binary (see
above). The distributions are nevertheless presented keeping in mind the
overall diverse demographics of field BH-star binary population from these
computed star clusters (c.f. Fig. 10), that may address such binaries
discovered in the future.

As evident from Figs. 10 and 11, the vast majority of dynamically assembled
field BH-star systems have $M_{\rm comp}<3\mathrm{M_{\odot}}$. Nearly half of
them have zero or small eccentricity, an MS companion, and $P_{\rm orb}<1$ day
— these very tight binaries are outcomes of tidal capture between the BH and
the companion star inside the parent cluster and the binary’s subsequent tidal
circularization333In the Ba21 runs, tidal energy loss and analytical tidal
circularization were enabled (option $27=1$; Aarseth 2003).. Inside clusters,
such close binaries behave like single, massive objects most of the time and
are ejected dynamically in ways similar to those of ejection of single BHs
(Kulkarni et al., 1993). These systems are potentially tidally interacting or
mass transferring and are not compatible with the Gaia detected, detached BH-
star binaries. There are also a few wide, $P_{\rm orb}\gtrsim 100$ day, MS-
companion binaries — these systems are assembled in the cluster via dynamical
exchange encounters and are also ejected via dynamical scatterings (see
above), making them eccentric in general.

Inside a cluster, dynamical formation of wide and eccentric BH-star binaries
with a low mass stellar companion is rather common (Banerjee, 2018b). After
formation, such a binary would undergo dynamical hardening whilst inside the
cluster. However, the BH-star pair is also vulnerable inside a cluster
containing a dynamically active BH subsystem, since the low mass stellar
companion can easily be displaced by a much more massive BH via a close
exchange encounter. In other words, the in-cluster BH-star phase is at best a
short-lived, intermediate phase (Banerjee, 2018a). This explains the general
dearth of field BH-low mass star systems from the computed models.

Dynamical ejection of wide BH-RG systems is further suppressed: owing to the
large envelope size of the RG, an in-cluster close interaction with such a
binary almost inevitably results in tidal interaction and/or geometrical
interception with the giant’s envelope, in turn resulting in a merger, CE, or
stable mass transfer with the RG. Indeed, for the Ba21 clusters, all
dynamically paired field BH-RG systems are stable-mass-transferring.
Inspection reveals that the interactions, that ejected these binaries, are
also responsible for putting them in a mass transferring state. However, it is
possible that an escaped, wide enough ($P_{\rm orb}\sim 10^{3}$ day),
eccentric BH-MS binary evolves into a detached, Gaia BH2-like BH-RG binary.

While no such case is recovered from the Ba21 set, one such field BH-MS binary
has been dynamically produced by a newly computed model cluster set. This set
comprises a grid of model clusters with homogeneous initial conditions and
properties: mass $M_{\rm cl}=3\times 10^{4}\mathrm{M_{\odot}}$, size $2$ pc,
“delayed” remnant mass model (Fryer et al., 2012), metallicities $Z=0.0002$,
0.001, 0.005, 0.01, 0.02, and overall and massive-star primordial binary
fractions (see above) of 10% and 100%, respectively. The distributions of the
primordial binaries are similar to those of Ba21. For each $Z$, model clusters
are evolved by placing them on circular orbits on the equatorial plane of an
axisymmetric, static, Milky Way-like, central mass-disk-halo galactic
potential (Allen & Santillan, 1991), at different galactocentric distances,
namely, $r_{\rm g}=1.0$, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 6.0, and 8.5 kpc
(hence, the grid contains 45 evolutionary models in total). All of these
clusters are evolved (using NBODY7) until they are completely dissolved into
the external galactic field (otherwise, at least until all BHs are depleted
from them). A detailed account of this model set (hereafter nicknamed as
“$r_{\rm g}$-grid”) and its outcomes will be described in a forthcoming
dedicated paper.

As for the BH-star yield, $r_{\rm g}$-grid has produced one field BH-MS binary
with $M_{\rm BH}=14.8\mathrm{M_{\odot}}$, $M_{\rm
comp}=1.1\mathrm{M_{\odot}}$, $P_{\rm orb}=1200$ day, $e=0.56$. This binary is
qualitatively similar to the binary demonstrated in Fig. 2 in the latter’s BH-
MS phase, and would evolve into a Gaia BH2-type detached (non-interacting)
binary. The total simulated initial mass of $r_{\rm g}$-grid being $M_{\rm
cl,tot}=1.35\times 10^{6}\mathrm{M_{\odot}}$, the unweighted formation rate of
field Gaia BH2-type systems from the set is $7.4\times
10^{-7}\mathrm{M_{\odot}}^{-1}$. These clusters have not produced any field
BH-star system that is compatible with Gaia BH1.

Notably, $v_{\rm ej}$s of the field BH-star binaries with low mass companion
are typically lower than the observed peculiar velocities of the Gaia BH
binaries (see Fig. 11, lowermost panel). However, the $v_{\rm ej}$
distributions have long tails that extend well beyond the velocities of both
Gaia BH1 and Gaia BH2. In the scenario of dynamical formation and ejection
from moderately massive clusters such as those in Ba21, a generally low
$v_{\rm ej}$ is expected due to the clusters’ low velocity dispersion and
escape velocity. Still, high velocity ejections can rarely occur due to
chaotic interactions inside the clusters, invoving compact subsystems. The
cluster Gaia BH1 and Gaia BH2 rates quoted in Table 5 do not include any
velocity criterion.

###  3.3 Discussions: on spin-orbit inclinations of dynamical BH-star binaries

As for the spin-orbit alignment, it can be expected that wide, dynamically
assembled field BH-star binaries (resembling Gaia BH1 and BH2; see above)
would have an isotropic spin-orbit misalignment w.r.t. the BH member, since
the exchanging BH is generally uncorrelated to the stellar binary (in a
similar way as for dynamical binary black hole mergers; see, e.g., Banerjee et
al. 2023). In case the pre-exchange (in-cluster) star-star binary is a
primordial binary, it may be spin-orbit aligned to some extent. However, after
the exchange, the exchanged BH’s spin would still be uncorrelated to and hence
isotropically misaligned w.r.t. the binary’s orbital angular momentum and the
stellar companion’s spin. Inside the cluster, close dynamical encounters would
impart additional, random, common spin-orbit tilts to the (BH-star) binary’s
members (Trani et al., 2021), so that the isotropic distribution is maintained
until the binary’s escape from the cluster.

On the other hand, for the very tight field BH-star systems that are formed
via in-cluster, dynamically mediated tidal interactions (see above; which are
unlike Gaia BH systems), the tidal interaction and circularization may align
the stellar member’s spin w.r.t. the orbital angular momentum. However, the BH
member would still remain misaligned, presumably isotropically, at least
initially. In the case of a mass transfer forming an accretion disc around the
BH in such a binary, the system can exhibit low frequency quasi-periodic
oscillations (QPO) related to Lense–Thirring precession of the disc (Ingram et
al., 2009, and references therein). In other words, stellar clusters may
dynamically produce a population of field stellar mass BH candidates that
exhibit potential misalignment signatures such as low frequency QPOs. Although
interesting, a detailed study of this is beyond the scope of this work.

A drawback of the present estimates from clusters is that the criteria for
selecting Gaia BH1 and Gaia BH2 systems could not be matched with those for
the isolated-binary counterpart. This is simply due to the sparsity of escaped
Gaia BH systems (see above), which makes applying well defined selection
criteria nearly impossible. Note that, in this work, we have simply compared
the two formation channels and have not attempted to combine them, unlike in
some recent works (Niccolò Di Carlo et al., 2023).

Table 5: Formation rates of field BH-star binaries from young massive and open star clusters. The rates in the first three rows are based on the evolutionary model star cluster set of Ba21. The Gaia BH2-like systems’ rate in the final row is based on a newly computed cluster model grid (see text). Case | Rate (plain) $[\mathrm{M_{\odot}}^{-1}]$ | Rate (weighted) $[\mathrm{M_{\odot}}^{-1}]$
---|---|---
Low mass companion | $5.8\times 10^{-6}$ | $5.3\times 10^{-6}$
Solar mass companion | $1.4\times 10^{-6}$ | $1.4\times 10^{-6}$
Gaia BH1-like | $2.7\times 10^{-7}$ | $3.5\times 10^{-7}$
Gaia BH2-like | $7.4\times 10^{-7}$ | —

Figure 10: Demographics of the BH-star binaries that have escaped into the
galactic field from the evolutionary model star clusters of Ba21. The filled
circles in the panels show the mass of the BH member ($M_{\rm BH}$; X-axes)
against the orbital period ($P_{\rm orb}$; left column), eccentricity (middle
column), and companion-star mass ($M_{\rm comp}$; right column) of these model
BH-star binaries. The binaries are plotted separately, based on the
evolutionary stage of the companion star (MS or evolved) and the assembly
channel (primordial or dynamical), as indicated in each row’s title. The data
points are colour-coded according to the velocity of ejection of the binary
($v_{\rm ej}$; colour bar). All of the plotted values correspond to the event
of the BH-star binary crossing the instantaneous tidal radius of its parent
cluster. For comparison, the observed BH-star binary candidates Gaia BH1, that
contains a main-sequence companion, and Gaia BH2, that contains an evolved
companion, are shown in the upper and lower pair of rows, respectively (empty
squares). For these observed data points, the colour coding is not followed
and their error bars are generally invisible due to the scale of the figure
axes. (The size of the squares are chosen for legibility and does not
represent measurement uncertainties.)  Figure 11: Normalized distributions of
BH mass ($M_{\rm BH}$), orbital period ($P_{\rm orb}$), eccentricity,
companion star mass ($M_{\rm comp}$), and ejection velocity ($v_{\rm ej}$) of
BH-star binaries that have escaped into the galactic field from the
evolutionary model star clusters of Ba21 (upper to lower panels; all values
correspond to the event of the BH-star binary crossing the instantaneous tidal
radius of its parent cluster). The distributions of binaries with a low mass
companion ($M_{\rm comp}<=3\mathrm{M_{\odot}}$; orange line) and with a solar
mass companion ($0.8\mathrm{M_{\odot}}<=M_{\rm comp}<=1.1\mathrm{M_{\odot}}$;
blue line) are shown separately. These distributions incorporate a weighted
distribution of cluster birth mass, $M_{\rm cl}$, according to the power law
$\propto M_{\rm cl}^{-2}$. Furthermore, the two distributions on each panel
are normalized according to the relative number of BH-low mass and BH-solar
mass binaries ejected from the weighted cluster population. The observed Gaia
BH1 and Gaia BH2 binaries are indicated in the panels; in the final panel,
peculiar velocities of the Gaia BH binaries are indicated.

##  4 Conclusions

In this paper we have investigated two possible channels of formation of two
dormant BH binaries Gaia BH1 and Gaia BH2: the isolated binary evolution (Sec.
2) and the dynamical interactions in young massive and open star clusters
(Sec. 3). We used the population synthesis code StarTrack for IBE calculations
(Sec. 2.1) and updated star-by-star N-body integration code NBODY7 for cluster
models (Sec. 3.1). The two systems from IBE Gaia-BH like binaries population,
that match well the parameters of Gaia BH1 and Gaia BH2, were evolved till the
Hubble time to find the characteristics of low-mass X-ray binaries that
possibly form from such systems (Sec.2.3.4). We have also analysed what are
the constrains on magnitude and direction of the natal kick imparted to BH set
by the orbital parameters and systemic velocities of Gaia BH1-like and Gaia
BH2-like binaries (Sec. 2.2.3 and Sec. 2.3.2).

As seen in Sec. 3.2, medium mass ($\sim 10^{4}\,\mathrm{M_{\odot}}$), pc-scale
stellar clusters struggle to produce dynamically paired Gaia BH-type binaries
in the field. Nevertheless, the overall population of field BH-star binaries
from such clusters does present interesting and diverse demographics. In the
case when no other constrains ($P_{\rm orb}$, $e$, $M_{\mathrm{BH}}$) were
imposed on the population of binaries except that they should consist of BH
and a low-mass star, the formation rate of field BH-star binaries is
$1.4-5.8\times 10^{-6}\,\mathrm{M_{\odot}}^{-1}$ depending on the mass of the
companion. This implies that further types of BH-star systems can potentially
be discovered in the future. Of course, internally, a cluster continues to
assemble BH-star systems (as long as its BH population is not fully depleted),
potentially at a much higher rate. The differences between internal and field
BH-star population from clusters will be explored in a future work.
In contrast, the IBE simulations result in a large Gaia BH-like binaries
population, but this is the effect of pre-choosing possible Gaia BH-like
progenitors as described in Sec. 2.3.1. Among Gaia BH-like binaries that are
produced through IBE there are systems with the parameters matching closely
those inferred from observations of Gaia BH1 and Gaia BH2.
The results from IBE simulations were obtained using the standard value of
parameter $\alpha_{\mathrm{ce}}$ which describes the fraction of the orbital
energy that is used to unbind the envelope in CE phase, that is
$\alpha_{\mathrm{ce}}=1.0$ (see Sec. 2.3.3).
While we acknowledge the uncertainties that the stellar and binary evolution
modelling is subject to, our results indicate that isolated binary evolution
should be considered as the scenario equivalent to dynamical interactions in
open clusters for Gaia BH1 and Gaia BH2 binaries formation.

The formation rate of Gaia BH-like binaries inferred from the simulations is
sensitive to the criteria adopted to classify a binary as Gaia BH-like (Sec.
2.2.1). In IBE channel it can differ by three orders of magnitude depending on
what we consider a Gaia BH-like binary (see Table 3). Ensuring that the common
criteria are set and the same total stellar mass per simulation is adopted for
IBE and dynamical interactions in young clusters, we find that the formation
rate of Gaia BH-like binaries for both channels is equal $2.7\times
10^{-7}\,\mathrm{M_{\odot}}^{-1}$. If we consider the general population of
field BH-low mass star ($<3\,\mathrm{M_{\odot}}$) binaries from young massive
and open clusters and corresponding BH-low mass star population from IBE the
formation rate of such binaries in both channels is on order of $\sim
10^{-6}\,\mathrm{M_{\odot}}^{-1}$.

Considering the IBE channel, the formation of Gaia BH1- and Gaia BH2-like
binaries depends not only on the natal kick magnitude but also on two angles:
the angle between the binary angular momentum and the natal kick vectors
($\gamma$), and the misalignment between the binary angular momentum and the
BH-spin vectors ($\theta$). We summarize the results that we inferred from the
model Gaia BH1- and Gaia BH2-like population as follows:

  * •

All Gaia BH1- and Gaia BH2-like binaries have systemic velocities $<80$ km/s
and the peculiar velocities at birth (which correspond to $v_{\mathrm{sys}}$
in the model) of the observed systems fit well into the $v_{\mathrm{sys}}$
distributions of Gaia BH1- and Gaia BH2-like binaries.

  * •

The natal kick velocity has to be lower than $100$ km/s for Gaia BH1- and Gaia
BH2-like binaries to form.

  * •

The population of Gaia BH1-like binaries is not vulnerable to the natal kick
velocity distribution adopted but a sevenfold increase in the number of Gaia
BH2-like binaries in model V2 comparing to model V1 favours the two-component
Maxwellian distribution with low- and high-velocity component ($\sigma\approx
21$ km/s and $\sigma\approx 107$ km/s) suggested by Zhao et al. (2023).

  * •

The median value of natal kick velocity imparted to the systems that become
Gaia BH1-like and Gaia BH2-like binaries is $\sim 39$ km/s for the former and
$\sim 19$ km/s for the latter (assuming the natal kick distribution from Zhao
et al. (2023)).

  * •

$\sim 94\%$ of Gaia BH1-like binaries that form in IBE channel have
$\theta<40^{\circ}$ and the median value is $\theta_{\mathrm{med}}\sim
10^{\circ}$.

  * •

For Gaia BH1-like binary to be formed the kicks should be directed out of the
binary orbital plane: the higher is the natal kick velocity the further from
orbital plane should the kick be directed.

  * •

$\sim 95\%$ of the binaries that become Gaia BH2-like experience the natal
kick directed within $\pm 15^{\circ}$ from the orbital plane.

Tracking the future evolution of Gaia BH-like binaries shows that those
systems may be the progenitors of LMXBs. Those LMXBs would be transient
systems with the recurrence times of order of tens to over hundred years
according to our calculations.

## Acknowledgements

The co-author, Prof. Krzysztof Belczynski, passed away on January 13th 2024.
He contributed to this paper with his ideas and discussions.
I.K. and K.B are supported by the Polish National Science Center (NCN) grant
Maestro (2018/30/A/ST9/00050). S.B. acknowledges funding for this work by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the
project “The dynamics of stellar-mass black holes in dense stellar systems and
their role in gravitational wave generation” (project number 405620641; PI: S.
Banerjee). Special thanks go to tens of thousands of citizen-science project
"Universe@home" (universeathome.pl) enthusiasts that help to develop the
StarTrack population synthesis code used in this study.

## Data Availability

The data presented in this work can be made available based on the request to
the corresponding author.

## References

  * Aarseth (2003) Aarseth S. J., 2003, Gravitational N-Body Simulations, Cambridge University Press, Cambridge, UK, pp. 430. ISBN 0521432723
  * Aarseth (2012) Aarseth S. J., 2012, MNRAS, 422, 841
  * Allen & Santillan (1991) Allen C., Santillan A., 1991, Rev. Mex. Astron. Astrofis., 22, 255
  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
  * Atri et al. (2019) Atri P., et al., 2019, MNRAS, 489, 3116
  * Banerjee (2018a) Banerjee S., 2018a, MNRAS, 473, 909
  * Banerjee (2018b) Banerjee S., 2018b, MNRAS, 481, 5123
  * Banerjee (2020) Banerjee S., 2020, Phys. Rev. D, 102, 103002
  * Banerjee (2021a) Banerjee S., 2021a, MNRAS, 500, 3002
  * Banerjee (2021b) Banerjee S., 2021b, MNRAS, 503, 3371
  * Banerjee et al. (2020) Banerjee S., Belczynski K., Fryer C. L., Berczik P., Hurley J. R., Spurzem R., Wang L., 2020, A&A, 639, A41
  * Banerjee et al. (2023) Banerjee S., Olejak A., Belczynski K., 2023, ApJ, 953, 80
  * Barack & Cutler (2004) Barack L., Cutler C., 2004, Phys. Rev. D, 70, 122002
  * Bastian et al. (2012) Bastian N., et al., 2012, MNRAS, 419, 2606
  * Belczynski et al. (2002) Belczynski K., Kalogera V., Bulik T., 2002, ApJ, 572, 407
  * Belczynski et al. (2008) Belczynski K., Kalogera V., Rasio F. A., Taam R. E., Zezas A., Bulik T., Maccarone T. J., Ivanova N., 2008, ApJS, 174, 223
  * Belczynski et al. (2010) Belczynski K., Bulik T., Fryer C. L., Ruiter A., Valsecchi F., Vink J. S., Hurley J. R., 2010, ApJ, 714, 1217
  * Belczynski et al. (2016) Belczynski K., Repetto S., Holz D. E., O’Shaughnessy R., Bulik T., Berti E., Fryer C., Dominik M., 2016, ApJ, 819, 108
  * Belczynski et al. (2022) Belczynski K., Doctor Z., Zevin M., Olejak A., Banerje S., Chattopadhyay D., 2022, ApJ, 935, 126
  * Bernardini et al. (2016) Bernardini F., Russell D. M., Shaw A. W., Lewis F., Charles P. A., Koljonen K. I. I., Lasota J. P., Casares J., 2016, ApJ, 818, L5
  * Blaauw (1961) Blaauw A., 1961, Bull. Astron. Inst. Netherlands, 15, 265
  * Chugai (1984) Chugai N. N., 1984, Soviet Astronomy Letters, 10, 87
  * De Marco et al. (2021) De Marco B., Zdziarski A. A., Ponti G., Migliori G., Belloni T. M., Segovia Otero A., Dziełak M. A., Lai E. V., 2021, A&A, 654, A14
  * El-Badry et al. (2023a) El-Badry K., et al., 2023a, MNRAS, 518, 1057
  * El-Badry et al. (2023b) El-Badry K., et al., 2023b, MNRAS, 521, 4323
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91
  * Giesers et al. (2018) Giesers B., et al., 2018, MNRAS, 475, L15
  * Guseinov & Zel’dovich (1966) Guseinov O. K., Zel’dovich Y. B., 1966, Soviet Ast., 10, 251
  * Hamann & Koesterke (1998) Hamann W. R., Koesterke L., 1998, A&A, 335, 1003
  * Hameury (2020) Hameury J. M., 2020, Advances in Space Research, 66, 1004
  * Hameury & Lasota (2020) Hameury J. M., Lasota J. P., 2020, A&A, 643, A171
  * Hobbs et al. (2005) Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS, 360, 974
  * Hurley et al. (2000) Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543
  * Hurley et al. (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, Monthly Notices of the Royal Astronomical Society, 329, 897
  * Ingram et al. (2009) Ingram A., Done C., Fragile P. C., 2009, MNRAS, 397, L101
  * Iorio et al. (2023) Iorio G., et al., 2023, MNRAS, 524, 426
  * Janka (2013) Janka H.-T., 2013, MNRAS, 434, 1355
  * Janka (2017) Janka H.-T., 2017, ApJ, 837, 84
  * Kalogera (1996) Kalogera V., 1996, ApJ, 471, 352
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kroupa et al. (1993) Kroupa P., Tout C. A., Gilmore G., 1993, MNRAS, 262, 545
  * Kulkarni et al. (1993) Kulkarni S. R., Hut P., McMillan S., 1993, Nature, 364, 421
  * Lasota (2001) Lasota J.-P., 2001, New Astron. Rev., 45, 449
  * Lasota et al. (1996) Lasota J. P., Narayan R., Yi I., 1996, A&A, 314, 813
  * Lasota et al. (2008) Lasota J. P., Dubus G., Kruk K., 2008, A&A, 486, 523
  * Mahy et al. (2022) Mahy L., et al., 2022, A&A, 664, A159
  * Nelemans et al. (2001) Nelemans G., Yungelson L. R., Portegies Zwart S. F., 2001, A&A, 375, 890
  * Niccolò Di Carlo et al. (2023) Niccolò Di Carlo U., Agrawal P., Rodriguez C. L., Breivik K., 2023, arXiv e-prints, p. arXiv:2306.13121
  * O’Doherty et al. (2023) O’Doherty T. N., Bahramian A., Miller-Jones J. C. A., Goodwin A. J., Mandel I., Willcox R., Atri P., Strader J., 2023, MNRAS, 521, 2504
  * Olejak & Belczynski (2021) Olejak A., Belczynski K., 2021, ApJ, 921, L2
  * Olejak et al. (2020) Olejak A., Belczynski K., Bulik T., Sobolewska M., 2020, A&A, 638, A94
  * Pavlovskii et al. (2017) Pavlovskii K., Ivanova N., Belczynski K., Van K. X., 2017, MNRAS, 465, 2092
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Portegies Zwart et al. (2010) Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A, 48, 431
  * Rastello et al. (2023) Rastello S., Iorio G., Mapelli M., Arca-Sedda M., Di Carlo U. N., Escobar G. J., Torniamenti S., Shenar T., 2023, arXiv e-prints, p. arXiv:2306.14679
  * Robson et al. (2019) Robson T., Cornish N. J., Liu C., 2019, Classical and Quantum Gravity, 36, 105011
  * Rodriguez et al. (2023) Rodriguez A. C., Cendes Y., El-Badry K., Berger E., 2023, arXiv e-prints, p. arXiv:2311.05685
  * Sana et al. (2012) Sana H., et al., 2012, Science, 337, 444
  * Sana et al. (2013) Sana H., et al., 2013, A&A, 550, A107
  * Shenar et al. (2022) Shenar T., et al., 2022, Nature Astronomy, 6, 1085
  * Tanikawa et al. (2023a) Tanikawa A., Cary S., Shikauchi M., Wang L., Fujii M. S., 2023a, arXiv e-prints, p. arXiv:2303.05743
  * Tanikawa et al. (2023b) Tanikawa A., Hattori K., Kawanaka N., Kinugawa T., Shikauchi M., Tsuna D., 2023b, ApJ, 946, 79
  * Trani et al. (2021) Trani A. A., Tanikawa A., Fujii M. S., Leigh N. W. C., Kumamoto J., 2021, MNRAS, 504, 910
  * Vink & de Koter (2005) Vink J. S., de Koter A., 2005, A&A, 442, 587
  * Vink et al. (2001) Vink J. S., de Koter A., Lamers H. J. G. L. M., 2001, A&A, 369, 574
  * Wang et al. (2016) Wang C., Jia K., Li X.-D., 2016, Research in Astronomy and Astrophysics, 16, 126
  * Webbink (1984) Webbink R. F., 1984, ApJ, 277, 355
  * Xu & Li (2010a) Xu X.-J., Li X.-D., 2010a, ApJ, 716, 114
  * Xu & Li (2010b) Xu X.-J., Li X.-D., 2010b, ApJ, 722, 1985
  * Zhao et al. (2023) Zhao Y., Gandhi P., Dashwood Brown C., Knigge C., Charles P. A., Maccarone T. J., Nuchvanichakul P., 2023, MNRAS, 525, 1498
  * de Kool (1989) de Kool M., 1989, in Bulletin of the American Astronomical Society. p. 1081

## Appendix A

]  onset finish $M_{\mathrm{BH}}$ [$\mathrm{M_{\odot}}$] 9.7 10.2
$\dot{M}_{\mathrm{tr}}$ [g/s] $4.45\times 10^{18}$ $6.07\times 10^{18}$
$R_{\mathrm{d}}$ [cm] $9.61\times 10^{12}$ $2.73\times 10^{13}$

Table 6: The parameters used to calculate the outburst recurrence times at the
onset and finish of LMXB phase.

##  Appendix B Gaia BH1-like binaries as LISA sources

We checked if Gaia BH1-like binaries become LISA sources at some point of
their late evolution using the parameters of a binary described in Sec. 2.2.2.
We calculate the gravitational waves (GW) strain amplitude $h$ taking into
account that Gaia BH1 will become BH-WD binary according to our results. We
start from the formula for the gravitational wave strain derived by Nelemans
et al. (2001):

$h=1.0\times
10^{-21}\frac{\sqrt{g(n,e)}}{n}\left(\frac{M_{\mathrm{Ch}}}{\mathrm{M}_{\odot}}\right)^{5/3}\left(\frac{P_{\mathrm{orb}}}{1\,\mathrm{hr}}\right)^{-2/3}\left(\frac{d}{1\,\mathrm{kpc}}\right)^{-1},$
(5)

where $g(n,e)$ is the Fourier decomposition of GW signal to $n$ harmonics,
given by the combination of Bessel functions of first kind (see Barack &
Cutler, 2004, Eq.10), $M_{\mathrm{Ch}}$ is a chirp mass in solar masses
defined as
$M_{\mathrm{Ch}}=(M_{\mathrm{BH}}M_{2})^{3/5}(M_{\mathrm{BH}}+M_{2})^{-1/5}$,
$P_{\mathrm{orb}}$ is the orbital period in hours and $d$ is the distance to
the binary in kilo-parsecs.
For circular binaries the factor $g(n,e)$ is non-zero only when $n=2$
($g(2,0)=1.0$), therefore, the frequency of GW emitted by circular binaries is
always twice their orbital frequency. We calculate the GW strain $h(n,e)$ for
$n=2$, $e=0.0$, $M_{\mathrm{BH}}=10.2\,$ $\mathrm{M_{\odot}}$, $M_{2}=0.5$
$\mathrm{M_{\odot}}$, $P_{\mathrm{orb}}=979.6$ days and $d=0.48$ kpc, that is
when Gaia BH1 becomes BH-WD binary:

$h(2,0)=2.03\times 10^{-23}\,\,\,\,\,\,\,f_{\mathrm{GW}}=4.99\times
10^{-7}\mathrm{Hz}$ (6)

Comparing the result with the expected LISA sensitivity (Robson et al., 2019)
it is clear that Gaia BH1 will not contribute to LISA sources.
The Gaia BH2 parameters, when it becomes BH-WD binary, are comparable to Gaia
BH1 except for the orbital separation which is $\sim 2.5$ times wider. It is
apparent that GW emitted from Gaia BH2 will be undetectable by LISA likewise,
therefore we skip the calculations of GW strain for Gaia BH2.