R. Spurzem, F. Rizzuto, M. Arca Sedda, A. Kamlah, P. Berczik, et al.

11institutetext: Astronomisches Rechen-Institut, Zentrum für Astronomie,
University of Heidelberg, Mönchhofstrasse 12-14, 69120 Heidelberg, Germany
11email<EMAIL_ADDRESS>22institutetext: Kavli Institute for
Astronomy and Astrophysics, Peking University,
Yiheyuan Lu 5, Haidian Qu, 100871, Beijing, China 33institutetext: Department
of Physics, University of Helsinki, Gustaf Hallstromin katu 2, 00014 Helsinki,
Finland 33email<EMAIL_ADDRESS>44institutetext: Dep. of
Physics and Astronomy ”Galileo Galilei” University of Padova, Via Marzolo 8,
35131 Padova, Italy 44email<EMAIL_ADDRESS>55institutetext: Max-
Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
55email<EMAIL_ADDRESS>66institutetext: Main Astronomical Observatory,
National Academy of Sciences of Ukraine, 27 Akademika Zabolotnoho St., 03143
Kyiv, Ukraine 66email<EMAIL_ADDRESS>77institutetext: Konkoly
Observatory, Research Centre for Astronomy and Earth Sciences, Eötvös Loránd
Research Network (ELKH), Konkoly Thege Miklós út 15-17, 1121 Budapest, Hungary
88institutetext: Department of Astronomy, School of Physics, Peking
University, Yiheyuan Lu 5, Haidian Qu, 100871, Beijing, China 88email:
<EMAIL_ADDRESS>99institutetext: Department of Earth Science and Astronomy,
College of Arts and Sciences, The University of Tokyo, 3-8-1 Komaba, Meguro-
ku, Tokyo 153-8902, Japan 99email<EMAIL_ADDRESS>1010institutetext: Max-Planck Institute for Astrophysics, Karl-Schwarzschild-
Str. 1, 85741 Garching, Germany 1010email<EMAIL_ADDRESS>
# Growth of Seed Black Holes in Galactic Nuclei

Rainer Spurzem 1122 Francesco Rizzuto 33 Manuel Arca Sedda 44 Albrecht Kamlah
5511 Peter Berczik 116677 Qi Shu (舒琦) 8822 Ataru Tanikawa (谷 川 衝) 99
Thorsten Naab 1010

###### Abstract

The evolution of dense star clusters is followed by direct high-accuracy
N-body simulation. The problem is to first order a gravitational N-body
problem, but stars evolve due to astrophysics and the more massive ones form
black holes or neutron stars as compact remnants at the end of their life.
After including updates of stellar evolution of massive stars and for the
relativistic treatment of black hole binaries we find the growth of
intermediate mass black holes and we show that in star clusters binary black
hole mergers in the so-called pair creation supernova (PSN) gap occur easily.
Such black hole mergers have been recently observed by the LIGO-Virgo-KAGRA
(LVK) collaboration, a network of ground based gravitational wave detectors.

##  1 Introduction

It is one of the grand challenges of theoretical astrophysics to understand
the dynamics of dense star clusters, both in their form as galactic globular
clusters orbiting in the Milky Way halo as well as nuclear star clusters,
surrounding the central supermassive black holes (SMBH) in our Galaxy and in
other galaxies. High precision dynamical simulations of star clusters use
direct orbit integration under the influence of (in principle) all other stars
allowing precise modeling of diffusive transport processes of mass, energy,
and angular momentum in the star cluster. This is important to understand
relaxation processes in the system. The physical and astrophysical challenge
is not only the gravitational million-body problem, but also the presence of a
large number of very tight binaries and stellar evolution with black holes,
neutron stars and white dwarfs forming. This becomes a strong multi-scale
problem, with orbital time scales of days coupling to overall crossing times
of the cluster of million years and ages of billion years.

The evolution of dense star clusters is not only governed by the aging of
their stellar populations and simple Newtonian dynamics. The stellar densities
become so high that stars can interact and collide, stellar evolution and
binary stars change the dynamical evolution, black holes can accumulate in
their centers and merge with relativistic effects becoming important. Recent
high-resolution imaging has revealed even more complex structural properties
with respect to stellar populations, binary fractions and compact objects as
well as – the still controversial – existence of intermediate mass black holes
in clusters of intermediate mass. Dense star clusters therefore are the ideal
laboratory for the concomitant study of stellar evolution and Newtonian as
well as relativistic dynamics. Last but not least black holes forming and
evolving in dense star clusters are one of the prominent sources of
gravitational waves (GW) across all frequency windows.

##  2 Current Astrophysical Updates

###  2.1 Stellar Evolution of Massive Stars

Before the first LIGO-Virgo-KAGRA (LVK) gravitational wave (GW) detection,
many theoretical models of stellar evolution predicted stellar black holes
(BHs) masses to be lower than 30 ${\rm M}_{\odot}$. These models remained
unchallenged for several years because all stellar BHs observed at the time
had masses $\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}\\!20\,{\rm
M}_{\odot}$ (Ziółkowski 2008; Özel et al. 2010). Surprisingly, the first LVK
detection, GW150914, revealed components more massive than 30 ${\rm
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 up-to-date stellar evolution models for the
correct interpretation and 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 [13, 49, 48]. Due to (P)PSN, isolated massive stars are not
supposed to collapse into BHs in the mass range of approximately
$50\\!-\\!130\,{\rm 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\,{\rm
M}_{\odot}$.

Figure 1:  Initial-Final mass relation (IFMR) for the escaping compact objects
of the MOCCA and Nbody6++GPU (Nbody) simulations. The keys refer to the Nbody-
delayedSNe-Uniform, Nbody-rapidSNe-Sana, MOCCA-delayedSNe-Uniform, and MOCCA-
rapidSNe-Sana simulations, respectively. They differ by the use of the code
(Nbody or MOCCA), and the remnant mass prescription for the stellar mass black
holes up to approximately the pair-instability mass gap (see Fryer et al.
2012). The black points show BH masses from another $N$-body simulation with
Level A parameters[10] (Plot taken from Kamlah2022).

Most of current updates of stellar evolution have been published for NBODY7[6]
and for NBODY6++GPU[27]. The major issues are summarized here:

  * (i)

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

  * (ii)

Pair-instability supernova and pulsation pair-instability supernova models
(according to Belczynski2016 incorporated in the remnant formation and
supernovae models as described in Fryer2012. Stars with helium core with
masses $>40{\rm M}_{\odot}$ undergo a violent phase of mass loss. For helium
cores in the range between $60\\!-\\!135/,{\rm M}_{\odot}$ the star is
completely destroyed.

  * (iii)

New prescription of BHs and NSs natal kick velocities that explicitly depend
on the fallback fraction[6].

  * (iv)

A model for electron capture supernovae (ECSN) that produces neutron stars
with low-velocity kicks that are therefore likely retained in medium-size star
clusters[36, 16].

Fig. 1 shows an example of how important the new updates of stellar evolution
are. It shows the initial-final mass relation (IFMR), which tells what is the
mass of a stellar evolution remnant as a function of the initial mass of a
star. On the y-axis all objects more massive than a few solar masses are black
holes; the figure compares the old IFMR with two variants of the new one
(delayed and rapid supernovae explosions), for our direct N-body simulations
as well as for an approximate Monte Carlo model (MOCCA)[18, 17].

###  2.2 Preparation for PopIII stars

Fitting formulae have been worked out for evolution tracks of massive stars
with $\lower 2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}M\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}160\,{\rm M}_{\odot}$ under extreme
metal poor (EMP) environments for $\log(Z/{\rm Z}_{\odot})=-2$, -4, -5, -6,
and -8, where ${\rm M}_{\odot}$ and ${\rm Z}_{\odot}$ are the solar mass and
metallicity, respectively[43]. They are based on reference stellar models,
newly obtained by simulating the time evolutions of EMP stars. The fitting
formulae take into account stars ending as blue supergiant (BSG) stars, and
stars skipping the Hertzsprung gap (HG) phases and blue loops, which are
characteristic of massive EMP stars. Here stars may remain BSG stars when they
finish their core Helium burning (CHeB) phase. The fitting formulae are in
good agreement with the stellar evolution models, and are now used for our
NBODY6++GPU code within the SSE/BSE packages; they are also used for other
codes such as PeTar and NBODY7 (see discussion in Section 4). The algorithms
should be useful to generate theoretical predictions for black holes and black
hole mergers under EMP environments. The work with NBODY6++GPU is currently in
progress on the Juwels-Booster system (Kamlah, Tanikawa, et al., in
preparation).

###  2.3 Relativistic Dynamics of Compact Objects

The gravitational energy loss and resulting merger of compact objects is
computed following the orbit-averaged approach[35], allowing for the final
coalescence if the orbit shrinking time due to gravitational wave emission
becomes shorter than an orbital time. We have added the following further
updates, which affect the formation and evolution of black holes[40]:

  * (i)

For collisions between a compact remnant and a main sequence star or red giant
a free parameter $f_{c}$ is introduced, which describes the mass loss from the
system in the process. The previous NBODY6 versions used only $f_{c}=1$, i.e.
no mass loss in the process.

  * (ii)

Simultaneous treatment of classical tidal interactions (Roche lobe overflow)
and Post-Newtonian orbit-averaged orbit shrinking due to gravitational wave
emission has been made possible. Both are treated technically in a similar
way, and can now be switched on together.

  * (iii)

Strongly bound binaries of two compact objects, which are subject to Post-
Newtonian relativistic energy loss are prevented from unperturbed two-body
integration, and defined as a new type of binary in the code.

Figure 2:  Visualization of the formation path towards a ”mass gap” merger
(grey region, third generation) of two black holes with mass $70{\rm
M}_{\odot}$ and $68{\rm M}_{\odot}$ developed in one of our N-body
simulations. The more massive BH (third generation) grew by two preceding
mergers involving black holes (first and second generation). The lower mass BH
(third generation) was created in a stellar merger of a red giant with a main
sequence star followed by the collision with a stellar mass BH. The masses of
the components (in solar masses) and the orbital periods (in days) are
indicated at the respective times of the merger (black horizontal lines) after
the start of the simulation (Figure from ArcaSedda2021).

Fig. 2 shows an example, and illustrates, how an intermediate mass black hole
is formed through several steps, involving collisions of two massive stars, a
collision between a massive star and a black hole and relativistic mergers of
black hole binaries. The procedure and parameters described above will affect
rates and timescales of the black hole formation process.

###  2.4 Relativistic Recoils at Coalescence and Spins

Another important physical process that has not yet been used for the
published papers so far is relativistic recoil for compact object
coalescences. Its absence might artificially enhance the probability of
forming massive black holes[5, 4]. We have implemented now the relativistic
recoils following Morawski2018, and currently tests are ongoing on the Juwels-
Booster system (Arca Sedda et al., in prep.).

Regarding spins (have been neglected so far, and are important to determine
proper values of recoil at coalescence) we have now included in our simulation
code a new spin treatment[7], which provides three alternative models for spin
evolution of stars, with variable efficiency of angular momentum transport, as
well as further options to choose different spin configurations for single and
binary black holes (Arca Sedda et al., in prep.). The one with highest
efficiency by Fuller2019 delivers at the end point of massive star evolution
black holes with low spins, consistent with current LVK observations.

Last, but not least the ultimate goal is to use a full Post-Newtonian dynamics
inside the regularized binary motion. In such approach relativistic periastron
precession and energy loss due to gravitational radiation can be followed at
each point of the orbit using a generalized quasi-classical equation of
motion[28]. Recently also the inclusion of spin dynamics and spin-spin and
spin-orbit interactions are included into the Post-Newtonian approach in our
codes[41]. Our cited papers describe the method to integrate this into our
codes; any reader interested in the relativistic theory for Post-Newtonian
dynamics please refer to papers cited therein.

##  3 Results

###  3.1 Intermediate Mass Black Hole Formation

Young dense massive star clusters are a promising environment for the
formation of intermediate mass black holes (IMBHs) through collisions (of
massive stars) and coalescences (of smaller mass black holes). We have
published a set of 80 simulations carried out with NBODY6++GPU using 10
different initial conditions, and shown that an IMBH can form in some cases
with – so far – up to $350\,{\rm M}_{\odot}$. We simulated compact star
clusters with $1.1\cdot 10^{5}$ particles ($\sim 7\cdot 10^{4}\,{\rm
M}_{\odot}$, core density of $\sim 10^{5}\,{\rm M}_{\odot}{\rm pc}^{-3}$) with
a resolved stellar population with 10% initial (primordial) hard binaries, and
find that very massive stars with masses up to $\sim 400\,{\rm M}_{\odot}$
grow rapidly by binary exchange and three-body scattering events with main
sequences stars in hard binaries. From them IMBHs with masses up to $350\,{\rm
M}_{\odot}$ form on timescales of order 15 Myr; the final mass depends
critically on an unknown parameter describing how much mass is accreted, if a
black hole collides with a main sequence star (see description above of
$f_{c}$ parameter).

Figure 3:  Top panel: GW strain evolution as a function of the frequency for
the mergers beyond the mass-gap (coloured lines). Simulated tracks are
overlapped to the sensitivity curves of LISA and LIGO-Virgo-KAGRA (LVK, solid
black lines), and DECIGO and Einstein Telescope (dashed black lines). The
white boxes indicate the time to merger for the heaviest merger. Bottom panel:
eccentricity evolution as a function of the frequency. The vertical lines
identify the moment in which the IMBH-BH mergers enter and exit the LISA
sensitivity window, whereas horizontal lines identify the eccentricity values
$e=0.1,0.5$. All the mergers are assumed to happen at a redshift $z=0.05$,
corresponding to a luminosity distance $D_{L}=230$ Mpc. The figure title
report the typical signal-to-noise ratio (S/N) $\,\sim\\!20\\!-\\!26$,
assuming for LISA a $T_{\rm obs}=4$ yr long mission. Figure from
ArcaSedda2021.

It is the first time that such IMBH formation process has been shown in a
direct N-body simulation of that quality and particle resolution. This process
was qualitatively predicted from Monte Carlo MOCCA simulations[18, 17]. After
formation, the IMBHs can experience occasional mergers with stellar mass black
holes in intermediate mass-ratio inspiral events on a 100 Myr timescale. For
more details please compare Rizzuto2021.

###  3.2 Black Hole Merger in the Forbidden Zone

The LIGO-Virgo-KAGRA Collaboration (LVC) discovered GW190521, a gravitational
wave (GW) source associated with the merger between two black holes (BHs) with
masses of 66 and $>\\!85\,{\rm M}_{\odot}$. GW190521 represents the first BH
binary merger with a primary mass falling in the PSN mass gap (see explanation
in Section 2.1) and leaving behind a $\sim 150\,{\rm M}_{\odot}$ remnant. So
far, the LVC has reported the discovery of four further mergers having a total
mass $>\\!100\,{\rm M}_{\odot}$, i.e., in the intermediate-mass black hole
(IMBH) mass range. In our simulations we discover the development of a
GW190521-like system as the result of a third-generation merger, and
furthermore four IMBH-BH mergers with total mass $(300-350)\,{\rm M}_{\odot}$.
We show that these IMBH-BH mergers are low-frequency GW sources detectable
with LISA and Decihertz Interferometer Gravitational wave Observatory (DECIGO)
out to redshift $z=0.01-0.1$ and $z>100$, and we discuss how their detection
could help unraveling IMBH natal spins. For the GW190521 test case, we show
that the third-generation merger remnant has a spin and effective spin
parameter that matches the 90% credible interval measured for GW190521 better
than a simpler double merger and comparable to a single merger. Due to GW
recoil kicks, we show that retaining the products of these mergers require
birth sites with escape velocities $\,\simeq\\!50\\!-\\!100$ km/s, values
typically attained in galactic nuclei and massive clusters with steep density
profiles. So, this is an explanation why observed black hole mergers do occur
in the PSN mass gap – they are the result of several generation mergers in
star clusters[3, 39].

Figure 4:  The panels show the primary ($m_{1}$) and secondary ($m_{2}$)
masses of all BH mergers in the simulations for an accretion fraction of
$f_{\rm c}=1.0$ (left, grey circles), $f_{\rm c}=0.5,$ (center, blue circles)
and $f_{\rm c}=0.0,$ (bottom, green circles). The currently available LIGO-
Virgo-KAGRA 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_{\rm
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. One of the $f_{\rm c}=0.0,$ realizations generated an
intermediate-mass ratio inspiral of two black holes with 31 and 181 ${\rm
M}_{\odot}$), respectively, as shown in the bottom panel (Figure from
Rizzuto2022).

We have presented more direct N-body simulations, carried out with
NBODY6++GPU, of young and compact low-metallicity (Z = 0.0002) star
clusters[39]. Early on, after tens of Myrs, every simulated cluster hosts
several black hole merger events which nearly cover the complete mass range of
primary and secondary black hole masses for current LIGO-Virgo-KAGRA
gravitational wave detections. The importance of gravitational recoil is
estimated statistically during post-processing analysis. We presented possible
formation paths of massive black holes above the assumed lower PSN mass-gap
limit ($45\,{\rm M}_{\odot}$) into the intermediate mass black hole (IMBH)
regime ($>100\,{\rm M}_{\odot}$) which include collisions of stars, black
holes and the direct collapse of stellar merger remnants with low core masses.
Fig. 3 shows how the black hole binaries found in our simulations would show
up in the detection sensitivity diagrams of current ground and future space
based gravitational wave detectors. Finally, Fig. 4 compares in a statistical
way current LIGO-Virgo-KAGRA detections of binary black hole mergers with a
collection of such mergers from our simulations, using in three different
panels a different mass loss factor $f_{c}$ for star-black hole collisions.
While we need more observational and simulation data to improve the
statistical quality it can be seen already that such comparisons would allow
to constrain the physics of black hole star collisions. In a similar way in
the future modelling the spin of black holes could and comparing with
observations could tell us something about initial and final spins of black
holes (which is quite difficult to observe directly via gravitational waves,
at least currently).

##  4 Initial Models and Codes

###  4.1 Initial Models

The simulations take into account full stellar evolution as well as the
formation and evolution of binary stars. We set the number of primordial
binaries typically to be 5 or 10 % of the total number of systems (these are
persistent hard binaries with binding energy higher than the r.m.s. random
kinetic energy of a star; a much larger number of soft binaries present at the
time of formation of the star cluster is being disrupted early on).

###  4.2 NBODY6++GPU

NBODY6++GPU111Link to repositories: tarballs/git service at
https://zenodo.org/record/6511341/files/Nbody6%2B%2BGPU-Jan2022.tgz
https://github.com/kaiwu-astro/Nbody6PPGPU-beijing is a high-precision direct
N-body simulation code based on the earlier N-body codes NBODY1-6[1] and
NBODY6++[42]. It uses for time integration Taylor series up to $4^{\rm th}$
order; due to the Hermite scheme it can be based on two time points only. This
together with the hierarchically blocked variable time step scheme allows an
efficient parallelization of the code for massively parallel supercomputers
(since NBODY6++); gravitational forces between particles are offloaded to
graphics processing units (GPUs), used for high-performance general purpose
computing (NBODY6++GPU[47]). The parallelisation is achieved via MPI and
OpenMP on the top level, distributing work within a group of particles due for
time integration, and efficient parallel use of GPU cores at the base level
(every MPI process using a GPU), for computing the gravitational forces
between particles. The GPU implementation in NBODY6++GPU provides a
significant performance improvement, especially for the long-range (regular)
gravitational forces (see Nitadori2012,Wang2015). Benchmarks and profiling are
published in detail in Wang2016,Huang2016.

Recipes to simulate single, binary, and multiple objects stellar evolution are
based on the SSE and BSE programs by Hurley[25, 24, 23] (see also code
references Hurley2013a,Hurley2013b). It includes rapid tidal circularization
for binaries with small pericenters and tidal captures[30]. The integrator
fully resolves orbits and dynamical evolution of binaries, even during phases
of mass loss or when one of the two stars undergoes a supernova explosion. The
binary orbit is adjusted to the corresponding loss of mass, energy and angular
momentum with appropriate time stepping; in case of a supernova explosion it
is always ensuring that the remnant and its companion leave the explosion with
the corrected orbital positions and velocities. Further recent improvements
are metallicity dependent winds, delayed or rapid white dwarf kicks[14],
electron-capture and pair instability supernova events. The reader interested
in more details please refer to our papers Banerjee2020,Kamlah2022, which
summarize the code updates and give all necessary further references.

###  4.3 New developments

A completely new code, called PeTaR[45] has been introduced. This code is
written in a coherent, modular way. It is supposed to contain all the physics
of NBODY6++GPU, but since it is a new code a lot of testing and comparison is
still required. Algorithmically it has two advantages over NBODY6++GPU
currently, one is the efficient parallelization of hard binaries, and the
other is the replacement of distant regular gravitational forces from
particles by a TREE based scheme, which makes it possible to take into account
also small perturbations on binaries from distant objects without too much
computational cost. It is possible to upgrade NBODY6++GPU accordingly; we
expect that both codes will co-exist for foreseeable time. Nevertheless we
have intensively tested PeTaR, see Fig. 5and may use it also in the next
computing period in this project.

A novel hierarchical $4^{\rm th}$ fourth-order forward symplectic integrator
and its numerical implementation has been shown in a new GPU-accelerated
direct-summation N-body code named FROST[37]. The new integrator also uses an
innovative MSTAR chain scheme[38] instead of the classical and algorithmic
chains by Mikkola1998. The integrator claims to be especially suitable for
simulations with a large dynamical range due to its hierarchical nature, and
for direct-summation N-body simulations beyond $N=10^{6}$ particles on systems
with several hundred and more GPUs. In that respect it is very similar to the
classic $4^{\rm th}$ order Hermite codes such as $\varphi$GRAPE[19, 11] or
HiGPU[12]. These codes, lacking the Ahmad-Cohen neighbour scheme, can easily
used hundreds of GPUs, because they compute full long-range gravitational
forces even for the smallest time steps. So, they gain parallelism by
introducing order N unnecessary computations, which can be avoided by using
NBODY6++GPU. In its AC neighbour scheme full force calculations (done on the
GPUs) are only done in order ten times larger time intervals than the smalles
steps. For the smallest steps only order 50-200 neighbour particle forces are
required for force calculation, and this is done efficiently by using OpenMP
on the multi-core host CPU[47, 46]. Therefore, NBODY6++GPU often appears to be
not to scale well to large GPU numbers, because it already obtains physical
results comparable to $\varphi$GRAPE or HiGPU (with many GPUs) by using only
few GPUs and OpenMP with CPU cores efficiently.

It should be noted, however, that the new code PeTar uses a TREE scheme for
distant forces, which in principle can be competitive with the AC scheme. A
full quantitative profiling analysis of these issues is still missing.
However, Fig. 5 shows first benchmarks done by our team using the PeTar code,
which are very promising and open up the path to $10^{7}$ particles in direct
N-body simulation, especially with large binary fractions (up to 50%, a region
in which NBODY6++GPU currently lags behind).

Figure 5:  Strong scaling of the novel PeTar code, showing wall clock times
obtained on the Juwels-Booster as a function of number of compute nodes.

## Acknowledgements

The authors gratefully acknowledge the Gauss Centre for Supercomputing 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
Supercomputer JUWELS (2021) at Jülich Supercomputing Centre (JSC). MAS is a
Marie Curie Skłodowska Fellow and Alexander von Humboldt Fellow. AK is a
fellow of the International Max Planck Research School for Astronomy and
Cosmic Physics at the University of Heidelberg (IMPRS-HD). PB and RS
acknowledge the support of the Volkswagen Foundation under the Trilateral
Partnerships 626 grant No. 97778 and the Strategic Priority Research Program
(Pilot B) Multi-wavelength gravitational wave universe of Chinese Academy of
Sciences (No. XDB23040100). PB acknowledges support from Volkswagen Foundation
under the special stipend No. 9B870 (2022), from a President’s International
Fellowship for Visiting Scientists of Chinese Academy of Sciences, from the
Science Committee of the Ministry of Education and Science of the Republic of
Kazakhstan (Grants No. AP08856184 and AP08856149), from the National Academy
of Sciences of Ukraine under the Main Astronomical Observatory GPU computing
cluster project No. 13.2021.MM, and by the special program of the National
Research Foundation of the Ukraine ‘Leading and Young Scientists Research
Support’ - “Astrophysical Relativistic Galactic Objects (ARGO): life cycle of
active nucleus”, No. 2020.02/0346. We thank Sambaran Banerjee, Mirek Giersz,
Jarrod Hurley, Arek Hypki, Nadine Neumayer, Long Wang, Kai Wu, Roberto
Capuzzo-Dolcetta, Andreas Just, M.B.N. (Thijs) Kouwenhoven, Xiaoying Pang for
helpful discussions, collaboration, and hospitality during visits. MAS
acknowledges funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement No.
101025436 (project GRACE-BH, PI Manuel Arca Sedda).

## References

  * [1] Aarseth, S. J., From NBODY1 to NBODY6: The Growth of an Industry, PASP1111333–13461999.
  * [2] Abbott, B. P., et al., Observation of Gravitational Waves from a Binary Black Hole Merger, PhRvL116061102-2016.
  * [3] Arca-Sedda, M., Rizzuto, F. P., Naab, T., Ostriker, J., Giersz, M., Spurzem, R., Breaching the Limit: Formation of GW190521-like and IMBH Mergers in Young Massive Clusters, ApJ920128-2021.
  * [4] Banerjee, S., Merger rate density of stellar-mass binary black holes from young massive clusters, open clusters, and isolated binaries: Comparisons with LIGO-Virgo-KAGRA results, PhRvD105023004-2022.
  * [5] Banerjee, S., Stellar-mass black holes in young massive and open stellar clusters - IV. Updated stellar-evolutionary and black hole spin models and comparisons with the LIGO-Virgo O1/O2 merger-event data, MNRAS5003002–30262021.
  * [6] Banerjee, S., Belczynski, K., Fryer, C. L., Berczik, P., Hurley, J. R., Spurzem, R., Wang, L., BSE versus StarTrack: Implementations of new wind, remnant-formation, and natal-kick schemes in NBODY7 and their astrophysical consequences, A&A639A41-2020.
  * [7] Belczynski, K., Klencki, J., Fields, C. E., et al., Evolutionary roads leading to low effective spins, high black hole masses, and O1/O2 rates for LIGO/Virgo binary black holes, A&A636A104-2020.
  * [8] Belczynski, K., Heger, A., Gladysz, W., Ruiter, A. J., Woosley, S., Wiktorowicz, G., Chen, H.-Y., Bulik, T., O’Shaughnessy, R., Holz, D. E., Fryer, C. L., Berti, E., The effect of pair-instability mass loss on black-hole mergers, A&A594A97-2016.
  * [9] Belczynski, K., Bulik, T., Fryer, C. L., Ruiter, A., Valsecchi, F., Vink, J. S., Hurley, J. R., On the Maximum Mass of Stellar Black Holes, ApJ7141217–12262010.
  * [10] Belczynski, K., Kalogera, V., Bulik, T., A Comprehensive Study of Binary Compact Objects as Gravitational Wave Sources: Evolutionary Channels, Rates, and Physical Properties, ApJ572407–4312002.
  * [11] Berczik, P., Spurzem, R., Wang, L., Zhong, S., Huang, S., Up to 700k GPU cores, Kepler, and the Exascale future for simulations of star clusters around black holes., hpc..conf52–592013.
  * [12] Capuzzo-Dolcetta, R., Spera, M., Punzo, D., A fully parallel, high precision, N-body code running on hybrid computing platforms, JCoPh236580–5932013.
  * [13] Fowler, W. A., Hoyle, F., Neutrino Processes and Pair Formation in Massive Stars and Supernovae., ApJS9201-1964.
  * [14] Fryer, C. L., Belczynski, K., Wiktorowicz, G., Dominik, M., Kalogera, V., Holz, D. E., Compact Remnant Mass Function: Dependence on the Explosion Mechanism and Metallicity, ApJ74991-2012.
  * [15] Fuller, J., Ma, L., Most Black Holes Are Born Very Slowly Rotating, ApJL881L1-2019.
  * [16] Gessner, A., Janka, H.-T., Hydrodynamical Neutron-star Kicks in Electron-capture Supernovae and Implications for the CRAB Supernova, ApJ86561-2018.
  * [17] Giersz, M., Askar, A., Wang, L., Hypki, A., Leveque, A., Spurzem, R., MOCCA survey data base- I. Dissolution of tidally filling star clusters harbouring black hole subsystems, MNRAS4872412–24232019.
  * [18] Giersz, M., Leigh, N., Hypki, A., Lützgendorf, N., Askar, A., MOCCA code for star cluster simulations - IV. A new scenario for intermediate mass black hole formation in globular clusters, MNRAS4543150–31652015.
  * [19] Harfst, S., Gualandris, A., Merritt, D., Spurzem, R., Portegies Zwart, S., Berczik, P., Performance analysis of direct N-body algorithms on special-purpose supercomputers, NewA12357–3772007.
  * [20] Huang, S.-Y., Spurzem, R., Berczik, P., Performance analysis of parallel gravitational N-body codes on large GPU clusters, Research in Astronomy and Astrophysics1611-2016.
  * [21] Hurley, J. R., Pols, O. R., Tout, C. A., SSE: Single Star Evolution, ascl.softascl:1303.015-2013.
  * [22] Hurley, J. R., Tout, C. A., Pols, O. R., BSE: Binary Star Evolution, ascl.softascl:1303.014-2013.
  * [23] Hurley, J. R., Pols, O. R., Aarseth, S. J., Tout, C. A., A complete N-body model of the old open cluster M67, MNRAS363293–3142005.
  * [24] Hurley, J. R., Tout, C. A., Pols, O. R., Evolution of binary stars and the effect of tides on binary populations, MNRAS329897–9282002.
  * [25] Hurley, J. R., Pols, O. R., Tout, C. A., Comprehensive analytic formulae for stellar evolution as a function of mass and metallicity, MNRAS315543–5692000.
  * [26] Jülich Supercomputing Centre, JUWELS Cluster and Booster: Exascale Pathfinder with Modular Supercomputing Architecture at Juelich Supercomputing Centre, Journal of large-scale research facilities A13872021.
  * [27] Kamlah, A. W. H., Leveque, A., Spurzem, R., Arca Sedda, M., Askar, A., Banerjee, S., Berczik, P., Giersz, M., Hurley, J., Belloni, D., Kühmichel, L., Wang, L., Preparing the next gravitational million-body simulations: evolution of single and binary stars in NBODY6++GPU, MOCCA, and MCLUSTER, MNRAS5114060–40892022.
  * [28] Kupi, G., Amaro-Seoane, P., Spurzem, R., Dynamics of compact object clusters: a post-Newtonian study, MNRAS371L45–L492006.
  * [29] Kustaanheimo P., Stiefel E.L., Perturbation theory of Kepler motion based on spinor regularization, J. Reine Angew. Math.2182041965
  * [30] Mardling, R. A., Aarseth, S. J., Tidal interactions in star cluster simulations, MNRAS321398–4202001.
  * [31] Mikkola, S., Aarseth, S. J., An efficient integration method for binaries in N-body simulations, NewA3309–3201998.
  * [32] Morawski, J., Giersz, M., Askar, A., Belczynski, K., MOCCA-SURVEY Database I: Assessing GW kick retention fractions for BH-BH mergers in globular clusters, MNRAS4812168–21792018.
  * [33] Nitadori, K., Aarseth, S. J., Accelerating NBODY6 with graphics processing units, Monthly Notices of the Royal Astronomical Society424545–5522012.
  * [34] ”Ozel, F., Psaltis, D., Narayan, R., McClintock, J. E., The Black Hole Mass Distribution in the Galaxy, ApJ7251918–19272010.
  * [35] Peters, P. C., Mathews, J., Gravitational Radiation from Point Masses in a Keplerian Orbit, PhRv131435–4401963.
  * [36] Podsiadlowski, P., Langer, N., Poelarends, A. J. T., Rappaport, S., Heger, A., Pfahl, E., The Effects of Binary Evolution on the Dynamics of Core Collapse and Neutron Star Kicks, ApJ6121044–10512004.
  * [37] Rantala, A., Naab, T., Springel, V., frost: a momentum-conserving CUDA implementation of a hierarchical fourth-order forward symplectic integrator, MNRAS5025546–55622021.
  * [38] Rantala, A., Pihajoki, P., Mannerkoski, M., Johansson, P. H., Naab, T., MSTAR - a fast parallelized algorithmically regularized integrator with minimum spanning tree coordinates, MNRAS4924131–41482020.
  * [39] Rizzuto, F. P., Naab, T., Spurzem, R., Arca-Sedda, M., Giersz, M., Ostriker, J. P., Banerjee, S., Black hole mergers in compact star clusters and massive black hole formation beyond the mass gap, MNRAS512884–8982022.
  * [40] Rizzuto, F. P., Naab, T., Spurzem, R., Giersz, M., Ostriker, J. P., Stone, N. C., Wang, L., Berczik, P., Rampp, M., Intermediate mass black hole formation in compact young massive star clusters, MNRAS5015257–52732021.
  * [41] Sobolenko, M., Berczik, P., Spurzem, R., Merging timescale for the supermassive black hole binary in interacting galaxy NGC 6240, A&A652A134-2021.
  * [42] Spurzem, R., Direct N-body Simulations, JCoAM109407–4321999.
  * [43] Tanikawa, A., Yoshida, T., Kinugawa, T., Takahashi, K., Umeda, H., Fitting formulae for evolution tracks of massive stars under extreme metal-poor environments for population synthesis calculations and star cluster simulations, MNRAS4954170–41912020.
  * [44] Vink, J. S., de Koter, A., Lamers, H. J. G. L. M., Mass-loss predictions for O and B stars as a function of metallicity, A&A369574–5882001.
  * [45] Wang, L., Iwasawa, M., Nitadori, K., Makino, J., PETAR: a high-performance N-body code for modelling massive collisional stellar systems, MNRAS497536–5552020.
  * [46] Wang, L., Spurzem, R., Aarseth, S., Giersz, M., Askar, A., Berczik, P., Naab, T., Schadow, R., Kouwenhoven, M. B. N., The DRAGON simulations: globular cluster evolution with a million stars, MNRAS4581450–14652016.
  * [47] Wang, L., Spurzem, R., Aarseth, S., Nitadori, K., Berczik, P., Kouwenhoven, M. B. N., Naab, T., NBODY6++GPU: ready for the gravitational million-body problem, MNRAS4504070–40802015.
  * [48] Woosley, S. E., Pulsational Pair-instability Supernovae, ApJ836244-2017.
  * [49] Woosley, S. E., Blinnikov, S., Heger, A., Pulsational pair instability as an explanation for the most luminous supernovae, Natur450390–3922007.
  * [50] Woosley, S. E., Heger, A., Weaver, T. A., The evolution and explosion of massive stars, RvMP741015–10712002.
  * [51] Ziółkowski, J., Masses of Black Holes in the Universe, Chinese Journal of Astronomy and Astrophysics Supplement8273–2802008.