# Compact Binary Formation in Open Star Clusters I: High Formation Efficiency
of Gaia BHs and Their Multiplicities

Ataru Tanikawa1, Savannah Cary2 , Minori Shikauchi3,4,5, Long Wang6,7, Michiko
S. Fujii2
1Department of Earth Science and Astronomy, College of Arts and Sciences, The
University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8902, Japan
2Department of Astronomy, Graduate School of Science, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
3Department of Physics, the University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo
113-0033, Japan
4Research Center for the Early Universe (RESCEU), the University of Tokyo,
7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan
5Department of Physics and Astronomy, the University of British Columbia, 6224
Agricultural Road, Vancouver, BC, V6T 1Z1, Canada
6School of Physics and Astronomy, Sun Yat-sen University, Daxue Road, Zhuhai,
519082, China
7CSST Science Center for the Guangdong-Hong Kong-Macau Greater Bay Area,
Zhuhai, 519082, China
E-mail<EMAIL_ADDRESS>0000-0002-8461-5517 0000-0003-1860-1632
0000-0002-3561-8658 0000-0001-8713-0366 0000-0002-6465-2978

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Gaia BHs, black hole (BH) binaries discovered from database of an astrometric
telescope Gaia, pose a question to the standard binary evolution model. We
have assessed if Gaia BHs can be formed through dynamical capture in open
clusters rather than through isolated binary evolution. We have performed
gravitational $N$-body simulations of $100$ open clusters with
$10^{5}M_{\odot}$ in total for each metallicity $Z=0.02$, $0.01$, and $0.005$.
We have discovered one Gaia BH-like binary escaping from an open cluster, and
found that the formation efficiency of Gaia BHs in open clusters ($\sim
10^{-5}M_{\odot}^{-1}$) is larger than in isolated binaries ($\sim
10^{-8}M_{\odot}^{-1}$) by 3 orders of magnitude. The Gaia BH-like binary is
the inner binary of a triple star system. Gaia BHs can have tertiary stars
frequently, if they are formed in open clusters. Combining additional $N$-body
simulations with 8000 open clusters with $8\times 10^{6}M_{\odot}$, we have
estimated the number of Gaia BHs in the Milky Way disk to $10^{4}-10^{5}$
(depending on the definitions of Gaia BHs), large enough for the number of
Gaia BHs discovered so far. Our results indicate that the discoveries of Gaia
BHs do not request the reconstruction of the standard binary evolution model,
and that Gaia BHs are a probe for the dynamics of open clusters already
evaporated.

###### keywords:

stars: black holes – galaxies: star clusters: general – Astrometry and
celestial mechanics – (stars:) binaries (including multiple): close

††pubyear: 2017††pagerange: Compact Binary Formation in Open Star Clusters I:
High Formation Efficiency of Gaia BHs and Their Multiplicities–References

##  1 Introduction

Black holes (BHs) are the final state of massive stars. Albeit of their
darkness, they have been discovered in binary systems, such as X-ray binaries
(Casares et al., 2017), spectroscopic binaries (Shenar et al., 2022),
astrometric binaries (El-Badry et al., 2023a, b; Tanikawa et al., 2023;
Chakrabarti et al., 2023), and gravitational wave sources (The LIGO Scientific
Collaboration et al., 2021). The discoveries of astrometric BH binaries, also
known as Gaia BH1 and BH2 (El-Badry et al., 2023a, b, respectively) in Gaia
DR3 (Gaia Collaboration et al., 2022), challenge the standard binary evolution
model. Gaia BH1 and BH2 (hereafter, Gaia BHs) have $\sim 1M_{\odot}$ visible
stars as companions, orbital periods of $\sim 10^{2}$-$10^{3}$ days, and
moderately high eccentricities ($\sim 0.5$). Isolated binaries can form such
configurations only if common envelope ejection is $\sim 10$ times more
efficient than previously expected (El-Badry et al., 2023a).

Another possible formation channel of Gaia BHs is that a BH dynamically
captures a $\sim 1M_{\odot}$ visible star in open clusters. This is similar to
double BH formation in open clusters (Ziosi et al., 2014; Kumamoto et al.,
2019, 2020; Di Carlo et al., 2019; Banerjee, 2021; Rastello et al., 2021). As
for double BHs, such dynamical capture can happen in globular clusters
(Portegies Zwart & McMillan, 2000; Tanikawa, 2013; Rodriguez et al., 2016;
Askar et al., 2017; Kremer et al., 2022), first star clusters (Wang et al.,
2022), and galactic centers (Tagawa et al., 2020). However, Gaia BHs should
not originate from globular clusters, first star clusters, and galactic
centers, since they are Milky Way disk components, and have near solar
metallicities. They are likely to be formed in open clusters if their
formation channel is dynamical capture.

Gaia BH1 and BH2 have different orbital periods by an order of magnitude, and
different types of visible stars: a main-sequence (MS) star for Gaia BH1 and a
red giant star for Gaia BH2. They may have different origins. What can be said
with certainty is that Gaia BHs have common features in which they are
difficult to be formed through isolated binary evolution, if common envelope
is ejected efficiently (El-Badry et al., 2023a, b). Thus, we explore
possibility that Gaia BHs are formed through dynamical capture in open
clusters.

In this paper, we report an extremely high formation efficiency of Gaia BHs in
open clusters. Moreover, our simulations show that they frequently host
tertiary stars. The structure of this paper is as follows. In section 2, we
describe our numerical method to follow Gaia BH formation in open clusters and
isolated binaries. In section 3, we show our simulation results. In section 4,
we conclude our results, and discuss the implication of our results for Gaia
BHs discovered recently. We focus on a Gaia BH-like binary formed in 100 open
clusters with $10^{5}M_{\odot}$ and metallicity of $Z=0.005$, while we discuss
the dependence of the formation efficiency of Gaia BHs on initial conditions
of open clusters in Appendix.

##  2 Method

We follow the dynamical evolution of open clusters by means of a gravitational
$N$-body code PeTar (Wang et al., 2020b). The PeTar code is based on the
particle-tree and particle-particle algorithm (Oshino et al., 2011; Iwasawa et
al., 2017) with the aid of FDPS (Iwasawa et al., 2016; Iwasawa et al., 2020).
Binary orbital evolution and close encounters are treated by the slow-down
algorithmic regularization method (SDAR: Wang et al., 2020a). The external
potential is calculated by GALPY (Bovy, 2015). We can see comparison between
PeTar and NBODY6++GPU (Nitadori & Aarseth, 2012; Wang et al., 2015; Kamlah et
al., 2022) for mass scale of small globular clusters in Wang et al. (2020b),
and open clusters in Appendix A.

The PeTar code is coupled with the BSE code (Hurley et al., 2002; Banerjee et
al., 2020) to solve single and binary star evolutions. We overview the model
of single and binary star evolution. Single stars evolve on a track developed
by Hurley et al. (2000) with stellar winds formulated by Belczynski et al.
(2010). Massive stars either cause supernovae or collapse to BHs at the end of
their lives, and leave behind NSs and BHs with masses modeled by the rapid
model (Fryer et al., 2012) with modification of pair instability supernovae
modeled by Belczynski et al. (2016b). Binary evolution model is the same as
Banerjee et al. (2020). It contains various binary evolution processes, such
as tidal interaction, wind accretion, stable mass transfer, and common
envelope. Here, we only describe parameters of common envelope which Gaia BHs
formation is the most sensitive to (El-Badry et al., 2023a). A binary can
experience common envelope evolution modeled as the $\alpha$ formalism
(Webbink, 1984), where we adopt $\alpha=3$ and $\lambda$ of Claeys et al.
(2014).

We generate $100$ open clusters for each metallicity of $Z=0.02$, $0.01$, and
$0.005$. The total mass of each open cluster is $10^{3}M_{\odot}$, and then
the total mass of open clusters is $10^{5}M_{\odot}$ for each metallicity. We
adopt the Plummer model for the initial phase space distribution of stars. We
set the half-mass radius to $0.8$ pc. Each cluster orbits around Milky Way at
$8$ kpc with $220$ km s-1. The initial binary fraction is $100$ %. Although
the binary fraction appears too high, it is not sensitive to the formation
efficiency and multiplicity of Gaia BHs as seen in Appendix B. The primary
stars of binaries have Kroupa’s initial mass function (Kroupa, 2001) with the
minimum and maximum masses of $0.08$ and $150M_{\odot}$, respectively. The
mass ratio, orbital period, and eccentricity distribution of binaries follow
initial conditions from Sana et al. (2012). The initial conditions are
generated using MCLUSTER (Küpper et al., 2011). We follow their evolutions
over $500$ Myr. Although the initial central number density seems to be high
($\sim 10^{3}$ pc-3), the central number density of visible stars (defined as
stars except BHs, neutron stars, and white dwarfs) drop to $\sim 3$ pc-3 at
$500$ Myr.

For comparison, we also calculated isolated binary evolution. We prepare two
types of initial conditions for calculations of isolated binary evolution.
Like the binaries in our open clusters, the 1st type of initial conditions for
isolated binary evolution are also adopted from Sana et al. (2012). We find
that the 1st type of initial conditions never form Gaia BHs, because the
minimum mass ratio is $0.1$. In our single star evolution model, $\gtrsim
20M_{\odot}$ stars can leave behind BHs. BHs cannot have $\lesssim 1M_{\odot}$
stars as companions. Thus, we prepare the 2nd type of initial conditions in
which the minimum mass ratio is reduced to $0.0005$, such that even
$150M_{\odot}$ stars can have $0.08M_{\odot}$ stars as companions. We would
like to remark that the 2nd type of initial conditions is much more likely to
form Gaia BHs than the initial conditions of open clusters unless anything
prohibits Gaia BHs from forming. We also note that we do not adopt such a
small mass ratio for initial binaries in open clusters. We follow their
evolution over $30$ Myr. Until that, all BH progenitors evolve to BHs.

Throughout this paper, we compare properties of BH binaries for open clusters
at $500$ Myr, and for isolated binaries at $30$ Myr. We underestimate the
formation efficiency of BH binaries with secondary masses of $\gtrsim
2.5M_{\odot}$, because $\gtrsim 2.5M_{\odot}$ stars evolve to remnant objects
until $500$ Myr. However, this does not affect our purpose. Gaia BHs
discovered so far have secondary stars of $\sim 1M_{\odot}$.

Our simulations generate $59$ BH binaries. Such BH binaries have various
secondary masses, periods, and eccentricities. In this paper, we define Gaia
BHs as BH binaries with secondary masses of $\leq 1.1M_{\odot}$, periods of
$100$-$2000$ days, and eccentricities of $0.3$-$0.9$ unless stated.

##  3 Results

We find that one open cluster with $Z=0.005$ contains a BH binary similar to
Gaia BHs, hereafter called “Gaia BH-like binary”. We summarize the properties
of the Gaia BH-like binary in Table 1. Its BH mass is $21.4M_{\odot}$. Its
companion (hereafter, secondary) star is a MS star and has a mass of
$0.82M_{\odot}$. The binary period is $8.3\times 10^{2}$ days, and its
eccentricity oscillates from $0.3$ to $0.8$ due to von Zeipel-Lidov-Kozai
(ZLK) mechanism (von Zeipel, 1910; Lidov, 1962; Kozai, 1962). We find that the
Gaia BH-like binary is the inner binary of a triple star system, and it is
because of a tertiary star that we see ZLK oscillation of the Gaia BH-like
binary. The tertiary star is a MS star and has $1.59M_{\odot}$ at $500$ Myr.
The outer period and eccentricity are $1.2\times 10^{6}$ days and $0.689$,
respectively.

Table 1: Summary of the Gaia BH-like binaryat $500$ Myr. Parameters | Values | Remarks
---|---|---
Metallicity (Z) | $0.005$ |
BH mass | $21.4M_{\odot}$ |
Secondary mass | $0.82M_{\odot}$ | MS star
Period |  $8.3\times 10^{2}$ days |
Eccentricity |  $0.3$-$0.8$ | Oscillating
Tertiary mass | $1.59M_{\odot}$ | MS star
Outer period |  $1.2\times 10^{6}$ days |
Outer eccentricity | $0.689$ |
Mutual inclination |  $34$-$59$ deg | Oscillating

Figure 1 shows the position and velocity of the Gaia BH-like binary at $500$
Myr. We can see that the Gaia BH-like binary has escaped from the open cluster
located at the coordinate origin. The Gaia BH-like binary will be long-lived
after $500$ Myr as long as it is not perturbed by the tertiary star (described
later). This Gaia BH-like binary escapes from the open cluster through either
of two-body relaxation or close binary interaction111Such a binary with
massive objects can be ejected from an open cluster (Fujii & Portegies Zwart,
2011). Looking at their position and velocity in Figure 1, the Gaia BH-like
binary belongs to a tidal tail of the open cluster. This means that it will be
identified as a Galactic disk component long after the open cluster has
evaporated.

Figure 1: Positions (top) and velocities (bottom) of the Gaia BH-like binary
(star points) and other stars (dots) at $500$ Myr. All the stars including
BHs, neutron stars, and white dwarfs, are plotted. The central number density
of visible stars is $\sim 3.4$ pc-3.

We show the time evolution of the Gaia BH-like binary and its progenitor in
Figure 2, and illustrate the formation process of the Gaia BH-like binary in
Figure 3. At the initial time, its BH progenitor has $22.6M_{\odot}$ and a
companion with $9.32M_{\odot}$. The binary period and eccentricity are
$2.8\times 10^{3}$ days and $0.167$, respectively. At $9.49$ Myr, its BH
progenitor evolves to a $16.7M_{\odot}$ BH. At $21.5$ Myr, this binary is
perturbed by an open cluster star, excites its eccentricity to nearly 1, and
merges. This merger leaves behind a $21.4M_{\odot}$ BH. Until $188$ Myr, this
BH does not form any hard binary. Although soft binaries are occasionally
formed, their binary periods are at least $\sim 10^{7}$ days. Such soft
binaries are quickly disrupted by perturbations of other stars. At $188$ Myr,
the BH captures an open cluster star, and form a relatively hard binary with a
period of $5.3\times 10^{5}$ days. After this, the BH companion is replaced
with another star several times (although omitted in Figure 3). At $244$ Myr
the BH companion becomes a $1.59M_{\odot}$ MS star to be finally the tertiary
star. At $260$ Myr, the BH companion is exchanged with a $0.82M_{\odot}$ MS
star to be finally the secondary star. The superseded companion is not fully
ejected from the Gaia BH-like binary, and stays as the tertiary star. Soon
after this interaction, this Gaia BH-like binary with the tertiary star
escapes from the open cluster.

Figure 2: Time evolution of the Gaia BH-like binary and its progenitor.
Secondary and tertiary stars are defined as temporary companions of the BH and
its progenitor at each time. We indicate the distance from the center,
component masses, binary periods, binary eccentricities, and mutual
inclination from top to bottom. These quantities are plotted for each $4$ Myr.
The shaded regions show binary periods of $100-2000$ days, which tends to be
detected by Gaia. Figure 3: Illustration of the Gaia BH-like binary formation.
Objects indicated by “BH”, “2” and “3” are the BH, secondary and tertiary
stars in the Gaia BH-like binary, respectively. Object “A” is the initial
companion of the BH progenitor. Objects “B”, “C”, “D” are perturbers for the
BH and its companions. Although the Gaia BH-like binary experience many
interactions during $188$-$244$ Myr, we omit them. Thus, an ejected object at
$244$ Myr is not actually object “C”.

We can see in Figure 2 that the Gaia BH-like binary eccentricity and mutual
inclination between the orbital planes of the inner and outer binaries
(hereafter, mutual inclination) are oscillating after the Gaia BH-like binary
escapes due to ZLK mechanism. The ZLK oscillation is modulated, and the
modulation amplitude appears inconstant. This is because the cadence of the
snapshot ($4$ Myr) is longer than the timescale of the ZLK oscillation. In
order to catch the ZLK oscillation by sufficiently high cadence of the
snapshot, we recalculate the orbital evolution of the triple system. For this
purpose, we extract the masses, positions, and velocities of the triple-system
components at $400$ Myr, and follow their evolution over $2$ Gyr by means of
the SDAR code (Wang et al., 2020a), where the SDAR code is incorporated into
the PeTar code to treat close binary and multiple star systems as described
above. We approximate that the triple system is completely isolated from other
stars, since the triple system escapes from the open cluster at $400$ Myr. We
do not take into account the tertiary star evolution to a red giant and white
dwarf phases.

Figure 4: Time evolution of the Gaia BH-like binary eccentricity (top) and
mutual inclination (bottom) over $2$ Gyr from $400$ Myr. The quantities are
plotted for each $0.1$ Myr.

Figure 4 shows the eccentricity and mutual inclination of the Gaia BH-like
binary. The cadence of the snapshot is $0.1$ Myr, which is sufficiently high
to show the ZLK oscillation. The oscillation on several Myr and modulation on
several $100$ Myr are due to the quadrupole-level and octupole-level
interactions, respectively (Antognini, 2015). The modulation amplitude is
constant, although it appears not to be constant due to the low cadence of the
snapshot ($4$ Myr) in Figure 2. We can see that the Gaia BH-like binary
eccentricity rises to at most $0.8$, and its pericenter distance is reduced
down to at least $0.97$ au. Since this pericenter distance is much larger than
the secondary radius, the Gaia BH-like binary will not merge. Note that the
secondary star, a $0.82$ MS star, will not evolve to a giant star within the
Hubble time. The tertiary star evolves to a red giant star and then a white
dwarf over $\sim 2$ Gyr. The red giant star will not fill the Roche lobe of
the outer binary, because the pericenter distance of the outer binary is
$2.0\times 10^{2}$ au much larger than the asymptotic giant branch radius
$\sim 1$ au. The tertiary star gradually decreases its mass and influence on
the Gaia BH-like binary. The tertiary star evolution is expected not to
trigger merger or disruption of the Gaia BH-like binary.

Figure 5: Formation efficiency of BH binaries per $M_{\odot}$ for each bin in
open clusters (top), the 1st type of isolated binaries (middle), and the 2nd
type of isolated binaries (bottom). The stellar type of the companions are
limited to MS stars. The grey, blue, and red histograms indicate All BH
binaries, BH binaries with $\leq 1.1M_{\odot}$ companions, and those with
$\leq 1.1M_{\odot}$ companions and eccentricities from $0.3$ to $0.9$ (i.e.
Gaia BHs). The shaded regions show binary periods of $100-2000$ days, which
tends to be detected by Gaia.

In Figure 5, we compare the formation efficiency among BH binaries in open
clusters to the 1st and 2nd types of isolated binary evolution. We focus on
only binary periods of $100-2000$ days, which tends to be detected by
Gaia222The binary periods of Gaia BH1 and BH2 are $\sim 180$ day (El-Badry et
al., 2023a) and $\sim 1300$ day (El-Badry et al., 2023b; Tanikawa et al.,
2023).. The formation efficiency of all BH binaries in open clusters is much
smaller than both the types of isolated binaries. However, when we limit BH
binaries to Gaia BHs, i.e. those with light companions ($1.1M_{\odot}$) and
moderately high eccentricities ($0.3-0.9$), the formation efficiency is
reversed. It is $\sim 10^{-5}M_{\odot}^{-1}$ in open clusters, while it is
$\sim 10^{-8}M_{\odot}^{-1}$ even in the 2nd type of isolated binaries. Thus,
the formation efficiency is larger by $3$ orders of magnitudes.

We again emphasize that the formation of Gaia BHs in open clusters is much
higher than in isolated binaries with the respect of initial conditions. The
numbers of binaries with $>20M_{\odot}$ and $<1.1M_{\odot}$ stars are $0$,
$0$, and $5.4\times 10^{-5}$ per $1M_{\odot}$ in open clusters and the 1st and
2nd types of isolated binaries, respectively. If we limit their orbital
periods to $100-2000$ days, the number is $9.2\times 10^{-6}$ per $1M_{\odot}$
in the 2nd types of isolated binaries. Note that $>20M_{\odot}$ stars evolve
to BHs long before $500$ Myr, and $<1.1M_{\odot}$ stars are still MS stars at
$500$ Myr if they evolve as single stars; such binaries can be potentially
progenitors of Gaia BHs. On the other hand, the formation efficiencies of Gaia
BHs are $\sim 10^{-5}$, $0$, and $\sim 10^{-8}M_{\odot}^{-1}$ in open clusters
and the 1st and 2nd types of isolated binaries, respectively. Open clusters
are more disadvantageous to form Gaia BHs than the 2nd type of isolated
binaries. However, open clusters can form many more Gaia BHs.

Figure 6: Formation efficiency of BH binaries per $M_{\odot}$ for each bin in
open clusters, and the 1st and 2nd types of isolated binaries as a function of
secondary masses. The top, middle, and bottom panels indicate $Z=0.005$,
$0.01$, and $0.02$, respectively. All the secondary stars are MS stars, and
all the binary periods and eccentricities are $100$-$2000$ days and
$0.3$-$0.9$, respectively.

Figure 6 shows the formation efficiency of BH binaries as a function of
secondary masses. From the top panel, we can confirm that, for the secondary
mass of less than $1.1M_{\odot}$, the formation efficiency in open clusters is
larger than in isolated binaries by 3 orders of magnitude for $Z=0.005$. If BH
binaries with the secondary mass of $\sim 2M_{\odot}$ can be regarded as Gaia
BHs, there is another Gaia BH in a $Z=0.01$ open cluster (see the middle panel
in Figure 6). Such a BH binary is also formed in open clusters more easily
than in isolated binaries by 3 orders of magnitude. Note that we find no
tertiary star for this BH binary.

The mass distribution in Figure 6 may have the sharp peak at $\sim 1M_{\odot}$
for the open clusters with $Z=0.005$ due to the small sampling. If we increase
the number of open cluster simulations to much more than $100$, the sharp peak
is likely to be smeared out around $\sim 1M_{\odot}$. Thus, some of them may
exceed $1M_{\odot}$, and may not be categorized as Gaia BHs. Nevertheless, it
should have a minor effect on the formation efficiency of Gaia BHs in open
clusters.

The low formation efficiency in isolated binaries is consistent with the
argument of El-Badry et al. (2023a). Here, we review its reason. Let’s
consider a binary with a BH progenitor ($\gtrsim 20M_{\odot}$) and a $\sim
1M_{\odot}$ star. Its period is $\sim 10^{3}$ days, and then its separation is
$\sim 10^{3}R_{\odot}$. When the BH progenitor evolves to a giant star, it
fills its Roche lobe, and starts Roche lobe overflow. The Roche lobe overflow
is typically unstable, since the donor star (here, the BH progenitor) is much
more massive than the accretor (here, the $1M_{\odot}$ star). Thus, the binary
experiences common envelope evolution. The common envelope evolution ends with
a merger of the BH progenitor core and $\sim 1M_{\odot}$ star. This is because
the gravitational binding energy between the core and star is too small to
fully eject the common envelope. Even if the binary survives against the
common envelope evolution, its orbital period is much shorter than the orbital
periods of Gaia BHs, and its eccentricity is reduced to $\lesssim 0.2$ (e.g.
Kruckow et al., 2021; Trani et al., 2022). The BH natal kick could excite its
eccentricity if its velocity is sufficiently high. However, the barycenter
velocities of the Gaia BHs are too small to suppose such high BH natal kick
velocities.

In contrast to isolated binaries, the formation of Gaia BHs in open clusters
avoids all the difficulties described above. As illustrated in Figure 3, Gaia
BHs can be formed through dynamical capture after BH progenitors evolve to BHs
in open clusters. If we dare to mention the difficulty of Gaia BH formation in
open clusters, BHs may rarely capture $\sim 1M_{\odot}$ stars. Generally, BHs
capture more massive stars more easily, and there are many $\gtrsim
1M_{\odot}$ MS and white dwarfs in open clusters.

Taking into account only the Gaia BH-like binary, we estimate the number of
Gaia BHs formed in open clusters over Milky Way as follows:

$\displaystyle N_{\rm GaiaBH,MW}$ $\displaystyle\sim 1.6\times
10^{4}\left(\frac{\eta}{10^{-5}M_{\odot}}\right)\left(\frac{M_{\rm
MW}}{6.1\times 10^{10}M_{\odot}}\right)$
$\displaystyle\times\left(\frac{f_{Z=0.005}}{0.26}\right)\left(\frac{f_{\rm
cluster}}{0.1}\right),$ (1)

where $\eta$ is the formation efficiency obtained from our simulation results,
$M_{\rm MW}$ is Milky Way stellar mass (e.g. Licquia & Newman, 2015), and
$f_{Z=0.005}$ is the mass fraction of $Z=0.005$ in the Milky Way, which is
calculated by integrating stellar mass between $Z=0.003$ and $0.007$ in a
Milky Way model of Wagg et al. (2022) and Shikauchi et al. (2023). The
fraction $f_{\rm cluster}$ is the mass fraction of stars formed in open
clusters (Misiriotis et al., 2006; Piskunov et al., 2007). We assume that the
lifetime of Gaia BHs is larger than the Hubble time because of the small mass
of the Gaia BH-like binary. Since the size of Milky Way disk is $\sim 10$ kpc,
the number of Gaia BHs can be $\sim 160$ within $\sim 1$ kpc. We expect that
the formation efficiency in open clusters will explain the number of Gaia BHs,
even if the number of Gaia BHs (currently, $2$) will grow in the future as
seen in the lists of BH binary candidates (Andrews et al., 2022; Shahaf et
al., 2023).

##  4 Conclusion and Discussion

Our main conclusions in this paper are as follows.

  1. 1.

Formation efficiency of Gaia BHs in open clusters is higher than in isolated
binaries by $3$ orders of magnitude.

  2. 2.

Gaia BHs are in multiple star systems with a high probability.

We find a Gaia BH-like binary in our simulations albeit of the small total
mass of open clusters. The Gaia BH-like binary has many points similar to Gaia
BHs. The visible star has $\sim 1M_{\odot}$, and the binary period and
eccentricity are $10^{2}$-$10^{3}$ days and $\sim 0.5$, respectively. Its
barycenter orbit indicates that it is a Galactic disk component. Since the
visible star is dynamically captured by the BH after the BH formation, it
should not be polluted by supernova ejecta from the BH, which is consistent
with the chemical abundance pattern of Gaia BH2 (El-Badry et al., 2023b). The
visible star has subsolar metallicity similar to Gaia BH1, whose visible star
has ${\rm[Fe/H]}=-0.2\pm 0.05$ (El-Badry et al., 2023a). The number of such BH
binaries is large enough to explain the presence of Gaia BHs. The Gaia BH-like
binary contains different points from Gaia BHs, especially Gaia BH2. The
visible star of Gaia BH-like binary has subsolar metallicity and does not
evolve to a red giant star within the Hubble time, while Gaia BH2’s visible
star has solar metallicity, and already evolves to a red giant star. However,
if we simulate the dynamical evolutions of many more open clusters, we are
likely to generate BH binaries similar to Gaia BH2.

We expect that Gaia BHs frequently have tertiary stars if they are formed in
open clusters. The frequency of the presence of tertiary stars is highly
uncertain. If we estimate the frequency from the small samples, the frequency
is $100$ %, because the Gaia BH-like binary, which is just one Gaia BH in our
simulation, has a tertiary star. Even if we include the BH binary with the
secondary mass of $\sim 2M_{\odot}$, this frequency is not changed much
because of its short lifetime. Even if we take into account more many runs in
Appendix B, the frequency is still high $\sim 50$ %. Nevertheless, we
recognize that this is a rough estimate. In order to determine the frequency
accurately, we need to perform many more cluster simulations, and acquire a
larger number of samples than now. On the other hand, there are no reports for
such tertiary stars in Gaia BHs. Its reason may be just that no one searches
for such tertiary stars (but see El-Badry et al., 2021). Moreover, the
tertiary star of the Gaia BH-like binary becomes a white dwarf after $\sim 2$
Gyr, and may be too faint to be detected by Gaia at the present day. Thus, we
are not always able to discover such tertiary stars if any.

As for multiplicity, our results show that a visible binary can orbit a BH. In
other words, the BH itself is a tertiary star of the multiple star system. We
observe several binaries orbiting BHs in our simulations. However, their outer
binary periods are $>10^{5}$ days. They cannot be discovered by Gaia, since
their periods are much larger than Gaia’s operation duration, a few $10^{3}$
days. We may generate such systems detectable by Gaia if we increase the
number of open cluster simulations again.

It appears that there are large uncertainties in our results, in other words,
the formation efficiency of Gaia BHs and multiplicity frequency. This is
derived from only one Gaia BHs formed from open clusters with $3\times
10^{5}M_{\odot}$ in total. However, we expect that our results would be
robust. This is because our results are not different from those in Appendix B
by an order of magnitude. In Appendix B, we investigate open clusters with
$8\times 10^{6}M_{\odot}$ in total and the number of Gaia BH-like binaries is
$11$.

At the moment when we posted this paper to arXiv, only Shikauchi et al. (2020)
have previously investigated Gaia BHs formed in open clusters, although many
studies have studied Gaia BHs formed in isolated binaries (Mashian & Loeb,
2017; Breivik et al., 2017; Yamaguchi et al., 2018; Kinugawa & Yamaguchi,
2018; Yalinewich et al., 2018; Shao & Li, 2019; Andrews et al., 2019;
Wiktorowicz et al., 2020; Chawla et al., 2022; Shikauchi et al., 2022, 2023).
They have not found BH binaries with $\sim 1M_{\odot}$ MS stars despite that
they treat open clusters with $9\times 10^{5}M_{\odot}$ in total, $3$ times
larger than those in this paper. No Gaia BH-like binary in Shikauchi et al.
(2020) can be explained by small statistics. No Gaia BH-like binary can be
formed even in some of model grids in which open clusters have
$10^{6}M_{\odot}$ in total according to Appendix B regardless of binary
fractions.

We should make a caveat that the presence of tertiary stars of Gaia BHs does
not mean the evidence of Gaia BH formation in open clusters. A few 10 % of
galactic field stars are in triple and quadruple systems (Tokovinin, 2014; Moe
& Di Stefano, 2017). Such multiple systems may form Gaia BHs with tertiary
stars. We expect that this possibility will be investigated owing to rapid
progress for studying evolution of multiple systems (Hamers et al., 2021;
Toonen et al., 2022; Stegmann et al., 2022; Dorozsmai et al., 2023), although
it has not yet.333El-Badry et al. (2023a, b) form Gaia BHs from triple
systems. However these Gaia BHs do not have tertiary stars, since the inner
binaries merge in the halfway of their evolutions. Comparison between tertiary
stars of Gaia BHs formed in open clusters and multiple systems will helpful
for identifying the origin of Gaia BHs. In any case, what we can say now is
that the open cluster scenario will be ruled out if no Gaia BHs have tertiary
stars at all.

In the near future, we will perform a larger number of open cluster
simulations than in this paper. These open clusters will have different mass
and density, and will be in different tidal fields. It will reveal the
formation efficiency of Gaia BHs in open clusters more accurately. We will
analyze the distribution of BH and visible star masses, binary periods,
eccentricities, and frequency of the presence of tertiary stars. It is also
important to obtain the distribution of the tertiary star masses, and outer
binary periods and eccentricities. If an outer binary has a short enough
period to cause the ZLK oscillation of its inner BH binary, we may observe an
eccentricity variability of the BH binary (Hayashi & Suto, 2020). We may also
see visible binaries orbiting BHs with periods of $100$-$2000$ days. The
presence of such binaries may be indicated by differences between
spectroscopic and astrometric mass functions for each binary (Tanikawa et al.,
2023). These distributions should be a clue to elucidate the formation process
of Gaia BHs.

This paper indicates two important possibilities. First, we do not need to
rebuild the theory of common envelope evolution, which is carefully
constructed to explain several types of compact binaries, such as double BHs
(e.g. Belczynski et al., 2016a). Second, if Gaia BHs are dominantly formed in
open clusters, they can be a useful probe to study open cluster dynamics.

## Acknowledgments

This research could not be accomplished without the support by Grants-in-Aid
for Scientific Research (17H06360, 19K03907) from the Japan Society for the
Promotion of Science. M.F. is supported by The University of Tokyo Excellent
Young Researcher Program. S.C. is supported by the Fulbright U.S. Student
Program, which is sponsored by the U.S. Department of State and Japan-U.S.
Educational Commission. Its contents are solely the responsibility of the
author and do not necessarily represent the official views of the Fulbright
Program, the Government of the United States, or the Japan-U.S. Educational
Commission. M.S. thanks a support by Research Fellowships of Japan Society for
the Promotion of Science for Young Scientists, by Forefront Physics and
Mathematics Program to Drive Transformation (FoPM), a World-leading Innovative
Graduate Study (WINGS) Program, the University of Tokyo, and by JSPS Overseas
Challenge Program for Young Researchers. L.W. thanks the support from the one-
hundred-talent project of Sun Yat-sen University, the Fundamental Research
Funds for the Central Universities, Sun Yat-sen University (22hytd09) and the
National Natural Science Foundation of China through grant 12073090 and
12233013. Numerical simulations are carried out on Small Parallel Computers at
Center for Computational Astrophysics, National Astronomical Observatory of
Japan, the Yukawa Institute Computer Facility, and Cygnus/Pegasus at the CCS,
University of Tsukuba.

## Data availability

Results will be shared on reasonable request to authors. The data is generated
by the software PeTar and SDAR, which are available in GitHub at
https://github.com/lwang-astro/PeTar and https://github.com/lwang-astro/SDAR,
respectively. The initial conditions of star cluster models are generated by
the software MCLUSTER (Küpper et al., 2011), which is available in GitHub at
https://github.com/lwang-astro/mcluster.

##  Appendix A Comparison between PeTar and NBODY6++GPU

We make comparison between simulation results of PeTar and NBODY6++GPU
(Nitadori & Aarseth, 2012; Wang et al., 2015; Kamlah et al., 2022). Indeed,
Wang et al. (2020b) have already done in detail. However, they have treated
small globular clusters with $\sim 10^{5}$ $M_{\odot}$, while we treat open
clusters with $\sim 10^{3}$ $M_{\odot}$. Thus, it might be worth comparing
simulation results between PeTar and NBODY6++GPU with respect of cluster mass
scale of $\sim 10^{3}$ $M_{\odot}$.

We generate 10 open clusters for each code. The metallicity is $Z=0.02$, and
initial binary fractions are 20 %. We do not take into account galactic tides,
since NBODY6++GPU does not support combination with GALPY. The other setup is
the same as section 2. Note that we do not remove any stars from simulations
even when stars go away from open clusters. This is true for simulations in
section 2.

Figure 7: (Top) Time evolution of the total masses of single and binary stars.
(Bottom) Time evolution of 10, 50, and 90 % Lagrangian radii.

Figure 8: (Top) Distribution of binary semi-major axes and eccentricity at 500
Myr. The point sizes show binary masses on linear scale. (Bottom) Distribution
of binary mass and mass ratio at 500 Myr. The point size indicates binary
semi-major axes on logarithmic scale.

Figures 7 and 8 shows the cluster evolution and binary properties at 500 Myr.
We do not find any different features between PeTar and NBODY6++GPU results.
This indicates that PeTar works well even for mass scale of open clusters.

##  Appendix B Effects of initial binary fractions

We investigate how the formation efficiency of Gaia BHs depend on initial
binary fractions. We adopt the fractions of 0, 20, 50, and 100 % for open
clusters with $Z=0.02$ and $0.005$. We generate 1000 open clusters for each
fraction and each metallicity in order to reduce statistical uncertainties.

Figure 9: (Top) Formation efficiency of Gaia BHs in open clusters (solid
curves) as a function of initial binary fractions. Error bars indicate
standard errors under assumption that Gaia BH formation is a Poisson process.
Dashed lines show the formation efficiency in isolated binaries. In order to
be easy to see, the initial binary fractions are shifted by $0.125$ leftward
and rightward for $Z=0.02$ and $0.05$, respectively. The definition of Gaia
BHs is the same as section 2. The black dotted curve indicates the formation
efficiency proportional to the initial binary fraction. (Bottom) Fraction of
Gaia BHs in multiple star systems as a function of initial binary fractions.
We calculate standard errors, assuming that Gaia BHs obtain tertiary stars in
a Poisson process. The initial binary fractions are shifted in the same way as
the top panel.

The top panel of Figure 10 shows that the formation efficiency of Gaia BHs in
open clusters is not sensitive to initial binary fractions in open clusters.
Although the error bars are large, the average values are not dependent on the
initial binary fractions. The bottom panel of Figure 10 indicates little
dependence of multiplicity on initial binary fractions. Although the results
suffer from small statistics, the multiplicity does not monotonically increase
with initial binary fractions. These results mean that our choice of initial
binary fractions does not affect our results in section 3.

The top panel of Figure 10 shows $\eta\sim 2\times 10^{-6}M_{\odot}^{-1}$ for
all binary fractions and metallicities. Then, we can rewrite Equation (1) as

$\displaystyle N_{\rm GaiaBH,MW}$ $\displaystyle\sim 1.2\times
10^{4}\left(\frac{\eta}{2\times 10^{-6}M_{\odot}}\right)\left(\frac{M_{\rm
MW}}{6.1\times 10^{10}M_{\odot}}\right)$
$\displaystyle\times\left(\frac{f_{\rm cluster}}{0.1}\right).$ (2)

This is similar to the estimate in Equation (1).

Figure 10: The same as Figure 10, except for the definition of Gaia BHs.
Period and eccentricity criteria are relaxed, such that the former and latter
are $0-4000$ days, and $0.0-1.0$, respectively. The companion criteria become
more strict, such that $<0.8M_{\odot}$. Note that $<0.8M_{\odot}$ stars do not
evolve to giant stars within the Hubble time.

We show the formation efficiency and multiplicity fraction of Gaia BHs. Here,
we change the definition of Gaia BHs to mimic the definition in Di Carlo et
al. (2023), such that period is $0-4000$ days, eccentricity is $0.0-1.0$, and
the companion mass is $<0.8M_{\odot}$. Compared to the definition in section
2, the formation efficiency is increased to $\sim 10^{-5}M_{\odot}^{-1}$ by 10
times. This is in good agreement with the results of Di Carlo et al. (2023).
The multiplicity of Gaia BHs in open clusters is $\sim 0.5$. As seen in the
top panel, the formation efficiency may increase with metallicity decreasing.
However, metallicity dependence is beyond the scope of this paper.

In summary, the formation efficiency of Gaia BHs is quite consistent with the
results of Rastello et al. (2023) and Di Carlo et al. (2023). Our formation
efficiency looks different from the formation efficiency of Di Carlo et al.
(2023). However, the difference comes from the difference of the definitions
of Gaia BHs.

## References

  * Andrews et al. (2019) Andrews J. J., Breivik K., Chatterjee S., 2019, ApJ, 886, 68
  * Andrews et al. (2022) Andrews J. J., Taggart K., Foley R., 2022, arXiv e-prints, p. arXiv:2207.00680
  * Antognini (2015) Antognini J. M. O., 2015, MNRAS, 452, 3610
  * Askar et al. (2017) Askar A., Szkudlarek M., Gondek-Rosińska D., Giersz M., Bulik T., 2017, MNRAS, 464, L36
  * Banerjee (2021) Banerjee S., 2021, 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
  * 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. (2016a) Belczynski K., Holz D. E., Bulik T., O’Shaughnessy R., 2016a, Nature, 534, 512
  * Belczynski et al. (2016b) Belczynski K., et al., 2016b, A&A, 594, A97
  * Bovy (2015) Bovy J., 2015, ApJS, 216, 29
  * Breivik et al. (2017) Breivik K., Chatterjee S., Larson S. L., 2017, ApJ, 850, L13
  * Casares et al. (2017) Casares J., Jonker P. G., Israelian G., 2017, X-Ray Binaries. p. 1499, doi:10.1007/978-3-319-21846-5_111
  * Chakrabarti et al. (2023) Chakrabarti S., et al., 2023, AJ, 166, 6
  * Chawla et al. (2022) Chawla C., Chatterjee S., Breivik K., Moorthy C. K., Andrews J. J., Sanderson R. E., 2022, ApJ, 931, 107
  * 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
  * Di Carlo et al. (2019) Di Carlo U. N., Giacobbo N., Mapelli M., Pasquato M., Spera M., Wang L., Haardt F., 2019, MNRAS, 487, 2947
  * Di Carlo et al. (2023) Di Carlo U. N., Agrawal P., Rodriguez C. L., Breivik K., 2023, arXiv e-prints, p. arXiv:2306.13121
  * Dorozsmai et al. (2023) Dorozsmai A., Toonen S., Vigna-Gómez A., de Mink S. E., Kummer F., 2023, arXiv e-prints, p. arXiv:2307.04793
  * El-Badry et al. (2021) El-Badry K., Rix H.-W., Heintz T. M., 2021, MNRAS, 506, 2269
  * 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
  * Fujii & Portegies Zwart (2011) Fujii M. S., Portegies Zwart S., 2011, Science, 334, 1380
  * Gaia Collaboration et al. (2022) Gaia Collaboration et al., 2022, arXiv e-prints, p. arXiv:2206.05595
  * Hamers et al. (2021) Hamers A. S., Rantala A., Neunteufel P., Preece H., Vynatheya P., 2021, MNRAS, 502, 4479
  * Hayashi & Suto (2020) Hayashi T., Suto Y., 2020, ApJ, 897, 29
  * 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, MNRAS, 329, 897
  * Iwasawa et al. (2016) Iwasawa M., Tanikawa A., Hosono N., Nitadori K., Muranushi T., Makino J., 2016, PASJ, 68, 54
  * Iwasawa et al. (2017) Iwasawa M., Oshino S., Fujii M. S., Hori Y., 2017, PASJ, 69, 81
  * Iwasawa et al. (2020) Iwasawa M., Namekata D., Nitadori K., Nomura K., Wang L., Tsubouchi M., Makino J., 2020, PASJ, 72, 13
  * Kamlah et al. (2022) Kamlah A. W. H., et al., 2022, MNRAS, 511, 4060
  * Kinugawa & Yamaguchi (2018) Kinugawa T., Yamaguchi M. S., 2018, arXiv e-prints,
  * Kozai (1962) Kozai Y., 1962, AJ, 67, 591
  * Kremer et al. (2022) Kremer K., Ye C. S., Kıroğlu F., Lombardi J. C., Ransom S. M., Rasio F. A., 2022, ApJ, 934, L1
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kruckow et al. (2021) Kruckow M. U., Neunteufel P. G., Di Stefano R., Gao Y., Kobayashi C., 2021, ApJ, 920, 86
  * Kumamoto et al. (2019) Kumamoto J., Fujii M. S., Tanikawa A., 2019, MNRAS, 486, 3942
  * Kumamoto et al. (2020) Kumamoto J., Fujii M. S., Tanikawa A., 2020, MNRAS, 495, 4268
  * Küpper et al. (2011) Küpper A. H. W., Maschberger T., Kroupa P., Baumgardt H., 2011, MNRAS, 417, 2300
  * Licquia & Newman (2015) Licquia T. C., Newman J. A., 2015, ApJ, 806, 96
  * Lidov (1962) Lidov M. L., 1962, Planet. Space Sci., 9, 719
  * Mashian & Loeb (2017) Mashian N., Loeb A., 2017, MNRAS, 470, 2611
  * Misiriotis et al. (2006) Misiriotis A., Xilouris E. M., Papamastorakis J., Boumis P., Goudis C. D., 2006, A&A, 459, 113
  * Moe & Di Stefano (2017) Moe M., Di Stefano R., 2017, ApJS, 230, 15
  * Nitadori & Aarseth (2012) Nitadori K., Aarseth S. J., 2012, MNRAS, 424, 545
  * Oshino et al. (2011) Oshino S., Funato Y., Makino J., 2011, PASJ, 63, 881
  * Piskunov et al. (2007) Piskunov A. E., Schilbach E., Kharchenko N. V., Röser S., Scholz R. D., 2007, A&A, 468, 151
  * Portegies Zwart & McMillan (2000) Portegies Zwart S. F., McMillan S. L. W., 2000, ApJ, 528, L17
  * Rastello et al. (2021) Rastello S., Mapelli M., Di Carlo U. N., Iorio G., Ballone A., Giacobbo N., Santoliquido F., Torniamenti S., 2021, MNRAS, 507, 3612
  * 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
  * Rodriguez et al. (2016) Rodriguez C. L., Chatterjee S., Rasio F. A., 2016, Phys. Rev. D, 93, 084029
  * Sana et al. (2012) Sana H., et al., 2012, Science, 337, 444
  * Shahaf et al. (2023) Shahaf S., Bashi D., Mazeh T., Faigler S., Arenou F., El-Badry K., Rix H. W., 2023, MNRAS, 518, 2991
  * Shao & Li (2019) Shao Y., Li X.-D., 2019, ApJ, 885, 151
  * Shenar et al. (2022) Shenar T., et al., 2022, Nature Astronomy, 6, 1085
  * Shikauchi et al. (2020) Shikauchi M., Kumamoto J., Tanikawa A., Fujii M. S., 2020, PASJ, 72, 45
  * Shikauchi et al. (2022) Shikauchi M., Tanikawa A., Kawanaka N., 2022, ApJ, 928, 13
  * Shikauchi et al. (2023) Shikauchi M., Tsuna D., Tanikawa A., Kawanaka N., 2023, arXiv e-prints, p. arXiv:2301.07207
  * Stegmann et al. (2022) Stegmann J., Antonini F., Moe M., 2022, MNRAS, 516, 1406
  * Tagawa et al. (2020) Tagawa H., Haiman Z., Kocsis B., 2020, ApJ, 898, 25
  * Tanikawa (2013) Tanikawa A., 2013, MNRAS, 435, 1358
  * Tanikawa et al. (2023) Tanikawa A., Hattori K., Kawanaka N., Kinugawa T., Shikauchi M., Tsuna D., 2023, ApJ, 946, 79
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration et al., 2021, arXiv e-prints, p. arXiv:2111.03606
  * Tokovinin (2014) Tokovinin A., 2014, AJ, 147, 87
  * Toonen et al. (2022) Toonen S., Boekholt T. C. N., Portegies Zwart S., 2022, A&A, 661, A61
  * Trani et al. (2022) Trani A. A., Rieder S., Tanikawa A., Iorio G., Martini R., Karelin G., Glanz H., Portegies Zwart S., 2022, Phys. Rev. D, 106, 043014
  * Wagg et al. (2022) Wagg T., Broekgaarden F. S., de Mink S. E., Frankel N., van Son L. A. C., Justham S., 2022, ApJ, 937, 118
  * Wang et al. (2015) Wang L., Spurzem R., Aarseth S., Nitadori K., Berczik P., Kouwenhoven M. B. N., Naab T., 2015, MNRAS, 450, 4070
  * Wang et al. (2020a) Wang L., Nitadori K., Makino J., 2020a, MNRAS, 493, 3398
  * Wang et al. (2020b) Wang L., Iwasawa M., Nitadori K., Makino J., 2020b, MNRAS, 497, 536
  * Wang et al. (2022) Wang L., Tanikawa A., Fujii M., 2022, MNRAS, 515, 5106
  * Webbink (1984) Webbink R. F., 1984, ApJ, 277, 355
  * Wiktorowicz et al. (2020) Wiktorowicz G., Lu Y., Wyrzykowski Ł., Zhang H., Liu J., Justham S., Belczynski K., 2020, ApJ, 905, 134
  * Yalinewich et al. (2018) Yalinewich A., Beniamini P., Hotokezaka K., Zhu W., 2018, MNRAS, 481, 930
  * Yamaguchi et al. (2018) Yamaguchi M. S., Kawanaka N., Bulik T., Piran T., 2018, ApJ, 861, 21
  * Ziosi et al. (2014) Ziosi B. M., Mapelli M., Branchesi M., Tormen G., 2014, MNRAS, 441, 3703
  * von Zeipel (1910) von Zeipel H., 1910, Astronomische Nachrichten, 183, 345


# Automating Microservices Test Failure Analysis using Kubernetes Cluster Logs

Pawan Kumar Sarika  Ericsson ABStockholmSweden<EMAIL_ADDRESS>,
Deepika Badampudi  Blekinge Institute of TechnologyKarlskronaSweden
<EMAIL_ADDRESS>,  Sai Prashanth Josyula  Blekinge Institute of
TechnologyKarlskronaSweden<EMAIL_ADDRESS>and  Muhammad Usman
Blekinge Institute of TechnologyKarlskronaSweden<EMAIL_ADDRESS>
(2023)

###### Abstract.

Kubernetes is a free, open-source container orchestration system for deploying
and managing Docker containers that host microservices. Kubernetes cluster
logs help in determining the reason for the failure. However, as systems
become more complex, identifying failure reasons manually becomes more
difficult and time-consuming. This study aims to identify effective and
efficient classification algorithms to automatically determine the failure
reason. We compare five classification algorithms, Support Vector Machines,
K-Nearest Neighbors, Random Forest, Gradient Boosting Classifier, and
Multilayer Perceptron. Our results indicate that Random Forest produces good
accuracy while requiring fewer computational resources than other algorithms.

Kubernetes cluster logs, microservices, machine learning

††journalyear: 2023††copyright: rightsretained††conference: Proceedings of the
International Conference on Evaluation and Assessment in Software Engineering;
June 14–16, 2023; Oulu, Finland††booktitle: Proceedings of the International
Conference on Evaluation and Assessment in Software Engineering (EASE ’23),
June 14–16, 2023, Oulu, Finland††doi: 10.1145/3593434.3593472††isbn:
979-8-4007-0044-6/23/06

##  1\. Introduction

Ericsson, a leading Information and Communication Technology (ICT) service
provider, is at the forefront of pioneering cloud
RAN111https://www.ericsson.com/en/ran/cloud, a technology that handles network
traffic in the cloud. To increase development efficiency, Ericsson developed
the Application Development Platform (ADP) ecosystem (Usman et al., 2022),
which among other things, includes a marketplace that hosts over 280
microservices. Among these, 50+ microservices are common microservices that
can be reused and integrated across all applications within the organization.
The microservices in the ADP ecosystem are hosted in Docker containers and
managed by Kubernetes.

Some applications need a specific set of microservices that may differ from
others. In addition, the applications may deploy microservices in different
environments, such as Amazon Web Service (AWS), or Microsoft Azure Cloud. To
ensure high quality, the ADP ecosystem includes a Continuous Integration and
Continuous Deployment (CICD) team, which simulates the different application
environments and performs various tests, including testing scalability,
robustness, and resilience.

Whenever a test fails, the developers must classify the failure into the
following categories - an issue in the cluster, artifactory, microservice,
CICD tests itself, or environment. The classification is done to report the
failure to the relevant team. For example, microservice bugs are assigned to
the relevant microservice team, while CICD bugs are reported to the CICD team
for resolution. In addition, knowing the failure reason can provide Ericsson
insights on what they should improve. These tests generate an average of 450
MB of log data. Developers use their knowledge and experience to classify the
failure, which may take up to two hours for initial analysis. As the number of
microservices and teams using them grows, it becomes increasingly challenging
and time-consuming for developers to classify these failures manually.

This paper reports the experience of automating the classification of failures
using machine learning techniques at Ericsson. The goal of automation is to
identify an efficient and cost-effective method to decrease the time required
for manual analysis. To achieve this goal, we utilized developers’ knowledge
and experiences to better understand the workflows and information necessary
for classifying the failure. Based on this understanding, we extracted
important log data to train a machine-learning model capable of accurately
predicting the failure’s cause.

##  2\. Study Design

Figure 1 represents the steps followed in this study. In this section, we
describe the process followed in each step.

Figure 1. Overall research process.

Sequence Flow diagram

###  2.1. Data Selection

As depicted in Figure 1, expert knowledge was utilized to identify the most
informative sections of Kubernetes cluster logs, which reduced the overall
data processing requirements. Moreover, before automation, it is important to
comprehend the manual process of failure classification. The first author, a
team member responsible for testing and deploying microservices, noted that
each team member had their own approach to analyzing failures. To establish a
consensus, we conducted four in-depth interviews to better understand the
manual failure classification process. The interviewees had varying levels of
experience in identifying reasons for failure, ranging from 1.5 to 4 years.
Furthermore, the interviewees have specialized expertise in several areas
related to manual failure analysis, such as test environments, end-to-end
testing, and coordination with microservice development teams. As a result,
the participants have diverse skills and knowledge that were valuable in the
failure classification process.

We started the interview by explaining the study objective. We then presented
a scenario in which a failure occurred and asked the interviewees to classify
the failure using the same steps they would use in a real-world situation. The
interviewee is free to use any log as an example. An experienced developer
typically examines the sections of the log that are likely to contain the
failure reason. This step allowed us to identify the most frequently visited
Kubernetes cluster logs.

The interviewees were asked to walk through the different steps they followed,
along with the inputs considered and the outputs of each step. Through this
process, we identified the crucial log parts for failure classification and
how developers identified failure reasons.

###  2.2. Data Cleaning

In this study, we utilized important data identified from the interview for
classification. Following data selection, data cleaning was conducted (Figure
1). We ensured data cleanliness by removing inaccurate, incomplete, or
irrelevant data. Human-readable elements (e.g., punctuation, newline
characters) were removed to reduce computational costs. We converted all log
elements to lowercase letters to avoid ambiguities, and unique elements (e.g.,
timestamps, IP addresses, line numbers) were removed. Stopwords and stemming
methods were not used as they rely on the English dictionary and may not work
for non-English words like ’requesthandlerclass.’

###  2.3. Feature Extraction

After cleaning the Kubernetes cluster logs data, the text format must be
converted into a numerical form using Feature Extraction (FE) for algorithms
to process it. Although word-to-vector (word2vec) is a popular method for
feature extraction, it was not considered appropriate due to the wide
variation in the size of each log. Instead, we utilized the TF-IDF method, a
common feature extraction procedure for sentences and documents. TF-IDF was
more suitable for logs since the words related to the error generally occur
once or twice in the entire log. This method highlights words distinctive from
the text, allowing the extraction of essential data features required for
classification. Using TF-IDF, we generate vocabulary, which is used to extract
the feature vectors of the log data.

###  2.4. Classification

Finally, the most crucial step, i.e., log classification, involves using
classification algorithms to construct models that can make accurate
predictions. This study employs supervised learning, as it is necessary to
make predictions for specified classes. Non-parametric algorithms are selected
for this study because they do not require prior knowledge and are
particularly suitable when working with large datasets (J Russell and Norvig,
2010). The chosen algorithms for this study include Support Vector Machines
(SVM), K-Nearest Neighbors (KNN), Random Forest Classifier, Gradient Boosting
Classifier, and Multilayer Perceptron (MLP).

###  2.5. Training

We proceeded with training the selected classification algorithms using the
feature vectors of the training data. We used the scikit-learn library to
obtain the trained models. The SVC() is trained with a linear kernel (SVM).
The RandomForestClassifier() is trained with 1000 trees. All other function
parameters of KNeighborsClassifier() and GradientBoostingClassifier() are set
to default values used by scikit-learn. MLP was built using keras sequential
model consisting of three dense layers and two Dropout layers. The F1 Score
weighted metric is used for evaluation along with the Adam optimizer set to a
learning rate of 0.001, and categorical cross-entropy is used as the loss
metric.

Following the training of the algorithms, we utilized them to predict class
labels for every log in the test data. We collected various metrics, such as
accuracy, F1-score, training time, and prediction time, to compare the
classification algorithms’ performances for predicting failure reasons. We
measured training and prediction times to assess computational cost, which can
be a better metric for real-world performance than time complexity (Lemire,
2021). Stratified 10-fold cross-validation was employed as the sampling
strategy. This approach helped identify the algorithm with the best predictive
accuracy and computational cost for a given hardware and software
configuration.

##  3\. Expert opinion findings on analyzing failure reasons

During the interviews, developers scrutinized a failed test case and analyzed
the Kubernetes cluster logs to classify the failure. Specifically, they
searched for logs associated with the microservice under test within the
pods/container or pods/describe directories. If no relevant logs were found,
the logs of dependent microservices were consulted. Our analysis revealed that
the error or relevant information could usually be found in these directories.
Even when it was not, there were clues indicating potential failure type. As a
result, we concluded that logs located in pods/containers and pods/describe
directories are crucial in most instances. To further strengthen this
conclusion, we also examined the tickets submitted by the developers. We found
that the root directories for 43 of the 50 reported log snippets were located
in pods/container or pods/describe, or both directories.

The pods/containers and pods/describe directories were determined to be
important during the interview phase. The logs were pre-processed, and the
logs’ size was significantly reduced (up to 96.08%) through data selection and
cleaning, reducing computational resources.

##  4\. Developing prediction models

The classification algorithms were used to process the training data and
generate a classification model, which was then used to predict the test
data’s class labels. The performance of the algorithms was evaluated based on
their accuracy, F1-score, training time, and prediction time.

Table 1. Model accuracy and F1-score for each algorithm Algorithm | Accuracy | F1-score
---|---|---
SVM | 0.6905 | 0.6580
KNN | 0.6576 | 0.6580
Random Forest | 0.7279 | 0.7106
Gradient Boosting | 0.7302 | 0.7071
MLP | 0.7086 | 0.6911

We compared the chosen classification algorithms statistically based on the
recommendations by Demšar (2006). We applied the Friedman test to test the
null-hypothesis that all algorithms are equivalent. If this hypothesis was
rejected, we conducted a post-hoc Nemenyi test to compare the algorithms
pairwise.

Table 1 presents the accuracy and F1-scores of each algorithm. The Friedman
test was performed to evaluate the statistical significance of the differences
between the algorithms, resulting in a p-value of 0.00074 and a Q value of
19.12, confirming significant differences in accuracy. We subsequently
conducted a Nemenyi test, and the p-values of the combinations of algorithms
are shown in Table 2 with P-values less than $\alpha$ (0.05) in bold. Random
Forest and Gradient Boosting were statistically better than KNN.

Table 2. Pairwise P-values for accuracy

| SVM | KNN |  | Random
---
Forest
| Gradient
---
Boosting
MLP
SVM | 1.0000 | 0.3517 | 0.4374 | 0.2418 | 0.9000
KNN |  | 1.0000 | 0.0037 | 0.0010 | 0.1143
| Random
---
Forest
|  | 1.0000 | 0.9000 | 0.7633
| Gradient
---
Boosting
|  |  | 1.0000 | 0.5626
MLP |  |  |  |  | 1.0000

Similarly, for F1-scores (Table 3), the Friedman test yielded a p-value of
0.00056 and a Q value of 19.76, indicating significant differences. The
Nemenyi test showed that Random Forest outperformed SVM and KNN with
statistically significant differences. Furthermore, KNN’s performance was
significantly lower than Gradient Boosting’s. Regarding both F1-score and
accuracy, Random Forest, Gradient Boosting, and MLP performed similarly.

We perform the Nemenyi test to determine which algorithms are statistically
different based on F1-scores. Table 2 represents the p-value of F1-score for
each combination of algorithms. The p-value in bold is less than $\alpha$
(0.05).

Table 3. Pairwise P-values for F1-score

| SVM | KNN |  | Random
---
Forest
| Gradient
---
Boosting
MLP
SVM | 1.0000 | 0.9000 | 0.0248 | 0.1143 | 0.4374
KNN |  | 1.0000 | 0.0022 | 0.0160 | 0.1143
| Random
---
Forest
|  | 1.0000 | 0.9000 | 0.6830
| Gradient
---
Boosting
|  |  | 1.0000 | 0.9000
MLP |  |  |  |  | 1.0000

Finally, the algorithms were evaluated based on their training and prediction
times. The average training time in minutes for all ten folds is shown in
Table 4. It was observed that Gradient Boosting took a significantly longer
time for training when compared to other algorithms. On the other hand, MLP
took more than 2x longer for training compared to the Random Forest. Regarding
prediction time, all five algorithms finished predictions in less than a
second, indicating no significant difference between them in real-world
performance. Based on the analysis, it can be concluded that the Random Forest
algorithm performed better in terms of both accuracy and training time.

Table 4. Training and prediction time in minutes Algorithm | Training time | Prediction time
---|---|---
SVM | $0.845$ | $0.082$
KNN | $0.004$ | $0.015$
Random Forest | $1.702$ | $0.006$
Gradient Boosting | $110.773$ | $0.002$
MLP | $4.111$ | $0.005$

##  5\. Piloting the prediction models

To test the model in a production environment, we developed an end-to-end
prototype. The results of the model were added to a dashboard for developer
analysis.

We received feedback from developers, who indicated that while the failure
classification helps them understand the reason for the failure and where they
should look in the logs, they would appreciate additional information beyond
the failure reason. As part of future work, we will devise an additional
solution to include the most similar log previously reported along with its
diagnosis report. Reviewing diagnosis reports of similar logs will help
developers quickly understand the cause of failure. This way, developers can
precisely look at the specific logs increasing their efficiency and helping
them close the cases faster.

In addition, developers expressed concerns about the trustworthiness of the
models’ predictions due to a lack of understanding of how the model arrived at
its conclusions. Further effort is required to improve the models’
interpretability and explainability.

Although Random Forest performed better, we need to tune the hyperparameters
to truly capture the effectiveness of Random Forest compared to other
algorithms.

##  6\. Conclusion

Identifying the reason for a failed CICD test is important to assign the
failure to the relevant team for resolution. The automatic identification of
failure reasons minimizes the manual effort needed to analyze the failure. In
this study, we report experiences of classifying large-scale Kubernetes
cluster logs using machine learning classification algorithms: Support Vector
Machines (SVM), K-Nearest Neighbors (KNN), Random Forest Classifier, Gradient
Boosting Classifier, and Multilayer Perceptron (MLP). Each log generates 450
MB of data on average, which was reduced to 91% by identifying relevant parts
of the logs using expert opinion. We used the TF-IDF feature extraction
method. Then we trained and evaluated the classification algorithms in terms
of accuracy, F1 scores, training, and prediction time.

The results indicated that Random Forest outperforms SVM, KNN, Gradient
Boosting, and MLP in terms of accuracy and computational cost. When evaluating
the proof-of-concept, the developers expect additional information to
strengthen the confidence in the prediction. For example, providing similar
logs and diagnosis reports can help the developers to understand how
previously failed tests were classified manually and arrived at the same
classification predicted by the model. As part of future work, we aim to make
the comparison of the algorithms more rigorous by applying hyperparameter
tuning and further improving the model’s interpretability and explainability.

As a prospect for future research, we intend to enhance the dimensionality of
the input data. Currently, we concatenate logs from different services into a
single file, creating a two-dimensional input. Our plan is to aggregate logs
from different pods within each service, producing a three-dimensional input
where each service has its own file. This approach would enable us to leverage
deep learning algorithms, such as CNNs, which are proficient in processing 3D
data. By doing so, our model would be able to identify the interdependencies
among the services, leading to improved results.

###### Acknowledgements.

The Knowledge Foundation supports this work through the OSIR project
(reference number 20190081) at Blekinge Institute of Technology, Sweden.

## References

  * (1)
  * Demšar (2006) Janez Demšar. 2006\.  Statistical comparisons of classifiers over multiple data sets.  _The Journal of Machine learning research_ 7 (2006), 1–30.
  * J Russell and Norvig (2010) Stuart J Russell and Peter Norvig. 2010.  Artificial Intelligence A Modern Approach Third Edition.
  * Lemire (2021) Author Daniel Lemire. 2021\.  Big-O notation and real-world performance.  Retrieved February 20, 2023 from https://lemire.me/blog/2013/07/11/big-o-notation-and-real-world-performance/
  * Usman et al. (2022) Muhammad Usman, Deepika Badampudi, Chris Smith, and Himansu Nayak. 2022.  An Ecosystem for the Large-Scale Reuse of Microservices in a Cloud-Native Context.  _IEEE Software_ 39, 05 (2022), 68–75.