# Does Gaia Play Dice? : Simple Models of non-Darwinian Selection

Rudy Arthur,1 Arwen Nicholson,2
1University of Exeter, Department of Computer Science
2University of Exeter, Department of Physics and Astronomy
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper we introduce some simple models, based on rolling dice, to
explore mechanisms proposed to explain planetary habitability. The idea is to
study these selection mechanisms in an analytically tractable setting,
isolating their consequences from other details which can confound or obscure
their effect in more realistic models. We find that while the observable of
interest, the face value shown on the die, ‘improves’ over time in all models,
for two of the more popular ideas: Selection by Survival and Sequential
Selection, this is down to sampling effects. A modified version of Sequential
Selection, Sequential Selection with Memory, implies a statistical tendency
for systems to improve over time. We discuss the implications of this and its
relationship to the ideas of the ‘inhabitance paradox’ and the ‘Gaian
bottleneck’.

##  1 Introduction

Relatively recent discussion about the persistence of life over long periods
has brought to the fore various selection principles [1, 2, 3, 4]. With the
recent launch of the James Webb Space Telescope [5, 6] these questions not
only have important implications for our understanding of Earth history, but
also for the search for other inhabited planets. In the future, with a large
enough catalogue of inhabited planets, it may be possible to experimentally
investigate alternative trajectories for life [7]. Until then, understanding
general principles behind planetary habitability and inhabitance is a way to
provide working hypotheses that can explain both how ‘lucky’ the Earth is to
be inhabited and what we might expect on planets orbiting other stars.

The most discussed of these selection principles is called ‘Selection by
Survival’ (SBS) in [3] though the idea has many names, see [8]. The essential
idea is that a population where the entities have different rates of survival
will be ‘purified’ so that, in the long run, surviving entities have must have
properties conducive to survival. Several works e.g. [1, 2, 9], attempt to
disentangle Darwinian selection from this ‘differential persistence’
(essentially a synonym for SBS). [2] and [9] emphasise the importance of this
selection principle acting on higher order phenomena like whole ecosystems,
planetary scale biogeochemical cycles and the entire life-Earth coupled system
i.e. Gaia [10]. They argue that in systems of hereditary replicators Darwinian
selection is more powerful, however for entities like populations, ecosystems
or bio-geochemical cycles [11], which do not have strict heredity and
reproduction, SBS will operate to favour certain macroscopic features.
Specific examples suggested by e.g. [8] are sexual reproduction and
macroevolutionary freezing.

[3] defines another, related, selection principle called Sequential Selection
(SS) [12]. This is a similar idea to SBS, but motivated by the frequent
upheavals in the history of life on Earth and meant to account for life’s
apparent stabilising effect on Earth’s habitability. [3] propose a simple
algorithm - evolutionary innovations have a stabilizing or destabilizing
effect on the environment. If they have a destabilizing effect, habitability
is reduced, eventually eliminating the destabilizing innovation. In this way
destabilizing effects are eliminated by ‘near fatal resets’ while stabilizing
innovations persist and accumulate.

In [13, 4] we argue for a refinement of this algorithm, emphasising that the
resets are‘near fatal i.e. the evolutionary innovations developed during the
previous stable period are not completely lost. The algorithm of [3] applies:
destabilizing innovations lead to resets which greatly reduce species
abundance but have a lesser effect on species diversity. The life-earth system
which arises after the reset is selected from a larger ‘pool’, which has the
potential to generate better, more stable ecosystems. Higher species and
functional diversity give Gaia more tools to generate stability. The process
is completely blind, so unstable states can also be selected, however, by
definition, these are short lived and eventually a long-lived stable state
will arise. During this stable period species diversity can increase again
leading to a kind of ratcheting effect. To emphasise this cumulative process,
in contrast to the sequential selection algorithm of [3], we call this
‘Sequential Selection with Memory’ (SSM).

A variety of abstract models of varying complexity have been proposed to
explore these selection principles e.g. [1, 14, 15, 4, 16]. Though these
models have great value, it can sometimes be unclear which of their features
are programmed in (as alleged by [11] of the famous Daisyworld [17] model) and
which are emergent. It must also be said that the mathematical or
computational complexity of these models can give them an air of mystery -
especially to biologists not well versed in these methods. Indeed, the very
fact that the key model features are often emergent means that understanding
how they emerge requires a detailed understanding of each model’s dynamics.

Here we propose an extremely simple probability model as a setting to study
selection principles. The aim is to strip out as much complexity as possible
to understand the core meaning of these principles and their consequences. A
very loose analogy would be trying to understand the approximately Gaussian
distribution of, say, human height. This has some genetic and environmental
causes which, with great difficultly, could be experimentally isolated and
formulated into a mechanistic model of height, simulated and shown to result
in a Gaussian distribution. However, a much simpler, and in many ways more
satisfactory, explanation is that a Gaussian distribution is the expected
outcome for an observable which is a sum of independent effects.

Continuing the height analogy, by simulating sums of random variables and
showing this results in a Gaussian distribution we might start to suspect that
a more general principle is operating, one which isn’t affected by the
particular details of our model, in this case, the Central Limit Theorem. This
paper doesn’t propose anything as general as a statistical convergence
theorem, what we do propose are models simple enough to be analytically solved
but complex enough to see selection principles operating. Theses models will
be shown to exhibit interesting behaviour which is also observed in more
complex models. The aim is to provide some clarity on exactly what non-
Darwinian selection principles can do in a clear and tractable setting.

##  2 Introducing the Model

Consider an $M$ sided die with the rule that, once rolled, whatever number is
showing on the top face gives the number of steps to wait before rolling again
or finishing the game. For $r\geq 1$ dice we roll each one independently to
get $x_{1},x_{2},\ldots,x_{r}$, and take the highest face value:
$\max(x_{1},x_{2},\ldots,x_{r})$. Based on this consider the following dice
games:

  1. 1.

Selection By Survival(SBS): Roll $N$ (where $N$ is a very large number)
independent dice once each.

  2. 2.

Sequential Selection (SS): Roll one die repeatedly for $T$ time steps.

  3. 3.

Sequential Selection with Memory (SSM):

    1. (a)

Starting with $r=1$, roll $r$ dice repeatedly for $T$ time steps. Add a new
die every time the top face shows the maximum value, $M$.

    2. (b)

Starting with $M=1$, roll an $M$ sided die repeatedly for $T$ times steps.
Every time the top face shows the maximum value $M$, increase $M$ by 1.

The SSM games are reminiscent of the Polya Urn model, though have not been
studied before to our knowledge. The quantity of interest will be the expected
face value at time $t$. The names chosen are based on the discussion in the
Introduction and follow the conventions of [3] and [4]. Our version of
Selection By Survival is much simpler than the (mostly verbal) models proposed
by others e.g. [9] and most closely follows the graphical model from [2].

As a rough mapping to reality - a die represents an inhabited ‘planet’. Each
roll is a period of stability for the planet’s biosphere. The face value
represents something akin to the ‘fitness’ of the biosphere on that planet,
i.e. how long it will persist. If we observe an inhabited planet at some
random point in its history we may see a biosphere with properties conducive
to long term stability (high face value) or only short term stability (low
face value). The question of interest for astrobiology is, if we were to
survey a large catalogue of inhabited planets, what would be the average
‘fitness’? For Earth history (or for the history of any inhabited planet) the
equivalent question is, if we were to observe a planet at a random point in
its history, what should we expect about the habitability properties of that
planet?

##  3 Selection by Survival

$t=1$

$t=2$

$t=3$

$t=4$

$t=5$

$t=6$

Figure 1: One possible unfolding of the SBS game with $N=25$ and $M=6$. At
$t=1$ we have our initial ensemble, at $t=2$ we have removed all the 1s, at
$t=3$ we remove all the 2s etc.

Figure 1 shows one realisation of the SBS game. At $t=1$ all of the dice are
in play and the average face value (over very large $N$ or many different
realisations of the same game) is

$(1+2+\ldots+M)/M$

at $t=2$ all of the dice showing 1 on the top face are removed. Restricting
our survey to inhabited planets, the average face value is now

$(2+3+\ldots+M)/(M-1)$

At time $t\leq M$ the average face value is

$\frac{(t+(t+1)+\ldots+M)}{(M-t)}=\frac{M+t}{2}$ (1)

So that average face value increases linearly with time.

Figure 2: Average face value in the SBS game as a function of $t$ for an
$M=10$ sided dice over $N=1000$ dice rolls.

Figure 2 shows the result of simulations of the game compared to equation 1.
In terms of ‘planets’ this model is simply stating the (obvious) fact that
planets which survive have properties (high face value) which allow them to
survive! Looking at the catalogue of inhabited planets will necessarily yield
planets with properties conducive to maintaining life, without the need for
any additional mechanism.

This realisation has all planets are seeded with life at the same time. More
complex games could be devised (say a constant rate of habitable planet
generation) to study how the generateion rate interacts with this simple
selection mechanism. For this paper, SBS represents a basic null model - older
inhabited planets must have features which have enabled them to remain
inhabited. The growth in fitness of the ‘surviving’ planets is simply a
sampling artefact, the average fitness of an inhabited planet increases
because we throw away more and more of the unfit planets from our average.
Considering our solar system according to SBS, the single inhabited planet we
see is habitable because if it wasn’t, we wouldn’t be looking at it, or living
on it. Thus in this context, SBS is nothing more than an observer effect or
anthropic principle.

##  4 Sequential Selection

The Earth has experienced numerous mass extinction events, had very different
planetary regulation mechanisms, atmospheric composition, levels of volcanic
activity and life has persisted the entire time [18]. We seek to model these
sequential resets with another simple game: repeatedly rolling a single die.

Figure 3: One possible unfolding of the SS game with $T=20$ and $M=6$. At
$t=1$ we roll 2 which shows for 2 steps, we roll 1 which shows for 1 step,
then 3 for 3 steps etc. The game is played a large number, $N$, of times as in
figure 1 so we are interested in average behaviour.

When observing the die at a random time $t$, what should we expect the face of
the die to show, on average? The chance of the die showing $k$ is proportional
to the probability of rolling a $k$, $p(k)$, times the number of ‘slots’ where
the observation could occur e.g. if the die is showing 3 this could be an
observation of the die on the first, second or third step where it is face up.
We normalise this probability and compute the expected value for the top face
as

$\sum_{k=1}^{M}k\frac{kp(k)}{\sum_{k=1}^{m}kp(k)}$ (2)

For one die $p(k)=1/M$ and this simplifies to

$\frac{2M+1}{3}$ (3)

Note this is larger than the expected value of a single dice roll,
$\frac{M+1}{2}$ for $M>1$.

Figure 4: Average face value in the SS game as a function of $t$ for an $M=10$
sided dice, averaged over $N=1000$ independent instances. The expected value
of a single $10$ sided dice roll is $5.5$ which is less than the expected face
value in the sequential sampling game, $7$. Note the logged x-axis.

Figure 4 shows the results of simulations of the game compared to equation 3.
The results here imply that when observing a ‘planet’ at a random time, it is
more likely to be in a state with stability enhancing properties (high face
value). Again this reflects the obvious fact that if we depict Earth’s history
as a time line and pick a random point on the line we are more likely to pick
a point in a long stable period. In particular, our present time is most
likely to be a stable period, without the need for any additional mechanism.
Like with SBS, Earth’s current stability is simply an observer effect.

One thing missing from this game (and from the algorithm of [3]) is the
possibility of total extinction, that is, finishing the game early. We could
implement an additional rule, say when we roll a 1, stop the game. This would
give a model where SBS and SS are both operating simultaneously. Here, for
simplicity and clarity, we don’t account for total extinction, so as not to
mix the mechanisms. All of our SS games persist for the same amount of time.
More complex models e.g. [4, 19] do have the possibility of early stopping and
come to very similar conclusions.

##  5 Sequential Selection with Memory

The continued inhabitance of Earth and the fact that biodiversity has
increased over time motivates the final games. Each reset does not start from
scratch, but builds on evolutionary and ecological innovations that came
before. We propose 2 models with an extremely simple ‘memory’. This memory is
implemented in two ways, first by adding extra dice at fixed $M$, second by
increasing $M$.

###  5.1 Game A: Adding dice

The face value in this game is determined by rolling multiple dice and
choosing the one with the maximum face value. The idea is that stable
biospheres outlast unstable ones. One could imagine independent ecosystems co-
existing with the final ‘reset’ only occurring when the most stable subsystem
collapses. If this seems contrived, in more complex models e.g. [4], a similar
feature emerges as a consequence of model dynamics rather than being enforced.

Figure 5: One possible unfolding of SSM game A with $T=20$ and $M=6$. At $t=3$
we roll 6 which adds an extra die. At $t=13$ we roll six again which adds a
third die. The bottom row is the observable, the other rows show dice rolls
which are not observed.

The analysis is a little more complex that the previous two games. The first
thing we need is the probability to get the face value $f$ when rolling $r$
dice and applying the rule $f=\max(x_{1},x_{2},\ldots,x_{r})$. This is

$p_{r}(f)=\sum_{i=1}^{r}\binom{r}{i}p(f)^{i}p(x<f)^{r-i}$ (4)

where $p(f)=1/M$ and $p(x<f)=\frac{f-1}{M}$. This is just the probability to
get at least one $f$ and nothing higher, multiplied by a combinatoric factor.
To simplify this, consider arranging all the possible outcomes of $r$ rolls in
an r-dimensional hypercube. The number of ways to obtain $f$ is given by the
difference in volumes between an $f$ and $f-1$ sided hypercube so

$p_{r}(f)=\frac{f^{r}-(f-1)^{r}}{M^{r}}$ (5)

It is shown in Appendix A that the expected face value for large $M$ is

$E[f|r,M\gg 1]=M\frac{r}{r+1}+\frac{1}{2}.$ (6)

In a game with $r$ dice, the expected number of dice rolls before hitting the
value $M$, where we add an extra die, is $1/p_{r}(M)$. Therefore, the expected
time spent playing with exactly $r$ dice is

$T(r)=\frac{\sum_{i=1}^{M}ip_{r}(i)}{p_{r}(M)}.$ (7)

For large $M$ (using the summation result from Appendix A) this is

$T(r)\simeq\frac{M^{2}}{r+1}$ (8)

To calculate the expected face value at $t$, we first compute the expected
number of dice at time $t$ by solving

$\sum_{k}^{r}T(k)=t$ (9)

for $r$. Using the large $M$ approximation

$\sum_{k}^{r}T(k)\simeq M^{2}\sum_{k}^{r}\frac{1}{i+1}=M^{2}(H_{r}-1)$ (10)

where $H_{r}$ is the $r^{th}$ harmonic number. This has a standard
approximation, valid for large $r$ but quite accurate even at $r=1$:
$H_{r}\simeq\gamma+\ln(r)$, where $\gamma$ is the Euler-Mascheroni constant.
Substituting and solving for $r$ gives

$r(t)=\exp\left(\frac{t}{M^{2}}+1-\gamma\right)=A\exp\left(\frac{t}{M^{2}}\right)$
(11)

where $A=\exp(1-\gamma)$. The number of dice grows exponentially, with growth
rate $1/M^{2}$. The more faces a die has, the longer we have to wait to land
on a specific one, for example the time to go from $r$ dice to $2r$ dice is
$\simeq M^{2}\ln 2$.

The expected face value at time $t$ is the expected face value with $r(t)$
dice. Still working in the large $M$ limit this is

$E[f:M\gg 1](t)\simeq
M\frac{r}{r+1}=M\frac{A\exp\left(\frac{t}{M^{2}}\right)}{1+A\exp\left(\frac{t}{M^{2}}\right)}$
(12)

This is a sigmoid function in the variable $t/M^{2}$. At large values of $t$
the value is $M$, as expected, we have so many dice we are virtually
guaranteed to roll at least one $M$. The small $t$ behaviour is interesting,
the function is roughly linear which means, despite the exponential growth in
the number of dice suggested by equation 11, the expected face value grows
much more slowly.

Figure 6: Average face value in the SSM game as a function of $t$ for an
$M=10$ sided dice, averaged over $N=1000$ independent instances.

Figure 6 shows the results of 1000 simulations of the game with $M=10$
compared to the ‘exact’ answer, equation 12. We observe convergence to the
upper bound $M$ at a rate that is roughly linear in log time. Such slow
convergence is seen in more complex evolutionary models, especially the
Tangled Nature Model [20, 21] and its variants [13, 4]. There, it arises from
the simulated ecosystem successively crossing ‘entropic barriers’ [22]. Each
time a barrier is crossed the system is likely to be in a more stable
configuration with a higher barrier. This behaviour has been discussed before
in the language of record statistics and is also observed in physical systems
like spin glasses, colloids and high temperature superconductors [23].

This model is simple enough for an approximate analytic solution. This shows
that there is competition between the growth in the number of dice over time
against the growth in the time taken between trials. There is also a trade-off
between large values of $M$, leading to higher average rolls, versus time
taken to add a new die. What this model suggests is that selection plus
accumulation leads to slow growth in stability. This model implies that older
inhabited planets should be more habitable, so our presence on Earth is not
just an observer effect but a statistically more likely outcome.

###  5.2 Game B: Increasing $M$

This game similar to the previous one, except instead of adding extra dice we
have just one die and add a extra faces to it, which makes this harder to play
with real dice! A rough analogy to a real ecosystem is to assume that species
diversity is not lost after each collapse (dice roll) and that species persist
at low abundance, in dormant states or isolated refugia. Reaching a fitness
peak (hitting the max value of $M$) generates more latent diversity and allows
ecosystems to explore more of the so-called fitness landscape [24]. Thus each
reset has the potential to find a more stable state because the space of
possibilities is wider.

Figure 7: One possible unfolding of SSM game B with $T=20$. The bottom row
shows the actual face values and the top row shows the number of sides of the
die. For example at $t=11$ we roll a 4 on a 4-sided die, increasing the number
of sides to 5 for the next roll.

If the die has $M$ sides, the expected number of rolls required to hit the $M$
face is just $1/M$. Each roll is expected to last $\frac{M+1}{2}$ steps so the
expected waiting time before increasing the number of faces is

$T(M)=M\frac{M+1}{2}$ (13)

Summing up the wait times from each $M$ gives the total duration of the
experiment

$t=\sum_{i=1}^{M}\frac{i(i+1)}{2}=\frac{1}{2}\left(\frac{M(M+1)(2M+1)}{6}+\frac{M(M+1)}{2}\right)$
(14)

Keeping only the terms of highest order in $M$ and solving for $t$ gives

$M=\sqrt[3]{6t}$ (15)

Substituting into equation 3 gives

$E[f](t)=\frac{2\sqrt[3]{6t}+1}{3}$ (16)

Figure 8: Average face value in SSM game B as a function of $t$, averaged over
$N=1000$ independent instances.

Figure 8 shows the results of simulations of the game compared to the exact
answer, equation 16. Unlike game A there is no convergence and the expected
face value grows without bound, though fairly slowly. Again there is a trade
off between increasing $M$ by performing a large number of trials and the time
it takes to complete those trials. This is again reminiscent of Tangled Nature
Model dynamics [22, 13, 4] and other physical systems which cross energetic or
entropic barriers [23].

##  6 Discussion

These three mechanisms give three reasonable ideas about what to expect when
surveying large catalogues of inhabited planets, or looking at an inhabited
one at a random point in its history. The first two (SBS, SS) have no role for
life. The stability properties of inhabited planets are down to observer
effects - unless they had these properties we wouldn’t be looking at them. The
third mechanism is more interesting. Once a planet is inhabited life can have
a positive effect on habitability. In particular - inhabited planets have
properties conducive to stability because of their history of inhabitance.

This idea has appeared previously as ‘The inhabitance paradox’ in [25] and is
closely related to the idea of the Gaian bottleneck [26]. This paradox says
that for a planet to be habitable, it must be inhabited. This means life must
seize the reins and exert a stabilising effect early in a planet’s history or
go extinct due to deteriorating geophysical conditions - an effect dubbed the
Gaian bottleneck. The SSM game shows a very simple mechanism by which this
could occur, combining the sequential selection algorithm of [3] with some
method of making cumulative improvements will tend to generate more stable
systems. We have argued previously [4] that such cumulative processes are
widespread on Earth, for example: microbial seed banks, dormancy [27] and
lateral gene transfer [28] all contribute to the maintenance of microbial
diversity and therefore the stabilising effect of functional redundancy.

We hope that this model and its analysis provides some clarity on selection
principles as well as providing a sandbox for studying selection effects. In
particular, we believe that Sequential Selection with Memory provides a
plausible way for a complex system, like an inhabited planet, to become more
stable over time. We propose that Gaia - the stabilising and symbiotic
feedback of life and the environment - can be born through this kind of
natural, but non-Darwinian, selection.

##  Appendix A Expected value for the max of $r$, $M$-sided dice.

Rolling $r$, $M$ sided dice gives the face value $f$ with probability

$p(f)=\frac{k^{r}-(k-1)^{r}}{M^{r}}$

as discussed in the text. The expectation for the face value is therefore

$E[f]=\sum_{k=1}^{M}k\frac{k^{r}-(k-1)^{r}}{M^{r}}$

Writing out the sum explicitly

$\displaystyle 1.1^{r}+2.2^{r}+3.3^{r}+\ldots+M.M^{r}$ $\displaystyle-($
$\displaystyle 1.0^{r}+2.1^{r}+3.2^{r}+\ldots+M.(M-1)^{r})$

Shows that we can regroup and rewrite as

$E[f]=\frac{1}{M^{r}}\left(M^{r+1}-\sum_{k=1}^{M-1}k^{r}\right)$

The sum can be simplified using Faulhaber’s formula [29]

$\sum_{k=1}^{M-1}k^{r}=\frac{(M-1)^{r+1}}{r+1}+\frac{(M-1)^{r}}{2}+O(M^{r-1})$

where the lower order terms are fairly complex coefficients involving the
Bernoulli numbers. Substituting and taking the limit of large $M$ we get

$E[f]=M\frac{r}{r+1}+\frac{1}{2}$

as stated in the text.

## References

  * [1] Pierrick Bourrat.  From survivors to replicators: evolution by natural selection revisited.  Biology & Philosophy, 29(4):517–538, 2014.
  * [2] W Ford Doolittle.  Natural selection through survival alone, and the possibility of gaia.  Biology & Philosophy, 29(3):415–423, 2014.
  * [3] Timothy M Lenton, Stuart J Daines, James G Dyke, Arwen E Nicholson, David M Wilkinson, and Hywel TP Williams.  Selection for gaia across multiple scales.  Trends in Ecology & Evolution, 33(8):633–645, 2018.
  * [4] Rudy Arthur and Arwen Nicholson.  Selection principles for gaia.  Journal of Theoretical Biology, 533:110940, 2022.
  * [5] Ignas A. G. Snellen, F. Snik, M. Kenworthy, S. Albrecht, G. Anglada-Escudé, I. Baraffe, P. Baudoz, W. Benz, J. L. Beuzit, B. Biller, J. L. Birkby, A. Boccaletti, R. van Boekel, J. de Boer, Matteo Brogi, L. Buchhave, L. Carone, M. Claire, R. Claudi, J. M. Demory, B. O. Désert, S. Desidera, B. S. Gaudi, R. Gratton, M. Gillon, J. L. Grenfell, O. Guyon, T. Henning, S. Hinkley, E. Huby, M. Janson, C. Helling, K. Heng, M. Kasper, C. U. Keller, O. Krause, L. Kreidberg, N. Madhusudhan, A. M. Lagrange, R. Launhardt, T. M. Lenton, M. Lopez-Puertas, A. L. Maire, N. Mayne, V. Meadows, B. Mennesson, G. Micela, Y. Miguel, J. Milli, M. Min, E. de Mooij, D. Mouillet, M. N’Diaye, V. D’Orazi, E. Palle, I. Pagano, G. Piotto, D. Queloz, H. Rauer, I. Ribas, G. Ruane, F. Selsis, A. Sozzetti, D. Stam, C. C. Stark, A. Vigan, and Pieter de Visser.  Detecting life outside our solar system with a large high-contrast-imaging mission.  Experimental Astronomy, October 2021.
  * [6] Sascha P. Quanz, Olivier Absil, Willy Benz, Xavier Bonfils, Jean-Philippe Berger, Denis Defrère, Ewine van Dishoeck, David Ehrenreich, Jonathan Fortney, Adrian Glauser, John Lee Grenfell, Markus Janson, Stefan Kraus, Oliver Krause, Lucas Labadie, Sylvestre Lacour, Michael Line, Hendrik Linz, Jérôme Loicq, Yamila Miguel, Enric Pallé, Didier Queloz, Heike Rauer, Ignasi Ribas, Sarah Rugheimer, Franck Selsis, Ignas Snellen, Alessandro Sozzetti, Karl R. Stapelfeldt, Stephane Udry, and Mark Wyatt.  Atmospheric characterization of terrestrial exoplanets in the mid-infrared: biosignatures, habitability, and diversity.  Experimental Astronomy, pages 1–25, September 2021.
  * [7] A Lenardic, JW Crowley, AM Jellinek, and M Weller.  The solar system of forking paths: bifurcations in planetary evolution and the search for life-bearing planets in our galaxy.  Astrobiology, 16(7):551–559, 2016.
  * [8] Jan Toman and Jaroslav Flegr.  Stability-based sorting: The forgotten process behind (not only) biological evolution.  Journal of theoretical biology, 435:29–41, 2017.
  * [9] Frédéric Bouchard.  Ecosystem evolution is about variation and persistence, not populations and reproduction.  Biological Theory, 9(4):382–391, 2014.
  * [10] James E Lovelock and Lynn Margulis.  Atmospheric homeostasis by and for the biosphere: the gaia hypothesis.  Tellus, 26(1-2):2–10, 1974.
  * [11] W Ford Doolittle.  Darwinizing gaia.  Journal of Theoretical Biology, 434:11–19, 2017.
  * [12] Richard A Betts and Timothy M Lenton.  Second chances for lucky gaia: a hypothesis of sequential selection.  Met Office, 2008.
  * [13] Rudy Arthur and Arwen Nicholson.  An entropic model of gaia.  Journal of theoretical biology, 430:177–184, 2017.
  * [14] Arwen E Nicholson, David M Wilkinson, Hywel TP Williams, and Timothy M Lenton.  Alternative mechanisms for gaia.  Journal of theoretical biology, 457:249–257, 2018.
  * [15] Timothy M Lenton, Timothy A Kohler, Pablo A Marquet, Richard A Boyle, Michel Crucifix, David M Wilkinson, and Marten Scheffer.  Survival of the systems.  Trends in Ecology & Evolution, 36(4):333–344, 2021.
  * [16] Richard A Boyle and Timothy M Lenton.  The evolution of biogeochemical recycling by persistence-based selection.  Communications Earth & Environment, 3(1):1–14, 2022.
  * [17] Andrew J Watson and James E Lovelock.  Biological homeostasis of the global environment: the parable of daisyworld.  Tellus B: Chemical and Physical Meteorology, 35(4):284–289, 1983\.
  * [18] Tim Lenton and Andrew Watson.  Revolutions that made the Earth.  OUP Oxford, 2013.
  * [19] Rudy Arthur and Arwen Nicholson.  A gaian habitable zone, 2023.
  * [20] Kim Christensen, Simone A Di Collobiano, Matt Hall, and Henrik J Jensen.  Tangled nature: a model of evolutionary ecology.  Journal of theoretical Biology, 216(1):73–84, 2002.
  * [21] Rudy Arthur, Arwen Nicholson, Paolo Sibani, and Michael Christensen.  The tangled nature model for organizational ecology.  Computational and Mathematical Organization Theory, 23(1):1–31, 2017.
  * [22] Nikolaj Becker and Paolo Sibani.  Evolution and non-equilibrium physics: A study of the tangled nature model.  EPL (Europhysics Letters), 105(1):18005, 2014.
  * [23] Paolo Sibani, Stefan Boettcher, and Henrik Jeldtoft Jensen.  Record dynamics of evolving metastable systems: theory and applications.  The European Physical Journal B, 94(1):1–23, 2021.
  * [24] Rudy Arthur and Paolo Sibani.  Decision making on fitness landscapes.  Physica A: Statistical Mechanics and its Applications, 471:696–704, 2017.
  * [25] Colin Goldblatt.  The inhabitance paradox: How habitability and inhabitancy are inseparable.  arXiv preprint arXiv:1603.00950, 2016.
  * [26] Aditya Chopra and Charles H Lineweaver.  The case for a gaian bottleneck: the biology of habitability.  Astrobiology, 16(1):7–22, 2016.
  * [27] Jay T Lennon and Stuart E Jones.  Microbial seed banks: the ecological and evolutionary implications of dormancy.  Nature reviews microbiology, 9(2):119–130, 2011.
  * [28] Nigel Goldenfeld and Carl Woese.  Biology’s next revolution.  Nature, 445(7126):369–369, 2007.
  * [29] Eric W. Weisstein.  Faulhaber’s formula. From MathWorld—A Wolfram Web Resource.  Last visited on 1/1/2023.


# MOCCA-SURVEY Database II – Properties of Intermediate Mass Black Holes
escaping from star clusters

Konrad Maliszewski,1 Mirek Giersz,2 Dorota Gondek-Rosińska,1 Abbas Askar,3
Arkadiusz Hypki2
1Astronomical Observatory Warsaw University, 00-478 Warsaw, Poland
2Nicolaus Copernicus Astronomical Centre, Polish Academy of Sciences, Warsaw,
Poland
3Lund Observatory, Department of Astronomy, and Theoretical Physics, Lund
University, Box 43, SE-221 00 Lund, Sweden E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In this work we investigate properties of intermediate-mass black holes
(IMBHs) that escape from star clusters due to dynamical interactions. The
studied models were simulated as part of the preliminary second survey carried
out using the MOCCA code (MOCCA-SURVEY Database II), which is based on the
Monte Carlo N-body method and does not include gravitational wave recoil kick
prescriptions of the binary black hole merger product. We have found that
IMBHs are more likely to be formed and ejected in models where both initial
central density and central escape velocities have high values. Most of our
studied objects escape in a binary with another black hole (BH) as their
companion and have masses between $100$ and $\rm 140{\rm~{}M}_{\odot}$.
Escaping IMBHs tend to build-up mass most effectively through repeated mergers
in a binary with BHs due to gravitational wave emission. Binaries play a key
role in their ejection from the system as they allow these massive objects to
gather energy needed for escape. The binaries in which IMBHs escape tend to
have very high binding energy at the time of escape and the last interaction
is strong but does not involve a massive intruder. These IMBHs gain energy
needed to escape the cluster gradually in successive dynamical interactions.
We present specific examples of the history of IMBH formation and escape from
star cluster models. We also discuss the observational implications of our
findings as well as the potential influence of the gravitational wave recoil
kicks on the process.

###### keywords:

stars: black holes – methods: numerical – globular clusters: general

††pubyear: 2021††pagerange: MOCCA-SURVEY Database II – Properties of
Intermediate Mass Black Holes escaping from star clusters–MOCCA-SURVEY
Database II – Properties of Intermediate Mass Black Holes escaping from star
clusters

##  1 Introduction

Black holes (BHs) with masses in the range $\rm
10^{2}-10^{5}{\rm~{}M}_{\odot}$ are known as intermediate-mass black holes
(IMBHs). This range defines BHs that are more massive than stellar-mass BHs
and less massive than supermassive black holes (SMBHs). The LIGO/Virgo
collaboration recently released the catalogue containing the detection of 8
binary black hole (BBH) mergers in which the remnant is an IMBH (Abbott et
al., 2021a). One such event, also known as GW190521, provided first definite
proof for the existence of low-mass IMBHs (Abbott et al., 2020). Until its
discovery, the existence of IMBHs had been heavily debated and concrete
evidence for their presence was lacking although several promising candidates
(see Greene et al., 2020, and references therein) had been observationally
identified through different methods (e.g., Lin et al., 2018a; Chilingarian et
al., 2018; Tanaka, 2018; Takekawa et al., 2020).

The two BHs that merged in the GW190521 event had masses of $\sim
85\rm{\rm~{}M}_{\odot}$ and $\sim 66\rm{\rm~{}M}_{\odot}$. BHs with these
masses are not expected to form through the isolated evolution of single
massive stars. This is because the progenitor stars that should form BHs with
these masses are likely to undergo pair or pulsational pair-instability
supernova (Fryer et al., 2012; Belczynski et al., 2016b). These supernova
events can lead to either enhanced mass loss or the complete explosion of the
star leaving behind no remnant. Therefore, they can prevent the formation of
BHs with masses in the range $\sim 50\ {\rm~{}M}_{\odot}$ to $120\
{\rm~{}M}_{\odot}$, which is referred to as the upper-mass gap (Woosley, 2017;
Spera & Mapelli, 2017). Different formation channels have been proposed to
explain the astrophysical origin of the merging BBHs that are being observed
by the LIGO/Virgo collaboration, including events like GW190521. Among others,
these channels include:

  1. 1.

dynamical BBH formation in dense stellar environments (e.g., Portegies Zwart &
McMillan, 2000; Askar et al., 2017; Banerjee, 2018; Di Carlo et al., 2021).
Repeated BH mergers (e.g., Rodriguez et al., 2019; Antonini et al., 2019;
Samsing & Hotokezaka, 2020; Fragione et al., 2021; Arca Sedda et al., 2021),
BHs forming from dynamical stellar mergers or binary evolution mergers have
been invoked to explain the formation of GW190521 in star clusters (e.g., Di
Carlo et al., 2020; Kremer et al., 2020; Banerjee, 2021; González et al.,
2021).

  2. 2.

isolated binary star evolution Dominik et al. (2013); Belczynski et al.
(2016a); Olejak et al. (2021)

  3. 3.

mergers in triple and quadruple systems (e.g., Antonini & Perets, 2012;
Martinez et al., 2020; Fragione et al., 2020; Vigna-Gómez et al., 2021)

  4. 4.

formation in AGN disks (e.g., McKernan et al., 2012; Bartos et al., 2017; Yang
et al., 2019; Tagawa et al., 2020; Samsing et al., 2020; Tagawa et al., 2021)

The high stellar densities in environments like globular clusters (GCs) can
lead to frequent encounters between stars and binary systems which can result
in the dynamical formation of BBHs. Additionally, dynamical interactions in
these dense environments may also lead to the formation of IMBHs. GCs have
been the subject of extensive theoretical research examining possible IMBH
formation mechanisms and the consequences that its presence can have on the
evolution of the cluster (e.g. Giersz et al., 2015; Morawski et al., 2018;
Arca Sedda et al., 2019; Hong et al., 2020a). In the cores of initially dense
clusters, massive stars may undergo runaway collisions (Portegies Zwart &
McMillan, 2002; Portegies Zwart et al., 2004; Sakurai et al., 2017) that could
lead to the formation of a very massive star larger than $\gtrsim 150\
{\rm~{}M}_{\odot}$ (Mapelli, 2016; Reinoso et al., 2018). In low metallicity
environments, where mass loss due to stellar winds is low, these massive stars
may evolve into IMBH seeds (Spera & Mapelli, 2017; Kremer et al., 2020). It is
also possible to form IMBH seeds in GCs through the gradual growth of stellar-
mass BHs via mergers with other BHs or stars (Miller & Hamilton, 2002; Giersz
et al., 2015; Rizzuto et al., 2021; Di Carlo et al., 2021). It has been shown
that seed IMBHs of $\rm 10^{2}-10^{3}\,\rm M_{\odot}$ can grow to larger
masses through tidal capture and disruption of stars in extremely dense
clusters (Alexander & Bar-Or, 2017; Stone et al., 2017; Sakurai et al., 2019)
with stellar densities $\rm\gtrsim 10^{6}\ {\rm~{}M}_{\odot}\ pc^{-3}$. It has
also been suggested that stellar-mass BHs can grow into IMBHs by accreting gas
in primordial, gas-rich, massive star clusters (Vesperini et al., 2010; Leigh
et al., 2013).

Observationally determining the presence of IMBHs in the crowded centres of
GCs is difficult (Mezcua, 2017; de Vita et al., 2017; Zocchi et al., 2019;
Aros et al., 2020; Häberle et al., 2021). Kinematic observations of the
central regions of many Galactic GCs have identified upper limits on central
IMBH masses (Noyola et al., 2008; Lützgendorf et al., 2013a; Lanzoni et al.,
2013; Häberle et al., 2021). However, the results of these works remain
inconclusive and there is no conclusive evidence for the presence of an IMBH
in the centre of any of the Galactic GCs (Mann et al., 2019; Baumgardt et al.,
2019). Observers have also looked for observational signatures of accreting
IMBHs at centre of Galactic GCs using X-ray and radio observations (Strader et
al., 2012). These studies show that is unlikely for Galactic GCs to be
harbouring IMBHs more massive than a $\sim 1000\ {\rm~{}M}_{\odot}$ (Tremou et
al., 2018). The best observational evidence for IMBHs in GCs come from
extragalactic observations. For instance, Lin et al. (2018b) identified an
X-ray flare in an extragalactic star cluster that was caused by a tidal
disruption event associated with a $\rm 1.75\times 10^{4}{\rm~{}M}_{\odot}\
pc^{-3}$ IMBH (Lin et al., 2020; Wen et al., 2021). Several IMBH candidates
have also been found in central regions of dwarf and low-mass galaxies through
dynamical mass measurements (Nguyen et al., 2019) and electromagnetic
observations (Baldassare et al., 2015). Several IMBH candidates have also been
identified close to the Galactic center through observations of molecular gas
streams (Takekawa et al., 2019; Takekawa et al., 2020). It is possible that
such IMBHs may have originated in disrupted star clusters that sink towards
the Galactic center due to dynamical friction (Arca-Sedda & Gualandris, 2018;
Askar et al., 2021a).

If an IMBH forms early during the cluster evolution, it can co-exist in the
cluster with several hundreds of stellar-mass BHs (Arca Sedda et al., 2019).
It is highly probable that these stellar-mass BHs will undergo strong
dynamical encounters with the IMBH in the cluster center. It is expected that
these encounters are likely to eject our stellar-mass BHs from the cluster
(Lützgendorf et al., 2013b; Hong et al., 2020b). However, it may also be
possible to eject out an IMBH in such interactions, particularly when a star
cluster hosts multiple BHs with masses larger than $100\ {\rm~{}M}_{\odot}$.
Although the possibility for an IMBH to escape from its host star cluster as a
result of dynamical interactions has not been intensively studied, numerical
studies indicate this can occur in realistic star clusters (Di Carlo et al.,
2021). Several works have investigated the ejection of IMBH and their seeds
from star clusters due to gravitational wave recoil kicks following the merger
of two BHs (e.g., Morawski et al., 2018; Rodriguez et al., 2019; Askar et al.,
2021b; Arca Sedda et al., 2021; Rizzuto et al., 2021). The magnitude of these
gravitational wave recoil kicks depends on the anisotropy in the masses and
spins of the merging BHs (González et al., 2007; Holley-Bockelmann et al.,
2008). If these kicks exceed the escape velocity of the cluster than this may
result in the ejection of an IMBH seed and hence lower its retention
probability (Gerosa & Berti, 2019). Such recoil kicks are are more likely to
eject potential IMBH seeds that are less massive than a $100\
{\rm~{}M}_{\odot}$.

In this paper, we investigate the dynamical ejection of IMBHs ($>100\
{\rm~{}M}_{\odot}$) in dense and massive GCs. We present results from star
cluster models simulated as part of the MOCCA-SURVEY Database II. In the
evolution of some of these models, we find that IMBHs can be ejected from the
cluster solely due to dynamical interactions with other objects in the core.
The paper is organized as follows. In the section 2 we will present the
methods used as well as models that were part of the Survey 2. In section 3 we
will discuss if there is any dependency on the initial parameters that were
used in the simulations. The physical parameters of escaping IMBHs will be
shown in section 4. In section 5 we will present what were the mass build-up
events that contributed to the mass of IMBH that was ejected compared to the
ones that stay in the system until the Hubble time. The physical processes
that are responsible for the IMBHs escapes are discussed in section 6. The
most peculiar cases of such escapers are shown in section 7. In section 8 we
try to quantify the effect that the gravitational recoil can have on the
ejection process of the IMBHs. Lastly, in section 9, we sum up the most
important results of the work while discussing main results and drawing
conclusions.

##  2 Models

The numerical simulations of GCs for the purpose of this work were done with
the mocca111http://moccacode.net code (Giersz et al., 2014; Hypki & Giersz,
2013; Giersz, 1998). mocca is one of the most feature rich codes which is able
to simulate GCs with realistic sizes. It is much faster than direct N-body
codes, it follows the cluster evolution very precisely and at the same time it
provides almost the same amount of details of the dynamical and stellar
evolution of all stars in the cluster (Wang et al., 2016). mocca code already
proved its value and versatility in a number of projects to simulate the
properties of realistic GCs (Giersz & Heggie, 2011), the evolution of compact
binaries e.g. BHs binaries (Hong et al., 2020b). mocca code uses sse/bse codes
to advance the stellar evolution of stars and binaries (Hurley et al., 2000,
2002). The close dynamical interactions, performed in N-body way, are done
with fewbody code (Fregeau et al., 2004).

Recently, mocca code was used to create a large number of numerical
simulations for various initial conditions – called mocca surveys. The newest
version of the surveys is MOCCA-SURVEY Database II which an upgrade of the
already published and widely used MOCCA-SURVEY Database I (Askar et al.,
2017).

The newest feature of the mocca code is the ability to follow the dynamical
evolution of the multiple stellar populations (Piotto et al., 2015). The
initial conditions for MOCCA-SURVEY Database II were chosen to test the code
and to test the dynamical mixing between two stellar populations for mostly
very dense models, up to central densities larger than $\rm\sim
10^{7}\>M_{\sun}\ pc^{-3}$. MOCCA-SURVEY Database II currently contains 212
very detailed models which cover a broad range of initial conditions (for
details see Hypki et al. (2021, in prep.)). In this study, the models from the
preliminary MOCCA-SURVEY Database II are used. These models will be extended
in the near future.

Among the many initial parameters for mocca models, the one which has
importance for the topic of this paper is the flag that determines the remnant
mass following neutron star (NS) and BH formation. For MOCCA-SURVEY Database
II the rapid and the delayed supernova models from Fryer et al. (2012) were
applied in order to have some BH population within the first mass gap range
($2-5$ ${\rm~{}M}_{\odot}$). The metallicity for all models is 0.001 (1/20 of
the solar metallicity). Models with various metallicities are planned for the
next stage of the MOCCA-SURVEY Database II. Wind prescriptions were chosen to
be similar to Belczynski et al. (2010). The final mass after mergers between
NSs/BHs with other stars was modified to be more conservative and now only 25
percent of the merging star’s mass is accreted by the compact object. In these
models, BH natal kicks are drawn from the same distribution that is used for
NSs (Hobbs et al. (2005); Maxwellian distribution with $\rm\sigma=265\ km\
s^{-1}$). However, BH kicks are corrected for linear momentum conservation and
thus depend on the final BH mass and the fallback prescription defined by the
supernova model. Fryer et al. (2012); Belczynski et al. (2016c); Banerjee
(2020). Pulsational pair instability (PPISN) and full-fledged pair-instability
supernova (PISN) were chosen to follow Belczynski et al. (2016c). Possibility
for the electron capture and accretion induced SNs are also switched on. The
effect of the tidal disruptions is not included in the simulated models. The
overall comparison of the code with the state-of-the-art N-body code
NBODY6++GPU and the treatment of the stellar winds can be found in Kamlah et
al. (2022).

The best indication on how well does the code work is to compare its output
models with the observed GCs. To do that we used the Baumgardt & Hilker (2018)
catalogue to compare the mass and half-mass radius of real GCs to the ones
produced by the mocca code. The results are presented in Fig. 1 and clearly
show that our models, which are evolved up to 12 Gyrs, have properties that
match the observations. However, it is worth noticing that our models do not
cover the whole range of the observed parameters and thus cannot be treated as
representative of all the clusters in the Universe. In fact, they represent
clusters which are rather dense.

Figure 1: Comparison of the global parameters of GCs produces with the mocca
code (green filled diamonds) and the observed ones from the Baumgardt & Hilker
(2018) catalogue (red empty circles). The MOCCA-SURVEY II models are evolved
to 12 Gyrs.

##  3 Initial model parameters

Firstly, in Table 1 we present initial conditions for the models that
constituted the MOCCA-SURVEY Database II. In order to determine, on which
initial parameters the IMBH formation and its ejection depends, we studied all
212 models available in the MOCCA-SURVEY Database II. We have distinguished
models in which:

  1. 1.

there was at least one IMBH escaper in the model and no IMBH was present in
the cluster at 12 Gyr (49 models)

  2. 2.

an IMBH was present in the cluster at 12 Gyr and there were no escaping IMBHs
during the evolution of the cluster (66 models)

  3. 3.

there was at least one IMBH escaper and an IMBH was present in the cluster at
12 Gyr (20 models).

This way we could differentiate between the IMBH formation inside the cluster
and processes that lead to their escapes. In Table 2, we present the initial
parameters that differ between the simulated star cluster models. For some of
the initial parameters, value ranges have been used to group models and Table
2 provides the total number of models belonging to each group.

We have found that for our models there was no dependency on the initial
number of objects in the cluster for IMBH formation. This is due to the fact
that models in the MOCCA-SURVEY Database II took rather extreme initial
parameters values, their initial central densities were high and clusters were
strongly tidally-underfilling, however, these conditions may correspond to the
initial values of the observed clusters (see Marks et al., 2022).

Both initial central density and initial central escape velocity, as expected
from theoretical considerations, seem to influence the process of IMBH
formation as well as those leading to its ejection from the cluster (see Fig.
2). The higher the values of those parameters are, the higher the probability
for an IMBH to form and stay in the cluster until 12 Gyr. For IMBHs that are
escaping there is a peak for models with the values of logarithm of the
initial central density (${\rm~{}M}_{\odot}\>\textrm{pc}^{-3}$) between 6.5
and 7. This is caused by the fact that for larger densities the pace of
dynamical interactions is higher (due to smaller relaxation time) and thus
IMBHs acquire their mass faster and are becoming less susceptible to be
removed.

Table 1: Initial conditions for the MOCCA-SURVEY Database II models. A two segmented initial mass function (IMF) as given by Kroupa (2001) was used for all models, $f_{\textrm{b}}$ \- binary fraction. If the binary fraction is equal to 0.95 then binary parameters are chosen according to Kroupa (1995), otherwise eccentricity distribution is thermal, mass ratio distribution is uniform and semi-major distribution is uniform in logarithm, between $2(R_{1}+R_{2})$ and 100 AU. Most of the models (147 out of 212) have the $f_{\textrm{b}}=0.95$, the rest are models with $f_{\textrm{b}}=0.1$ (46) or $f_{\textrm{b}}=0.0$ (12). $N$ \- initial number of stars and binaries, $R_{\textrm{t}}$ \- tidal radius, $R_{\textrm{h}}$ \- half-mass radius, $W_{0}$ \- King model parameter, $Z$ \- cluster metallicity. The upper mass limit for second generation stars was set to either $50$ or $150{\rm~{}M}_{\odot}$, the lower mass limit is $0.08{\rm~{}M}_{\odot}$. For each initial number of objects different combinations of parameters are used to generate the initial model. $N$ | $R_{t}\>(\textrm{pc})$ | $R_{\textrm{h}}\>(\textrm{pc})$ | $W_{0}$ | $Z$ | $f_{\textrm{b}}$ | Natal Kicks | # of models
---|---|---|---|---|---|---|---
$6\cdot 10^{5}$ | 60, 120 | 0.6, 1.2, 2.0, 4.0, 6.0 | 6, 8 | 0.001 | 0.1, 0.95 | Fallback | 66
$8\cdot 10^{5}$ | 60, 120 | 0.6, 1.2 | 6, 8 | 0.001 | 0, 0.1, 0.95 | Fallback | 71
$10^{6}$ | 60, 120 | 0.6, 1.2 | 6, 8 | 0.001 | 0, 0.1, 0.95 | Fallback | 75
Table 2: Studied initial parameters of the models which constituted the
survey. $N$ is the initial number of objects in the simulation,
$\rho_{\textrm{c0}}$ is the initial central density expressed in $\rm
M_{\sun}\>\textrm{pc}^{-3}$, $r_{\textrm{t0}}/r_{\textrm{c0}}$ is the
concentration parameter (initial tidal radius to core radius),
$r_{\textrm{t0}}$ is the initial tidal radius in pc, $v_{\textrm{ce0}}$ is the
initial central escape velocity in $\textrm{km}\>\textrm{s}^{-1}$. For every
value/range of the parameter the total number of models belonging to it is
given.

| $N$ | $\log_{10}(\rho_{\textrm{c0}})$
---|---|---
| $6\cdot 10^{5}$ | $8\cdot 10^{5}$ | $10^{6}$ | $[4.2,6.5)$ | $[6.5,7.0)$ | $[7.0,8.4)$
total # of models | 66 | 71 | 75 | 87 | 68 | 57

| log${}_{10}(r_{\textrm{t0}}/r_{\textrm{c0}})$ | $r_{\textrm{t0}}$
---|---|---
| $[1.75,2.5)$ | $[2.5,2.75)$ | $[2.75,3.35)$ | $60$ | $120$
total # of models | 68 | 81 | 63 | 140 | 72

| $v_{\textrm{ce0}}$
---|---
| $[40,90)$ | $[90,110)$ | $[110,140)$ | $[140,195)$
total # of models | 48 | 62 | 54 | 48

Figure 2: Initial central density (as defined in Casertano & Hut, 1985) and
initial central escape velocity plot. Side histograms show the share in the
total number of models belonging to the same value range. Models where only
IMBH escape have been observed are marked with red color (square marks and
dotted histogram). Models with an IMBH present in the cluster at 12 Gyr and no
escaping IMBHs during the evolution of the cluster are marked with green color
(asterisk marks and oblique lines histogram). Those where both escape of an
IMBH and its presence in the system at the Hubble time were observed are
marked with blue color (diamond marks and horizontal lines histogram).

The concentration parameter comes from the King (1966) model and it reflects
the same dependency as the initial central escape velocity and initial density
in the centre, which means that for an IMBH the higher the value of
concentration parameter the more likely it will form. Our models are initially
very dense (see Table 1 and 2) and the process of the core collapse due to the
mass segregation starts within the first few or few hundred Myr. After that
the energy that is generated in dynamical interactions with the massive
binaries starts to control the evolution of the system and the phase of
evolution after the core collapse.

Additionally, we have checked whether IMBH formation or its ejection is
dependent on the initial tidal radius and initial half-mass radius. We found
out that both parameters should be checked on simultaneously and that IMBH
ejection is not dependant on those parameters, but the formation is more
likely to occur in models with higher tidal radius and lower half-mass radius.
Such models can be treated as strongly tidally-underfilling and are well
approximated by isolated system for the extended period of time, which means
that they won’t lose substantial part of their mass over given period of time.

Lastly, we found that our models match the theoretical predictions regarding
the correlation between the central velocity dispersion and the possibility of
IMBH formation due to mergers during few-body interactions (see Miller &
Davies, 2012). Within $66$ models in which an IMBH was present in the system
at Hubble time and no previous IMBH have been observed only $5$ of them
constituted cases where the velocity dispersion values were less than
$40\>\textrm{km}\>\textrm{s}^{-1}$ (minimal value was $27\
\textrm{km}\>\textrm{s}^{-1}$). The rest of them belonged to the range between
$40$ and $87\>\textrm{km}\>\textrm{s}^{-1}$ which is in agreement with
theoretical predictions for the velocity dispersion values in the centre of
the cluster that is required to form a massive BH. The
$40\>\textrm{km}\>\textrm{s}^{-1}$ velocity dispersion threshold is a value
below which heating from binaries prevents full core collapse. It is a result
of a theoretical study which bases on assumptions and approximations regarding
physical processes happening inside the cluster and thus cannot be treated as
exact.

##  4 Escapers parameters

There were in total $66$ escape events which included an IMBH and occurred in
49 models. The maximum of 3 IMBHs escaped from the single cluster. Most of all
the ejection events have happened to the IMBH in a binary with another BH (60
out of 66). Two escapes concerned case in which both binary escapers were
IMBHs.

In Fig. 3 we present the mass ratio of the escaping binaries as a function of
the escaping IMBH mass. Single IMBH escapes and are not present in the plot.
Masses of single escapers do no exceed $150{\rm~{}M}_{\odot}$ whereas those
that correspond to IMBH escapers in a binary can have masses that are twice as
large. The mass ratio distribution of IMBH binary escapers with merger time
less than 12 Gyr has two peaks: at around $0.1$ and at values between
$0.3-0.4$ which corresponds to the companion mass of about $40\
{\rm~{}M}_{\odot}$ and is similar to the mass of the most massive BH that can
be formed via stellar evolution in the MOCCA-SURVEY Database II. Mergers
occurring after 12 Gyr will most likely involve binary systems that have mass
ratios between $0.2-0.4$. For most cases, masses of the escaping IMBHs do not
exceed $300\>{\rm~{}M}_{\odot}$, however, a BH of nearly
$600{\rm~{}M}_{\odot}$ was ejected from the simulation and we present the
exact history of this object in section 7.

Figure 3: Mass ratio as a function of escaping IMBH mass. The histograms (to
the right and above the scatter plot), show the frequency distribution of the
mass ratio and escaping IMBH mass. There were 6 single esacapers present in
our models which are not plotted in this graph. Their masses were between
$102$ and $115{\rm~{}M}_{\odot}$. The red diamond mark represents binaries
which will merge after 12 Gyr since the beginning of simulation, the blue
asterisk mark is for the binaries that are going to merge after 12 Gyr.

Fig. 4 shows the semi-major axis (taken in logarithm) and squared eccentricity
of the IMBHs escaping in a binary. The mean value of eccentricity is equal to
$0.508$ with the dispersion of the mean value of $0.277$. After sufficient
number of interactions, that the binary in a star cluster experiences, the
eccentricity distribution should be thermal (see Heggie, 1975; Geller et al.,
2019, and references therein). Our data seems to be consistent with those
theoretical predictions (mean value of eccentricity for thermal distribution
is $2/3$) assuming the sigma test. This means that binaries that contain an
IMBH have had a large number of interactions in the cluster before they
managed to escape from it and are not primordial. It is also worth to note
that there are no binaries with the low semi-major axis values and high
eccentricities. It comes from the fact that for such objects the time scale
between mergers due to gravitational wave emission is shorter than the mean
time between subsequent dynamical interactions. Most of the studied objects
have the semi-major axis between $30$ to $100\>R_{\sun}$ and their
distribution is not dependant on eccentricity for $\log{a}>1.5$.

Figure 4: Semi-major axis logarithm and squared eccentricity plot of the IMBHs
that are escaping in a binary. There are in total 60 of such cases. On the
sides are plotted histograms with the number of occurrence as the $y$-value.
Mean value of eccentricity equals $0.508$ and its dispersion is $0.277$.

Lastly, we present the velocity at infinity of the escaped object as a
function of the time of the escape (taken in logarithm) in Fig. 5. The
velocity at infinity tends to get smaller as the time of the escape
progresses. For a more dynamically evolved system there are objects with
smaller average masses and thus the dynamical interactions become weaker.
Additionally the half-mass radius is bigger as well as the central density is
smaller, which lowers down the escape velocity. The values of the velocity are
not high and in most cases they don’t exceed
$100\>\textrm{km}\>\textrm{s}^{-1}$. Taking into account the proper motion of
the real GCs (e.g. Dinescu et al., 1999; Bajkova & Bobylev, 2021) we can
conclude that substantial fraction of IMBHs that have escaped from the cluster
is likely to be found in the vicinity of GC tidal tails.

Figure 5: The velocity at infinity of the escaping IMBH (single or in a
binary) as a function of the escape time in Myr taken in logarithm. On the
sides are histograms that show the number of occurrence on the $y$-axis that
correspond to different parameters.

##  5 mass build-up events

In order to determine how IMBHs that are escaping grew in mass we studied all
mass build-up events that such BHs experienced from the time the BH seed have
been formed up to its ejection as IMBH. We have put them into different
categories to identify processes that lead to their mass increase. The result
of the study is presented in Table 3.

In Table 4 we have presented the same results for the IMBH that remained in
the system and evolved until the Hubble time, in order to establish whether
different processes are responsible for their mass build-up.

Table 3: Table showing mass build-up events contribution to the total mass of escaping IMBHs. Dynamical events include: 3b,4b BH:BH – three- or four-body encounter in which merger of two BHs occur, 3b,4b other – three- or four-body encounter in which merger of a BH with object other than a BH (e.g. with a MS star, WD, NS) occurs and direct collisions (only with objects that are not BH). There were no collisions with other BHs. Evolutionary events consist of: BH:BH merger – BBH merger, mass transfer in a binary and BH:other merger – binary merger with object that is not a BH. Both total mass and mean mass are expressed in $M_{\sun}$. |  | total mass | # of events | mean mass
---|---|---|---|---
{314.5mm[  dynamical ] | 3b,4b BH:BH | 484.7 | 11 | 44.07
| 3b,4b other | 115.9 | 48 | 2.41
| other collision | 18.27 | 6 | 3.04
{314.5mm[  evolution ] | BH:BH merger | 4801 | 114 | 42.11
| mass transfer | 38.55 | 10 | 3.85
| BH:other merger | 245.2 | 12 | 20.44
Table 4: Table showing mass buildup events contribution to the total mass of IMBHs that are staying inside the cluster at 12 Gyr. Dynamical events include: 3b,4b BH:BH – three- or four-body encounter in which merger of two BHs occur, 3b,4b other – three- or four-body encounter in which merger of a BH with object other than a BH (e.g. with a MS star, WD, NS) occurs and direct collisions (both with BHs and objects other than a BH). Evolutionary events consist of: BH:BH merger – BBH merger, mass transfer in a binary and BH:other merger – binary merger with object that is not a BH. Both total mass and mean mass are expressed in $M_{\sun}$. |  | total mass | # of events | mean mass
---|---|---|---|---
{414.5mm[  dynamical ] | 3b,4b BH:BH | 1435 | 39 | 36.8
| 3b,4b other | 19881 | 15626 | 1.27
| BH collision | 129751 | 19162 | 6.77
| other collision | 19611 | 28798 | 0.68
{314.5mm[  evolution ] | BH:BH merger | 21049 | 617 | 34.1
| mass transfer | 41 | 906 | 0.05
| BH:other merger | 226 | 308 | 0.73

We have found that IMBHs that are escaping build-up their mass most
effectively through reapeted BBH mergers due to gravitational wave emission.
In contrast those that are present in the system at 12 Gyr acquire their mass
substantially through collisions (mostly with BHs). This is with agreement
with the results from the MOCCA-SURVEY Database I (see Giersz et al., 2015)
that for a massive BH (over $\sim 1000{\rm~{}M}_{\odot}$) the physical
collisions with other objects dominate the mass build-up process of an IMBH.

Firstly, IMBH seed will acquire its mass most effectively by the repeated BBH
mergers. It happens in the very begining of the simulation when there is
sufficient amount of BHs that constitute subsystem. While comparing the mean
masses of the dynamical mergers between two BHs in Table 3 and 4, it is clear
that IMBH escapes happen at the very beginning of the simulation when
evolutional mergers take place with most massive BHs present in the system.
Mean increase of mass through BH mergers is also in agreement with the
theoretical predictions, that firstly more massive BHs will be depleted from
the BH subsystem. In this paper, however, we didn’t study the influence of the
gravitational recoil on the process. Due to anisotropy of masses and spins
such mergers will gain velocity that can be greater than the escape velocity
and in result can be ejected from the cluster. The code itself didn’t account
for the recoil so the situation in real cluster may be different. We studied
the effect in the Section 8.

If the BH is not ejected through dynamical interactions it will continue to
acquire mass through collisions with other objects. It is then harder or even
impossible to remove it from system both due to dynamical interactions or
gravitational recoil kick as it will quickly deplete the population of
stellar-mass BHs in the cluster (see Leigh et al., 2014; Giersz et al., 2015).

##  6 Cause of the escape

The purpose of the paper is to show the possibility of dynamical ejections of
IMBHs from globular clusters as our models did not account for the
gravitational wave recoil kick caused by the asymmetry of masses and spins of
the merging BHs. Studying IMBHs that are escaping from the cluster, we
investigated the masses of interacting objects in the last encounter. In most
cases the last interaction will be the one responsible for the ejection of our
studied objects, which is connected with the fact that after a strong
interaction binding energy will raise to the value which is sufficient for the
object to escape from the cluster (it has to be greater than the tidal energy
$-1.5GM/R_{\textrm{t}}$, where $M$ is the current cluster mass). In most cases
the binding energy inside the cluster is positive, which means that an object
is not bound to it.

Subsequently, while studying the processes that can lead to an IMBH escape the
key information is what are the masses of objects involved in the last
encounter. In Fig. 6 we present both masses of an intruder and escaping
object. In this case we define escaping object as a single IMBH or a binary
involving an IMBH. Intruder is a single object or a binary that interacts with
the escaping object in the last encounter. Two peaks of intruder mass at
around $40$ and $80$ solar masses is caused by the fact that the most massive
BH that can be formed in single star evolution is around $40{\rm~{}M}_{\odot}$
in MOCCA-SURVEY Database II. This means that the intruder, in most cases, is a
massive BH that did not further grow its mass, a binary consisting of two such
BHs or their merger product.

Figure 6: Escaping IMBH mass and intruder mass plot. Both masses are expressed
in $M_{\sun}$. Escaping object is defined as single IMBH or a binary that
involves an IMBH. Intruder is a single object (red color and diamond mark) or
a binary (blue color or square mark) that interacts with the escaping object
in the last encounter.

It is worth noticing that the intruder mass does not have to be high for an
IMBH to escape from the cluster. In fact, in most cases the mass of the
intruder is 5-10 times smaller than the escaping object mass. This suggests
that the binary increases its energy gradually and the final interaction is
strong not due to the fact that the intruders mass is high but due to its
really close approach to the very hard IMBH binary.

To support this claim, we have plotted the binding energy as the function of
the binding energy change to mean kinetic energy of objects in the core for
all 60 IMBHs that were escaping in a binary. Plots were obtained while
studying the history of the binary that have escaped from the time it was
formed up to the last interaction it took part in. We have distinguished the
last interaction that they have experienced in the cluster and those that have
happened before to compare them. Results are presented in Fig. 7 and clearly
show that the energy changes in last interaction were comparable to the
binding energy of the IMBH binary. Studied binaries are also very hard (their
binding energy is over thousand times bigger than the mean kinetic energy in
the core), which means that they must have encountered hundreds different
interactions before being ejected from the cluster. This seems to support the
picture that was resulting from the study of the masses of intruders. One
interaction seems to stand out of all others and has the binding energy change
at the same level as the mean kinetic energy in the core. This peculiar
example is presented in section 7.1. We finally studied the distance of
closest approach for every last interaction that was a fly-by (28 cases) or
exchange (15 cases) in a three-body encounter. We have found that in most
cases, this distance is about hundredth part of a solar radius (calculated
from the definition of the impact parameter:
$b^{2}=r_{p}^{2}(1+2G(m_{1}+m_{2}+m_{3})/(u^{2}r_{p}))$, where $r_{p}$ is the
pericentre with respect to the centre of mass of the binary, $m_{i}$ are the
masses of interacting objects and $u$ is the intruders relative velocity with
respect to the centre of mass of the binary) and the intruder is always a BH.
It means that these are very strong interactions that can significantly affect
the orbit of the IMBH companion.

Figure 7: Binding energy of the escaping binary as the function of the binding
energy change (both scaled to the mean kinetic energy) in the interaction.
Plot was made for objects that were escaping in a binary and included only
history of such objects from the time they were formed up to the last
interaction. We have distinguished the last interaction with purple color
(marked with dot), all previous are marked with green color (or square mark).

##  7 IMBH escape examples

###  7.1 600 $M_{\sun}$ IMBH

Fig. 8 represents the history of nearly $600{\rm~{}M}_{\odot}$ IMBH that
escapes from the cluster. The initial parameters of the model are:

  * •

central density – $6.94\cdot 10^{6}\>\rm M_{\sun}\>\textrm{pc}^{-3}$

  * •

central escape velocity – $122\>\textrm{km}\>\textrm{s}^{-1}$

  * •

number of objects – $6\cdot 10^{5}$

  * •

tidal radius – $60\>\rm{pc}$

The initial parameters in this specific case correspond to the values which
are on lower end of the values of both central escape velocity and central
density.

Our given BH forms at $3.7$ Myr with $201{\rm~{}M}_{\odot}$ as a result of a
single star evolution of $232{\rm~{}M}_{\odot}$ that has experienced multiple
stellar mergers before. Such massive stars may form through stellar mergers or
collisions as described in González et al. (2021). In the next phase it builds
the mass through 19 different events that were presented graphically and are
not in a chronological order. The biggest mass contribution in the final mass
of a BH hold the BBH mergers that in 7 mergers acquired mass of nearly
$400{\rm~{}M}_{\odot}$. Our IMBH forms its final binary at $881$ Myr with
semi-major axis of $7296\>R_{\sun}$ and $0.99$ eccentricity. Ensuing, it
starts to harden due to interactions with other objects in the centre of GC
and prior to the interaction responsible for its ejection, semi-major axis is
lowered by the factor of $4$ which means that the binary binding energy
increase was of the same level. It is worth mentioning that the interaction
responsible for the ejection was not strong enough to remove the
$600{\rm~{}M}_{\odot}$ IMBH from the simulation immediately. The energy after
the interaction was not sufficient for the object to escape, however, due to a
large number of previous dynamical interactions the orbit of the object was
highly eccentric. As a result, there was no need for strong interaction to
eject it from the cluster because its binding energy was already greater than
the tidal energy for the cluster.

Figure 8: Diagram representing history of two BHs (blue color – IMBH and red –
regular BH) escaping in a binary. BHs are marked with a dot whereas star
represents any other object. A direct transition from one state of the system
to another is marked with a solid line while dashed line is a non direct
transition. Abbreviations used in the description of events: 3b-m –
coalescence due to three-body interaction, 2b-gw – BBH merger due to the
emission of gravitational waves, 4b-F – four-body flyby, 4b-FD – four-body
flyby with a disruption of a binary, 3b-F – three-body flyby.

###  7.2 Two IMBHs binary

Another interesting example of the escaping IMBHs are the ones that are
escaping together as a binary. In Fig. 9, we present one out of two such cases
that occured in our simulated models. The initial parameters of the model are:

  * •

central density – $2.41\cdot 10^{7}\>\rm M_{\sun}\>\textrm{pc}^{-3}$

  * •

central escape velocity – $169\>\textrm{km}\>\textrm{s}^{-1}$

  * •

number of objects – $10^{6}$

  * •

tidal radius – $120\>\rm{pc}$

Both of the BHs form roughly at the same time in the evolution of the binary
star system. One of the BHs forms with $89{\rm~{}M}_{\odot}$ and builds its
mass in two events. First it merges with another BH and a star in a three-body
interaction, after which it forms a binary again and due to three-body
interaction a star merges directly into it. The second BH forms with
$80{\rm~{}M}_{\odot}$ and gains its mass solely due to BBH merger. They both
form a binary together at $146$ Myr and acquire energy needed to escape by the
interactions in the centre of GC with other objects. The one responsible for
the ejection of a binary is a flyby of $53{\rm~{}M}_{\odot}$ intruder.

Figure 9: Diagram representing history of two IMBHs escaping in a binary. BHs
are marked with a dot whereas star represents any other object. A direct
transition from one state of the system to another is marked with a solid line
while dashed line is a non direct transition. Abbreviations used in the
description of events: 3b-m – coalescence due to three-body interaction, 2b-gw
– BBH merger due to the emission of gravitational waves, 4b-F – four-body
flyby, 4b-FD – four-body flyby with a disruption of a binary, 3b-F – three-
body flyby.

##  8 Gravitational kicks

In this section we will present our predictions regarding the influence of the
gravitational kicks on our results. We tried to quantify it by calculating the
retention probability for each BBH merger that involved our studied IMBHs. The
gravitational kick prescriptions were taken from Morawski et al. (2018). These
are the steps we used to calculate the retention probability for each merger:

  1. 1.

First the initial spin model had to be chosen between:

    1. (a)

Constant spin model with the initial value of spin $0.1$. The assumed value is
motivated by the gravitational wave observations (Abbott et al., 2021b; Nitz
et al., 2021b) and is thought to represent the small value of the observed BBH
mergers components spin.

    2. (b)

Random spin model – Gaussian distribution of the spin centered around $0.5$
and a cutoff at $0.5/3$. The purpose of this model is to study the broader
range of the initial spin values and it covers nearly the whole range of the
spin magnitude.

  2. 2.

Next we drew both spin orientations $10^{5}$ times by drawing them uniformly
on a sphere of a radius corresponding to the magnitude of the spin. For each
spin configuration we calculated the recoil velocity of the merger product and
compared it with the GCs central escape velocity at that time. Nearly all of
our studied mergers were the case of the dynamically formed BBH systems (with
only 1 exception of the primordial one) and is the reason for the assumed spin
orientation treatment. The number of cases in which the recoil velocity was
smaller than the central escape velocities divided by the total number of
configurations is the retention probability for the BBH merger product.

  3. 3.

Previous point was repeated 10 times and the average as well as dispersion of
the results had been calculated. From now on we will refer to the mean of the
retention probabilities as being the retention probability of the merger
product.

  4. 4.

For each merger chain in the studied IMBHs history, we calculated the overall
retention probability of the whole chain of mergers by multiplying
probabilities for each node (merger event). These overall probabilities,
marked with a number corresponding to the length of the merger chain, are
shown in Fig. 10 (constant low spin model) and Fig. 11 (random spin model).

For the constant (and low) initial spin model the overall cumulative retention
probability is higher for IMBHs that have encountered 1 BBH merger. This comes
from the fact that the first BBH merger involves BHs of a similar mass and
their mass ratio is close to 1. The recoil velocity of the merger product for
low spins (both spins having same magnitude) and high mass ratio is highly
dependant on the spin orientations and its dispersion is large. If on the
contrary the spins differ and not necessarily take low values the recoil
velocity is very high compared to the central escape velocity and thus the
retention probability will be low.

It is also worth noticing that not the first merger but the second one is the
one that will most likely break the merger chain before the BH will gather
enough mass to become IMBH. It comes from the fact that the first merger
involves components of similar mass and thus their mass ratio is close to 1.
After merging the outcome mass will be more less twice higher and it will
experience higher recoil velocity due to higher mass asymmetry (there is a
peak in the mean recoil velocity for $q\approx 0.6$ – see Morawski et al.,
2018) during the next merger, with different BH.

Nevertheless, our post processing study shows that there are still some BHs
that will form IMBHs, which will later on be ejected from the cluster solely
due to dynamical interactions inside it. Depending on the assumed initial spin
model there are about 7 – 15 BHs (10 – 22% of all IMBH escapers) that will be
removed from the cluster in the way presented in this paper.

Figure 10: The overall (cumulative) retention probability for each of the IMBH
escapers. The probabilities are marked with a number corresponding to the
length of the merger chain that given IMBH was part of. The initial spin model
is the constant spin model. In the plot are present only merger chains for
which the overall retention probability was greater than 0. These include: 9
chains with length of 1, 15 chains with length of 2, 4 chains with length of
3, 6 chains with length of 4, 2 chains with length of 5, 1 chain with length
of 7. Figure 11: Similar to Fig. 10, we show the overall (cumulative)
retention probability for each of the IMBH escapers. The probabilities are
marked with a number corresponding to the length of the merger chain that
given IMBH was part of. The initial spin model is the random spin model. In
the plot are present only merger chains for which the overall retention
probability was greater than 0. These include: 2 chains with length of 1, 14
chains with length of 2, 4 chains with length of 3, 6 chains with length of 4,
2 chains with length of 5, 1 chain with length of 7.

##  9 Discussion and conclusions

In this work we explored the possibility of IMBHs being ejected from the star
cluster solely due to dynamical interactions. Using a set of 212 dense and
strongly tidally-underfilling star cluster simulated as part of the
preliminary MOCCA-SURVEY Database II, we studied how different initial
parameters affect the process. We found that IMBHs are most likely to be
ejected in models with initial central densities between $\sim 3\cdot
10^{6}-10^{7}{\rm~{}M}_{\odot}\ \rm{pc}^{-3}$. For larger densities, the rate
of dynamical interactions is higher and thus IMBHs acquire their mass faster
and are less susceptible to be removed from the cluster. On the other hand,
for lower densities, dynamical encounters are less frequent and such models
are likely to never form an IMBH (Hong et al., 2020b).

In total we observed 68 IMBH escapers coming from our models and most of them
have a mass lower than $200{\rm~{}M}_{\odot}$. Due to the fact that a
significant number of our models had large initial central densities and
additionally the code itself did not account for the gravitational recoil
kicks of the merging BBHs, our results must be treated as the upper limit
estimates of the process. Moreover, our models represent GCs which are rather
dense and in many observed clusters the initial central densities may be
lower, which will result in a lower probability of forming an IMBH and a
smaller number of escapers. We estimated that around 10 - 22 % of the escapers
present in the Survey will be subject to this dynamical ejection scenario,
whereas the rest will be removed from the cluster due to gravitational kicks
from the merger with another BH. Escaping IMBHs acquire their mass most
effectively through BBH mergers and thus are likely to be removed before they
reach $100\>M_{\sun}$ due to gravitational wave recoil kicks (Gerosa & Berti,
2019). One of the concerns regarding gravitational kicks and their influence
on the cluster evolution could be that the BH population is depleted and thus
will not allow the IMBHs formation or that there will not be enough dynamical
interactions with massive objects to eject them from GCs as we presented here.
The number of BH:BH mergers is however small (order of tens) compared to the
number of unique BHs in our models (order of thousands) and will not strongly
affect the clusters evolution. Even when the merger product remains in the
cluster its orbit can be changed into more elongated one. This does not
necessary imply that it will no longer be subject to dynamical interactions in
the core as its barycenter will be positioned there and every orbital period
the BH will pass through the dense core where the probability of dynamical
interactions is substantial due to its large mass.

The tidal disruptions were not considered in the code, however, the effect of
the mass increase for dynamically escaping IMBHs can be neglected. Roche
distance is (for mean mass escaping IMBH and MS star of 1 stellar mass):
$d=R_{\sun}(M_{IMBH}/M_{\sun})^{1/3}\approx 4.9R_{\sun}$. In the code the
threshold for calculating the collision probability is $R_{\sun}$. This means
that the mass increase effect would be even smaller than that and combined
with the picture resulting from Table 3, we conclude that this additional
contribution to the total mass would still be insignificant compared to the
mass gained through BH:BH mergers.

In the post processing of the data from our simulated stellar cluster models,
we calculated the retention probability of merged BHs assuming different
values of initial spin magnitudes and orientations. The spins were drawn from
a random uniform distribution and we took $10^{5}$ different orientations. For
every set of spin orientations, kick velocities were calculated from
prescriptions taken from Morawski et al. (2018). We found that if we assume
the initial spin magnitudes to be low (we assumed it to be 0.1, which is in
fact consistent with observational evidence - Abbott et al. (2021c); Nitz et
al. (2021a)), the retention probability becomes relatively high (around 0.4)
for mass ratios close to 1, which corresponds to the first merger. This is
caused by the fact that for low spin model, the kick velocity has its peak at
around $q=0.4$ and is approximately constant for higher $q$ values. In
addition to this, for higher $q$ the spin orientation starts to play key role
for the kick velocity as the mean kick velocity dispersion grows. The problem
of retaining the merger product is more complicated for repeated mergers. On
one hand, each merger will gradually increase the spin magnitude and kick
velocities will become really high. On the other, the build-up of IMBH mass
will gradually lower down the value of $q$. It is also worth noting that the
mass build-up of an IMBH is connected with mass accretion in collisions and
with mergers with stars, which will result in an increase of the spin of a
given IMBH.

Binaries play key role in dynamical processes that lead to an IMBH ejection as
they are able to acquire substantial amount of energy needed for those massive
objects to escape from star cluster. The reason of escape is very close
approach (hundredth parts of solar radii) of the intruder that can
significantly affect the orbit of the IMBH companion.

Finally, we calculated the merger rate for IMBHs that have escaped in a binary
(merger rate outside the cluster for IMBH binary)
$\mathcal{R}=0.7\>\textrm{Gpc}^{-3}\>\textrm{yr}^{-1}$. This value is the
upper limit for merger rate estimation as our models are very dense and the
effect of the gravitational kick due to BBH merger has not been accounted for.

## Acknowledgements

MG, AH and AA were partially supported by the Polish National Science Center
(NCN) through the grant UMO-2016/23/B/ST9/02732. AH is also supported by the
Polish National Science Center grant Maestro 2018/30/A/ST9/00050. DGR and AA
acknowledge support through the COST action CA16104. DGR was partially
supported by the NCN grant UMO-2017/26/M/ST9/00978 and by POMOST/2012-6/11
from the Program of Foundation for Polish Science co-financed by the European
Union within the European Regional Development Fund. AA is supported by the
Swedish Research Council through the grant 2017-04217.

KM would like to express his gratitude towards MG whose patient supervision
has made publishing this work possible.

## Data Availability

The data underlying this article will be shared on reasonable request to MG or
AH.

## References

  * Abbott et al. (2020) Abbott R., et al., 2020, Phys. Rev. Lett., 125, 101102
  * Abbott et al. (2021a) Abbott R., et al., 2021a, arXiv e-prints, p. arXiv:2111.03606
  * Abbott et al. (2021b) Abbott R., et al., 2021b, Physical Review X, 11, 021053
  * Abbott et al. (2021c) Abbott R., et al., 2021c, Physical Review X, 11, 021053
  * Alexander & Bar-Or (2017) Alexander T., Bar-Or B., 2017, Nature Astronomy, 1, 0147
  * Antonini & Perets (2012) Antonini F., Perets H. B., 2012, ApJ, 757, 27
  * Antonini et al. (2019) Antonini F., Gieles M., Gualandris A., 2019, MNRAS, 486, 5008
  * Arca-Sedda & Gualandris (2018) Arca-Sedda M., Gualandris A., 2018, MNRAS, 477, 4423
  * Arca Sedda et al. (2019) Arca Sedda M., Askar A., Giersz M., 2019, arXiv e-prints, p. arXiv:1905.00902
  * 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
  * Aros et al. (2020) Aros F. I., Sippel A. C., Mastrobuono-Battisti A., Askar A., Bianchini P., van de Ven G., 2020, MNRAS, 499, 4646
  * Askar et al. (2017) Askar A., Szkudlarek M., Gondek-Rosińska D., Giersz M., Bulik T., 2017, MNRAS, 464, L36
  * Askar et al. (2021a) Askar A., Davies M. B., Church R. P., 2021a, arXiv e-prints, p. arXiv:2107.10862
  * Askar et al. (2021b) Askar A., Davies M. B., Church R. P., 2021b, MNRAS, 502, 2682
  * Bajkova & Bobylev (2021) Bajkova A. T., Bobylev V. V., 2021, arXiv e-prints, p. arXiv:2108.08507
  * Baldassare et al. (2015) Baldassare V. F., Reines A. E., Gallo E., Greene J. E., 2015, ApJ, 809, L14
  * Banerjee (2018) Banerjee S., 2018, MNRAS, 473, 909
  * Banerjee (2020) Banerjee S., 2020, Monthly Notices of the Royal Astronomical Society, 500, 3002–3026
  * Banerjee (2021) Banerjee S., 2021, arXiv e-prints, p. arXiv:2109.14612
  * Bartos et al. (2017) Bartos I., Kocsis B., Haiman Z., Márka S., 2017, ApJ, 835, 165
  * Baumgardt & Hilker (2018) Baumgardt H., Hilker M., 2018, MNRAS, 478, 1520
  * Baumgardt et al. (2019) Baumgardt H., et al., 2019, MNRAS, 488, 5340
  * Belczynski et al. (2010) Belczynski K., Bulik T., Fryer C. L., Ruiter A., Valsecchi F., Vink J. S., Hurley J. R., 2010, The Astrophysical Journal, 714, 1217–1226
  * Belczynski et al. (2016a) Belczynski K., Holz D. E., Bulik T., O’Shaughnessy R., 2016a, Nature, 534, 512
  * Belczynski et al. (2016b) Belczynski K., et al., 2016b, A&A, 594, A97
  * Belczynski et al. (2016c) Belczynski K., et al., 2016c, A&A, 594, A97
  * Casertano & Hut (1985) Casertano S., Hut P., 1985, ApJ, 298, 80
  * Chilingarian et al. (2018) Chilingarian I. V., Katkov I. Y., Zolotukhin I. Y., Grishin K. A., Beletsky Y., Boutsia K., Osip D. J., 2018, ApJ, 863, 1
  * Di Carlo et al. (2020) Di Carlo U. N., Mapelli M., Bouffanais Y., Giacobbo N., Santoliquido F., Bressan A., Spera M., Haardt F., 2020, MNRAS, 497, 1043
  * Di Carlo et al. (2021) Di Carlo U. N., et al., 2021, MNRAS, 507, 5132
  * Dinescu et al. (1999) Dinescu D. I., Girard T. M., van Altena W. F., 1999, The Astronomical Journal, 117, 1792
  * Dominik et al. (2013) Dominik M., Belczynski K., Fryer C., Holz D. E., Berti E., Bulik T., Mandel I., O’Shaughnessy R., 2013, ApJ, 779, 72
  * Fragione et al. (2020) Fragione G., et al., 2020, ApJ, 900, 16
  * Fragione et al. (2021) Fragione G., Kocsis B., Rasio F. A., Silk J., 2021, arXiv e-prints, p. arXiv:2107.04639
  * Fregeau et al. (2004) Fregeau J. M., Cheung P., Portegies Zwart S. F., Rasio F. A., 2004, MNRAS, 352, 1
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91
  * Geller et al. (2019) Geller A. M., Leigh N. W. C., Giersz M., Kremer K., Rasio F. A., 2019, The Astrophysical Journal, 872, 165
  * Gerosa & Berti (2019) Gerosa D., Berti E., 2019, Phys. Rev. D, 100, 041301
  * Giersz (1998) Giersz M., 1998, MNRAS, 298, 1239
  * Giersz & Heggie (2011) Giersz M., Heggie D. C., 2011, MNRAS, 410, 2698
  * Giersz et al. (2014) Giersz M., Leigh N., Marks M., Hypki A., Askar A., 2014, preprint,  (arXiv:1411.7603)
  * Giersz et al. (2015) Giersz M., Leigh N., Hypki A., Lützgendorf N., Askar A., 2015, MNRAS, 454, 3150
  * González et al. (2007) González J. A., Sperhake U., Brügmann B., Hannam M., Husa S., 2007, Phys. Rev. Lett., 98, 091101
  * González et al. (2021) González E., Kremer K., Chatterjee S., Fragione G., Rodriguez C. L., Weatherford N. C., Ye C. S., Rasio F. A., 2021, ApJ, 908, L29
  * Greene et al. (2020) Greene J. E., Strader J., Ho L. C., 2020, ARA&A, 58, 257
  * Häberle et al. (2021) Häberle M., et al., 2021, MNRAS, 503, 1490
  * Heggie (1975) Heggie D. C., 1975, MNRAS, 173, 729
  * Hobbs et al. (2005) Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS, 360, 974
  * Holley-Bockelmann et al. (2008) Holley-Bockelmann K., Gültekin K., Shoemaker D., Yunes N., 2008, ApJ, 686, 829
  * Hong et al. (2020a) Hong J., Askar A., Giersz M., Hypki A., Yoon S.-J., 2020a, MNRAS, 498, 4287
  * Hong et al. (2020b) Hong J., Askar A., Giersz M., Hypki A., Yoon S.-J., 2020b, MNRAS, 498, 4287
  * Hurley et al. (2000) Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543
  * Hurley et al. (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897
  * Hypki & Giersz (2013) Hypki A., Giersz M., 2013, MNRAS, 429, 1221
  * Hypki et al. (2021) Hypki A., Giersz M., Vesperini E., Hong J., Leveque A., Askar A., Belloni D., Otulakowska-Hypka M., 2021
  * Kamlah et al. (2022) Kamlah A. W. H., et al., 2022, MNRAS, 511, 4060
  * King (1966) King I. R., 1966, AJ, 71, 64
  * Kremer et al. (2020) Kremer K., et al., 2020, ApJ, 903, 45
  * Kroupa (1995) Kroupa P., 1995, MNRAS, 277, 1507
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Lanzoni et al. (2013) Lanzoni B., et al., 2013, ApJ, 769, 107
  * Leigh et al. (2013) Leigh N. W. C., Böker T., Maccarone T. J., Perets H. B., 2013, MNRAS, 429, 2997
  * Leigh et al. (2014) Leigh N. W. C., Lützgendorf N., Geller A. M., Maccarone T. J., Heinke C., Sesana A., 2014, MNRAS, 444, 29
  * Lin et al. (2018a) Lin D., et al., 2018a, Nature Astronomy, 2, 656
  * Lin et al. (2018b) Lin D., et al., 2018b, Nature Astronomy, 2, 656
  * Lin et al. (2020) Lin D., et al., 2020, ApJ, 892, L25
  * Lützgendorf et al. (2013a) Lützgendorf N., et al., 2013a, A&A, 552, A49
  * Lützgendorf et al. (2013b) Lützgendorf N., Baumgardt H., Kruijssen J. M. D., 2013b, A&A, 558, A117
  * Mann et al. (2019) Mann C. R., et al., 2019, ApJ, 875, 1
  * Mapelli (2016) Mapelli M., 2016, MNRAS, 459, 3432
  * Marks et al. (2022) Marks M., Kroupa P., Dabringhausen J., 2022, A&A, 659, A96
  * Martinez et al. (2020) Martinez M. A. S., et al., 2020, ApJ, 903, 67
  * McKernan et al. (2012) McKernan B., Ford K. E. S., Lyra W., Perets H. B., 2012, MNRAS, 425, 460
  * Mezcua (2017) Mezcua M., 2017, International Journal of Modern Physics D, 26, 1730021
  * Miller & Davies (2012) Miller M. C., Davies M. B., 2012, ApJ, 755, 81
  * Miller & Hamilton (2002) Miller M. C., Hamilton D. P., 2002, MNRAS, 330, 232
  * Morawski et al. (2018) Morawski J., Giersz M., Askar A., Belczynski K., 2018, MNRAS, 481, 2168
  * Nguyen et al. (2019) Nguyen D. D., et al., 2019, ApJ, 872, 104
  * Nitz et al. (2021a) Nitz A. H., Capano C. D., Kumar S., Wang Y.-F., Kastha S., Schäfer M., Dhurkunde R., Cabero M., 2021a, arXiv e-prints, p. arXiv:2105.09151
  * Nitz et al. (2021b) Nitz A. H., Capano C. D., Kumar S., Wang Y.-F., Kastha S., Schäfer M., Dhurkunde R., Cabero M., 2021b, ApJ, 922, 76
  * Noyola et al. (2008) Noyola E., Gebhardt K., Bergmann M., 2008, ApJ, 676, 1008
  * Olejak et al. (2021) Olejak A., Belczynski K., Ivanova N., 2021, A&A, 651, A100
  * Piotto et al. (2015) Piotto G., et al., 2015, AJ, 149, 91
  * Portegies Zwart & McMillan (2000) Portegies Zwart S. F., McMillan S. L. W., 2000, ApJ, 528, L17
  * Portegies Zwart & McMillan (2002) Portegies Zwart S. F., McMillan S. L. W., 2002, ApJ, 576, 899
  * Portegies Zwart et al. (2004) Portegies Zwart S. F., Baumgardt H., Hut P., Makino J., McMillan S. L. W., 2004, Nature, 428, 724
  * Reinoso et al. (2018) Reinoso B., Schleicher D. R. G., Fellhauer M., Klessen R. S., Boekholt T. C. N., 2018, A&A, 614, A14
  * Rizzuto et al. (2021) Rizzuto F. P., et al., 2021, MNRAS, 501, 5257
  * Rodriguez et al. (2019) Rodriguez C. L., Zevin M., Amaro-Seoane P., Chatterjee S., Kremer K., Rasio F. A., Ye C. S., 2019, Phys. Rev. D, 100, 043027
  * Sakurai et al. (2017) Sakurai Y., Yoshida N., Fujii M. S., Hirano S., 2017, MNRAS, 472, 1677
  * Sakurai et al. (2019) Sakurai Y., Yoshida N., Fujii M. S., 2019, MNRAS, 484, 4665
  * Samsing & Hotokezaka (2020) Samsing J., Hotokezaka K., 2020, arXiv e-prints, p. arXiv:2006.09744
  * Samsing et al. (2020) Samsing J., et al., 2020, arXiv e-prints, p. arXiv:2010.09765
  * Spera & Mapelli (2017) Spera M., Mapelli M., 2017, MNRAS, 470, 4739
  * Stone et al. (2017) Stone N. C., Küpper A. H. W., Ostriker J. P., 2017, MNRAS, 467, 4180
  * Strader et al. (2012) Strader J., Chomiuk L., Maccarone T. J., Miller-Jones J. C. A., Seth A. C., Heinke C. O., Sivakoff G. R., 2012, ApJ, 750, L27
  * Tagawa et al. (2020) Tagawa H., Haiman Z., Kocsis B., 2020, ApJ, 898, 25
  * Tagawa et al. (2021) Tagawa H., Kocsis B., Haiman Z., Bartos I., Omukai K., Samsing J., 2021, ApJ, 908, 194
  * Takekawa et al. (2019) Takekawa S., Oka T., Tokuyama S., Tanabe K., Iwata Y., Tsujimoto S., Nomura M., Shibuya Y., 2019, PASJ, 71, S21
  * Takekawa et al. (2020) Takekawa S., Oka T., Iwata Y., Tsujimoto S., Nomura M., 2020, ApJ, 890, 167
  * Tanaka (2018) Tanaka K., 2018, ApJ, 859, 86
  * Tremou et al. (2018) Tremou E., et al., 2018, ApJ, 862, 16
  * Vesperini et al. (2010) Vesperini E., McMillan S. L. W., D’Ercole A., D’Antona F., 2010, ApJ, 713, L41
  * Vigna-Gómez et al. (2021) Vigna-Gómez A., Toonen S., Ramirez-Ruiz E., Leigh N. W. C., Riley J., Haster C.-J., 2021, ApJ, 907, L19
  * Wang et al. (2016) Wang L., et al., 2016, MNRAS, 458, 1450
  * Wen et al. (2021) Wen S., Jonker P. G., Stone N. C., Zabludoff A. I., 2021, ApJ, 918, 46
  * Woosley (2017) Woosley S. E., 2017, ApJ, 836, 244
  * Yang et al. (2019) Yang Y., Bartos I., Haiman Z., Kocsis B., Márka Z., Stone N. C., Márka S., 2019, ApJ, 876, 122
  * Zocchi et al. (2019) Zocchi A., Gieles M., Hénault-Brunet V., 2019, MNRAS, 482, 4713
  * de Vita et al. (2017) de Vita R., Trenti M., Bianchini P., Askar A., Giersz M., van de Ven G., 2017, MNRAS, 467, 4057