# Black hole mergers in compact star clusters and massive black hole formation
beyond the mass-gap

Francesco Paolo Rizzuto1, Thorsten Naab1, Rainer Spurzem2,3,4, Manuel Arca-
Sedda2 , Mirek Giersz5, Jeremiah Paul Ostriker6,7, Sambaran Banerjee8,9
1Max-Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1,D-85741,
Garching, Germany
2Astronomisches Rechen-Institut, Zentrum für Astronomie, University of
Heidelberg, Mönchhofstrasse 12-14,69120, Heidelberg, Germany
3National Astronomical Observatories and Key Laboratory of Computational
Astrophysics, Chinese Academy of Sciences, 20A Datun Rd.,
Chaoyang District, 100101, Beijing, China
4Kavli Institute for Astronomy and Astrophysics, Peking University, Yiheyuan
Lu 5, Haidian Qu, 100871, Beijing, China
5Nicolaus Copernicus Astronomical Centre, Polish Academy of Sciences, ul.
Bartycka 18, 00-716 Warsaw, Poland
6Department of Astronomy, Columbia University, New York, NY 10027, USA
7Department of Astrophysical Sciences, Princeton University, Princeton, NJ
08544, USA
8Helmholtz-Institut fur Strahlen- und Kernphysik (HISKP), Nussallee 14-16,
D-53115 Bonn, Germany
9 Argelander Institut für Astronomie (AIfA), Auf dem Hügel 71, D-53121, Bonn,
Germany Contact e-mail<EMAIL_ADDRESS>Fellow at
Frankfurt Institute for Advanced Studies (FIAS)

###### Abstract

We present direct N-body simulations, carried out with nbody6++gpu, of young
and compact low metallicity star clusters with $1.1\times 10^{5}$ stars, a
velocity dispersion of $\sim$ 10 $\mathrm{km/s^{-1}}$, a half mass radius
$R_{h}=0.6$ pc, and a binary fraction of $10\%$ including updated evolution
models for stellar winds and pair-instability supernovae (PISNe). Within the
first tens of megayears of evolution, each cluster hosts several black hole
(BH) merger events which nearly cover the complete mass range of primary and
secondary BH masses for current LIGO/Virgo/Kagra gravitational wave
detections. The importance of gravitational recoil is estimated statistically.
We present several possible formation paths of massive BHs above the assumed
lower PISNe mass-gap limit ($45M_{\odot}$) into the intermediate-mass BH
(IMBH) regime ($>100M_{\odot}$) which include collisions of stars and BHs as
well as the direct collapse of stellar merger remnants with low mass cores.
The stellar evolution updates result in the early formation of higher mass
stellar BHs than for the previous model. The resulting higher collision rates
with massive stars support the rapid formation of massive BHs. For models
assuming a high accretion efficiency for star-BH mergers, we present a first-
generation formation scenario for GW190521-like events – a merger of two BHs
in the PISN mass-gap – which is dominated by star-BH mergers. This IMBH
formation path is independent of gravitational recoil and therefore
conceivable in dense stellar systems with low escape velocities. One simulated
cluster even forms an IMBH binary (153$M_{\odot}$,173$M_{\odot}$) which is
expected to merge within a Hubble time.

###### keywords:

gravitational waves – methods: numerical – stars: black holes – stars:
dynamics – stars: mass-loss – galaxies: star clusters: general.

††pubyear: 2020††pagerange: Black hole mergers in compact star clusters and
massive black hole formation beyond the mass-gap–LABEL:lastpage

##  1 Introduction

Before the first LIGO gravitational-wave (GW) detection, many theoretical
models of stellar evolution predicted stellar black holes (BHs) masses to be
lower than $30$ $\mathrm{M_{\odot}}$. These models remained unchallenged for
several years because all stellar BHs observed at the time had masses
$\lesssim 20$ $\mathrm{M_{\odot}}\,$(Ziółkowski, 2008; Özel et al., 2010).
Surprisingly, the first LIGO detection, GW150914, revealed components more
massive than $30$ $\mathrm{M_{\odot}}$(Abbott et al., 2016). Such masses had
been predicted by stellar evolution models at low metallicity introducing a
dependence between stellar winds mass loss and metallicity (see Woosley et
al., 2002; Vink et al., 2001, and references therein). This highlights the
importance of stellar evolution models for the correct
interpretation/prediction of GW events. An accurate theory for the evolution
of massive stars is particularly important to predict the mass distribution of
stellar BHs at their formation. For this, precise models of stellar winds and
a correct description of the last stages of the stellar evolution before the
collapse are required. At the onset of stellar collapse, stars with
sufficiently large helium cores undergo a phase of electron-positron pair
production that in turn leads to one or more violent explosions. Depending on
the initial mass of the core, the star can experience pulsation pair-
instability supernovae (PPSN) getting partially destroyed or it can experience
the more violent pair-instability supernovae (PSN) and is destroyed completely
(Fowler & Hoyle, 1964; Woosley et al., 2007; Woosley, 2017). Due to (P)PSN,
isolated massive stars are not supposed to collapse into BHs in the mass range
of approximately $50-130$ $\mathrm{M_{\odot}}$. This gap in the stellar BH
mass distribution is known as the (P)PSN mass gap. The mass limits of this gap
are affected by various uncertainties and therefore they depend on the details
of the stellar evolution adopted. In this study, the assumed mass gap is
$45-195$ $\mathrm{M_{\odot}}$.

The 66 - 85 $\mathrm{M_{\odot}}\,$BH merger (GW190521) detected by the
LIGO/Virgo/Kagra collaboration (Abbott et al., 2020) has started a debate
about the origin of the two objects involved in the collision as both of BHs
fall in the (P)PSN mass gap. In particular, the mass of the primary in
GW190521 ($85$ $\mathrm{M_{\odot}}$) poses a challenge to theoretical models.
In a recent paper, Vink et al. (2021) has proposed a modified stellar
evolution model such that $\sim 100$ $\mathrm{M_{\odot}}\,$stars lose only
little mass via stellar winds, small enough for the stars to collapse into BHs
in the range of $\sim 85$ $\mathrm{M_{\odot}}$. A recent study, leveraging
newly estimated uncertainties on (P)PSN mass loss, affirm that GW190521 could
have formed through the classical isolated binary evolution (Belczynski,
2020). Other potential formation channels for BHs of such masses might involve
population III stars (Liu & Bromm, 2020; Kinugawa et al., 2021) or might
depends on details of the rates of carbon-oxygen nuclear reactions (Farmer et
al., 2019).

As both observed merging BHs in question are above the conventional mass gap
lower limit of $\sim 50$ $\mathrm{M_{\odot}}\,$(Farmer et al., 2019; Woosley
et al., 2020), and because there is evidence for residual eccentricity
(Romero-Shaw et al., 2020; Gayathri et al., 2020) and/or spin precession
(Abbott et al., 2020), many authors have argued that GW190521 could have been
produced dynamically through repeated mergers of smaller BHs (Abbott et al.,
2020). The most popular models involve multiple mergers of low mass BHs in
dense and massive stellar environments with an escape velocity of a few
hundred km/s (Arca Sedda & Benacquista, 2019; Fragione et al., 2020; Romero-
Shaw et al., 2020; Kimball et al., 2020; Dall’Amico et al., 2021; Mapelli et
al., 2021). According to these studies, young star clusters and globular
clusters are less favoured because the recoil kick imparted to the merger
remnant could eject the final product of a BH - BH merger and halt its mass
growth too early.

However, many of these studies might have underestimated the importance of
star - BHs collisions for the formation of GW events mass range of GW190521.
In agreement with previous studies (Giersz et al., 2015; Di Carlo et al.,
2019; Kremer et al., 2020; Becker et al., 2020; Rizzuto et al., 2021; Arca
Sedda et al., 2021), the star cluster simulations we present in this paper
show that low mass BHs can reach the (P)PSN mass gap or even the IMBH mass
range through stellar mergers. A similar outcome could be produced invoking
tidal capture events (Stone et al., 2017). Moreover, they illustrate that
GW190521-like events could be hosted by $\sim 10^{5}$
$\mathrm{M_{\odot}}\,$compact star clusters if the BHs involved in the
collisions also grow by the accretion of stellar material.

This work is an extension of Rizzuto et al. (2021) where we have shown that a
few hundred solar masses IMBHs can form in low metallicity compact star
cluster environments via multiple collisions111As pointed out by Mark Gieles
in private communication, the process that leads to the formation of massive
objects in our simulation is not a runaway growth because the growth rate does
not increase with growing mass.. We use updated stellar evolution models with
pulsation pair instability prescriptions and new stellar winds recipes
(following Banerjee et al., 2020) to evolve the most promising initial
conditions presented in Rizzuto et al. (2021). We describe the details of the
stellar evolution models in the next section (Section 2) and present the
initial conditions in Section 3. In Section 4 we discuss our results in
comparison to our previous work. In the final section, we discuss and
summarise the main results.

##  2 The Method

Model Name | $r\mathrm{{}_{c}}$ | $\rho_{c}$ | $R\mathrm{{}_{h}}\,$ | $\sigma$ | $t_{\mathrm{rh}}$ | $t_{\mathrm{s}}$ | $f\mathrm{{}_{c}}\,$ | IMBH | $M_{\mathrm{IMBH}}$ | $t_{\mathrm{form}}$
---|---|---|---|---|---|---|---|---|---|---
| [pc] | [$\mathrm{M_{\odot}}\,$/ pc 3] | [pc] | [km/s] | [Myr] | [Myr] |  |  | [$\mathrm{M_{\odot}}$] | [Myr]
R06W61.0 | $0.19$ | $1.1\times 10^{5}$ | $0.6$ | 15 | 56 | 1.4 | $1.0$ | $6/8$ | 102, 108, 153, 156, 172, 293 | 76, 3.9, 29, 21, 6.9, 44
R06W60.5 | $0.19$ | $1.1\times 10^{5}$ | $0.6$ | 15 | 56 | 1.4 | $0.5$ | $4/8$ | 98, 99, 109, 113 | 34, 12, 6.5, 4.3
Table 1: Model parameters of the cluster simulations: $r\mathrm{{}_{c}}\,$:
initial core radius; $\rho\mathrm{{}_{c}}\,$: initial central density;
$W\mathrm{{}_{0}}\,$: central potential parameter for the King density profile
(King, 1966); $R\mathrm{{}_{h}}\,$: half mass radius; $\sigma$ : dispersion
velocity; $t_{\mathrm{rh}}$: half mass relaxation time; $t_{\mathrm{s}}$:
segregation time scale for $100$ $\mathrm{M_{\odot}}$; $f\mathrm{{}_{c}}\,$:
fraction of mass absorbed by a compact object during a direct collision with a
star; # IMBH : Number of BHs with masses $\gtrsim 100$
$\mathrm{M_{\odot}}\,$formed out of 8 realisations; $M_{\mathrm{IMBH}}$: IMBH
masses; $t_{\mathrm{form}}\,$: IMBHs formation times. Figure 1: The final
black hole mass as a function of initial stellar mass for stellar evolution
model bse (red) used in Rizzuto et al. (2021) and the new model bse-psn
(black) used here with updated stellar wind and pair instability supernova
models. This is for single stars in isolation at metallicity $Z=0.0002$. For
bse-psn the mass gap with with no stellar mass black hole remnant ranges from
$45$ $\mathrm{M_{\odot}}\,$to $195$ $\mathrm{M_{\odot}}$.

We present 16 new simulations of compact star clusters adopting the same
initial conditions as in Rizzuto et al. (2021) using updated the stellar
evolution models in the nbody6++gpu (Wang et al., 2015) code. nbody6++gpu is
capable of evolving star clusters with a realistic number of stars following
the dynamical evolution of the system accurately as well as the stellar
evolution of single and binary stars (see also Wang et al., 2016). The
integrator follows the dynamical evolution of binary stars, even in phases of
mass loss or when one of the two stars experiences a supernova explosion. In
this case, the new trajectories of the remnant and its companion are computed
accurately. (for more details see e.g. Rizzuto et al., 2021, and references
therein). In extension to Rizzuto et al. (2021), the stellar evolution models
of nbody6++gpu have been significantly updated and improved (see Kamlah et
al., 2021, for more details). Most of the updates have been published for
nbody7 already in Banerjee et al. (2020). The major updates are:

  1. 1.

New stellar wind models following Belczynski et al. (2010). which in turn
follow the wind mass-loss rates given by Vink et al. (2001). With these
models, BHs masses that originate from single stars depend strongly on the
metallicity. For instance, a $100$ $\mathrm{M_{\odot}}\,$isolated main-
sequence star would leave a $15$ $\mathrm{M_{\odot}}\,$BH at solar metallicity
($Z=0.02$). At very low metallicity ($Z=0.0002$), however, it can form a BH of
about $60$ $\mathrm{M_{\odot}}$, in the absence of pair-instability models.

  2. 2.

Pair-instability supernova and pulsation pair-instability supernova models
(according to Belczynski et al., 2016) incorporated in the remnant formation
and supernovae models as described in Fryer et al. (2012). Stars with helium
core with masses $>40$ $\mathrm{M_{\odot}}\,$undergo a violent phase of mass
loss. For helium cores in the range between $60-135$ $\mathrm{M_{\odot}}\,$the
star is completely destroyed.

  3. 3.

New prescription of BHs and NSs natal kick velocities that explicitly depend
on the fallback fraction (Banerjee et al., 2020).

  4. 4.

A model for electron capture supernovae (ECSN) following Podsiadlowski et al.
(2004) and Gessner & Janka (2018) that produces neutron stars with low-
velocity kicks that are therefore likely retained in medium-size star
clusters.

The gravitational energy loss and resulting merger of compact objects is
computed following Peters & Mathews (1963) and the implementation is presented
in Rizzuto et al. (2021). An important physical process that is not yet
implemented in the code is relativistic recoil for compact object collisions.
As we will see in the next sections, its absence might artificially enhance
the probability of forming massive BHs (see Banerjee, 2021, for more details
on relativistic recoil implementation in nbody7).

Hereafter, We use the term bse to refer to the previous stellar evolution
prescription (Rizzuto et al., 2021) while for the updated recipes we use the
term bse-psn. The updates in bse-psn have a strong impact on the mass
distribution of the stars and their remnants in our simulations. With bse, an
isolated massive star at metallicity $Z=0.0002$, even as massive as $500$
$\mathrm{M_{\odot}}$, can never collapse into stellar-mass BH more massive
than $30$ $\mathrm{M_{\odot}}$, as shown in Fig. 1 (red line). At the same
metallicity, the updated model bse-psn with new stellar wind and the (P)PSN
models, results in more massive stellar BHs. Main-sequence stars with masses
in the range $80-100$ $\mathrm{M_{\odot}}\,$leave a remnant of about $40$
$\mathrm{M_{\odot}}\,$and very massive stars with masses of a few $\sim
10^{2}$ $\mathrm{M_{\odot}}\,$collapse directly into IMBHs, as illustrated by
the black dots in Fig. 1. The mass-gap, without any remnant black holes, is
thus limited to $45-195$ $\mathrm{M_{\odot}}\,$(Fig. 1) in our models.

Figure 2: The plot shows in black (red) the BH mass range, at different times,
produced by evolving heavy stars of $25-100$ $\mathrm{M_{\odot}}\,$in
isolation using bse-psn (bse). While bse-psn produces BHs of $\sim$ 40
$\mathrm{M_{\odot}}$, already after $4$ Myr, bse starts to generate more
massive BHs only after 6 Myr. Therefore, in the bse-psn (contrary to bse)
simulations, 40 $\mathrm{M_{\odot}}\,$BHs can coexist with massive stars of
$\sim 50-70$ $\mathrm{M_{\odot}}$. The BHs have a higher probability to
collide with more massive stars and grow more rapidly in mass.

The new stellar evolution model also affects the BH formation time. The black
bars in Fig. 2 indicate the stellar BH mass ranges as a function of time. To
estimate these mass ranges we evolved massive stars in isolation with bse-psn
setting $Z=0.0002$. Already after $4$ Myr BHs with $35-45$
$\mathrm{M_{\odot}}\,$form as remnants of the most massive stars ($\gtrsim 80$
$\mathrm{M_{\odot}}$). Lighter stars explode later, pushing the lower BH mass
limit down to $\sim 9$ $\mathrm{M_{\odot}}\,$after $10$ Myr. On the contrary,
with the bse prescriptions, heavy stars in the range between $80$
$\mathrm{M_{\odot}}\,$and $100$ $\mathrm{M_{\odot}}\,$lose $\sim 90\%$ of
their mass leading to early BHs with much lower masses of $\sim 9-12$
$\mathrm{M_{\odot}}\,$. With time the upper mass limit increases only to $\sim
25$ $\mathrm{M_{\odot}}$(red bars in Fig. 2). As we will see in Section 4,
this difference between the old and the new stellar evolution models impacts
the overall formation path of massive BHs: in the bse-psn simulations, stellar
BHs above $35$ $\mathrm{M_{\odot}}\,$appear in the system when $\sim 70$
$\mathrm{M_{\odot}}\,$stars are still alive. Therefore they have a higher
probability to increase their mass by accretion of stellar material through
close stellar interactions.

There are also updates concerning the neutron star (NS) modelling. While in
bse almost all NSs were ejected from the cluster due to the high natal kick
velocities, in bse-psn the lower velocity kicks associated with ECSNe enable
the formation of a NS population after 10 Myr. However, for the short cluster
evolution times considered in this study, none of our models shows any NS-NS
or NS-BH merger.

Tidal capture events might play an important role in the growth of massive
BHs. A recent analytical model, Stone et al. (2017) show that if every tidal
capture event of a BH results in a collision, stellar-mass BHs, located in
dense stellar environments, can grow up to $10^{6}$
$\mathrm{M_{\odot}}\,$within $\sim 1$ Gyr. Our simulations include tidal
capture prescriptions that lead to rapid circularization rather than a
collision. The radius of the final circular orbit is set assuming angular
momentum conservation. Therefore, in general, tidal capture events do not lead
to collisions and might not contribute directly to the BHs mass growth.
Nevertheless, the circular BH-star binary, once formed, have good chance to
experience mass transfer or even mergers. Moreover, our instant tidal capture
prescription requires knowledge of the orbital elements (pericenter and
eccentricity) of the interacting objects. It is not activated during few-body
chaotic encounters, where pericenter and eccentricity are not well-defined
quantities. Therefore in some close interactions, we might underestimate the
total number of tidal capture events.

##  3 Initial Conditions

Figure 3: The formation paths of all six (A - F) BHs more massive than $100$
$\mathrm{M_{\odot}}\,$formed in the $f\mathrm{{}_{c}}=1.0\,$model. Main
sequence stars and red giants are indicated by blue and red stars,
respectively. Black holes are shown by filled black circles. The symbol sizes
indicate the respective masses. The collisions on the right of the grey
vertical segment would be certainly suppressed by relativistic kicks. Time
goes from left to right and we give the masses and the eccentricities at the
time of the collision. Paths C and D involve BH mergers before a GW190521-like
mass gap merger which could have been prevented by recoil kicks. Path B offers
a GW190521-like formation path only involving stellar accretion. Paths A, E,
and F result in the formation of massive BHs in the IMBH range without stellar
black hole mergers. For all examples, except C, we show the full formation
path. The full path of C up to 100 Myr is shown in Fig. 4. Figure 4: The
formation path of the most massive BH formed in the 16 simulations. The final
IMBH of 293 $\mathrm{M_{\odot}}\,$has formed through a chain of eight
collisions. Panel C in Fig. 3 shows only the first 4. 6 collisions involve
only BHs. This IMBH formed only because we did not include relativistic recoil
kicks for BH - BH mergers. With kicks the mergers after the grey vertical line
would not have happened as one of the progenitors had left the system already.
Figure 5: The formation paths of the five (A-E) massive black holes in the
$f\mathrm{{}_{c}}=0.5\,$model. The notation is the same as in Fig. 3. While
this model also produces massive black holes, none of the formation paths
involves a merger of two mass-gap BHs, like GW190521.

Star clusters are expected to form from massive turbulent molecular clouds
where gas is rapidly converted into stars. The high level of turbulence in
these regions generates clumpy substructures in the whole star forming cloud.
Numerical simulations indicate that violent relaxation erases initial
substructure on a dynamical time scale (Scally & Clarke, 2002; Hurley & Bekki,
2008) and leaves a spherical density profile compatible with an Elson, Fall
and Freeman (Elson et al., 1987) or a King (King, 1966) density profile.
Moreover, in the early stages, gas and newly formed stars coexist in the same
system. For moderately dense environments, only a small fraction of the cloud
is converted into stars, the remaining gas is later ejected by stellar
feedback and supernovae explosions. Although star clusters at formation are
very compact, since they typically have sizes $\lesssim 1$ pc (Marks & Kroupa,
2012; Fujii et al., 2021), gas ejection can lead to a rapid, strong expansion
that can increase the size of the system to $\sim 5-10$ pc. This is especially
true for low-mass star clusters ($10^{3}-10^{4}$ $\mathrm{M_{\odot}}$). These
types of systems are expected to lose their primordial gas in an interval of
time of about $\sim 9$ Myr (Pelupessy & Portegies Zwart, 2012; Fujii et al.,
2021).

On the other hand, observations show that young massive clusters (YMCs) such
as NGC 3603, Trumpler 14, and Arches are devoid of dense gas despite being
about $\sim 2$ Myr old (sse Longmore et al., 2014, and citations therein). All
these systems are very compact, with a photo-metric mass $>10^{4}$
$\mathrm{M_{\odot}}\,$and and radii respectively equal to $0.7,0.5,0.4$ pc
(Portegies Zwart et al., 2010). They are likely to originate from very dense
gaseous environments, which thanks to their high density can maintain a high
stellar formation efficiency, rapidly converting most of the gas into stars.
In other words, these YMCs remove a big fraction of their primordial gas
through exhaustion rather than ejection (Longmore et al., 2014). Therefore
they can retain their initial compactness. The stellar formation phase for
these systems has a very short duration, leaving a star cluster with a very
narrow age spread (see e.g. Lahén et al., 2019; Lahén et al., 2020, for the
simulated formation of such systems in a galactic environment).

A recent three-dimensional radiation hydrodynamical simulations of YMCs
formation support the picture just presented, indicating that stellar feedback
suppresses star formation only in clouds with low surface density; When the
initial surface density exceeds a certain threshold most of the gas is
converted into stars (Fukushima & Yajima, 2021). In this last case, the star
formation phase can last as little as 1 Myr.

Therefore, compact YMCs are ideal systems to be followed with direct N-body
simulations, in spite of the fact that they are little affected by processes
associated with gas ejections. We adopted the same initial conditions of the
eight realizations of the model R06W6 presented in Rizzuto et al. (2021). The
model represents the state of young compact star clusters of $110.000$ stars
with $10\%$ in binaries after gas removal. Their initial densities follow a
King profile with initial half mass radius
$R\mathrm{{}_{h}}=0.6\,\mathrm{pc}\,$and central potential parameter
$W\mathrm{{}_{0}}\,$222The central potential parameter is defined as
$W\mathrm{{}_{0}}\,$$=\frac{\psi(0)}{\sigma^{2}}$ where $\psi(0)$ is the
potential at the centre of the cluster and $\sigma^{2}$is a parameter linked
to the velocity dispersion of the system. See (King, 1966) for more details.
$=6$. Again, we assume that the systems maintain their initial compactness as
the initial gas is removed by exhaustion from star formation. In addition, we
assume the star formation phase is short-lived (see e.g Lahén et al., 2019).
Hence we neglect the initial age spread and we initialise zero-age main-
sequence stars with a Kroupa initial mass function (Kroupa, 2001). All the
systems started with a primordial binary fraction of $10\%$, initialised with
a uniform semi-major axis distribution on a logarithmic scale from $0.001$ AU
to $100$ AU, a uniform distribution of mass ratios, and thermal distribution
of eccentricities.

In the simulations we allow BHs to accrete mass from a star in case of a
collision. The fraction of stellar mass absorbed by a BH during such a star -
BH collision, which we define as collision factor $f\mathrm{{}_{c}}\,$, has a
significant impact on the formation of IMBHs (see Rizzuto et al., 2021). This
mass fraction, in general, is difficult to estimate, as it depends on the
relative velocity of the two colliding objects as well as the internal
structure of the stars, etc. A quantitative assessment of this number requires
complex hydrodynamical simulations. In this work, we treat this fraction
$f\mathrm{{}_{c}}\,$as a free parameter. We run $8$ simulations setting
$f\mathrm{{}_{c}}=0.5\,$, and the other $8$ with $f\mathrm{{}_{c}}=1.0\,$.

##  4 Results

In Fig. 3, Fig. 4, and Fig. 5 we show the formation paths of the most massive
BHs in the simulations for $f\mathrm{{}_{c}}=1.0\,$and
$f\mathrm{{}_{c}}=0.5\,$, respectively. Overall, these simulations with
updated binary stellar evolution models confirm our previous result (Rizzuto
et al., 2021) that massive black holes in the mass-gap and even IMBHs more
massive than $100$ $\mathrm{M_{\odot}}\,$can form in compact star clusters
through multiple mergers. With the new stellar prescriptions, the runs
registered an even higher probability to form massive BHs. The latter build up
their mass through at least one of these three channels:

  1. 1.

Star - BH mergers: stellar BHs increase their mass by collisions with one ore
more massive stars (see for instance panels A and B of Fig. 3 and panel D of
Fig. 5).

  2. 2.

BH - BH mergers: BHs grow by mergers with other BHs (i.e. panel C of Fig. 3
and panel C of Fig. 5).

  3. 3.

The collapse of stars that are the product of previous collisions: a VMS
formed through a series of subsequent collisions could acquire a large
hydrogen-rich envelope maintaining a small helium core ($M_{c}<45$
$\mathrm{M_{\odot}}$), it is therefore little or not affected by (P)PSN and
can collapse directly into an IMBH (i.e. panel F of Fig. 3 and panel A of Fig.
5.

This last formation channel was also observed in a recent set of N-body
simulations by Di Carlo et al. (2021), who also assume no mass loss in stellar
merger events. The importance of stellar mergers, however, depends on the
prescription used to model collisions between stars. Our assumption of no mass
loss might overestimate the impact of stellar mergers. For example, Glebbeek
et al. (2009) use a detailed stellar evolution code and take into account the
internal structure of massive stars formed through multiple collisions. They
indicate that steller mergers likely lose most of their mass in the form of
stellar winds. Moreover, our method does not account for remnant rotation.
Even though the consequences of rotation have not yet been completely explored
(Burrows & Vartanyan, 2021), high stellar spin is expected to enhance the
mass-loss rate (Maeder & Meynet, 2000).

In general, at least two of the three processes mentioned above are always
involved in the formation of a massive BH. For instance, the IMBH shown in
panel C of Fig. 3 forms through the first and the second channel; while panel
E of Fig. 5 presents an IMBH generated via the second and the third channel.
The most massive BH in all simulations has a mass of $293$
$\mathrm{M_{\odot}}$. This object forms through $8$ collisions $6$ of them are
BH - BH collisions as shown in Fig. 4. This long chain of BH mergers happens
when, as in our simulations, relativistic recoil kicks are not considered.
When two BHs collide, the relativistic recoil kick velocity depends on the
mass ratio as well as the spin amplitude and direction of the two colliding
objects. It can reach a few thousand kilometres per second (Campanelli et al.,
2007; Baker et al., 2008; Lousto & Zlochower, 2009; Kulier et al., 2015;
Morawski et al., 2018; Zivancev et al., 2020). Since the typical escape
velocity in our clusters is 20-40km/s, it seems extremely unlikely that the
$293$ $\mathrm{M_{\odot}}\,$IMBH with this formation path can form in star
cluster environments.

The first four collisions on the chain of mergers that lead to the $293$
$\mathrm{M_{\odot}}\,$are reported in panel (C) of Fig. 3. Overall this panel
shows two BH - star mergers and two BH - BH collisions that in the end
generate a $146$ $\mathrm{M_{\odot}}\,$IMBH. Since we expect that most
stellar-mass BHs are born with a low initial spin, as indicated partially by
the current GW data (Abbott et al., 2019; Miller et al., 2020; Roulet et al.,
2020), and suggested by recent theoretical models (Fuller & Ma, 2019), the
$146$ $\mathrm{M_{\odot}}\,$IMBH has a non-negligible probability to form even
with the inclusion of relativistic recoil as we shall discuss below.

Therefore, the $f\mathrm{{}_{c}}=1.0\,$simulations provide two possible
formation paths for GW190521. Path C is a second-generation event and has a
low probability due to a first generation BH merger, as discussed above. Path
B is more likely as the BH merger event itself is a first generation BH
merger.

###  4.1 The impact of stellar evolution updates

Figure 6: Left panel: Mass distribution of all BHs formed in the bse model
(red) and the bse-psn model (black) with $f\mathrm{{}_{c}}=1.0\,$after $100$
Myr of evolution. Simulations with bse-psn form more and more massive BHs in
the (P)PSN mass gap and the IMBH mass regime. Right panel: Comparison of the
BH mass distribution of the bse-psn model with $f\mathrm{{}_{c}}=1.0\,$(black,
same as left panel) to the model with $f\mathrm{{}_{c}}=0.5\,$(blue). Lower
accretion fractions for star-BH collisions result in less massive BHs. Figure
7: The two left (right) panels show the BH (massive stars) mass distribution
after 5 Myr (top) and 10 Myr (bottom) of cluster evolution. In the bse-psn
simulations (black histograms) $40$ $\mathrm{M_{\odot}}\,$BHs form early in
the simulation while massive stars $\sim 60$ $\mathrm{M_{\odot}}\,$are still
present (top panels). In contrast, the bse clusters contains only $\sim 10$
$\mathrm{M_{\odot}}\,$BHs after $5$ Myr (red histograms). After 10 Myr the
bse-psn model formed more massive BHs (bottom left panel) with a comparable
stellar mass distribution (bottom right panel).

After 100 Myr of cluster evolution, the BH mass distribution of the updated
bse-psn model is very similar to the previous results as shown in the left
panel of Fig. 6). The typical formation paths for the most massive BHS,
however, differ significantly. In the Rizzuto et al. (2021) simulations (red
histogram of Fig. 6, left panel) the BHs in the IMBH mass regime manly through
mergers of BHs with VMSs: massive star experiences multiple stellar collisions
reaching masses of $200-400$ $\mathrm{M_{\odot}}$. A large fraction of this is
then accreted by a stellar BH in a collision. Typically, it is a single
collision between a stellar BH and a VMS that determines whether or not the
cluster forms an IMBH. While this formation path can also occur in bse-psn
(see e.g. panel D of Fig. 5) it is not the dominant one.

In the bse-psn simulations, stellar BHs, before reaching the IMBH mass range,
typically experience multiple collisions with massive stars and other BHs. In
several cases, IMBHs form entirely, or almost entirely, via repeated BH
mergers, as illustrated in panels B and C of Fig. 5 and panel D of Fig. 3.
These formation paths are extremely rare in bse; only 1 out of 80 simulations
formed an IMBH through multiple BH - BH collisions.

The origin of this difference is the collision rate. The updated bse-psn model
has about four times higher collision rate between massive stars and stellar
BHs and three times higher BH-BH collision rate than bse , and BH - BH. To
understand this difference we briefly review the early cluster evolution.
Within the first 4 Myr, the clusters undergo core-collapse driven by mass
segregation of the most massive objects. Right after this initial collapse
fast expansion is triggered by the mass loss of massive stars. For bse-psn
massive stellar BHs form rapidly (see Fig. 2 and the top left panel of Fig. 7)
and already populate the inner part of the cluster while the central density
is at its maximum and when massive stars with a large collision cross-sections
are still alive (top right panel of Fig. 7). High collision rates for star-BH
and BH-BH mergers are the results.

This is not the case for bse simulations. The first BHs that appear in the
core-collapse phase of the bse simulations have a mass of only $\sim 10$
$\mathrm{M_{\odot}}\,$(see top left panel of Fig. 7) and they need a
relatively long time to reach the inner part of the system, they typically
reach the centre. At this time the clusters are not very dense anymore and
$\gtrsim 40$ $\mathrm{M_{\odot}}\,$stars are already gone. More massive
stellar BHs ($\sim 30$ $\mathrm{M_{\odot}}$) start to populate the bse
simulations only after the most massive stars have already exploded (bottom
panels of Fig. 7).

In contrast, $\sim 40$ $\mathrm{M_{\odot}}\,$BHs populate the bse-psn clusters
while massive stars are still alive (top panels of Fig. 7). Due to their short
segregation time, these BHs can reach the core of the cluster, while it is
still dense, and collide with one or more massive stars. At the same time
while sinking into the centre they support and enhance the core collapse. To
summarize, in the bse clusters, since the first stellar BHs have low mass (10
$\mathrm{M_{\odot}}$), they have little chance to reach the centre when it is
still very dense; Instead, in bse-psn the first stellar BHs that appear in the
system are more massive (40 $\mathrm{M_{\odot}}$), therefore, they can reach
the centre while it is still compact and therefore they have a higher
probability to collide with other objects.

Overall, the simulations with bse-psn seem to generate a larger number of
massive BHs. The left panel of Fig. 6 shows that bse-psn formed about 10 BHs
more massive than 75 $\mathrm{M_{\odot}}$, while bse simulation only form 4
BHs. This difference is also connected with the fact that bse-psn can form
massive BHs, as we mentioned previously, through direct collisions of VMSs.
The direct collapse channel cannot occur in bse because the stellar winds
prescriptions are extreme and all VMSs lose most of their material before the
collapse.

###  4.2 Comparison between $f\mathrm{{}_{c}}=0.5\,$and
$f\mathrm{{}_{c}}=1.0\,$

Figure 8: The mass of the three most massive BHs in each of the 8 realisations
of the cluster simulations with $f\mathrm{{}_{c}}=1.0\,$(left panel) and
$f\mathrm{{}_{c}}=0.5\,$(right panel). The first massive black hole is
typically heavier in the $f\mathrm{{}_{c}}=1.0\,$simulations, while the second
and third have comparable masses to $f\mathrm{{}_{c}}=0.5\,$models.  Figure 9:
The Fig. shows in what percentage BH - BH mergers, BH - star coalescence, and
direct collapse of stars (produced by previous collisions) contribute to the
re-population of the (P)PSN mass gap in our simulations.

Fig. 8 compares the most massive BHs formed in the
$f\mathrm{{}_{c}}=1.0\,$simulations with the
$f\mathrm{{}_{c}}=0.5\,$simulations. As expected the former is more likely to
generate BHs with higher masses. $f\mathrm{{}_{c}}=1.0\,$models appear to be
two times more likely to form IMBHs since already a single star - BH collision
is enough to generate a BH heavier than $100$ $\mathrm{M_{\odot}}\,$(as shown
in panel (A) of Fig. 3). This can also be seen in the right panel of Fig. 6
showing the BH mass distributions generated by the two models: below $\sim 40$
$\mathrm{M_{\odot}}\,$the two distributions are very similar. However, only
the $f\mathrm{{}_{c}}=1.0\,$model results in the formation of IMBHs above 150.
The massive BHs of $f\mathrm{{}_{c}}=0.5\,$model are lower masses in the range
between $60$ and $110$ $\mathrm{M_{\odot}}$.

As we mentioned earlier, because of the absence of observational and
theoretical constraints, we choose $f\mathrm{{}_{c}}\,$to be a free parameter.
The constraints on $f\mathrm{{}_{c}}\,$presented in the literature refer to
tidal disruption events between a low mass star and a supermassive BH. They
estimate an $f\mathrm{{}_{c}}\,$in the range between $\sim 0.2-0.5$ (Shiokawa
et al., 2015; Lu & Bonnerot, 2020; Wen et al., 2020) because most of the mass
is lost in high-speed jets formed during mass accretion. However, when a low
mass BH collides with a massive star, the scenario might be very different.
Massive stars form very rapidly a compact core, so it might not be surprising
that the actual core of the star is entirely absorbed by the BH right after
the collision. Later on, the resulting BH, immersed in the envelope of the
star, might continue to absorb stellar material. Due to the large size of the
envelope of a massive star, the high-speed jets might be contained by the
envelope itself, keeping most of the gas in the vicinity of the black hole.
This situation resemble a similar scenario as described in Safarzadeh & Haiman
(2020) and Rice & Zhang (2021). Under these assumptions,
$f\mathrm{{}_{c}}\,$might be close to unity.

The charts in Fig. 9 visualise the formation paths of the BH more massive
lower limit of the (P)PSN mass gap ($>45$ $\mathrm{M_{\odot}}$) for the
$f\mathrm{{}_{c}}=1.0\,$and $f\mathrm{{}_{c}}=0.5\,$models. Most stellar BHs
reach the mass gap by a merger with another BH. The next most important path
is a collision with a star, followed by the direct collapse of a stellar
merger remnant. The relative importance of the formation paths does not change
significantly with the assumed accretion fraction - even setting
$f\mathrm{{}_{c}}=0.5\,$, the star - BH collision channel produced $11$ BHs in
the mass gap just $2$ less than the one produced with
$f\mathrm{{}_{c}}=1.0\,$.

BH merger |  $V_{\mathrm{esc}}$ (km/s) | P1 | P2 | P3
---|---|---|---|---
Fig. 3 panel B Gen-a | 38.1 | 3.05 % | 82.2 % | 100.0 %
Fig. 3 panel C Gen-a | 40.3 | 1.91 % | 0.00 % | 0.00 %
Fig. 3 panel C Gen-b | 38.6 | 0.00 % | 0.00 % | 0.00 %
Fig. 3 panel D Gen-a | 29.2 | 0.00 % | 0.00 % | 0.00 %
Fig. 3 panel D Gen-b | 25.8 | 0.00 % | 0.00 % | 0.00 %
Fig. 5 panel B Gen-a | 31.2 | 0.00 % | 0.00 % | 0.00 %
Fig. 5 panel B Gen-b | 27.1 | 0.00 % | 0.00 % | 0.00 %
Fig. 5 panel C Gen-a | 35.7 | 4.44 % | 99.0 % | 100.0 %
Fig. 5 panel E Gen-a | 41.5 | 0.00 % | 0.00 % | 0.00 %
Table 2: The table reports, for each BH - BH collision illustrated in Fig. 3
and 5, the retention probability of the BH merger remnant in the star cluster.
$V_{\mathrm{esc}}$: cluster escape velocity of the BH merger remnant. P1:
retention probability for an initial Gaussian spin distribution peaked on
$S_{1}=0.2$, with a $\sigma=0.2$ for the stellar BHs. P2: retention
probability for $S_{2}=10^{-2}$. P3: retention probability for
$S_{3}=10^{-3}$.

###  4.3 Comparison with LIGO/Virgo/Kagra gravitational wave detections

Figure 10: The panels show the primary ($m_{1}$) and secondary ($m_{2}$) BH
masses of all BH mergers in the simulations for an accretion fraction of
$f\mathrm{{}_{c}}=0.5\,$(left, grey circles) and
$f\mathrm{{}_{c}}=1.0\,$(right, blue circles). The letters indicate the
simulated formation paths highlighted in Fig. 3 and 5. The currently available
LIGO/Virgo/Karga gravitational wave detections including error bars are
indicated in orange. BH merger events that might be excluded due to
gravitational recoil kicks are indicated with open circles. In general, the
simulated events cover a similar parameter space as all currently available
observations. The $f\mathrm{{}_{c}}=1.0\,$simulations provide two possible
formation paths for GW190521. Path C is a second-generation event and has a
low probability due to a first generation BH merger. Path B is more likely as
the event itself is a first generation BH merger.

Our runs generated a large variety of BH mergers. The
$f\mathrm{{}_{c}}=1.0\,$simulations, for example, resulted in 27 BH - BH
collisions (just $2$ more than the $f\mathrm{{}_{c}}=0.5\,$runs). These BH
merger events cover almost the entire mass range of the LIGO/Virgo/Karga
detections as shown in Fig. 10. Di Carlo et al. (2020) presented comparable
results for a larger set of simulations of lower mass metal-poor star
clusters. Interestingly, two-star cluster realisations with
$f\mathrm{{}_{c}}=1.0\,$have produced events with primary and secondary masses
very similar to the GW190521 mass-gap detection. These BH merger events are
highlighted in panels (B) and (C) of Fig. 3. Path B is particularly
interesting as it indicates that BH mergers in the GW190521 mass range could
originate from BHs which have grown by BH - star collisions. This formation
path is not influenced by BH merger kicks. Therefore GW190521-like events
might naturally occur in compact star clusters and they do not necessarily
require very massive environments such as nuclear star clusters. The only
requirements are that the environment is compact enough to trigger several BH
- star mergers and that a large fraction of the stellar mass is accreted onto
the BH during the collision.

As we have mentioned above, the simulations presented here do not include
relativistic velocity recoil kicks for compact object merger remnants.
Therefore all massive BHs formed through two or more BH - BH collisions would
have a much lower formation probability in the presence of kicks. The
probability to retain a compact merger remnant in the cluster depends on its
escape velocity as well as the spin and the mass ratio of the colliding BHs.
Following Arca Sedda et al. (2021) we compute the kick velocity distributions
for different assumed natal spins and quantify the retention probability of
the BH merger remnants. We test three different initial spin distributions for
stellar-mass BH, centered around $S_{1}=0.1$ with a Gaussian standard
deviation of $\sigma=0.1$, and fixed spins of $S_{2}=0.01$ and $S_{3}=0.001$,
respectively.

We compute the retention probability for each BH - BH collision reported in
Fig. 3 and 5 for all three spin distribution and present the results in Tab.
2.

Under the assumption that stellar BHs are born with very low spins ($S_{2}$ or
$S_{3}$, Fuller & Ma, 2019), only the events (B) in Fig. 3 and (C) in Fig. 5
would survive the relativistic kick. In these cases, since the mass ratio is
close to unity, the recoil velocity is lower than the escape velocity. For all
the other cases $V_{\text{recoil}}\gtrsim 50$ km/s therefore merger remnants
have little to no chance to stay bound to the cluster.

If we assume the stellar BHs to be slowly rotating with distribution $S_{1}$
(as suggested by observations, e.g. Abbott et al., 2019), some BH - BH mergers
have a non-negligible probability to stay in the star cluster for certain spin
orientations. Under this assumption the first BH collision reported in the
panel (C) of Fig. 3 has a probability of 1.9 % to stay in the cluster because
of comparable masses. Therefore, similar events, which in turn resemble
GW190521, could have been hosted in compact but not very massive star clusters
as long as the BHs involved in the collisions are born with a low initial spin
similar to $S_{2}$.

Another interesting event, that has non negligible probability to occur (see
Tab. 2), is illustrated in panel (C) of Fig. 5. The plot shows the merger of a
57 $\mathrm{M_{\odot}}\,$and a 41 $\mathrm{M_{\odot}}\,$BH. Due to the similar
masses the BHs this collision might indicate a possible formation channel for
detected events like GW190701_203306 ($m_{1}\approx 53.6$
$\mathrm{M_{\odot}}$, $m2\approx 40.7$ $\mathrm{M_{\odot}}$ Abbott et al.,
2020).

To estimate the gravitational recoil velocity we implicitly assume that the
spin of a BH is not affected during a collision with a massive star. While the
BH spins of several low mass X-ray binaries are likely to be influenced by
mass transfer from the companion into the BH (Podsiadlowski et al., 2003;
Fragos & McClintock, 2015; Sørensen et al., 2017), all the three BH spins
measured in high-mass X-ray binaries appear to have been affected very little
by mass accretion (Qin et al., 2019). However, even allowing for larger BH
spins for star-BH collisions our results would not change.

Two IMBH binaries form at the centre of two different realizations. The first
one, a ($173$ $\mathrm{M_{\odot}}\,$, $153$ $\mathrm{M_{\odot}}\,$) BH binary,
whose formation paths are shown in panels (E) and (F) of Fig. 3, forms in the
second realisation of $f\mathrm{{}_{c}}=1.0\,$(see Fig.8). Right after
formation, the two IMBHs evolve into a hard binary. At the end of the
simulation ($t\approx 100$ Myr) this binary has an eccentricity of $0.81$ and
a semi-major axis of 1.88 AU, therefore, if evolved in isolation, will merge
in about 11.5 Gyr. Both IMBHs in this binary formed from the direct collapse
of a VMS. According to the stellar model we adopted, a star-forming from
multiple stellar mergers might retain a small helium core. This in turn leads
to a weaker (P)PSN allowing the star to collapse into a massive BH or even
into an IMBH. We note that this result is connected to significant
uncertainties of stellar evolution and we do not include the effect of
rotation in VMSs, which could increase stellar winds mass-loss rate (Maeder &
Meynet, 2000) and it could affect the final stage of the massive star Burrows
& Vartanyan (2021). The second IMBH binary forms in the sixth realization of
model $f\mathrm{{}_{c}}=0.5\,$(see right panel of Fig. 8) and it involves two
black holes of $113$ and $99$ $\mathrm{M_{\odot}}\,$(see panels (A) and (E) of
Fig. 5). However, since the $99$ $\mathrm{M_{\odot}}\,$object is the product
of a $70$ and $29$ $\mathrm{M_{\odot}}\,$BH merger, it is very likely to
escape the cluster right after formation (see Tab. 2). Even if it stays in the
cluster this IMBH binary would have an expected merging time in isolation
greater than the Hubble time (at the end of the run the binary has
eccentricity equal to 0.75 and semi-major axes of 3.25 AU). Nevertheless, the
merging time could be considerably reduced by chaotic and hierarchical
interactions.

##  5 Conclusions

We study a collisional formation scenario for IMBHs and BHs in the (P)PSN mass
gap in young compact star clusters utilizing direct N-body simulations carried
out with a version of nbody6++gpu that contains relativistic treatments to
include the GW energy loss in compact binaries. We rerun $16$ realizations of
the model R06W6 introduced in Rizzuto et al. (2021) with updated stellar
evolution recipes. The upgrades follow the implementations in Banerjee et al.
(2020) and include new metallicity-dependent stellar winds prescriptions,
recipes for the electron-capture supernovae, and treatments for (P)PSN. The
latter creates a mass gap in the BH mass spectrum between $45$ and $195$
$\mathrm{M_{\odot}}\,$(see Fig. 1).

Overall our simulations show that $\sim 10^{5}$ $\mathrm{M_{\odot}}\,$compact
star clusters have a good chance to generate BHs in the (P)PSN mass gap or
even in the IMBH mass range ($\gtrsim 100$) as shown in Fig. 6. Together with
BH - BH collisions, BH - star collisions seem to play an important role in the
formation of massive BHs. We parameterize the fraction of stellar mass
absorbed by the BH to be unity ($f\mathrm{{}_{c}}=1.0\,$) or 50%
($f\mathrm{{}_{c}}=0.5\,$). The Simulations with $f\mathrm{{}_{c}}=1.0\,$and
$f\mathrm{{}_{c}}=0.5\,$can generate BH merger events with primary and
secondary BH masses compatible with nearly the whole mass range of
LIGO/Virgo/Kagra gravitational wave detections (see Fig. 10). As we do not
follow gravitational recoil, we have, a posterior, estimated, the retention
probability for a BH - BH merger in the star cluster and have identified all
unlikely formation paths.

The model assuming efficient accretion ($f\mathrm{{}_{c}}=1.0\,$) of stellar
material in a merger with a BH generates two BH merger events with primary and
secondary masses compatible with the mass-gap merger event GW190521 (Abbott et
al., 2020). For one event, both colliding BHs grow entirely through BH - star
mergers (see panel (B) Fig. 3). Therefore, the final merger is independent of
gravitational recoil kicks and might not only be expected in high escape
velocity systems such as nuclear star clusters. GW190521 like mergers might
also form in lower mass compact star clusters at low initial metallicity. The
only requirement is that the system is compact enough to trigger many BH -
star collisions and that the amount of mass accreted by the BH during the
collision is close to unity.

The chain of collisions that led to the second GW event involves two BH - star
collisions and one BH - BH collision before the final merger (see panel (C) of
Fig. 3). When these two compact objects merge, the final product should
experience a recoil kick which depends on the spin orientation and magnitude
of the colliding objects. For low spins, $\sim 0.1$ the probability that the
product of this BH - BH collision stays in the cluster is low ($1.9\%$) but
non-zero, therefore this formation path is not excluded.

Contrary to the $f\mathrm{{}_{c}}=1.0\,$model, the $8$ realizations of the
model $f\mathrm{{}_{c}}=0.5\,$did not form BH mergers in the mass range of
GW190521 because $f\mathrm{{}_{c}}=0.5\,$simulations tend to form less massive
progenitor BHs (see Fig. 8, and right panel of Fig. 6). However, the total
number of BHs in the (P)PSN mass gap generated by the two models through BH -
star collision is comparable (see Fig. 9). Even imposing
$f\mathrm{{}_{c}}=0.5\,$, BH - star mergers can bring stellar BHs in the
(P)PSN mass range.

In one simulation we even form a hard IMBH binary of
153$\mathrm{M_{\odot}}\,$and 173 $\mathrm{M_{\odot}}\,$(see (E) and (F) of
Fig. 3) which is expected to merge within a Hubble time. Its formation,
however, depends on details of the evolution of very massive stars which are
only approximately captured by the stellar evolution models used in this
study.

## Acknowledgements

We thank Albrecht Kamlah, Peter Berczik, and Jarrod Hurley for their
cooperation by incorporating the new stellar evolution features into
NBODY6++GPU. TN acknowledges support from the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) under Germany’s Excellence Strategy -
EXC-2094 - 390783311 from the DFG Cluster of Excellence "ORIGINS". The authors
gratefully acknowledge the Gauss Centre for Supercomputing (GSC) e.V.
(www.gauss-centre.eu) for funding this project by providing computing time
through the John von Neumann Institute for Computing (NIC) on the GCS
Supercomputers JUWELS and JUWELS-BOOSTER at Jülich Supercomputing Centre
(JSC). MG was partially supported by the Polish National Science Center (NCN)
through the grant UMO-2016/23/B/ST9/02732. SB acknowledges support from the
Deutsche Forschungsgemeinschaft (DFG; German Research Foundation) through the
individual research grant "The dynamics of stellar-mass black holes in dense
stellar systems and their role in gravitational-wave generation" (BA 4281/6-1;
PI: S. Banerjee)

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2019) Abbott B. P., et al., 2019, ApJ, 882, L24
  * Abbott et al. (2020) Abbott R., et al., 2020, Phys. Rev. Lett., 125, 101102
  * Arca Sedda & Benacquista (2019) Arca Sedda M., Benacquista M., 2019, MNRAS, 482, 2991
  * Arca Sedda et al. (2021) Arca Sedda M., Rizzuto F. P., Naab T., Ostriker J., Giersz M., Spurzem R., 2021, arXiv e-prints, p. arXiv:2105.07003
  * Baker et al. (2008) Baker J. G., Boggs W. D., Centrella J., Kelly B. J., McWilliams S. T., Miller M. C., van Meter J. R., 2008, ApJ, 682, L29
  * Banerjee (2021) Banerjee S., 2021, MNRAS, 500, 3002
  * 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
  * Becker et al. (2020) Becker D., Spera M., Kremer K., Rasio F., 2020, in American Astronomical Society Meeting Abstracts #235. p. 369.10
  * Belczynski (2020) Belczynski K., 2020, ApJ, 905, L15
  * 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., et al., 2016, A&A, 594, A97
  * Burrows & Vartanyan (2021) Burrows A., Vartanyan D., 2021, Nature, 589, 29
  * Campanelli et al. (2007) Campanelli M., Lousto C. O., Zlochower Y., Merritt D., 2007, Physical Review Letters, 98, 231102
  * Dall’Amico et al. (2021) Dall’Amico M., Mapelli M., Niccolò Di Carlo U., Bouffanais Y., Rastello S., Santoliquido F., Ballone A., Arca Sedda M., 2021, arXiv e-prints, p. arXiv:2105.12757
  * 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. (2020) Di Carlo U. N., et al., 2020, MNRAS, 498, 495
  * Di Carlo et al. (2021) Di Carlo U. N., et al., 2021, arXiv e-prints, p. arXiv:2105.01085
  * Elson et al. (1987) Elson R. A. W., Fall S. M., Freeman K. C., 1987, ApJ, 323, 54
  * Farmer et al. (2019) Farmer R., Renzo M., de Mink S. E., Marchant P., Justham S., 2019, ApJ, 887, 53
  * Fowler & Hoyle (1964) Fowler W. A., Hoyle F., 1964, ApJS, 9, 201
  * Fragione et al. (2020) Fragione G., Loeb A., Rasio F. A., 2020, ApJ, 902, L26
  * Fragos & McClintock (2015) Fragos T., McClintock J. E., 2015, ApJ, 800, 17
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91
  * Fujii et al. (2021) Fujii M. S., Saitoh T. R., Hirai Y., Wang L., 2021, arXiv e-prints, p. arXiv:2103.02829
  * Fukushima & Yajima (2021) Fukushima H., Yajima H., 2021, arXiv e-prints, p. arXiv:2104.10892
  * Fuller & Ma (2019) Fuller J., Ma L., 2019, ApJ, 881, L1
  * Gayathri et al. (2020) Gayathri V., et al., 2020, arXiv e-prints, p. arXiv:2009.05461
  * Gessner & Janka (2018) Gessner A., Janka H.-T., 2018, ApJ, 865, 61
  * Giersz et al. (2015) Giersz M., Leigh N., Hypki A., Lützgendorf N., Askar A., 2015, MNRAS, 454, 3150
  * Glebbeek et al. (2009) Glebbeek E., Gaburov E., de Mink S. E., Pols O. R., Portegies Zwart S. F., 2009, A&A, 497, 255
  * Hurley & Bekki (2008) Hurley J. R., Bekki K., 2008, MNRAS, 389, L61
  * Kamlah et al. (2021) Kamlah A. W. H., et al., 2021, arXiv e-prints, p. arXiv:2105.08067
  * Kimball et al. (2020) Kimball C., et al., 2020, arXiv e-prints, p. arXiv:2011.05332
  * King (1966) King I. R., 1966, AJ, 71, 64
  * Kinugawa et al. (2021) Kinugawa T., Nakamura T., Nakano H., 2021, MNRAS, 501, L49
  * Kremer et al. (2020) Kremer K., et al., 2020, ApJ, 903, 45
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kulier et al. (2015) Kulier A., Ostriker J. P., Natarajan P., Lackner C. N., Cen R., 2015, ApJ, 799, 178
  * Lahén et al. (2019) Lahén N., Naab T., Johansson P. H., Elmegreen B., Hu C.-Y., Walch S., 2019, ApJ, 879, L18
  * Lahén et al. (2020) Lahén N., Naab T., Johansson P. H., Elmegreen B., Hu C.-Y., Walch S., Steinwandel U. P., Moster B. P., 2020, ApJ, 891, 2
  * Liu & Bromm (2020) Liu B., Bromm V., 2020, ApJ, 903, L40
  * Longmore et al. (2014) Longmore S. N., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 291 (arXiv:1401.4175), doi:10.2458/azu_uapress_9780816531240-ch013
  * Lousto & Zlochower (2009) Lousto C. O., Zlochower Y., 2009, Phys. Rev. D, 79, 064018
  * Lu & Bonnerot (2020) Lu W., Bonnerot C., 2020, MNRAS, 492, 686
  * Maeder & Meynet (2000) Maeder A., Meynet G., 2000, A&A, 361, 159
  * Mapelli et al. (2021) Mapelli M., et al., 2021, MNRAS, 505, 339
  * Marks & Kroupa (2012) Marks M., Kroupa P., 2012, A&A, 543, A8
  * Miller et al. (2020) Miller S., Callister T. A., Farr W. M., 2020, ApJ, 895, 128
  * Morawski et al. (2018) Morawski J., Giersz M., Askar A., Belczynski K., 2018, MNRAS, 481, 2168
  * Özel et al. (2010) Özel F., Psaltis D., Narayan R., McClintock J. E., 2010, ApJ, 725, 1918
  * Pelupessy & Portegies Zwart (2012) Pelupessy F. I., Portegies Zwart S., 2012, MNRAS, 420, 1503
  * Peters & Mathews (1963) Peters P. C., Mathews J., 1963, Physical Review, 131, 435
  * Podsiadlowski et al. (2003) Podsiadlowski P., Rappaport S., Han Z., 2003, MNRAS, 341, 385
  * Podsiadlowski et al. (2004) Podsiadlowski P., Langer N., Poelarends A. J. T., Rappaport S., Heger A., Pfahl E., 2004, ApJ, 612, 1044
  * Portegies Zwart et al. (2010) Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A, 48, 431
  * Qin et al. (2019) Qin Y., Marchant P., Fragos T., Meynet G., Kalogera V., 2019, ApJ, 870, L18
  * Rice & Zhang (2021) Rice J. R., Zhang B., 2021, ApJ, 908, 59
  * Rizzuto et al. (2021) Rizzuto F. P., et al., 2021, MNRAS, 501, 5257
  * Romero-Shaw et al. (2020) Romero-Shaw I., Lasky P. D., Thrane E., Calderón Bustillo J., 2020, ApJ, 903, L5
  * Roulet et al. (2020) Roulet J., Venumadhav T., Zackay B., Dai L., Zaldarriaga M., 2020, Phys. Rev. D, 102, 123022
  * Safarzadeh & Haiman (2020) Safarzadeh M., Haiman Z., 2020, ApJ, 903, L21
  * Scally & Clarke (2002) Scally A., Clarke C., 2002, MNRAS, 334, 156
  * Shiokawa et al. (2015) Shiokawa H., Krolik J. H., Cheng R. M., Piran T., Noble S. C., 2015, ApJ, 804, 85
  * Sørensen et al. (2017) Sørensen M., Fragos T., Steiner J. F., Antoniou V., Meynet G., Dosopoulou F., 2017, A&A, 597, A12
  * Stone et al. (2017) Stone N. C., Küpper A. H. W., Ostriker J. P., 2017, MNRAS, 467, 4180
  * Vink et al. (2001) Vink J. S., de Koter A., Lamers H. J. G. L. M., 2001, A&A, 369, 574
  * Vink et al. (2021) Vink J. S., Higgins E. R., Sander A. A. C., Sabhahit G. N., 2021, MNRAS, 504, 146
  * 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. (2016) Wang L., et al., 2016, MNRAS, 458, 1450
  * Wen et al. (2020) Wen S., Jonker P. G., Stone N. C., Zabludoff A. I., Psaltis D., 2020, ApJ, 897, 80
  * Woosley (2017) Woosley S. E., 2017, ApJ, 836, 244
  * Woosley et al. (2002) Woosley S. E., Heger A., Weaver T. A., 2002, Reviews of Modern Physics, 74, 1015
  * Woosley et al. (2007) Woosley S. E., Blinnikov S., Heger A., 2007, Nature, 450, 390
  * Woosley et al. (2020) Woosley S. E., Sukhbold T., Janka H. T., 2020, ApJ, 896, 56
  * Ziółkowski (2008) Ziółkowski J., 2008, Chinese Journal of Astronomy and Astrophysics Supplement, 8, 273
  * Zivancev et al. (2020) Zivancev C., Ostriker J., Kupper A. H. W., 2020, arXiv e-prints, p. arXiv:2004.06083


# Building Extractive Question Answering System to Support
Human-AI Health Coaching Model for Sleep Domain

Iva Bojic1 Qi Chwen Ong1 Shafiq Joty2,3 Josip Car1,4
1Centre for Population Health Sciences, Lee Kong Chian School of Medicine,
NTU, Singapore
1School of Computer Science and Engineering, NTU, Singapore
3Salesforce AI Research, USA
4Department of Primary Care and Public Health, School of Public Health,
Imperial College London, UK

##  1 Introduction

Non-communicable diseases are the leading cause of deaths globally, imposing a
huge burden on healthcare systems Organization et al. (2014). This leads to
the increasing emphasis on primary prevention, which seeks to avert the onset
of diseases by addressing their causes and risk factors, such as altering
unhealthy behaviors (Kisling and Das, 2021). Health coaching has emerged as an
effective strategy for primary prevention of chronic conditions by
facilitating lifestyle behavior change Yang et al. (2020). With assistance
from Question Answering (QA) system, health coaches could be empowered with
evidence-based information, streamlining the task of helping clients to adopt
new practices.

An innovative human-Artificial Intelligence (AI) health coaching model
incorporating QA systems holds the potential to revolutionize the health
coaching domain and the practice of preventive healthcare. This model involves
the interaction of three components: a health coach, a client, and a QA system
(see Fig. 1). The coaching process takes place online (e.g., through text
messaging). If needed, the coach can consult the QA system for assistance in
responding to factoid questions posed by the client Wang et al. (2006). In
this way, human coaching is supported by the evidence-based knowledge supplied
by the QA system. This not only accelerates the health coaches’ response time
but also ensures that the information provided by the health coach is grounded
in evidence.

Figure 1: Human-AI health coaching model.

This paper outlines the creation of a domain-specific extractive QA system
employed as the AI component within a Human-AI health coaching model.
Generally, an extractive QA system composes two phases: a passage retrieval
followed by a passage reader (Chen et al., 2017). For a given question, the
retrieval initially extracts a set of pertinent passages from a knowledge base
(e.g., a text corpus), after which the reader selects an answer (for instance,
a text span) from one of those passages (Zhu et al., 2021). Although
generative large language models like ChatGPT show great ability in QA, the
lack of references and hallucinations limit their utilizations in healthcare
at present Zhang et al. (2023). _Domain-specific_ systems typically address
questions confined to a single area of expertise, such as medicine (Kia et
al., 2022).

##  2 Methods

In the pursuit of developing a domain-specific extractive QA system to bolster
a Human-AI health coaching model, we manually assembled a dataset called
SleepQA, comprising 7,005 passages focused on the sleep domain Bojic et al.
(2022). This dataset featured 5,000 expert-annotated labels111Code and dataset
are available at:
https://github.com/IvaBojic/SleepQA. Subsequently, we fine-tuned various
domain-specific BERT models on our SleepQA dataset and conducted a
comprehensive assessment of the resulting end-to-end extractive QA system
using both automatic and human evaluation methods.

To further improve the QA system’s performance, we introduced a data-centric
framework aimed at refining domain-specific extractive QA datasets222Code and
datasets are available at:
https://github.com/IvaBojic/framework. Our proposed framework consists of
methods for (i) improving the selection of negative passages during the
retrieval fine-tuning, and (ii) question reformulation through paraphrasing,
word replacement, and back translation Bojic et al. (2023). Using the proposed
framework, we first produced new datasets by employing the four aforementioned
data quality enhancement methods. Subsequently, we fine-tuned both the
retrieval and reader models individually on each new dataset. Following this,
we fine-tuned the retrieval/reader models continuously, commencing from the
best individual checkpoint, using the enhanced datasets that demonstrated
progress in the initial step. Lastly, we generated new augmented datasets by
merging all datasets that exhibited fine-tuning enhancements in the initial
step.

##  3 Results

The initial fine-tuned QA system consisted of fine-tuned
PubMedBERT333https://huggingface.co/microsoft/BiomedNLP-PubMedBERT-base-
uncased-abstract-fulltext retrieval model and BioBERT
BioASQ444https://huggingface.co/gdario/biobert_bioasq reader model. Automatic
evaluation of the initial fine-tuned QA system on test labels demonstrated
that the proportion of accurate predictions, specifically those that precisely
align with the correct answers (known as an Exact Match or EM), reached 24%,
as opposed to the 30% achieved by the baseline QA system (see Table 1). The
baseline QA system refers to the combination of Lucene BM25 (Lin et al., 2021)
used as a baseline retrieval model and BERT SQuAD2 QA (Rajpurkar et al., 2018)
as a baseline reader model.

Table 1: Automatic evaluation of two extractive QA systems on test labels using EM Bojic et al. (2022). Extractive QA system name | EM
---|---
Lucene BM25 \+ BERT SQuAD2 QA | 0.30
PubMedBERT \+ BioBERT BioASQ | 0.24

Although our initial fine-tuned QA system did not outperform the baseline on
the test labels, results from human evaluation of the two QA systems on real-
world questions suggested otherwise. Real-world questions are sleep-related
questions posed by university students, which mimic the types of questions
that clients would ask their sleep health coaches in real life. The findings
presented in Bojic et al. (2022) showed that our initial fine-tuned QA system
outperformed the baseline QA system in almost 40% of cases in its intended
purpose, which is accurately providing sleep health coaches with the
appropriate responses to sleep-related inquiries they might receive from their
clients.

In order to increase the performance of our initial fine-tuned QA system on
the test labels as well, we introduced various data quality enhancement
methods. The best method to increase the performance of retrieval fine-tuning
was word substitution while fine-tuning of the reader model resulted in the
best results when using augmentation. The overall relative increase in
performance of the entire enhanced fine-tuned QA system, compared to the
initial one, was 17% (i.e. EM = 0.28). Although the relative increase was not
negligible, the enhanced fine-tuned QA system still did not outperform the
baseline QA system on the test labels.

##  4 Discussion and Conclusions

This paper outlines the creation of a domain-specific extractive QA system,
tailored for the sleep domain, utilizing a two-stage approach. First, we fine-
tuned our initial QA system using 5,000 experts-annotated labels from the
SleepQA dataset. Second, we enhanced the initial SleepQA dataset by employing
a data-centric approach and subsequently fine-tuned the QA system on the
enhanced dataset. The enhanced QA system was then integrated into a Human-AI
health coaching model.

Automatic evaluation showed that both fine-tuned QA systems did not surpass
the performance of the baseline QA system. Despite employing data-centric
methods to reformulate manually collected questions, it is likely that the
similarities remain between initial and reformulated questions. The
superiority of the Lucene BM25 approach over all fine-tuned retrievals could
be attributed to this.

Nevertheless, human evaluation indicated that the initial fine-tuned QA system
outperformed the baseline, displaying significant improvement in domain-
specific natural language processing on real-world sleep-related questions.
Human evaluation of enhanced fine-tuned QA system was done during a pilot
Randomized Controlled Trial (RCT). Participants in the intervention arm
interacted with health coaches who were supported by the enhanced fine-tuned
QA system, while in the control group, they directly interacted with the QA
system without the intermediate interaction with a coach. Our future work
includes analyzing the 5,000 questions and participants’ feedback on the
provided answers collected from this pilot RCT.

## Acknowledgment

The authors would like to acknowledge the Accelerating Creativity and
Excellence (ACE) Award (NTU-ACE2020-05) and center funding from Nanyang
Technological University, Singapore. Josip Car’s post at Imperial College
London is supported by the NIHR NW London Applied Research Collaboration.

## References

  * Bojic et al. (2023) Iva Bojic, Josef Halim, Verena Suharman, Sreeja Tar, Qi Chwen Ong, Duy Phung, Mathieu Ravaut, Shafiq Joty, and Josip Car. 2023.  A data-centric framework for improving domain-specific machine reading comprehension datasets.  _arXiv preprint arXiv:2304.00483_.
  * Bojic et al. (2022) Iva Bojic, Qi Chwen Ong, Megh Thakkar, Esha Kamran, Irving Yu Le Shua, Jaime Rei Ern Pang, Jessica Chen, Vaaruni Nayak, Shafiq Joty, and Josip Car. 2022.  Sleepqa: A health coaching dataset on sleep for extractive question answering.  In _Machine Learning for Health_ , pages 199–217. PMLR.
  * Chen et al. (2017) Danqi Chen, Adam Fisch, Jason Weston, and Antoine Bordes. 2017.  Reading Wikipedia to answer open-domain questions.  _arXiv preprint arXiv:1704.00051_.
  * Kia et al. (2022) Mahsa Abazari Kia, Aygul Garifullina, Mathias Kern, Jon Chamberlain, and Shoaib Jameel. 2022.  Adaptable closed-domain question answering using contextualized CNN-attention models and question expansion.  _IEEE Access_ , 10:45080–45092.
  * Kisling and Das (2021) Lisa A Kisling and Joe M Das. 2021.  Prevention strategies.  In _StatPearls [Internet]_. StatPearls Publishing.
  * Lin et al. (2021) Jimmy Lin, Xueguang Ma, Sheng-Chieh Lin, Jheng-Hong Yang, Ronak Pradeep, and Rodrigo Nogueira. 2021.  Pyserini: A Python toolkit for reproducible information retrieval research with sparse and dense representations.  In _Proceedings of the 44th Annual International ACM SIGIR Conference on Research and Development in Information Retrieval (SIGIR 2021)_ , pages 2356–2362.
  * Organization et al. (2014) World Health Organization et al. 2014.  _Global status report on noncommunicable diseases 2014_.  WHO/NMH/NVI/15.1. World Health Organization.
  * Rajpurkar et al. (2018) Pranav Rajpurkar, Robin Jia, and Percy Liang. 2018.  Know what you don’t know: Unanswerable questions for SQuAD.  _arXiv preprint arXiv:1806.03822_.
  * Wang et al. (2006) Mengqiu Wang et al. 2006.  A survey of answer extraction techniques in factoid question answering.  _Computational Linguistics_ , 1(1):1–14.
  * Yang et al. (2020) Juan Yang, Brent A Bauer, Stephanie A Lindeen, Adam I Perlman, Kasey R Boehmer, Manisha Salinas, Susanne M Cutshall, et al. 2020.  Current trends in health coaching for chronic conditions: A systematic review and meta-analysis of randomized controlled trials.  _Medicine_ , 99(30).
  * Zhang et al. (2023) Muru Zhang, Ofir Press, William Merrill, Alisa Liu, and Noah A Smith. 2023.  How language model hallucinations can snowball.  _arXiv preprint arXiv:2305.13534_.
  * Zhu et al. (2021) Fengbin Zhu, Wenqiang Lei, Chao Wang, Jianming Zheng, Soujanya Poria, and Tat-Seng Chua. 2021.  Retrieving and reading: A comprehensive survey on open-domain question answering.  _arXiv preprint arXiv:2101.00774_.