Computational methods for collisional stellar systems

Dense star clusters are spectacular self-gravitating stellar systems in our Galaxy and across the Universe—in many respects. They populate disks and spheroids of galaxies as well as almost every galactic center. In massive elliptical galaxies nuclear clusters harbor supermassive black holes, which might influence the evolution of their host galaxies as a whole. The evolution of dense star clusters is not only governed by the aging of their stellar populations and simple Newtonian dynamics. For increasing particle number, unique gravitational effects of collisional many-body systems begin to dominate the early cluster evolution. As a result, stellar densities become so high that stars can interact and collide, stellar evolution and binary stars change the dynamical evolution, black holes can accumulate in their centers and merge with relativistic effects becoming important. Recent high-resolution imaging has revealed even more complex structural properties with respect to stellar populations, binary fractions and compact objects as well as—the still controversial—existence of intermediate mass black holes in clusters of intermediate mass. Dense star clusters therefore are the ideal laboratory for the concomitant study of stellar evolution and Newtonian as well as relativistic dynamics. Not only the formation and disruption of dense star clusters has to be considered but also their galactic environments in terms of initial conditions as well as their impact on galactic evolution. This review deals with the specific computational challenges for modelling dense, gravothermal star clusters.

Abstract Dense star clusters are spectacular self-gravitating stellar systems in our Galaxy and across the Universe -in many respects. They populate disks and spheroids of galaxies as well as almost every galactic center. In massive elliptical galaxies nuclear clusters harbor supermassive black holes, which might influence the evolution of their host galaxies as a whole. The evolution of dense star clusters is not only governed by the aging of their stellar populations and simple Newtonian dynamics. For increasing particle number, unique gravitational effects of collisional many-body systems begin to dominate the early cluster evolution. As a result, stellar densities become so high that stars can interact and collide, stellar evolution and binary stars change the dynamical evolution, black holes can accumulate in their centers and merge with relativistic effects becoming important. Recent high-resolution imaging has revealed even more complex structural properties with respect to stellar populations, binary fractions and compact objects as well as -the still controversial -existence of intermediate mass black holes in clusters of intermediate mass. Dense star clusters therefore are the ideal laboratory for the concomitant study of stellar evolution and Newtonian as well as relativistic dynamics. Not only the formation and disruption of dense star clusters has to be considered but also their galactic environments in terms of initial conditions as well as their impact on galactic evolution. This review deals with the specific computational challenges for modelling dense, gravothermal star clusters.
Keywords Numerical Methods · Star Clusters · Stellar Evolution · Direct N -Body Simulation

Preface 1 Astrophysical Introduction
Stars play a fundamental role in astronomy -a large piece of information available to us about the Universe we inhabit comes from stars. Coeval associations of stars, also called star clusters, are the birth place of most if not all stars. Star clusters form, evolve and "die" by dissolution all across the cosmic time, which is covered in the excellent review "Star Clusters across Cosmic Time" that focuses on the interplay between observations and theoretical knowledge [1]; our review, on the other hand, focuses more on the computational challenges to model a special class of dense star clusters. The term "dense" is not very well defined. We use it here in the sense that both the system should be "gravothermal" and that during at least some phases of the evolution close encounters or even direct collisions between stars or binaries occur. This definition constrains us to globular and young dense star clusters, as well as nuclear star clusters (NSCs). On nuclear star clusters there is another nice review [2]; nuclear star clusters are not in the focus of this work; however, if special computational issues have to be taken into account for nuclear star clusters we may elaborate on them here. Globular star clusters (GCs) are thought to be the oldest objects in our Galaxy, their age covering a large fraction of the age of the Universe, and they are considered as fossil records of the time of early galaxy formation. GCs of variable age are found near all galaxies (except for the smallest dwarfs) and their specific frequency (number of clusters per galaxy mass unit, see e.g. [3]) differs as a function of galaxy type, highlighting the close relation between cluster and galaxy formation. The approximately 150 globular clusters of our own Milky Way have been studied in much more detail for their proximity. Today, star-by-star observations with the Hubble Space Telescope (HST), and proper motion studies using Gaia with high resolution spectroscopy to determine their stellar velocity dispersions [4] are possible. Small and big galaxies in the Local Group have systems of GCs, e.g. the Andromeda galaxy and the Magellanic clouds. Globular clusters in huge quantities have been detected around massive galaxies like M87 or other bright central cluster galaxies [5] or at sites of star formation near the Antenna galaxies. Still this is -in cosmological scales -our neighborhood. Do clusters form normally following the cosmic star formation history, which peaks at redshifts of around 2 [6]? Or do massive clusters form preferentially as special objects at much higher redshifts (e.g. z ∼ 6 [7])? Computer simulations of structure formation in the Universe begin to resolve GC scales [8], but they cannot compensate the current lack of deep observations. Only recently gravitational lensing from galaxy clusters has helped to identify candidates for proto-GCs at redshifts of z > 3 [9], and more recently even out to z = 6 [10,11,12]. Future instruments such as Extremely Large Telescope (ELT) and others, will improve the situation significantly, while novel instruments such as James Webb Space Telescope (JWST) are already producing ground-breaking results here [13,14]. We return to the "dense" and "gravothermal" nature of such star clusters. In gravothermal star clusters it is essential to consider the mutual gravitational interactions between many if not all of their stars. The cumulative effect of small angle gravitational deflections (encounters) between distant stars generates transport of heat and angular momentum (relaxation processes connected to these encounters are analogous to heat conduction and viscosity in a gaseous system), and this effect is dominant during certain stages of the star cluster evolution. In order to get the proper time scales connected to such relaxation processes in N -body simulations of star clusters typically a very large number of pairwise distant encounters needs to be followed (which asymptotically results in the computational complexity of all force calculations at a point of time ≈ N 2 , assuming a simple approach without parallelization and hybrid hardware). This physical constraint has led to the use of thermodynamics and statistical mechanics to model gravothermal star clusters [15,16,17,18]. The star cluster could be modelled using computer codes for the gas dynamical evolution of stars. The relaxation time in a star cluster is long compared to the dynamical timescale, while in stars the opposite is true (considering the photon diffusion time as relaxation time). This occurs because the mean free path in stellar systems is much longer than in stellar gaseous matter. Therefore, conductivity and viscosity need to be defined in a different functional form than for the interior of stars (cf. the more detailed discussions in [19,20]). An area where observations have picked up greatly through increased angular resolution and sensitivity of spectrographs is the identification of stellar binaries (e.g. [21,22,23]). Binary stars are an extremely important component of star clusters, because they form a dynamically active population which has a dramatic impact on the evolution of the host cluster (e.g. [24,25,26]): for instance stellar exotica observed in clusters (blue stragglers, fast rotating stars, and X-ray binaries) all originate from binary star evolution. The special role that binary stars play in the life cycle of a cluster requires that we pin down as accurately as possible what fraction of stars form in binaries if we are to make progress when predicting the statistics of stellar populations at the later stages of a cluster's evolution. Close stellar encounters, including direct collisions, become a reality in the densest region of dense star clusters. How often these take place, and what the outcomes of such events are remain a puzzle that is only now beginning to be solved. Last, but not least compact objects form in binaries and take part in few-body interactions and stellar evolution of binaries; ultimately, binaries consisting of compact objects only are possible sources of gravitational waves to be detected by ground-based gravitational wave dectectors that are operational today, such as (Advanced) Laser Interferometer Gravitational-Wave Observatory ((a)LIGO) [27,28,29], (Advanced) Virgo Interferometer ((a)Virgo) [30,28,29] and Kamioka Gravitational Wave Detector (KAGRA) (e.g. [28,31,32]). Note also the next (third) generation space based detectors in planning (Einstein Telescope [33] and Cosmic Explorer [34]). The numerical and computational tools to model such dense, massive star clusters are based on either approximate models from statistical physics, or the direct N -body simulation approach. Currently, the latter approach is very dominant since it allows to include details of astrophysics (binaries, stellar evolution, tidal fields) more easily. However, in order to establish the degree of reliability of N -body simulations approximate models have been very important. These are mostly based on the Fokker-Planck approximation in statistical mechanics, and the numerical solution of resulting kinetic equations by using, e.g., Monte Carlo techniques, direct numerical solution, or gaseous and moment models). Since they are important and fundamental to understand results N -body simulations this review contains a basic description of them, too. GCs and NSCs (dense and gravothermal) are ideal laboratories to examine the theoretical physical processes (heat conduction, angular momentum transport through viscosity) and their influence on the formation and evolution of extreme stellar populations like X-ray binaries or blue stragglers, compact objects (neutron stars, black holes), and are ideal test beds for stellar population synthesis models and stellar evolution.

Theoretical Foundations
Computational modelling of star clusters requires to follow the complex interplay of thermodynamic processes such as heat conduction and relaxation with the physics of self-gravitating systems, the stochastic nature of star clusters having finite particle number N , and the astrophysical knowledge and models for the evolution of single and binary stars and of external tidal forces. GCs are a very good laboratory for this, because their dynamical and relaxation timescales are well separated from each other and from the lifetime of the cluster and the Universe in its entirety. This article deals with "direct" N -body simulations, which are suitable for systems where the interaction between dynamics and relaxation is important (sometimes also denoted as "gravothermal" systems [15]). Other kinds of N -body simulations are useful for example for hydrodynamics ("smoothed particle hydrodynamics"), galaxy dynamics ("collisionless systems") or cosmological N -body simulations of structure formation in the Universe, and are not covered here. The main distinction of those from the models presented here, is that the dynamics of systems dominated by twobody relaxation ("collisional systems") require typically very high accuracy (typical energy error per crossing time ∆E/E < 10 −5 or smaller) over very long physical integration times (thousands of crossing times). The term "collisional" refers here to elastic gravitational encounters (relaxation encounters), which drive the "thermal" cluster evolution. Other processes, such as close encounters, encounters involving one or more binaries, and direct collisions also happen in the system. Let us begin with the definition of some useful time scales. A typical particle crossing time t cr in a star cluster is where r h is the radius containing 50 % of the (current) total mass and σ h is a typical velocity associated with the root mean square random motion (velocity dispersion) taken at r h . If virial equilibrium prevails, we have σ 2 h ≈ GM h /r h (where the sign ≈ here and henceforth means "approximately equal" or "equal within an order of magnitude"), thus This relation is equal to the dynamical timescale, which is also used for example in the theory of stellar structure and evolution. Global dynamical adjustments of the system, like oscillations, are connected with this timescale. Taking the square of Eq. 2 yields t 2 cr ≈ r 3 h /(GM h ) which is related to Kepler's third law, because the orbital velocity in a Keplerian point mass potential has the same order of magnitude as the velocity dispersion in virial equilibrium. Unlike most laboratory gases stellar systems are not usually in thermodynamic equilibrium, neither locally nor globally. Radii of stars are usually extremely small relative to the average inter-particle distances of stellar systems (e.g. the radius of the sun is r ≈ 10 11 cm, a typical distance between stars in our galactic neighbourhood is of the order of 10 18 cm). Only under rather special conditions in the centres of galactic nuclei and during the short high-density core collapse phase of a globular cluster, stellar densities might become large enough that stars come close enough to each other to collide, merge or disrupt each other. Therefore it is extremely unlikely under normal conditions that two stars touch each other during an encounter; encounters or collisions usually are elastic gravitational scatterings. The mean inter-particle distance is large compared to p 0 = 2Gm/σ 2 , which is the impact parameter for a 90 o deflection in a typical encounter of two stars of equal mass m, where the relative velocity at infinity is √ 2σ, with local 1D velocity dispersion σ. Thus most encounters are smallangle deflections. The relaxation time t rx is defined as the time after which the root mean square velocity increment due to such small angle gravitational deflections is of the same order as the initial velocity dispersion of the system. We use the local relaxation time as defined by [35]: G is the gravitational constant, ρ the mean stellar mass density, σ the 3D velocity dispersion, N the total particle number; this definition was used by [36,37], because it naturally occurs when computing collisional terms (as in Eq. 11), if the velocity distribution function is written as a series of Legendre polynomials [38], with numerical factors being unity (for equipartition terms of lowest order [38]) or only little different from unity, such as 9/10 for the collisional decay of anisotropy [37]). Other definitions of relaxation can be found frequently, for example in [39]. They differ only by numerical factors, except for the so-called Coulomb logarithm ln(γN ), which may take different functional forms. For common forms of the Coulomb logarithms only γ is of order unity, but may take different values (e.g. 0.11 [40], or 0.4 [39]). Assuming virial equilibrium a fundamental proportionality turns out: (cf. e.g. [39]). As a result, for very large N , dynamical equilibrium is attained much faster than thermodynamic equilibrium. Therefore, even if treated them as "gaseous" spheres, stellar systems evolve qualitatively different from stars; in stars the thermal timescale is short compared to the dynamical timescale [41]. Another interesting consequence of the long thermal timescale in star clusters is that anisotropy can prevail for many dynamical times. If one assumes a purely kinetic temperature definition, it ensues that in star clusters the temperatures (or velocity dispersions) can remain different for different coordinate directions over many dynamical times. For example, in a spherical system (using polar coordinates) the radial velocity dispersion of stars ("temperature") σ 2 r could be different from the tangential one σ 2 t . For the relaxation time above the 3D velocity dispersion σ 2 = σ 2 r + 2σ 2 t is used. If axisymmetric or triaxial the tangential velocity dispersion can be decomposed into two different dispersions 2σ 2 t = σ 2 θ + σ 2 φ . A full account on the relevance of anisotropy for star clusters is beyond the scope of this paper; exemplary we mention here that interest in anisotropy was recently sparked by anisotropic mass segregation in rotating star clusters, both globular and nuclear [42,43,44,45,46].

Direct Fokker-Planck and moment models
Models based on the Fokker-Planck approximation (also denoted as approximate or statistical models) have been designed and implemented in times when it was very difficult to simulate large star clusters directly by N -body simulations. Dramatic development in hardware and software has made now possible direct N -body simulations of up to a million bodies with realistic astrophysics and binaries [47,48] (Dragon simulations). However, the use of the approximate models is very useful to understand the nature of physical processes in star clusters (such as heat conduction or viscosity); by comparison with N -body models they can mutually support each other. The Fokker-Planck approximation to describe two-body relaxation in spherical star clusters is the foundation for all Monte Carlo models used nowadays (MOCCA and CMC, see Sect. 4 below). Therefore it is deemed useful to provide a deeper than usual insight into its theoretical foundations here. Statistical models have been employed to clarify the physical nature of relaxation processes in star clusters, such as core collapse [15], post-collapse evolution due to an energy source from binaries undergoing close encounters with single stars [49], gravothermal oscillations [50,41,51]. These methods have also been important for the study of anisotropy, mass segregation, and rotation later on. Comparison and mutual adjustment of parameters in order to get agreement between statistical models and direct N -body simulations was started by [52] for a first pre-collapse comparison, followed by an extended study using statistical averages of N -body simulations to match gaseous models [40,53,54], see also Fig. 1; in [55] it was shown that relaxation processes consistent with theory dominate core collapse in star clusters, and [56] showed for the first time a signature of gravothermal oscillations in a N -body simulation. Classical models based on the Fokker-Planck approximation use quite strong approximations, like spherical symmetry (in general, with some exceptions allowing axisymmetry), dominance of relaxation encounters, modelling all few-body effects (binary-single and binary-binary close encounters) in a statistical way. A mass spectrum would be modelled by discrete dynamical components with different masses (except Monte Carlo models, see below). With the increasing need to have star cluster models matching detailed observations of star clusters, the use of Fokker-Planck type models was no more practical. That hardware and software developments made more and more realistic particle numbers possible in direct N -body simulations has been another reason for the decline of  [20], isotropic gas model [57], an isotropic Fokker-Planck model, and direct N -body simulation (from [54], the figure is also contained in Chapter 8 of [58]. statistical models. There is one remarkable exception, the Monte Carlo technique -while it is also based on Fokker-Planck theory it uses a quasi N -body realization and allows state-of-the-art models up to the current time (see Sect. 4). In the following subsections, nevertheless, we will outline the basic theory beneath the statistical models, because in some areas (rotating, flattened star clusters, NSCs, anisotropy) they are still important today in order to analyse results of N -body simulations.

Fokker-Planck approximation
The Fokker-Planck approximation truncates the so-called BBGKY hierarchy (named after Bogoliubov-Born-Green-Kirkwood-Yvon) of kinetic equations at lowest order assuming that for most of the time all particles are uncorrelated with each other and only coupled via the smooth global gravitational potential (see Chapter 8.1, Sects. 2 and 3 of [58]; the following paragraph is a summary from their text). We start with most general N -particle distribution function f (N ) , which depends on positions and velocities r i , v i of a set of N particles and time t: It provides a probability to find all the particles at their given positions and velocities. In 6N dimensional phase space the particles are an incompressible fluid following Liouville's "big" equation where the derivative is the Lagrangian derivative. If all particles are uncorrelated f (N ) = (f (1) ) N , i.e. the N -particle distribution function is just the N -fold product of a single particle distribution function. That is the case for example in a collisionless stellar system, where all particles just follow their trajectories determined by a global smooth gravitational potential and any direct interaction between two or few particles (stars) is negligible. For collisional stellar systems, however, gravitational encounters (two-body relaxation) change the phase space distribution, particles are not uncorrelated anymore. The theoretical ansatz in that case would be to define a two-body correlation function g by So, from knowing f (1) and g we get f (2) ; by using higher order correlation functions one can get from f (n) to f (n+1) . Integration of Eq. 6 step-by-step over single particles provides a sequence of equations for f (N −1) to f (1) , which is the BBGKY hierarchy. However it is usually not very helpful, because all the correlation functions for 2, ...N−1 need to be known. For practical purposes in collisional stellar systems, where two-body relaxation is important (e.g. open, globular, and nuclear star clusters), it is sufficient to deal with the two-body correlation, which is done phenomenologically in two different ways for distant and close correlations (encounters and binaries), as described below. Higher than two-body correlations are rarely important. There could be a relation to Sundmann's famous theorems [59,60], which state that in the threebody problem direct three-body collisions occur only with a negligibly small probability; it means whenever three particles get close to each other, there will be always a sequence of separate close two-body encounters, practically never 1 the three bodies will simultaneously get extremely close together. Burrau's three body problem is a nice demonstration of that behaviour [61]; whether the fundamental assumption of dominance of two-body correlations is in fact realized or not can only be checked computationally by comparison of models based on the Fokker-Planck approximation (such as also Monte Carlo models) with direct N -body simulations. An example for a situation of a very high density in a collapsing core of a star cluster, where higher correlations become important, can be found in [62]. Instead of determining a general correlation function one resorts to a phenomenological description of the effects of collisions by computing local diffusion coefficients directly from the known solution of the two-body problems. Diffusion coefficients D(∆v i ) and D(∆v i v j ) denote the average rate of change of v i and v i v j due to the cumulative effect of many small angle deflections during two-body encounters, at a given radius r (from here assuming spherical symmetry). Let m, v and m f , v f be the mass and velocity of a star from a test and field star distribution, respectively (both distributions can but need not to be the same). In Cartesian geometry (for the local velocities) the diffusion coefficients are defined by Local means here that we do not explicitly consider the dependence of f on the spatial coordinate. g, h are the Rosenbluth potentials defined in [63] h Note that provided the distribution function f is given in terms of a convenient polynomial series as in Legendre polynomials the Rosenbluth potentials can be evaluated analytically to arbitrary order, as was seen already by [63], see for a modern re-derivation and its use for star cluster dynamics [54,38,64].
With these results we can finally write down the local Fokker-Planck equation in its standard form for the Cartesian coordinate system of the v i : The subscript "enc" should refer to encounters, which are the driving force of two-body relaxation. Still Eq. 10 is a six-dimensional integro-differential equation; its direct numerical simulation in stellar dynamics can presently only be done by further simplification. If the encounter term is zero, Eq. 10 is transformed into Liouville's equation for a collisionless system. For a selfgravitating system Eqs. 10 and 11 are not sufficient, since the knowledge of the gravitational potential Φ is necessary. This can be seen above from the v i term -its computation requires to know the gravitational force. In moment or gas models described below (for spherical symmetry) Poisson's equation takes the simple form Eq. 28; for orbit averaged Fokker-Planck or Monte Carlo models (see Sect. 3.3 and 4.2) the gravitational potential enters directly into the energy as constant of motion (cf. Eq. 30).

Moment or Gas Models
The local Fokker-Planck equation (Eq. 10) is utilized in another way for gaseous or conducting gas sphere models of star clusters. Integrating it over velocity space with varying powers of the velocity coordinates yields a system of equations in the spatial coordinates; the local approximation is used in the sense that the orbit structure of the system is not taken into account, diffusion coefficients and all other quantities are assumed to be well defined just as a function of the local quantities (density, velocity dispersions and so on). The system of moment equations is truncated in third order by a phenomenological equation of heat transfer. Such an approach has been suggested by [19,57] and generalized to anisotropic systems by [65,37,20]. In the following the derivation of the model equations is described.

The "Left Hand Sides"
In spherical symmetry, polar coordinates r θ, φ are used and t denotes the time. The vector v = (v i ), i = r, θ, φ, denotes the velocity in a local Cartesian coordinate system at the spatial point r, θ, φ. In the interest of brevity The distribution function f , which due to spherical symmetry is a function of r, t, u, v 2 + w 2 only, is normalized according to where ρ(r, t) is the mass density; if m denotes the stellar mass, we get the particle density n = ρ/m. Then is the bulk radial velocity of the stars. Note that for the analogously defined quantitiesv andw we have in spherical systemsv =w = 0 (rotating, axisymmetric systems:w = 0). In order to go ahead to the anisotropic gaseous model equations we now turn back to the left hand side of the Fokker-Planck equation Eq. 10, which is the collisionless Boltzmann or Vlasov operator. For practical reasons we prefer for the left hand side local Cartesian velocity coordinates, whose axes are oriented towards the r, θ, φ coordinate space directions. With the Lagrange function the Euler-Lagrange equations of motion for a star moving in the cluster potential Φ become: The complete local Fokker-Planck equation, derived from Eq. 10, attains the form where the term subscribed by "enc" denotes the terms involving diffusion coefficients as in Eq. 11. Moments i, j, k of f are defined in the following way (all integrations range from −∞ to ∞): 2, 0, 0 : 1, 2, 0 := F θ +ūp θ = uv 2 f dudvdw 1, 0, 2 := F φ +ūp φ = uw 2 f dudvdw .
Note that the definitions of p i and F i are such that they are proportional to the random motion of the stars. Due to spherical symmetry we have p θ = p φ =: p t and F θ = F φ =: F t /2. By p r = ρσ 2 r and p t = ρσ 2 t the random velocity dispersions are given, which are closely related to observables in GCs and galaxies. It is convenient to define velocities of energy transport by By multiplication of the Fokker-Planck equation 16 with various powers of u, v, w we get up to second order the following set of moment equations forū dropped in the following): The terms labeled with "enc" and "bin3" symbolically denote the collisional terms resulting from the moments of the right hand side of the Fokker-Planck equation (Eq. 11) and an energy generation by formation and hardening of three body encounters. Both will be discussed below. With the definition of the mass M r contained in a sphere of radius r ∂M r ∂r = 4πr 2 ρ the set of Eqs. 25-27 is equivalent to gas-dynamical equations coupled with Poisson's equation. Since moment equations of order n contain moments of order n + 1, it is necessary to close the system of the above equations by an independent closure relation. Here we choose the heat conduction closure, which consists of a phenomenological Ansatz in analogy to gas dynamics. It was first used (restricted to isotropy) by [19]. It is assumed that heat transport is proportional to the temperature gradient, where we use for the temperature gradient an average velocity dispersion σ 2 = (σ 2 r + 2σ 2 t )/3 and assume v r = v t (this latter closure was first introduced by [37]). Therefore, the last two equations to close our model are With Eqs. [25][26][27]28, and 29 we have now seven equations for our seven dependent variables M r , ρ, u, p r , p t , v r , v t .

The "Right Hand Sides"
All right hand sides of the moment equations 25-27 are calculated by multiplying the right hand side (the encounter term) of the Fokker-Planck equation as it occurs in Eq. 11 with the appropriate powers of u, v and w and integrating over velocity space. Since the diffusion coefficients in Eq. 8 also contain the distribution function f , the equation depends non-linearly on it. That has led in the early papers to a simplification by using an isotropized background distribution function f b inside the diffusion coefficients, different from the actual one [36,66]. In [37,67,64] there is always full consistency between the background and the actual distribution function.

Orbit averaged Fokker-Planck models and rotation
The direct solution of the six-dimensional integro-differential equations (Eq. 10 and Eq. 11) is generally not possible. To have numerical solutions of the Fokker-Planck equation directly one applies Jeans's theorem and transforms f into a function of the classical integrals of motion of a particle in a potential under the given symmetry, as e.g. energy E and modulus of the angular momentum J 2 in a spherical potential or E and z-component of angular momentum J z in axisymmetric coordinates. Thereafter, the Fokker-Planck equation can be integrated over the accessible coordinate space for any given combination of constants of motion and the orbit-averaged Fokker-Planck equation ensues. By transformation from v i to E and J and via the limits of the orbital integral ω 0 Trot /T kin (%) e dyn (0) r tid /rc(0) r h /rc(0) τrc(0) τ rh (0) 0.00 0.  Table 1: Initial conditions of rotating King models from [67] with W 0 = 6. T rot /T kin is the ratio of bulk rotational energy to total kinetic energy in percent, e dyn (0) is the dynamical ellipticity, r tid /r c (0) is the ratio of the tidal radius to core radius, r h /r c (0) is the ratio of the half-mass radius to core radius, τ rc (0) is the central relaxation timescale and τ rh (0) is the half-mass relaxation timescale. All of these quantities are shown for t = 0 of system time units. Table data are from [67], but note that in this paper the per cent sign (%) in the header for T rot /T kin had been omitted erroneously.
the potential enters both implicitly and explicitly. In a two-step scheme alternatively solving the Poisson-and Fokker-Planck equation a direct numerical solution is obtained for spherical systems in 1D (using E only [66]), or in 2D (using both E and J 2 [68,69,70,71]). One of the main uncertainties in this method is that for non-spherical mass distributions the orbit structure in the system may depend on unknown non-classical third integrals of motion which are neglected. First 2D models of axisymmetric, rotating globular star clusters [67] used initial models obtained published earlier [72,73,74], which are generalizations of the standard King models [75], using its dimensionless central potential value W 0 and a new dimensionless rotation parameter ω 0 . These models have a rigid rotation in the core, maximum of the rotation curve around the half-mass radius and a differentially decreasing rotation in the halo. They are still mainly supported by "thermal" pressure (velocity dispersions), the rotational energy provides a smaller part of the energy. Typical initial data for W 0 = 6 and different rotation parameters are seen in Table 1: the ratio of total rotational to total kinetic energy, dynamical ellipticity, ratios of tidal and half-mass radii to initial core radii, the central relaxation time and finally the half-mass relaxation time in system units are given (for definitions see [67]).
Fokker-Planck models showed that in presence of rotation there is an effective viscosity transporting angular momentum outwards and accelerating cluster evolution significantly as compared to a spherical cluster (see Fig. 2 and [67]). A series of follow-up papers include post-collapse and multi-mass models [76,77,78] and found an accelerated rotation in the core for heavy masses sinking to the core -as it was predicted by the combined gravogyro  Table 1) ω 0 = 0.0, 0.3, 0.6 as indicated. Shown are the radii for mass columns containing the indicated percentage of total initial mass in the direction of the θ = 54 • .74 angle. From [67]. and gravothermal "catastrophes" predicted by [17,79]. One rotating model included in the Dragon simulations [48], however, did not show accelerated evolution. Whether this is due to heavy mass loss by stellar evolution (not included in earlier papers) or due to a small deviation from the proper initial model is not clear. There is an urgent need for more coverage of rotating stellar clusters by direct N -body simulations, see some first progress [80,81,82]. The initial models of Table 1 are still in frequent use, in particular if realized as N -body configurations for N -body models [83,84,80,81,82]. Notice also the alternative rotating models of Varri [85,86], which are more suitable with regard to the outer cluster zones under influence of tidal fields.

Monte Carlo Models
Monte Carlo models of star clusters are the only ones which are still intensively used up to the present time, even though they are based on the Fokker-Planck approximation, in the same way as Fokker-Planck or gaseous/moment models. Sometimes this may not be clear to every reader of current papers using Monte Carlo models, because they provide data equivalent to N -body simulationsparticles with masses, positions and velocities at certain times. Astrophysics (stellar single and binary evolution, stellar collisions, relativistic binaries...) has been included very much like in N -body models. This review is not about Monte Carlo models, but a brief summary of their history and entry points to the current literature should be given.

Hénon and Spitzer type method
As the name suggests Monte Carlo models are based on the principle that stars have an orbit in a known self-consistent potential; random perturbations are applied, which model the effect of relaxation by distant gravitational encounters. Spitzer's method follows the orbits of stars in the global potential of the cluster and randomly applies kicks in velocity to the stars; at the end of a long series of papers they included binaries and a mass spectrum [87,88,89,90,91,92,93,94]. Hénon's method is using the phase space of constants of motion of a star in a spherically symmetric potential, energy and angular momentum. Deflections are selected randomly, and their effect on angular momentum and energy computed and applied [95]. The method was extended to include astrophysical effects, including binaries and stellar evolution [96,97]. These models still allowed for "superstars", i.e. one particle in the Monte Carlo model could represent many real stars. Current Monte Carlo models are based on Hénon's method, but restricted to star-by-star modelling (much like N -body), where every star is a particle in the Monte Carlo simulation. This only made it possible to include all astrophysical effects in the same way than it is done in N -body simulations. This new line of Monte Carlo models was initiated by Giersz [98] (code name MOCCA) and the Northwestern team [99] (code name CMC).

MOCCA and CMC
The Monte Carlo codes based on the Hénon scheme use constants of motion (specific energy E, specific angular momentum L) as basic variables, properties of stars in the simulation. If the spherically symmetric gravitational potential Φ(r) is known, the pericenter r min and apocenter r max of the orbit are known. At every point of the orbit r the radial velocity is known from The orbital integral defines the orbital time τ by With p(r) = (2/τ )·(dr/v r ) one gets a probability distribution function, used to randomly pick a radial position r i for the star on its orbit (which should be distributed according to p(r)). Let m i be the stellar mass of stars (i = 1, . . . , n), then the spherically symmetric gravitational potential can be computed according to [95] In addition to that two angles θ and φ are randomly picked, so as to have a three dimensional position of the star. Velocities are obtained from E, L, and U (r i ) (one more random number needed). In that way a model star cluster is produced whose data structure is three dimensional -equivalent to that of an N -body simulation. To model the relaxation effect, two neighbouring stars are selected and a mean squared deflection angle chosen, which is proportional to the time-step over the relaxation time. Using this angle changes in E and L are computed. Binaries and close encounters between them have been first modelled completely stochastically as well (using random impact parameters, and using random realization of known cross sections). More recently a fewbody integration is done in both MOCCA and CMC codes. This is a very rough account of Monte Carlo principles, the reader interested in more details is referred to the papers cited in the next paragraphs. An account of the state of the MOCCA code and comparison with N -body simulations is published here [100,101]. Recently, it has been used for a large number of simulations of Galactic and extragalactic clusters, the MOCCA Survey Database has been published [102,103,104,105]. The CMC code [106] has been developed in parallel, with matching models to observations [107], and an overview of the current state of the code [108]. Examples of current use of this code focus on compact remnants and their gravitational wave emission, such as e.g. [109,110,111]. Both Monte Carlo codes have been very successful in terms of generating a large amount of cluster simulations to be compared with observational data and also to follow the evolution of special objects. However, one should not forget their serious limitations: if we have a number of massive objects in a central high density regionthe assumption of a smooth spherically symmetric potential breaks down; at high densities and if many binaries are present, the assumption that there are uncorrelated two-body relaxation encounters and close few-body encounters, which can be clearly separated, breaks down. -Taking into account external tidal fields is quite difficult, though in simple cases not impossible, due to the strictly spherical cluster centered gravitational potential.
The bottom line is that Monte Carlo models have to be used in order to get an overview of large parameter ranges of star cluster evolution, but in many cases a check by comparison with direct N -body simulations is desirable. They do not suffer from all the problems mentioned above; however, also direct numerical solutions of the N -body problem have certain issues, see Sect. 5.4.
A nice overview of current Monte Carlo models is in [112], who also present a somewhat restricted Monte Carlo code for rotating systems (see Sect. 6.4.6).

Direct N -body Simulations -Methods and Algorithms
To integrate the orbits of particles in time under their mutual gravitational interaction the total gravitational potential at each particle's position is required. Poisson's equation in integral form gives the potential Φ generated at a point in coordinate space r due to a smooth mass distribution ρ( r) A discrete particle distribution in N -body simulations is given by with N particles of mass m i distributed at positions r j . Putting this into the integral Poisson equation (33) we get Newton's law for point masses: 5.1 Nbody -the growth of an industry It was already discovered by Sebastian von Hoerner in the earliest published N -body simulations that the relaxation time [35] is relevant for star cluster evolution and that the formation of close and eccentric binaries occurs as the rule rather than as an exception. It was particularly difficult to accurately integrate them, effectively the simulation would be stopped if close binaries demanded too small time-steps [152,153]. About at the same time a young postdoc -Sverre Aarseth -in Cambridge developed a direct N -body integrator for galaxy clusters with gravitational softening, thereby avoiding von Hoerner's problems with tight binaries [113]. His code was based on Taylor series evaluation of the gravitational force up to its second derivative. Eight years later regularization methods [115] were implemented in Aarseth's direct N -body code [154]. This allowed to proceed past the binary deadlock detected in von Hoerner's models. Another direct N -body code by Roland Wielen appeared on the market, and   [134], corrected in some places, but expanded to more recent developments. The abbreviations used are as follows: Data) Extension for vectorization in the CPU (central processing unit). -P 3 T / SDAR: particle-particle particle-tree/Slow-Down Agolrithmic chain -MSTAR / BiFrost: Minimum spanning tree + algorithmic chain / Binaries in Frost in a seminal paper [52] fair agreement was shown between Aarseth's and Wielen's codes and a Monte Carlo code by Lyman Spitzer (see above Sect. 4.1). However, only at the turn of the century Aarseth and von Hoerner could compare their codes, and von Hoerner published a remarkable account on "how it all started" [155]. Already in 1985 the code Nbody5 [120] had become a kind of "industry standard", attaining world wide use. It employed Taylor series using up to the third derivative of the gravitational force, in a divided difference scheme based on four time points, with individual particle time-steps. Also there were regularizations for more than two bodies, such as the classical chain regularization [123], and the Ahmad-Cohen [118] neighbour scheme already in Nbody5. The advent of vector and parallel computers demanded an optimization towards hierarchically blocked time-steps and the Hermite scheme (Sect. 5.2.1) [122,128], because it used only two time points, which made memory management easier. This became known as Nbody6. The growth of the "industry" [134] included further improvements in the regularization techniques [156,132,157] and a comprehensive book summary [158]. Table 2 summarizes the main algorithmic, hardware and software development stepping stones in the direct N -body community up until today.

The Nbody6 scheme
In the following, the Nbody6 integrator is described in some more detail (note that Nbody7 already contains parallelization through GPU acceleration and will be treated in the next section). Nbody6 and its parallelized and GPU accelerated offspring (Nbody6++, Nbody6GPU, Nbody6++GPU, Nbody7, see Table 3) is still the most widely used method for direct Nbody simulations, but recently also new approaches have been published (cf. Sect. 5.5).

The Hermite scheme
The Hermite scheme and Nbody6 go back to [128]; in conjunction with a hierarchically blocked time-step scheme (see below and [122]) it improved the performance on vector computers and turned out to be efficient for all of recent parallel and innovative hardware (general and special purpose parallel computers, GRAPE and GPU). Assume a set of N particles with positions r i (t 0 ) and velocities v i (t 0 ) (i = 1, . . . , N ) is given at time t = t 0 , and let us look at a selected test particle at r = r 0 = r(t 0 ) and v = v 0 = v(t 0 ). Note that here and in the following the index i for the test particle i and also occasionally the index 0 indicating the time t 0 will be dropped for brevity; sums over j are to be understood to include all j with j = i, since there should be no self-interaction. Accelerations a 0 and their time derivativesȧ 0 are calculated explicitly: where By low order predictions, new positions and velocities for all particles at t > t 0 are calculated and used to determine a new acceleration and its derivative directly according to Eq. 36 at t = t 1 , denoted by a 1 and ȧ 1 . On the other hand a 1 and ȧ 1 can also be obtained from a Taylor series using higher derivatives of a at t = t 0 : If a 1 and ȧ 1 are known from direct summation (from Eq. 36 using the predicted positions and velocities) one can invert the equations above to determine the unknown higher order derivatives of the acceleration at t = t 0 for the test particle: This is the Hermite interpolation, which finally allows to correct positions and velocities at t 1 to high order from Taking the time derivative of Eq. 44 it turns out that the error in the force calculation for this method is O(∆t 4 ), as opposed to standard leapfrog scheme, which has a force error of O(∆t 2 ) (but see new developments in Sect. 5.5). Additional errors induced by approximate potential calculations (particle mesh or Tree) create potentially even larger errors than that. In Fig.  3, however, it is shown that the above Hermite method used for a real N -body integration sustains generally an error of O(∆t 4 ) for the entire calculation.

Time-step Choice
Aarseth [120] provides an empirical time-step criterion The error is governed by the choice of η, which in most practical applications is taken to be η = 0.01 − 0.04. It is instructive to compare this with the inverse square of the curvature κ of the curve a(t) in coordinate space Clearly under certain conditions the time-step choice of Eq. 45 becomes similar to choosing the time-step according to the curvature of the acceleration curve; since it was determined just empirically, however, it cannot generally be related to the curvature expression above. In [127] a different time step criterion has been suggested, which appears simpler and more straightforwardly defined, and couples the time-step to the difference between predicted and corrected coordinates. The standard Aarseth time-step criterion from Eq. 45 has been used in most N -body simulations so far, because it seems to achieve an optimal step better than (on average) the mathematically more sound Makino step (see the time-step related discussion in [159]).
Since the position of all field particles can be determined at any time by the low-order prediction from Eq. 38, the time-step of each particle (which determines the time at which the corrector of Eq. 44 is applied) can be freely chosen according to the local requirements of the test particle; the additional error induced due to the use of only predicted data for the full N sums of Eq. 36 is negligibly small, for the benefit of not being forced to keep all particles in lockstep. Such an individual time-step scheme is in particular for non-homogeneous systems very advantageous, as was quantitatively pointed out by [160]. Particles in the high density core of a star clusters need to be updated much more often than particles on orbits very far from the centre. They show that the gain in computational speed due to the individual timestep scheme (as compared to a lockstep scheme where all particles share the minimum required time-step) is of the order N 1/3 for homogeneous and N 1 for strongly spatially structured systems [160].
For the purpose of vectorization and parallelization it is better not to have the particles continuously distributed on a time axis. Consequently, [161] uses a hierarchical scheme, still on the basis of Eq. 45; but a change of the time-step is considered only if that equation yields a variation of ∆t compared to the last step by more than a factor of 2 (increase or decrease). If this is the case a variation by 2 is applied only. Thus in model units all time-steps are selected from the set {2 −i |i = 0, ...i max } with k = i max determined by the condition that ∆t min > 2 −imax for the minimum time-step ∆t min determined from Eq. 45. For core collapse simulations of star clusters of a few ten thousand particles i max goes up to about 20; empirically and theoretically [160] ∆t min ∝ N −1/3 , so for large N i max becomes larger, however, on the other hand, how large i max grows for fixed N depends on the selected criteria for so-called KS regularisation of perturbed two-body motion (see below). The implementation of the block step scheme indeed uses an even stronger condition than the above described one, it is demanded that not only the time-steps, but also the individual accumulated times of each particles are commensurate with the time-step itself. This ensures that for any particle i and any time T i = t i + δt i all particles with δt j < δt i have for their own time T j = t j + δt j = T i , where the last equality is the non-trivial one. Such procedure is important for the parallelization of the algorithm. For example it has as a consequence that at the big timesteps always huge groups of particles are due for correction, sometimes even all particles (at the largest steps).

Ahmad-Cohen neighbour scheme
Another refinement of the Hermite or Aarseth "brute force" method is the twotime-step scheme, denoted as neighbour or Ahmad-Cohen scheme [118]. For each particle a neighbour radius is defined, and a and ȧ are computed due to neighbours and non-neighbours separately. Similar to the Hermite scheme the higher derivatives are computed separately for the neighbour force (irregular force) and non-neighbour force (regular force). Computing two time-steps, an irregular small ∆t irr and a regular large ∆t reg , from these two force components by Eq. 45 yields a time-step ratio of γ := ∆t reg /∆t irr being in a typical range of 5-20 for N of the order 10 3 to 10 4 . The reason is that the regular force has much less fluctuations than the irregular force. The Ahmad-Cohen neighbour scheme is implemented in a self-regulated way, where at each regular timestep a new neighbour list is determined using a given neighbour radius r si for each particle. If the neighbour number found is larger than the prescribed optimal neighbour number, the neighbour radius is increased or vice versa. In [120,160] more complicated algorithms to adjust the neighbour radius are described. On the contrary to [160], who find an optimal neighbour number of N n,opt ∝ N 3/4 we find that adopting a constant neighbour number of the order of 50 − 200 is sufficient at least up to N = 10 6 . The reason is that by using special purpose machines or parallelization for parts of the code, an optimal neighbour number is not well defined, so the neighbour number can be selected according to accuracy and efficiency requirements. After each regular time-step the new neighbour list is communicated along with the new particle positions to all processors of the parallel machine, thus making it possible to do the irregular time-step in parallel as well.
Using a two-time-step or neighbour scheme again increases the computational speed of the entire integration by a factor of at least proportional to N 1/4 [127]. Both the regular and irregular time-steps are arranged in the hierarchical, commensurable way.

Regularizations
As the relative distance r of the two bodies becomes small, their time-steps are reduced to prohibitively small values, and truncation errors grow due to the singularity in the gravitational potential. Such a close encounter is characterised by an impact parameter p and a relative velocity at "infinity" (in practice some distance inside the cluster) v ∞ . A close encounter is characterized by where p 90 is the impact parameter related to a 90 degree deflection in a twobody problem; G, m 1 , m 2 , v ∞ are the gravitational constant, the masses of the two particles and their relative velocity at infinity. In the cluster centre, it is very likely that two stars come very close together in a hyperbolic encounter. So, if the separation of two particles gets smaller than p 90 they are candidates for regularization. To be actually regularized, the two particles have to fulfil two more sufficient criteria: that they are approaching each other, and that their mutual force is dominant. These sufficient criteria are defined as Here, a pert is the vectorial differential force exerted by other perturbing particles onto the two candidates, R, R, V are scalar and vectorial distance and relative velocity vector between the two candidates, respectively. The factor 0.1 in the upper equation allows nearly circular orbits to be regularized; γ < 0.25 demands that the relative strength of the perturbing forces to the pairwise force is one quarter of the maximum. These conditions describe quantitatively that a two-body subsystem is dynamically separated from the rest of the system, but not necessarily unperturbed.
The idea is to take both stars out of the main integration cycle, replace them by their centre of mass (c.m.) and advance the usual integration with this composite particle instead of resolving the two components. The internal motion of the two members of the regularized pair (henceforth KS pair, for Kustaanheimo and Stiefel [115]) is done in a separate coordinate system. However, as was already noted by [154] there is no need for the perturbation of the KS pair from other stars to be small. The internal motion of a KS pair is integrated in a 4D vector space obtained from a combined canonical and time transformation of relative Cartesian positions and velocities. The coordinate transformation goes back to Levi-Civita in 2D [162]. A full generalization to higher dimensions is only possible over the mathematical object of a field, the next one to be quaternions in 4D. Kustaanheimo and Stiefel found a way to transform forward from 3D to 4D and back from 4D to 3D by working over a skewed field of quaternions (sacrificing some commutativity rules; their mathematical language was different though). A modern theoretical approach to this subject can be found e.g. in [163]; the complete formalism including also the time transformation can be found in [164]. Aarseth uses this method to integrate the KS pairs in 4D space, and when using the back-transformation automatically returning to Cartesian 3D space [154]. The KS transformation converts the motion in a singular Newtonian gravitational potential into a harmonic oscillator in 4D space, which has no singularity. Since the harmonic potential is regular, numerical integration with high accuracy can proceed with much better efficiency, and there is no danger of truncation errors for arbitrarily small separations. The internal time-step of such a KS-regularized pair is independent of the eccentricity and of the order of some 50-100 steps per orbit. While regularization can be used for any analytical two-body solution even across a mathematical singularity (collision), it is practically applied to perturbed pairs only. Once the perturbation γ falls below a critical value of ≈ 10 −6 , a KS pair is considered unperturbed, and the analytical solution for the Keplerian orbit is used instead of doing numerical integration. The two-body KS regularization occurs in the code either for short-lived hyperbolic encounters or for persistent binaries. Close encounters between single particles and binary stars are also a central feature of cluster dynamics. The chain regularization [132] is invoked if a KS pair has a close encounter with another single star or another pair. Especially, if systems start with a large number of primordial (initial) binaries, such encounters may lead to stable (or quasi-stable) hierarchical triples, quadruples, and higher multiples. They are treated by using special stability criteria [165]. Every subsystem -KS pair, chain or hierarchical subsystem -could be perturbed by other single stars. Perturbers are typically those objects that get closer to the object than R sep = R/γ 1/3 min , where R is the typical size of the subsystem; for perturbers, the components of the subsystem are resolved in their own force computation as well. Algorithmic regularization is a relatively recent method based on a time transformed leap-frog method [143]; it does not employ the KS transformation. See for its use and application the next subsection.

Parallel and GPU computing and Nbody
A fundamental problem was raised by Daiichiro Sugimoto about 30 years ago [166] -direct numerical simulations of globular star clusters -with order of a million stars in direct N -body -could not be completed for decades if extrapolating the standard evolution of computational hardware at that time (Moore's law) for the future. Therefore astronomers in the Department of Astronomy at Tokyo University started wire-wrapping and designing a new integrated circuit, a special purpose computer chip named GRAPE (=GRAvity PipE). The work was continued with great success by the team of Junichiro Makino, the GRAPE chips were finally assembled into GRAPE accelerator mainboards containing several chips (such as HARP, GRAPE-4, GRAPE-6) [167,168,169,140]. The GRAPE chip is an application specific integrated circuit (ASIC), which could only compute gravitational forces between particles (it also computed the time derivative of the force, to be directly applicable to the Hermite scheme of Nbody6). A GRAPE board is a multi-core (multi-chip) parallel computing device (e.g. GRAPE-4 board contained 48 chips with shared memory, each chip contained one pipeline for force calculation; the GRAPE-6 chip contained 6 force pipelines [140]). Custom built computing clusters using GRAPE were built outside of Japan e.g. in Rochester and Heidelberg [170]. In the following years, graphical processing units (GPU) widely replaced GRAPE; direct N -body implementations were done on GPU clusters [142,171]. Interfaces have been written such that GRAPE users could right away also use GPU with the newly invented programming language CUDA (Yebisu [172], Sapporo [173], Kirin [174,175]). Still somewhat state of the art is Nbody6GPU, which includes GPU acceleration of Nbody6 using CUDA kernels for single node servers [145]. Many of these kernels written by Keigo Nitadori are still in current use, even for the massively parallel programs such as Nbody6++GPU, see below. At the same time another development started, parallelization of Nbody6 with the (at that time) new standard MPI (message passing interface). Nbody6++ [137] uses the SPMD (Single Program Multiple Data) scheme to run many instances of the code in parallel, while distributing force computations for different particles to the processors of a massively parallel computer. From time to time data transfers using MPI communication routines are necessary, to make sure all processors are synchronous. Systems with hundreds of processing units were used at the time (e.g. CRAY T3E), which demanded efficient coding of the communication scheme. Copy and ring algorithms were developed, and asynchronous data transfer and computation implemented [176,177]. A copy algorithm keeps always a complete copy of all particle data on every parallel process; parallelization is over groups of particles due for the correction step; communication sends around all new particle positions and velocities in the Hermite scheme to all other processes. In contrast the ring algorithm uses a domain decomposition, every process has its specific set of particles (at least for some time), and instead of communicating particle positions and velocities partial gravitational forces and their time derivatives are communicated. A copy algorithm has been implemented by [137,178] for Nbody6++, and a ring algorithm is used in phiGRAPE [170]. All these communication algorithms have been implemented long time ago using the MPI SENDRECV routine in a cyclic fashion -for p processes p − 1 communication steps are needed 2 . Every process simultaneously sends data to its next neighbour and receives data from its other neighbour, in a ring structure. Therefore these algorithms are also denoted as systolic communication algorithms (both copy and ring). Nowadays MPI ALLGATHER or MPI ALLREDUCE may be used, but their implementations are not transparent and vary; the latter would normally use a Tree-based implementations (instead of systolic) -the number of communication steps is then only log 2 (p) (while our systolic algorithm needs O(p) steps). It can be shown that asymptotically (for large data chunks and low latency) both algorithms are equivalent with regard to the total time required, because the Tree based algorithm uses increasingly large data packages, while in systolic algorithms every step communicates the same amount of data. [177]. Hence, currently still the systolic communication with a copy algorithm is used in Nbody6++GPU. If going to ten or hundred million bodies the copy algorithm may become too large for system memories, and should be updated to the ring algorithm with domain decomposition, which is not a fundamental problem (already used in the phiGRAPE code, which is a simple variant of the Hermite scheme with blocked hierarchical time steps). While both ring and copy algorithms scale linearly with p a hypersystolic algorithm exists which scales only with √ n p [179,180]. For GRAPE a spe-cial implementation of a hypersystolic algorithm for 2D meshes of processing elements has been proposed [176]. Hypersystolic and other Tree based communication algorithms can play out their strengths in case of a huge number of processes (as in case of GRAPE for example) with relatively modest computational load on every process. On the contrary the current Nbody6++GPU algorithm requires modest numbers of processes (order 10-100 with current particle numbers of up to around a million, may be more in the future), which have big computing loads and very large data chunks to communicate. Nbody6++ [137] parallelizes both force loops with MPI, for the regular and the neighbour force in the Ahmad-Cohen scheme.
A ring communication algorithm with domain decomposition in the future would also help for situations when there are many small block time-steps with few particles to integrate. The current code Nbody6++ (and its successors Nbody6++GPU using GPU) only invoke parallel MPI execution if the number of particles in a time block is large enough (like e.g. 50-100, best value has to be tested for every hardware). For smaller blocks all processors are redundantly computing everything without communication, to avoid the overhead connected with MPI. Since the special hierarchical time-step scheme of Nbody6 favours time blocks with many particles this is for usual globular cluster simulations no bottleneck. However, in case of very high central density, like in nuclear star clusters with central supermassive black hole (see Sect. 7) the parallel performance gets degraded.
With the advent of clusters, where nodes would be running MPI, and each node having a GPU accelerator, Nbody6++GPU was created -on the top level MPI parallelization is done for the force loops (coarse grained parallelization) and at the bottom level each MPI process calls its own GPU to accelerate the force calculation [146], using Nitadori's Yebisu library [172] for the regular force only. Secondly, an AVX/SSE implementation accelerates prediction and neighbour (irregular) forces, and also a number of other features had been severely optimized and improved (such as particle selection at block times) [47]. This code currently keeps the record of the largest direct N -body simulation of a globular cluster with all required astrophysics (single and binary stellar evolution, stellar collisions, tidal field), simulated over 12 Gyrs [48].
In recent years, inspired also by LIGO/Virgo/KAGRA gravitational wave detections [181], numerous current updates have been made with regard to stellar evolution of massive stars in singles and binaries [46], and with regard to collisional build up of stars (mass loss at stellar collisions allowed) and intermediate mass black holes [182,183,184]. The current code is available via Github 3 . Note that a different service is provided by Long Wang, quoted here [86]. That alternative version of Nbody6++GPU has been recently used by the Padova team, replacing the older SSE/BSE stellar evolution by MOBSE (see e.g. [185]). Star clusters with primordial (initial) binaries inevitably lead to binaries of black holes. If two black holes get close enough to each other, either during a hyperbolic encounter or due to close Newtonian three-body or four body interactions, Post-Newtonian corrections have to be taken into account. They take the form of an expansion of the relative acceleration between the two bodies in terms of (v/c) 2i , denoted as PNi-terms. PN1, PN2, and PN3 are conservative, producing periastron shifts of orbits, while PN2.5 and PN3.5 provide the energy and angular momentum loss due to gravitational radiation. The first implementation was done in Nbody5 up to PN2.5 [141]; and for Nbody7 [186]. Also relativistic spin-spin and spin-orbit interactions of orders PN1.5, PN2.0, PN2.5 have been recently included [147]. The most recent version of Nbody7 [187] includes also the full PN terms by using the ARChain (algorithmic regularization chain) method [143]. Nbody7 is GPU accelerated, but has not yet the MPI parallelization of Nbody6++ and Nbody6++GPU. Generally binary evolution governed by Post-Newtonian terms has been compared with full numerical solutions of general relativity; deviations between fully relativistic and Post-Newtonian evolution only occur during the final period of merger, in a time span usually negligible for astrophysical purposes. The reader interested in the original citations with regard to the derivation and justification of Post-Newtonian terms is referred to [141,143,147] for further references therein. When black holes merge they experience a kick due to asymmetric gravitational wave emission, see e.g. the MOCCA implementation [188]; a similar model is already used in Nbody7 [187], and this is a field where Nbody6++GPU is currently lagging behind, current work is ongoing on it. The following Table 3 gives a summary of the features of different variants of the Nbody codes. Table 4 shows for Nbody6++GPU a model fit, obtained from a number of simulations using a range of particle numbers N and MPI process number N p , where each MPI process also uses a GPU [190]. Eight different pieces of the code have been profiled as indicated. The fit shows the following key informations: The timing model is already a few years old, the current code version has made progress in MPI parallelization of prediction (pos. 3). To improve the communication scaling faster MPI or NVLink 4 communication hardware will be beneficial (pos. 5,6). Note that all numerical factors in the fit dependent on the specific hardware used -CPUs, GPUs, communication lines between CPU nodes and between CPU and GPU.
: Included in standard version of that level ITS: Individual time-steps [120] ACS: Ahmad-Cohen neighbour scheme [118] KS: KS-regularization of few-body subsystems [115] HITS: Hermite scheme integration method combined with hierarchical block time-steps [128] PN: Post-Newtonian [141,143,186] r: restricted PN, only orbit-averaged energy loss by gravitational radiation [182,183,189] ( ): only included in special version of the code [141] AR: Algorithmic regularization [143] CC: Classical chain regularization [132] MPI: Message Passing Interface, multi-node multi-CPU parallelization [137] GPU: use of GPU acceleration [145] (if also MPI: multi-node many GPU [146])  4 shows a similar information in principle than Table 4, but here the eye should inspect the relative weight of the different components, when increasing the number of MPI processes. The coloured fields correspond to the code parts discussed above, but a little more segmented: a Reg. and Irr. correspond to regular and irregular force computation in Table 4; b Pred. is prediction; c Move is data moving; d Comm.R, Send.R., Comm.I. and Send.I is MPI communication (regular, irregular) e Barr. is synchronization f Init.B., Adjust, KS, refer to sequential parts on the host.
The bottom line to stress from these results is that even for one million bodies the bottleneck of the parallel code is NOT the regular force (which would be extremely dominant in a sequential processing), so it is NOT the stumbling block for going to much higher particle number, these are prediction and communication. There are also phiGPU [146], phiGRAPE [170], and HiGPU [191,192]. All of them are using the Hermite scheme with hierarchically blocked time-steps, and are fully parallelized and GPU accelerated. There is no Ahmad-Cohen neighbour scheme and no regularization, which means that on a serial computer they would be much slower than e.g. Nbody6++GPU. But with a very efficient parallelization and GPU acceleration this is partly compensated; they have been used for astrophysical problems where star-bystar modelling could be neglected, such as e.g. galactic nuclei and galaxy mergers with supermassive black holes (cf. e.g [193,193,194,195,196]). An interesting feature is that these codes implement higher order Hermite integration schemes (in phiGPU 4 th , 6 th , or 8 th order can be chosen, HiGPU uses 6 th order. There is another 6 th and 8 th order Hermite integrator [172]; so far these higher order integrators have seen relatively little use, consistent with the conclusion that the 4 th order integrator is an optimal choice for performance and accuracy [127].

Are N-Body simulations reliable?
At this point the reader may expect that direct N -body simulation turn out to be the most reliable (although computationally most expensive) way to simulate the dynamical evolution of a gravitating system consisting of N point masses. It does not involve any serious approximations and assumptions, as e.g. the Fokker-Planck approximation and the Monte Carlo codes. By reducing the η-values in the time-step (Eq. 45) any accuracy can be achieved in principle, as long as machine accuracy permits it. Usually for accuracy and time-step choice globally conserved quantities are used, such as energy and angular momentum, and center of mass conservation (position and velocity). However, for a system with N particles phase space has 6N dimensions, and a check of say energy, angular momentum, and center of mass alone only checks whether the numerically calculated system remains within an allowed 6N − 9 dimensional hypervolume. There is no a priori information how "exact" the "true" individual trajectories are reproduced in the simulation within this hypervolume. It was early pointed out that, due to repeated close encounters occurring between particles, initial configurations that are very close to each other, quickly diverge in their evolution from each other [197]. In that work it was shown that the separation in phase space of two trajectories increases exponentially with time, or with other words, the evolution of the configuration is extremely sensitive to initial conditions (particle positions and velocities). The timescale of exponential instability is as short as a fraction of a crossing time, and the accurate integration of a system to core collapse would require of order O(N ) decimal places [198,199]. These papers argue that the problem is caused by two-body encounters, but chaotic orbits in non-integrable potentials can be a source of exponential instability and thus cause unreliable numerical integrations as well. However, the situation is not as bad as it seems. N -body simulations of star clusters or galactic nuclei do not always exploit the detailed configuration space of all particles. Quantities of interest are global or somehow averaged quantities, like Lagrangian radii or velocity dispersions averaged in certain volumes. As it was nicely demonstrated in a pioneering series of papers [40,53,200,201] such results are not sensitive to small variations of initial parameters. They took statistically independent initial models (positions and velocities at the beginning selected by different random number sets) and showed that the ensemble average of the dynamical evolution of the system always evolved predictably and in remarkable accord with results obtained from the Fokker-Planck approximation. The method was also partly and successfully used in [54], which focused on the evolution of anisotropy and comparisons with the anisotropic gaseous models of the author of this paper, or in more recent examples [182,183] where the formation of intermediate mass black holes was analyzed over a large set of N -body simulations, using statistically independent initial models. As a consequence, it should be remembered, however, that great care has to be taken when interpreting results of N -body simulations on a particle by particle basis, for example determining rates of specific types of encounters, which could produce mergers in a large direct N -body model. The long-term behaviour of dynamical systems as the solar system are being studied by N -body simulations as well, but clearly there are much higher requirements on the accuracy of the individual orbits in contrast to the star cluster problem. Therefore for the solar system dynamics symplectic methods, using a generalized leap-frog, like the widely used Wisdom-Holman symplectic mapping method [202] are the standard integration method. Symplectic mapping methods do not show secular errors in energy and angular momentum. However, in their standard implementation they require a constant time-step (but see recent new developments described in the following subsection). A generalization using a time transformation simultaneously with the generalized leap-frog has been suggested which can cope with variable time-steps [203]. It has been proposed to reduce secular errors in Hermite schemes and direct N -body simulations to a level comparable with symplectic methods by using a time-symmetric scheme. A small variation in the Hermite corrector is needed, which allows to iterate to convergence (few iterations usually enough) and individual time-steps made reversible through another iteration [204,205,168]. How well this generally works and its relation to symplectic schemes is presently not clear. But it has been well used for direct N -body simulations of planet formation and planetary systems [206,168]. These codes though are still on the level of Nbody4, because they do not use the Ahmad-Cohen neighbour scheme -even in the smallest steps full force calculations over all N particles are needed. Nbody6++ has been similarly improved using an extended Hermite scheme to allow iteration, for a hybrid N -body and Fokker-Planck simulation of planetesimal growth in protoplanetary disks [207,208] (no GPU implementation with Nbody6++GPU yet). In [132] it is stressed that even with a newly applied classical method secular errors in the integration of close binaries can be strongly reduced. One should keep in mind though, that the N -body integration schemes discussed in this paper yield excellent results in the star cluster research (see Sect. 4) but are unsuitable for long-term solar system studies, because they generally have secular errors, although small. Due to the inherently physically chaotic nature of star clusters remaining small secular errors can usually be tolerated. It means that the solution found in the computer always stays near a permitted solution of the underlying Hamiltonian, even if it does not stay on the one trajectory which belongs to the initial conditions [209]. But a recent dynamical study has reiterated that it may not be sufficient just to check a few globally conserved quantities, because that could be dominated by a few high energy objects (binaries) and could cover up errors in other parts of the system [210]. As outlined above in star cluster simulations the secular errors are being kept small relative to typical values of energy and angular momentum and an accurate reproduction of all individual stellar orbits is not generally required.  [150], showing wall clock times obtained on the Juwels-Booster as a function of number of compute nodes, for particle numbers up to eight million and binary fractions of up to 50%. The benchmarks have been done by Qi Shu, and will be published in [219].

New Approaches
A completely new code, called PeTaR has been introduced [150]. It is a hybrid N -body code, which combining the P 3 T (particle-particle particle-tree) method [211,212,213,214] and a slow-down time transformed symplectic integrator (SDAR) [149]. The latter is mathematically similar to a KS [115] regularization, using a time transformation in a similar way (based on the Poincaré transform of the Hamiltonian [215,143]), but regarding the canonical coordinate transformation it uses an extended classical phase space rather than the 4D KS space. Both regularization methods also employ a slow-down procedure. PeTar uses the parallelization framework for developing particle simulation codes (FDPS [213,216]) to manage the particle-tree construction and long-range force calculation. For single and binary stellar evolution the standard SSE and BSE packages are used as in e.g. MOCCA [217,187,46]. The code is conceptually ahead of Nbody6++GPU in several respects; parallelization of a large number of hard binaries is included and a domain decomposition makes it easier to go to particle numbers much larger than 10 6 , as well as the use of the Tree scheme for distant groups of particles, rather than the Ahmad-Cohen neighbour scheme in Nbody6++GPU. We show in Fig. 5 its excellent strong scaling obained on the Juwels Booster supercomputer in Germany [218]. The paper presenting the PeTar code makes some much stronger statements about Nbody6++GPU, claiming it is impossible to implement binary parallelization and domain decomposition in it. We think the given arguments are not sound, a thorough discussion of these issues will follow here in the future depending on further work with Nbody6++GPU. Also the paper presents a couple of convincing comparisons between both codes, but a real long time model comparison with many binaries is still work in progress. The problem is, that the single and binary stellar evolution in Nbody6++GPU has a much stronger coupling to the dynamical evolution of binaries and close encounters -that is from a programmer's viewpoint a trouble (as stated in [150]), but it is important when following very hard and relativistic binaries [182,183]. Even the Tree algorithm is not really essential -inherently it can be much faster than the direct Hermite scheme (even with Ahmad-Cohen neighbour scheme), but it depends on the accuracy required. We have shown by comparison of Nbody6++GPU with the Bonsai Tree code [220] that if same accuracy is required, performance difference is not significant [190]. This can be also seen in Sect. 5.3 from the performance analysis of Nbody6++GPU -even for more than a million bodies the regular distant forces are not the bottleneck, they have been effectively parallelized. Another novel approach is based on forward symplectic integrators (FSI) [221,222,223,148]. It is called Frost [151], uses MPI parallelization and GPU acceleration and shows competitive benchmarks for very large N (million and more). The authors note that the FSI approach has hardly been used in the community, even though it generates very accurate orbital integrations. However, as discussed above, for most star cluster simulations orbit integration with maximum machine precision is not really needed; however, for the innermost regions of nuclear star clusters, with stellar orbits around supermassive black holes, this will be important and Frost may turn out to be excellently suited for such environment. On the other hand the planetary system and protoplanetary disk simulation community regularly uses either symplectic schemes or the improved iterated Hermite schemes (see above Sect. 5.3). A very significant innovation in Frost is Mstar -a fast parallelized algorithmically regularized integrator. Instead of classical and algorithmic chains, which assemble their particles linearly (in a chain), Mstar uses a minimum spanning Tree, so the chain can have branches. That leads to smaller average distances within chain particles as before, reduces computing time and errors. Such an algorithm could be as well used for the other chain regularizations. Finally, hybrid codes have been constantly designed and used for a while to amend the direct N -body for a large number of distant particles. It started already with Nbody5 which was coupled to a Tree scheme [125]; Nbody6++ has been hybridized with a series expansion code [178] (sometimes also called self-consistent field [224], SCF). Last, but not least there is a new hybrid code Etics [225], which has been coupled with phiGRAPE, and applied to the loss cone problem in a nuclear star cluster around a binary supermassive black hole [226].
6 Astrophysics in star clusters

Single Stellar Evolution
In realistic star cluster simulations all stars undergo stellar evolution as time proceeds, see e.g. [227] and therefore, a large array of stellar evolutionary processes must be considered. We briefly outline the fundamentals of single stellar evolution (Sect. 6.1) because it is essential to understand the complexities that need to be modelled before we move on to an area, in which collisional N -body simulations find some of their strongest applications, which is binary stellar evolution (Sect. 6.2) in dense star clusters. The discussion in this Sect. 6.1 is primarily based on [228], but a more recent review may also be found in [229]. A star is a self-gravitating object of a hot plasma, which emits energy at the surface in form of photons (and from the inner regions in the form of neutrinos). Furthermore, it is spherically symmetric in the absence of rotation, magnetic fields and a sufficiently close companion or multiple companion stars that induce interior oscillations and bulges through tidal interaction or deform the star through mass transfer. These are typical assumptions in 1D modelling of single stars and they yield four fundamental structure equations that govern stellar evolution under the assumption of hydrostatic equilibrium. Any deviation from hydrostatic equilibrium will become increasingly important in harder binary stars. Energy transport in a star is either radiative or convective (where convective transport can also include some conduction, which is not that important). Whether it is one or the other is given by the Schwarzschild stability criterion, which compares the temperature gradient in the radiative case with the temperature gradient by an adiabatic movement of matter elements: ∇ rad < ∇ ad . The less practical Ledoux criterion also takes into account a possible gradient in the density and chemical composition of a star. If some matter is unstable according to the Ledoux criterion, then convection will set in and will mix the material until stellar homogeneity. This process will diminish these gradients. Therefore, in practice the Schwarzschild criterion is more commonly used. Radiative energy transport is commonly described using a diffusion equation. For convection, on the other hand, there is a convective gradient and since no complete theory of convection exists, the problem is approximated using mixing-length theory (MLT). MLT describes the convective temperature gradient ∇ c surprisingly well despite a large number of unrealistic assumptions, e.g., all convective "bubbles" are assumed to travel the same distance due to buoyancy forces until they disappear leading to the dispersion of all their energy to the next level. MLT is parameterised globally by α MLT , which is the ratio of the mixing length to the pressure scale height. α MLT is around 2 for Solar models. Since in deep stellar interiors, convection is very efficient and thus the "blobs" move adiabatically, it is ∇ c ∇ ad . In the outermost layers (low densities), however, convection is not so efficient, a lot of energy is lost by a blob moving up and the energy transport is super-adiabatic even approaching the radiative gradient in the extreme: ∇ rad > ∇ c > ∇ ad . The chemical composition of a star changes with time due to nuclear reactions in its interior. It can also be subject to convective mixing, sedimentation, rotation (angular momentum transport) and hydrodynamical instabilities. The inclusion of all of these effects is difficult, because it requires 3D treatment; but most currently used stellar evolution codes, such as Modules for Experiments in Stellar Astrophysics (MESA) [230,231,232,233,234,235] or HOngo Stellar Hydrodynamics Investigator (HOSHI) [236,237,238,239] are 1D.

Two fundamental principles of stellar evolution
The general evolution of a star following the assumptions above is governed by two fundamental principles: the first one is the virial theorem (gravitational energy E g = −2E i total internal energy), which follows from the assumption of hydrostatic equilibrium in the star that is represented by a self-gravitating sphere (Sect. 6.1). The virial theorem implies that on contraction of a star that is modelled as an ideal gas, half of the liberated energy is radiated away and the other half is used to increase the internal energy, which means that the star is heating up. In other words, if stars lose energy from the surface, the star must contract overall and heat up, which is a consequence of its negative heat capacity. That does not mean that some parts like the stellar envelope are not expanding over the star's evolution, but what is certain that the largest part of the star is contracting over the life-time and heating up. Interestingly, massive stars, which are radiation pressure dominated, approach the limit of an unbound structure, which is one of the reasons why they lose mass much more easily. The second fundamental principle is the Coulomb repulsion, which determines the sequence of nuclear burning phases. Due to the virial theorem above that leads to a general increase in the interior stellar temperature, nuclear burning phases follow a sequence of light to heavier elements, i.e. they start with hydrogen (H) burning (the main sequence (MS) phase), followed by helium (He) burning (horizontal branch (HB) phase), the carbon (C) burning phase and so on. This burning sequence stops when an iron (Fe) core is reached, because any further nuclear fusion is endothermic. We reach an "onion-like" stellar structure: in the outer layers original stellar material is still processing (H fusing to He), while at the centre an Fe and Nickel (Ni) core forms simultaneously if the stellar mass is large enough.

Timescales, energy conservation and homology
The following timescales are extremely useful in characterising the evolution of stars: 1. hydrostatic timescale τ hydro : let us assume that the internal stellar forces are not balanced and the star is not in hydrostatic equilibrium anymore. The timescale to return to hydrostatic equilibrium is: τ hydro 1 2 (Gρ) −1/2 . It is extremely short, on the order of seconds for White Dwarfs (WDs), of minutes for the Sun and of the order of days for Red Giants (RGs).
2. Kelvin-Helmholtz (thermal) timescale τ KH : it is defined as the timescale during which the entire internal energy of star would be radiated away by its current luminosity. For the Sun it is on the order of 10 million years.
3. nuclear timescale τ nuc : let us assume that the whole luminosity comes only from the nuclear energy reservoir within the star and that the luminosity stays constant at the current state for the duration of the thought experiment. For the Sun this the emission of all nuclear energy as radiation is on the order of 70 billion years.
In most phases of stellar evolution, we have τ hydro τ KH τ nuc and mostly also even τ KH τ nuc for MS and core He burning (CHeB) stars. In late stellar evolution phases τ KH can approach τ nuc . If we look at the global energy conservation in stellar evolution, we arrive at the homology ("similarity") relations for stars. From these we can derive a mass-luminosity relation that is very fundamental in stellar physics. For MS stars, the homology analysis yields L µ 4 M 3 , where µ is the mean molecular weight (rT ∼ µm). This relation implies that the luminosity does not directly depend on energy generation; also the proportionality factor predominantly depends on the opacity of the stellar material, which in turn is determined by its chemical composition. If the energy generation in the star changes, it will adjust itself such that is has the same luminosity as before. Finally, we discuss the homologous contraction of a gaseous sphere. This analysis yields a relation between the central temperature and central density. For ideal gases, the contraction thereof leads to heating of the gas and for non-relativistic strongly degenerate gases, this contraction leads to cooling in a transition from non-degenerate to strongly degenerate region. This means that low mass stars will never ignite certain elements, because at some stage they become degenerate in the core and the central temperature drops upon further contraction.

Fundamental parameters -mass and composition
While they are incredibly useful to understand fundamental relations in stellar astrophysics, the homology relations (see Sect. 6.1.2) cannot be applied over the full evolution of the star and are typically only applied to MS stars. We need other ways to describe the full evolution of a star. In general, the fundamental parameters of stellar evolution are the zero-age MS (ZAMS) mass and the (homogeneous) chemical composition. Other very important parameters independent of mass and composition are rotation and magnetic fields. Rotation can lead to additional interior mixing, which changes the chemical composition of the star. Magnetic fields may influence the pressure balance and interact with convection and rotation, which is probably most important for massive stars.

Mass change of stars -stellar winds
The masses of all stars change throughout their lives through winds, parameterised by a stellar mass loss rateṀ . Stellar winds are the outflows of matter leaving the stellar surface with an energy sufficient to escape from the star's gravity. The main question is what the nature of the force is that is powerful enough to overcome the star's gravity. Different types of stars have different winds. Recently, excellent reviews of the winds of lower mass stars were written by [240] and similarly of high mass stars by [241]. Many stellar evolution models used inside N -body codes express wind acceleration by a Γ factor. Γ is defined as the ratio of radiative over gravitational acceleration. Radiative acceleration is due to radiative pressure and introduces an extra force acting on a spherically symmetric, isothermal wind. It is related to electron scattering Γ e or dust scattering Γ d , for example. These quantities are introduced into the momentum equation of an isothermal, spherically symmetric stellar wind, which leads to an effective gravitational acceleration g eff (r). Using g eff (r), we can calculate the escape velocities and these are lower by the introduction of the extra force. However, it depends very strongly on the distance to the stellar surface, where this additional force is introduced; the farther out it occurs, the less impactful it becomes on the overall stellar mass loss rate. Therefore, since dust grains form very close to the star (e.g. in CMSs), these are very impactful on the mass loss rate. In red supergiants (RSGs), on the other hand, these grains form much farther out and therefore, dust-driven winds are generally not relevant here. Moreover, radiation transport and the chemistry in the wind are both essential to a full modelling of a stellar wind. It is important to state that in general, there is no full theory of stellar winds available [240]. Furthermore, the layperson is overwhelmed by the large number of mass loss rate prescriptions derived predominantly from observations, which differ enormously in magnitude and slope [240]. The choice of mass loss recipe has an enormous impact on the outcome of realistic N -body simulations and the dynamics of the star cluster as described in this review. As an astrophysical community, we are just at the beginning of unravelling the complexities of specific stellar winds, such as Wolf-Rayet (WR) stars [242] or the impact of pulsations and variability on winds in AGB and post-AGB stars [243] before a fully self-consistent theory can be envisioned.

Formation of compact objects and their natal masses, kicks and spins
Depending on the progenitor star core mass, a compact object such as a white dwarf (WD), neutron star (NS) or black hole (BH) may form. Oftentimes binary processes are involved [250,251,252,253], but these are discussed in the next sub-chapter. The following processes apply to all single stars in the relevant mass ranges. The formation of a compact object is associated with a natal remnant mass, a natal kick and a natal spin, which are all subject to significant theoretical and observational uncertainty. Nevertheless, it is important to model these as accurately as possible, because the global dynamical evolution of a collisional stellar system critically depends on these. The natal mass depends on a number of factors and we will only focus now on the collapse mechanism and its associates fallback and not the mass loss in the progenitor star, although it is also instrumental. The impact of the mass loss has been discussed already in a previous section. Traditionally, the natal masses of the WDs (and their three main types HeWDs, COWDs, ONeWDs) and their dependence on the progenitor masses are modelled by [138,254]. For NSs a maximum mass of around 2.5 M [255,256] and the relationship follows typically [138], but the exact masses are unknown because of the large uncertainties mainly in the internal structure of a NS [257,258]. In addition to [138], the possibility of a so-called electron-capture SNe (ECSNe) that leads to the formation of a NS [259,260,261,262,263,264], which has very important properties that are discussed below, has been included in many codes [265,187,244]. Most attention has arguably been paid to the remnant BH masses [266,265,249] and a number of collapse mechanisms for certain mass ranges have been proposed.
In simulations the most popular prescriptions are the rapid or delayed corecollapse SNe models by [249] in combination with various (pulsational) pair instability SNe ((P)PISNe) stellar evolution recipes [267,268,269,270,271,245,246,247]. Fig. 6 shows a suite of small simulations when McLuster [272,244,273] is used as a population synthesis tool with level C stellar evolution [244]. It shows all relevant remnant mass phases, which can be sub-divided into a core-collapse SNe, PPISNe, PISNe and a direct collapse phase in increasing ZAMS mass (this is an extension of the core-collapse SNe models for ZAMS masses above which PISNe is ineffective; in our case an extension of the rapid SNe models by [249]). Two interesting conclusions can immediately be drawn here: first, the metallicity is incredibly important for the production of high mass BHs, because progenitor stars with high metallicities will contain more metal lines for radiative wind mass loss. Secondly, the (P)PISNe prescriptions available from theory can have an enormous impact on the abundance of BHs. This might particularly important in Population III (Pop-III) star clusters, where intermediate mass black hole (IMBH) progenitor stars are postulated to have large enough masses and crucially also low enough metallicities from birth to evolve by (P)PISNe from interior evolution alone (e.g. [82,274] for recent N -body simulations of these clusters; see Sect. 6.4.3 for a more general McLuster version uses level C stellar evolution [244]. Shown are the recipes for the "strong" (psflag=1) on top [245], "weak" (psflag=2) in the middle [246,247] and the "moderate" (psflag=3) (P)PISNe on the bottom [246,247]. The ZAMS stars suffer wind mass loss via metallicity-dependent winds (mdflag=4) (no bi-stability jump) from [248] and the core-collapse SNe are set to "rapid" [249] (Figure from [244]).

discussion of Pop-III stars in the initialisation of star cluster simulations).
The magnitude of natal kicks results, broadly speaking, from an inherent asymmetry in the SNe process and generally their magnitude is rather uncertain [276,275]. They affects the dynamical stability of a binary (if one of the binary stars is forming a compact object) and are even able to disrupt a binary completely. This also implies that a large amount of gravitational binding energy in binaries may be removed from the cluster in this way and this will consequently impact the global cluster evolution. WDs are associated with low velocity kicks of the order of 10 0 km s −1 [277], while neutron stars may reach kicks above even 10 3 km s −1 except in the case in which they form as a result of an electron-capture SNe (ECSNe). Here they receive then kicks of order of only 10 0 km s −1 [278] meaning that they can be retained in a star cluster (simulation) [279,280,281,282,244]. Natal kicks for BHs and NSs, which do not undergo an ECSNe (or AIC or MIC, see next chapter), receive kicks traditionally scaled by fallback during the SNe in simulations [265,249]. The more fallback of stellar material there is onto the proto-remnant core, the lower the resulting kick is. Furthermore, in simulations, it is typically assumed that the asymmetry is produced by a dominant process [187,283]: convectionasymmetry driven kicks [284,285,286], collapse-asymmetry driven kicks [287,288,289,290] or neutrino-driven natal kicks [291,292,187,283]. These lead to different retention fractions of BHs in star cluster simulations [187], which can be seen in Fig. 7 for a sample of Nbody7 simulations from [187]. It is apparent that for these settings the postulated collapse asymmetry driven kicks will produce most (stellar mass) BHs below v esc of the cluster. The natal spins of compact objects are important in general binary evolution and can also have significant impact on the mergers of compact objects, for example in a BH-BH merger [188,293]. In the following, we focus on BHs, but the same arguments can be extended to NSs and WDs and the discussion is largely taken from [244]. In general, the spin angular momentum of the parent star does not necessarily translate directly into the natal spin angular momentum of the BH upon collapse. To quantify the spin, [294] define a dimensionless parameter a spin that accounts for the natal spin angular momentum. [283] assumes that the magnitude of a spin for the BHs is set directly at the moment of birth without any related mass accretion of GR coalescence processes.
In the following we highlight three natal spin models that are available now in Nbody7, Nbody6++GPU, McLuster, PeTar and MOCCA, see also [244]. The simplest model of BH natal spins, the Fuller model, produces zero natal spins [283] as here the Tayler-Spruit magnetic dynamo can essentially extract all of the angular momentum of the proto-remnant core, leading to nearly non-spinning BHs [297, 298,299]. The second spin model is the Geneva model [300,301,283]. The basis for this model is the transport of  [187]. A metallicity of Z=0.0001 is assumed here. The models feature rapid core-collapse SNe from [249] and strong (P)PISNe from [245]. Due to the logarithmic vertical axis, direct-collapse BHs with a fallback fraction, f fb = 1 and v kick = 0 are not shown in these panels. The sharp drop in v kick with increasing m CO or m rem is the approach towards direct collapse. The typical v esc for the M cl (0) 5.0 × 10 4 M and r h (0) 2 pc clusters considered here (blue, solid line). The velocity dispersion of the Maxwell distribution from all the kick models are scaled is 265.0 km s −1 from [275]. It is apparent that for these settings the collapse asymmetry driven kicks will produce most (stellar mass) BHs below v esc of the cluster (Figure adapted from [187]). , as generated by Nbody7 in [283]. Top panels: the N -body models corresponding to these panels employ the "Geneva model" of [295] for BH spin and comprise only single stars initially, whose ZAMS masses range from (0.08 − 150.0) M and which are distributed according to a standard Kroupa IMF [296]. The models feature rapid corecollapse SNe from [249] and strong (P)PISNe from [245]. The models are shown for four metallicities, Z= 0.0002, 0.001, 0.01, and 0.02 as indicated in the legends. Bottom panels: the N -body models corresponding to these panels employ the "MESA model" of [295] for BH spin. The other model characteristics are the same as those in the top panels except that the "weak" PPISNe mass prescription [246] is utilized (resulting in the non-monotonic behaviour with respect to M BH , which, here, extends up to 50 M as opposed to the models in the top panels, where M BH is capped at 40.5 M due to the use of [245] (Figure adapted from [283]).
the angular momentum from the core to the envelope. This is only driven by convection, because the Geneva code does not have magnetic fields in the form of the Taylor-Spruit magnetic dynamo. This angular momentum transport is comparatively inefficient and leads to high natal spins for low to medium mass parent O-type stars, whereas for high mass parent O-type stars, the angular momentum of the parent star may already have been transported away in stellar winds and outflows and thus the natal BH spins may be low. The third and last spin model is the MESA model, which also accounts for magnetically driven outflows and thus angular momentum transport [297,230,232,299,283]. This generally produces BHs with much smaller natal spins than the Geneva model described above. The Geneva and the MESA models and their metallicity dependence are shown in Fig. 8.

Binary stellar evolution
In addition to the astrophysical processes that affect all stars in isolation, the proximity (orbital period P orb ≤ 10 4 days [302]) to another star or compact object through the frequent encounters in collisional stellar systems or through intrinsic binary evolution, can affect the individual stars or compact objects dramatically and we need to account for these in the simulations. A population synthesis code should include them all [303].

Stellar Spin and orbital changes due to mass loss or gain
If two stars are in a binary, they can transfer mass via stellar winds and therefore also transfer angular momentum even if they are not yet undergoing Roche-lobe overflow (RLOF) [139,303,304]. If a secondary star accretes mass by passing through the wind of the primary star, it is spun up intrinsically by a fraction of the spin angular momentum that is lost by the donor star. The accretion rate is traditionally modelled by [305], which depends on the wind velocity v W . This quantity is observationally difficult to determine [240] and should be proportional to the escape velocity from the stellar surface of the star [139]. The mass variations between companion stars also changes the orbital parameters of the binary star. In general, the eccentric orbit is circularised as a result of mass transfer being more effective at periastron than apastron. Additionally, the accretor star is slowed down by the drag induced by the wind it passes through and this dissipates angular momentum away from the system. The orbital circularisation timescale τ circ as result of mass transfer is orders of magnitudes larger than the equivalent timescale caused by tidal friction for the same binary star system (see Sect. 6.2.2).

Effects of tidal damping
Observations show that the rotation of close binary stars is synchronised with the orbital motion without any dynamical mass transfer having taken place [306,307,308]. Therefore, there must exist a torque that transfers angular momentum between the stellar spin and the orbit in such a way that the binary approaches the observed equilibrium state that is characterised by corotation (spin-orbit synchronisation timescale τ sync ) and a circular orbit (circularisation timescale τ circ ) [309,310,139,304]. Alternatively, dissipation of energy might also lead to an accelerated in-spiral of the binary stars [311,312,304]. When two binary star members are detached but sufficiently close, tidal interaction between them becomes important. The mere presence of a companion star causes a tidal force that elongates a star along the line between the centres of mass, thereby resulting in tidal bulges (see e.g. [139]). When the binary component rotates uniformly with a circular orbital motion, then the tidal bulges on its stellar surfaces are steady and the stars are in hydrostatic equilibrium. In such a scenario, we also speak of equilibrium tides. However, when this condition no longer holds, the hydrostatic equilibrium is disrupted and the star undergoes forced stellar oscillations. This scenario is described by a combination of equilibrium and now also dynamical tides, the latter of which produce much smaller tidal bulges than the former and they can also take any orientation [313,303,139,314,315,316,309,317]. Dissipative processes within a star cause the tides to be misaligned with the line of centres. This results in a torque that transfers angular momentum between the stellar spin and the orbit [139]. This dissipation is non-conservative and happens on relatively long timescales [303]. The dissipative processes within a star depend on the stellar structure. Typically, a distinction is made between stars with appreciably deep convective envelopes and stars with radiative envelopes. The tides dissipate energy and the binary system approaches an equilibrium state that is characterised by a circular orbit and corotation [309,310,139,318].
In stars with appreciably deep convective envelopes, turbulent viscosity that acts on equilibrium tides (the same effect on dynamical tides is negligible [316,309]) is the most efficient form of dissipation [319,310,139]. The dissipation takes shorter than the nuclear burning timescale τ nuc (see Sect. 6.1.2) [320,321,139].
In stars with radiative envelopes, radiative dissipation near the surface of the star causes an asymmetry in the internal stellar oscillations induced by tides and the tidal field itself. This leads to a torque that is necessary for the binary system to approach the equilibrium state [309,320,322,139] and in sufficiently close binaries this happens on shorter timescales than the nuclear burning timescale τ nuc [316]. This radiative damping on the dynamical tides is the most efficient process to achieve the equilibrium state in binary stars with member stars that do not have an outer convective zone. However, if they do then the aforementioned turbulent friction on the equilibrium tides provides the primary torquing [316,309,320]. τ sync and τ circ in binary stars with convective envelopes are typically orders of magnitude smaller than those with radiative envelopes [309,139]. τ sync and τ circ are generally not equal except in a limiting case [309].
If the stars are degenerate but have sufficient stellar structure, i.e. WDs and NSs, then the above two dissipative mechanisms cannot be used as the stellar structure is significantly different. WDs will have very low spins, because the progenitor AGB star has already spun down in its expansion. Furthermore, in WD-WD binaries, the orbit will already be circularised (in the absence of WD natal kicks [277]) due to the stellar wind mass and thus angular momentum loss. For this reason only the synchronisation timescale τ sync due to degenerate damping is of importance here and it is only applicable for close systems. τ sync in WD-WD, WD-NS and NS-NS binaries could exceed the age of the Universe [323].

Dynamical mass transfer and its stability
Apart from mass transfer through stellar winds, mass transfer can also happen via RLOF. This happens when the primary star fills it RL as a result of stellar expansion or in-spiral. The subsequent mass transfer then happens through the innermost Lagrange point. Typically, this process depends strongly on the mass ratio of the binary [324] and happens in corotating, circularised binaries but in some instances, it can also occur in highly eccentric binaries, that are a result of tidal capture.
In the RLOF mass transfer, also angular momentum is transferred. The stability of the mass transfer traditionally determined by three logarithmic derivatives of radii with respect to the mass of the RL-filling star following [325,326]: the derivative describing the rate of change of the RL radius R L for conservative mass transfer (total mass and angular momentum conservation) ζ L [303], the derivative at constant entropy s and composition of each isotope X i throughout the star ζ ad and a third derivative that describes the rate of change of the radius of the primary with mass in equilibrium ζ eq . The mass transfer rateṀ depends on the relative values of these derivatives [304]: 1. ζ ad < ζ L →Ṁ increases rapidly, there is positive feedback and the mass transfer is unstable, the secondary star cannot accrete at such a high rate and it expands → formation of a common envelope (CE) around the two stars [327,328,329]. 2. ζ eq < ζ L < ζ ad →Ṁ decreases in its immediate response, but then expands on a thermal timescale. 3. ζ L < ζ ad & ζ L < ζ eq →Ṁ decreases initially, because the stellar radius decreases. RLOF happens again, when the primary fills it RL again.
On these basis of these exponents alone, it is possible to make a number of arguments on the evolution of Cataclysmic Variables (CVs), Algols and other exotic binary stars.

Common-envelope evolution
The process of CE evolution (CEE) is instrumental in compact binary and close binary formation. [328,330,331,329,327]. A CE is the outcome when ζ ad < ζ L in RLOF or when two stars collide, where one of the stars has a dense core. Generally, CEE occurs when the primary star transfers more mass on dynamical timescales than secondary can accept. It strongly depends on the instabilities in the RLOF preceding the formation of a CE [332]. The CE expands and thus rotates more slowly than the orbit of the secondary and primary star. This causes friction, the binary spirals in and transfers orbital energy to the envelope. Either so much energy in this process is transferred that the envelope is expelled completely leaving behind a close binary in corotation or in the process of in-spiral the binaries coalesce [303,139,129]. The CE is traditionally modelled with the "αλ" energy-formalism [333,129,139], which assumes energy is conserved and where α (α < 1 if no other energy sources other than the binding and orbital energy are present; it can be as high as α = 5 otherwise [334]) is the "efficiency" of the energy re-use and λ is a measure of the binding energy between the envelope and the core of the donor star and should depend on the type of the star, its mass and its luminosity [335,336,329,332]. This picture is very simplistic and does not take into account the myriad of processes that go on during CEE, which are also not fully understood yet [328,337,329,338]. On the other hand, the αλ energy-formalism is computationally very efficient and therefore it is widely used in population synthesis codes that require fast and robust stellar evolution computations [139,265,336,339,340,341,244]. Some of these also allow for recombination energy of hydrogen in the cool outer layers of the CE being transferred back into the binding energy of the CE.
Recently, a new formalism has been developed by [342], which solves a binary orbit under gas friction with numerical integration. This means that the authors do not approximate CE as an instantaneous process, unlike in many binary population synthesis (BPS) codes around. The new formalism, which can be easily implemented in BPS codes, provides a significant upgrade, which can explain observations of post-CE binaries which non-zero eccentricities [343]. In a binary consisting of a NS or a BH and a giant star, after the CE has been ejected and if the binary survives this phase, the H-rich envelope of giant stars might be stripped completely off. Now, the binary consists of a BH or a NS orbiting a naked He star. There might now be subsequent mass transfer from the naked He star to the NS of BH. This post-CE RLOF mass transfer leaves behind a so-called "ultra-stripped" He star that explodes in an ultrastripped SNe [344,345,346]. This type of SNe is significantly different from the typical core-collapse SNe and the process of ultra-stripping leads to a significant decrease in BH-NS and BH-BH mergers and a slight increase in NS-NS mergers [347].

Mergers and general relativistic merger recoil kicks
An outcome of CEE may be the coalescence of the two stars. The subsequent merger product depends on the relative compactness of the two stars and Fig. 9: Left hand side: Conceptual picture presenting components of the recoil velocity. Dashed lines represent a Cartesian coordinate system in the orbital plane: e 1 and e 3 . A vertical dotted line is a line perpendicular to orbital plane (e 3 , parallel to orbital angular momentum). The red vector is the kick component related to spin asymmetry, and magenta vectors are its projections on the plane and parallel to e 3 . The blue vector represents the mass inequality contribution. The black filled circles represent a pair of BHs, their spins and orientation in a spherical coordinate system are illustrated. This drawing also reflects typical proportions between recoil velocity components. Right hand side: Example of how each recoil velocity component depends on mass ratio q for a metallicity dependent spin vector from [348]. q is the only variable for the determination of v m . Other components and the overall kick velocity depend also on spin magnitudes and orientations, in this case the mean value is plotted. As we can see, the major component is almost always v , others only play a role for low q. (Figures and captions adapted from [188] (combination of Fig. 2 and Fig. 3)).
thus it depends on the stellar evolutionary stage [129,139]. If similar in stellar type, then the two stars mix completely. If one is much more compact than the other, then the more compact core sinks to the centre and the other mixes with the envelope. An unstable Thorne-Żytkow object is created if the merger involves a NS or a BH [349]. Detailed calculations on the merger outcomes following coalescence and collisions, which are less likely than coalescence, but still relevant in star clusters (e.g. [182,183]), depending on the initial stellar types have been tabulated in [139]. A coalescence for our purposes here means that at least one of the members is a star with a core and that the binary has a circular orbit before merging, while a collision means an actual physical collision, where none of the binary members is an evolved stellar type, but the member can also be a compact object. Generally, the mixing and the final masses of the merger products are highly uncertain and only approximations can be made according to our current knowledge [350,244]. There are recent attempts to unravel the masses and compositions of merger products of mas-sive stars with hydrodynamical codes [351,352] and they can be used to give approximate formulae for N -body or BPS codes in the future. The merger of compact objects is associated with a general relativistic (GR) merger recoil kick due to the asymmetry in the GW (see also [244] for a brief discussion with respect to Nbody6++GPU and MOCCA). The recoil velocity in this process depends on the mass ratio of the two compact objects and their spin vectors [353] and can reach several hundreds km s −1 on average [188,293], which is much larger than typical star cluster escape speeds. Fig. 9 (from [188,293]) shows the a conceptual picture of the geometry of a GR merger recoil kick in a BH-BH merger and the dependence of the mean recoil velocity on the mass ratio q of the two BHs for a metallicity dependent spin model from [348]. It can be seen that q has a huge impact on whether a GR merger recoil kick velocity exceeds the escape speed of the surrounding stellar (and gaseous) material or not. Equal mass mergers might be retained in nuclear star clusters [354] and extreme mass ratio mergers might theoretically even be retained in open clusters (although IMBHs will probably not form there) [355,356,357,358]. For (nearly) non-spinning BHs (Fuller model [298]), the kick velocity is smaller than for high spins. For non-aligned natal spins and small mass ratios, the asymmetry in the GW may produce GR merger recoils that reach thousands of km s −1 [356,359]. The calculation of the mass ratio is straightforward and the spins may be calculated from e.g. [360] or [361]. Generally, the orbital angular momentum of the BH-BH dominates the angular momentum budget that contributes to the final spin vector of the post-merger BH and therefore, within limits, the final spin vector is mostly aligned with the orbital momentum vector [283]. In the case of physical collisions and mergers during binary-single interactions, the orbital angular momentum is not dominating the momentum budget and thus the BH spin can still be low. [283] also includes a treatment for random isotropic spin alignment of dynamically formed BHs. Additionally, [283] assumes that the GR merger recoil kick velocity of NS-NS and BH-NS mergers [362,363] to be zero but assigns merger recoil kick to BH-BH merger products from numerical-relativity fitting formulae of [359] (which is updated in [364]). The final spin of the merger product is then evaluated in the same way as a BH-BH merger. The inclusion of these kicks in direct N -body simulations is still unusual (e.g. [365,366,367,185,182,183,244,46] all do not include these in addition to missing PN terms), but it is worth mentioning [184] do include the GR merger recoil kicks by posterior analysis. Nbody7 [186,187,283] on the other hand does include GR merger recoil kicks based on [353,360]. In MOCCA numerical relativity (NR) models [368,369,370,359,361] have been used to formulate semi-analytic descriptions for MOCCA and Nbody codes [188,293,283,371,184,364]. Recently, GR merger recoil kicks have also been added to Nbody6++GPU [189] following [368,361] and with this code version, the whole kick process can be followed self-consistently.

Accretion or merger induced collapse
In sufficiently close double degenerate COWD-COWD, ONeWD-ONeWD or COWD-ONeWD binary stars, sufficiently high and dynamically stable RLOF mass accretion of hot CO-rich matter may lead to a heating of the outer layers of the secondary, which will result in the ignition of nuclear burning [372]. If carbon burning is ignited in the COWD envelope, the heat will be transported the stellar core by conduction and then the secondary will evolve into an ONeWD [373,374], which will eventually collapse into a NS if the critical mass of the ONe core is surpassed (M ecs =1.38 M ) [259,260,265]. This ONeWD collapse is referred to as accretion induced collapse (AIC). If, on the other hand, the ignition happens in the centre then the star will undergo a SN-Ia explosion, which leaves no remnant behind. Double degenerate COWD binaries may also coalesce without undergoing dynamically stable mass transfer. During this process the less massive forms a thick, turbulent accretion disk and the more massive COWD will accrete matter close to the Eddington limit. Here the carbon will be ignited on the envelope of the secondary and thus the outcome will be a ONeWD and no SN-Ia will happen [372]. Again, if the ONe core mass surpasses M esc , then the ONeWD will collapse into a NS and this is known as a merger-induced collapse (MIC). Other pathways for MIC are mergers of a ONeWD with any type of WD companion if the resulting merger mass surpasses the critical mass for NS formation [374,265]. The distinction between AIC and MIC is made, because the former may be observed already through their stable mass transfer phase or in low-mass Xray binary stars and the latter may be observed through gravitational waves observed with LISA (Laser Interferometer Space Antenna [375]).

Gravitational radiation and magnetic braking
Gravitational radiation emitted from sufficiently close binary stars (P ≤ 0.6 days) transports angular momentum away from the system and drives it to a mass transfer state that might result in coalescence [376,139,303]. The effect this radiation has on the orbit of the binary (excluding PN terms) may be obtained by averaging the rates of energy loss and angular momentum loss over an approximately Keplerian orbit [377,303]. Gravitational radiation will circularise the orbit on the same timescale as the orbit shrinks until coalescence. In co-rotating and sufficiently close binary stars, magnetic braking slows down the rotation of the individual star with a convective envelope, but also drains angular momentum from the orbit of the binary star, because tidal friction between the stars may conserve co-rotation [378,379,380,303]. As a result, this process will force a close binary to a state of RLOF within Hubble time.
In some situations, this process is dominating binary evolution, such as in CVs above the orbital period gap [381,382,383]. In spin-spin period evolu-tion (P −Ṗ ) of pulsars this process is also important (e.g. [384,385]). Both processes outlined above are non-conservative.

Combining stellar evolution with collisional N-body codes
There are two main methods that stand out in practice concerning the integration of the complicated stellar evolution into N -body codes. Both of these, interpolation between tables or approximation of stellar evolution data by some interpolation (fitting) formulae as functions of mass, age and metallicity, has unique advantages and disadvantages that have been known for a long time [302]. As it stands now, the two approaches are not in competition, but rather complement one another [138].

Interpolation between tables
This method calculates stellar parameters from detailed evolutionary tracks (e.g. [386]). These evolutionary tracks are derived from 1D stellar evolution codes and are in tabular format. They are necessarily rather large and therefore, this approach has historically been limited by memory availability on hardware [302,138,387]. Unlike fitting formulae, stellar parameters from the given set of detailed tracks are calculated in real time with this method. Hence, one just needs to change the input stellar tracks to generate a new set of stellar parameters. It has been claimed that this approach is the most flexible, robust and efficient today when combining detailed stellar evolution with stellar dynamics [387]. [388,389,390,391,392,393,394,395] constructed such tables, which were later then expanded upon and refined by [386]. In the aforementioned works, the convective mixing or overshooting length l OV presents another hurdle, which describes the average distance by which convective cells push into stable regions (or radiative regions from Schwarzschild condition [396,397]) beyond the convective boundary [389,386,398]. This treatment was modified by [386] and replaced with a "∇ prescription", which is based on the stability criterion itself (δ OV = 0.12 was found to best reproduce observations [399,400,386,138]). This new criterion avoids physical discontinuities for disappearing classical convective cores. Further quantities that will influence the calibration of the luminosity L of a stellar evolution model are the nuclear reaction rates and the core Helium abundance Y . Another source of large uncertainty was left largely unchanged by [386]. This uncertainty has been described by the [386] as the "Achilles heel" in stellar evolution codes. This uncertainty is in the mixing length of α MLT , which is derived from mixing-length theory [401] to describe heat transport in the convective regions of stars [398,402] (see also Sect. 6.1). [386] set α MLT =2.0 (based on the Solar model). But not all stars with convective regions exhibit identical convective properties and α MLT can show large variations from star to star [398].
Even today methods stellar evolution by interpolation between tables are being developed with increasing success as hardware memory capabilities also improve: -SEVN [403,269,404], which has been completed for binary evolution [405] has been used extensively to study the evolution gravitational wave source progenitor stars. Additionally, it is not available as SEVN2.0, which is integrated in PeTar [150]. and COMBINE [406] codes, which also has binary evolution implemented [407,343] has also been used extensively to study the evolution gravitational wave source progenitor stars. -METISSE code [387], which is based on the STARS [408,409,410,411,412,400,399], MESA [230,231,232,233,234,235] and BEC [413,414,415,416,417,418]. Unlike SEVN or COMBINE, this code does not yet account for binary stars. In general, METISSE will be another promising candidate for combining full stellar dynamics with detailed stellar evolution.

Interpolation/Fitting formulae
A first attempt to incorporate simple stellar evolution fitting formulae in a direct N -body code was done by [419] on the basis of [420]. Later, as a successor to [420] was created using the method developed by [386]. They based their code on the original Cambridge STARS stellar evolution program by [408,409,410,411,412,400,399]. The result are the famous Single Stellar Evolution (SSE) fitting formulae, which for the first time included metallicity as a free parameter [138,421,422]. Fig. 10 shows the complex discretization of stellar phases and the possible evolutionary pathways between them in the SSE package. The figure has been included, because this fundamental structure still remains in many stellar evolution production codes today (see below).
In general, such fitting formulae take much more care and thus time to set up than method of interpolating between tables [227], because the movement of a star in the Hertzsprung-Russell-Diagram (HRD) is highly non-uniform and erratic. Furthermore, they are also less adaptable to changes in stellar tracks, for example, when they need to be adjusted due to some new development in astrophysics. On the other hand, the SSE provides us with rapid, robust and analytic formulae, which can be easily modified and integrated into an N -body code along the lines of [419] and give stellar luminosity, radius and core mass of the stars as functions of mass, metallicity and age for all stellar evolutionary phases [138,423]. However, these formulae necessarily also discard a lot of crucial stellar evolution information [421]. For example, stellar mixing depends on several timescales and internal stellar structure parameters [350] and so these cannot be modelled directly by the fitting formulae. Only the outcomes can be parameterised for stellar types of the individual stars along the lines of [139]. Despite these fundamental complications in stellar evolution modelling that persist to this day (see e.g. [398,402,424,425]) and which translate directly into the continuous and differentiable fitting formulae (polynomial form from least square fitting [138]), the SSE code has successfully, for the first time, provided us with a method by which we can evolve stars from ZAMS masses (0.1-100) M (the models from [386] only reach 50 M , but the SSE formulae can be safely extrapolated to 100 M [421]) rapidly and accurately (within 5% of detailed stellar evolution models over all phases of the evolution [138]) in N -body simulations throughout all evolutionary phases taking into account all of the astrophysical processes outlined in Sect. 6.1 and offering a metallicity range from 0.0001 to 0.03 with Z 0.02 being Solar metallicity as an input parameter. However, for a complete picture we also need to model the binary evolution processes outlined in Sect. 6.2. For the fitting formulae this is provided by the Binary Stellar Evolution (BSE) code [139,426,427], which is an addon of the SSE package. This has been a huge success story and many full dynamical cluster simulations have utilised SSE & BSE to evolve the stars, e.g. [48,428,365,366,367,185,182,183,244]. The SSE & BSE codes have been the foundation for many other BPS codes: [431,432,433] and related code called ASPS also used in [434], -StarTrack [435,265,295] -COSMIC [341] and its implementations in CMC [436,108] -BSE-LevelC [244] and its implementation in McLuster [272,244].
The fitting formulae from the SSE code are also implemented in BPS code binary c [437,438,439]. New fitting formulae have recently been constructed, which are derived from fitting to 1D HOSHI stellar evolution models [236,237,238,239] to extremely massive low metallicity (EMP; Pop-III) stars [440,441,442,443]. These are constructed such that they can be implemented into any of the BSE variants mentioned above in a straightforward fashion and therefore also into stellar dynamics codes such as Nbody6++GPU [47].

Global star cluster initial conditions
Defining appropriate global initial conditions for star cluster simulations is highly non-trivial as the formation of a star cluster and the stars within it depend on a large number of parameters that are very uncertain due to a lack of better theoretical understanding and or observations. In the following, we give an overview of the most important parameters in this context for N -body simulation of star clusters.

Initial 6D phase space distribution
In order to initialise an N -body star cluster simulation, we need to distribute the N particles in 6D phase space. A statistical approach as described in Sect. 3.1 is taken to realize a star cluster, which follows the probability density distribution f ( r, v, t). The full 6D distribution function is rarely known explicitly; under the assumption of steady state, such that f does not depend on time, the Jeans's theorem [58] allows us to express f as a function of integrals of motion of a single star moving in the gravitational potential Φ(r). For now we assume spherical symmetry, so we have for example specific energy and specific angular momentum: f = f ( r, v) = f (E, L), which are defined as E = v 2 /2 + Φ(r) and L = | r × v| (cf. Sect. 4.2). Deviations from spherical symmetry can be taken into account as well, see for example Sect. 6.4.6 for the importance of initial bulk rotation. Examples of such self-consistent distribution functions are given by where n is an integer index and F n a normalization factor to make sure that f (E) is properly normalized as a probability density function. For n = 5 this results in the famous Plummer model [444], and for n = 7/4 another famous solution, a density cusp [445,446] around supermassive black holes is found [447]. In [58] these models are also called stellar polytropes, because their density distribution is the same as a gaseous polytrope [448] of the same index n. Analytical density distributions exist for n = 0, n = 1, and n = 5 [228], but for stellar systems only n = 5 is physically useful. The theory of gaseous spheres also knows the isothermal solution, which is obtained for n = ∞; in stellar dynamics the equivalent is the isothermal sphere Here σ 2 is the r.m.s. stellar velocity dispersion, analogous to the temperature in a gaseous sphere. These models have some problem, because their radial extent is unlimited (Plummer and isothermal) or even their mass is infinite (isothermal). Therefore, and since real star clusters are often subject to a tidal cutoff due to the host galaxy, a cutoff radius is introduced (connected to a cutoff energy). If at the cutoff radius the gravitational potential of an isolated star cluster would be Φ 0 , then a relative potential Ψ and a relative energy ε are defined by In that way the cluster extends from the center out to ε = 0 (and Ψ = 0), and lowered isothermal or Plummer distribution is defined as follows: Again f 5 and f ∞ have to be properly chosen normalization factors. The model Eq. 51 is the widely used King model [75] ("lowered isothermal"). Note that some papers and books also prefer to change the sign of E (or ε), such that bound objects have a positive value. We do not follow this here to avoid confusion. Even in spherical symmetry the distribution function could be 2D, since we have E and |L| as constants of motion; it corresponds to the possibility that in spherical star clusters still at any given radius r the radial and tangential velocity dispersion may be different. So, a more general approach for the distribution function in case of an isothermal is which is also known as Michie [449] distribution.  [451]. We note that [452] developed a new family of lowered isothermal models called the limepy models. Based on the 1D models of King a generalization in 2D for rotating star clusters is now being used and often described as rotating King models (see Sect. 6.4.6 and citations there).
Note that from f (E) or f (E, L) directly a numerical star cluster cannot be constructed, because E depends explicitly on r and implicitly through the gravitational potential (Eq. 30). Therefore, in order to be self-consistent, the gravitational potential has to be determined by a velocity space integration over the distribution function and then Poisson's equation solved to obtain the stellar density as function of radius (see e.g. the textbook [58] for examples). In a final step a random procedure has to be used to obtain stellar positions and velocities. If density or gravitational potential are analytically known functions (like in case of a Plummer model) the entire self-consistent model can be constructed in one loop using random numbers [52].

Initial stellar mass function
In order to initialise star cluster simulations, we need to draw the ZAMS masses from an assumed distribution. For this purpose, we use an initial stellar mass function (IMF), a "Hilfskonstrukt" [453,454], as a mathematical formulation of an idealised stellar population that has formed from a singular star formation event. We will discuss in Sect. 6.4.5 that this is not the case in nature. An excellent review on the IMF and its construction has been provided by [455] (and sources therein) and it also explores the universality of the IMF ("unchanging distribution regardless of environment and over the entirety of cosmic history") and concludes that general the IMF is not universal. This has consequences for the initial conditions of star cluster simulations across cosmic time and we need to model the IMF of different stellar populations individually. The IMF was established as a concept in a pioneering work by [456] as a quantisation of stellar masses in the Universe [457]. In general, the number of stars in the IMF is given by where dN is the number of stars formed in a small region, i.e. an embeddedcluster-forming molecular cloud core, in the mass interval m to m + dm [458]. Typically, we express the IMF as a (multi-)power law (powers are typically denoted by "α") depending on the stellar population that we want to model. For example, for Population I (Pop-I) stars we typically choose an IMF from [296,459], [460] or [461]; they are quite similar. The standard Nbody-codes, such as Nbody6++GPU or Nbody7 provide tools to initialize a star cluster model with generalized Salpeter or Kroupa [296] IMF's, in a mass range from 0.08 to 100 or 150 M . Note that also the initialization of lower mass objects has been prepared in the codes by Pavel Kroupa. Despite many observations over the last decades, the IMF for Pop-I stars is still quite uncertain, see [455] and sources therein. For Pop-III stars, the IMF is very different. It becomes increasingly top-heavy for decreasing metallicity [462,463,464,465,466,458,467]. However, we do not have observations of Pop-III stars (yet; although even with the JWST it will be a difficult or impossible undertaking [468]. On the other hand, [469] claim that some hundred SNe detections by JWST may be enough to constrain the IMF of Pop-III stars. See also [470] for a further discussion) and therefore, we do not have statistics from which to conclude an IMF. A flat IMF with α −1.0 between 8 M and 300 M for Pop-III stars has been proposed by [471]. However, [472] use a Salpeter IMF [456] of slope α −2.35 with a maximum mass of around 87 M instead. We will have to wait for observations of Pop-III stars or their remnants before we can reliably constrain their IMF.

Initial binary population
Almost all stars form in binary systems and some in higher order multiple systems [473,474,475]. As with the IMF in Sect. 6.4.3, there is some debate on the universality of an initial binary population (IBP). In other words, it is contentious that the IBP is independent of environment, in which binaries form [476,477]. This is typically quietly assumed in the initialisation of star cluster simulations, at least in simulations of clusters made up of Pop-I stars (e.g. [428,244]). Although, we would expect this to vary for decreasing metallicity and higher redshift, because the environments and also the primordial gas from which the stars form have very different properties from Pop-I stars (e.g. [478,479]). In any way, the IBP evolves on a cluster crossing timescale t cr . The widest binaries that form are dynamically disrupted, while in star clusters the hardest binaries harden further [25,480]. This leads to a pronounced SN-Ia rate in star clusters [481]. Binaries generally dominate the global, dynamical evolution of the star cluster by close Newtonian few-body interactions (binary-single and binary-binary encounters) [482,483]. A distinction is made following [484,474,477] between a birth stellar population, where all protostars are embedded in circum-protostellar material, and an initial stellar population, which consists of pre-main-sequence stars and which are not embedded in circum-protostellar material.
The process from an initial to a fully formed main-sequence binary star population is called pre-main-sequence eigenevolution. Eigenevolution is the sum of all dissipative physical processes that transfer mass, energy and angular momentum between the companions when they are still very young and accreting (sentence directly cited from [474]). Most obvious is the process of tidal circularization (e.g. [165]) of tight and initially highly eccentric binaries, which leads to a depletion of high eccentricities for small semi-major axes. In [474] there is a comprehensive description of eigenevolution and the relation between binary parameters at birth and after the pre-main sequence evolution, when typically N -body simulations start (see also further citations in [474] about observational binary data therein, and also [453,272,423,477]). We present this as an example here, but it is not clear whether under all astrophysical circumstances such universal binary parameters are realized. Dynamically, a binary star depends on four parameters: its system mass m sys = m 1 + m 2 , its period P and correspondingly its semi-major axis a (via Kepler's third law), its mass ratio q = m 2 /m 1 ≤ 1 and its eccentricity e=(r apo −r peri )/(r apo +r peri ) [474], where r apo and r peri are the apocentric and pericentric distances, respectively. Thus, a complete initial binary population in a star cluster depends on the stellar IMF ξ * (m) (see Sect. 6.4.3), the period distribution f P (logP ), the mass ratio distribution f q (q) and the eccentricity distribution f e (e) [474,485,477]: 1. f P (logP ): [484] show that where P min = 10 days, δ = 45, η = 2.5 and P max = 2.188 × 10 8 days, because the initial binary fractions f b is 100% [473]. Adjustments to this distribution were later made for high mass stars with m > 5 M following [486,487], where for these stars f P ∝ (log(P )) −0.55 with P min = 1.412 days and P max = 3.162 × 10 3 days.
2. f e : typically, we distribute the binaries thermally, meaning angular momenta are distributed equally f e = 2e [484], although this might greatly overpredict observed merger rates according to [488]. 3. f q : the binary stars with members below 5 M are distributed randomly and for masses above 5 M , the binary mass ratios are distributed uniformly (0.1 < q < 1.0) [489,490,486,491].
Many direct N -body simulations usually start at a time when the star cluster is gas-free, when eigenevolution has terminated. Also dynamical interactions of binaries with each other and with single stars cause a further dynamical evolution of binary properties during the cluster formation phase. The "Kroupa" [484] period distribution includes many wide binaries. They are dynamically ionized (disrupted) in a time scale very short compared to the relaxation time scale. Monte Carlo simulations (some of them start with binary fractions close to unity [492]) and direct N -body models, starting with lower binary fraction (but most of them tightly bound) converge quickly, as can be seen in [244]. Therefore one uses typically rather small binary fractions in initial models for long-term evolution of star clusters (order of 5% -20%; [48,493,494,495]. Also for the same reason simpler power laws for the initial distributions of logP , e, and q are used (see e.g. [274,496]). Note that in case of rotating star clusters the initial inclination of binaries with respect to the rotation axis of the cluster is a new possibly important parameter, which has not yet been examined. It is sometimes quoted that small binary fractions are chosen in direct N -body simulations, because tightly bound binaries are computationally expensive; this may not be a problem in the near future. Already the PeTar code [497,149] can use efficiently a very large number of binaries in the simulations (see Fig. 5). Also Nbody6++GPU potentially allows a good KS binary parallelization, which is work in progress. Since the treatment of interacting and relativistic binaries in PeTar is not equivalent to the one used in Nbody6++GPU, tests and comparisons of both codes with respect to binary evolution are ongoing.
In summary of what we discussed so far in Sect. 6.4.1 to Sect. 6.4.4, during the initialisation of N -body simulations we need to distribute the stars in 6D phase space, draw their masses from some IMF, and distribute primordial binaries according to IBP [474,272,498]. Furthermore, we generally need to make a decision if our star cluster is mass segregated at the beginning of the simulation [499,500,501], shows fractality [502] and if it is or is not in virial equilibrium. Below we highlight two more areas of active research when it comes to simulations of star clusters and their initialisation: multiple stellar populations in Sect. 6.4.5 and initial bulk rotation of the star cluster in Sect. 6.4.6, respectively.

Multiple stellar populations
Modern observational methods have made it possible to resolve multiple stellar populations (MSPs) in globular clusters, which can mostly be inferred from photometric diagrams such as colour-magnitude diagrams (CMDs) from multiband HST photometry (e.g. [503,504,505,475,506,507,508]), as well as from abundance variations in integrated light [509,510]. Nowadays, MSPs have been confirmed in around 70 Galactic and extragalactic clusters [511,512,513,514,515,516,517,507,508]. MSPs can have very diverse origins, see e.g. the review [518]. Recently, kinematic differences between multiple populations have been found in Galactic globular clusters [519]. Dynamical simulations of different populations can be useful to unravel the origin of MSP and the relation between initial states and currently observed dynamics. Our direct Nbody6++GPU code allows for the distinction between different populations by defining a corresponding label for each star; this has been used by [520] to constrain the dynamical origin of multiple populations in intermediate-age clusters in the Magellanic Clouds.
Another rarely used feature of Nbody6++GPU (but see [521]) is that it can start with individual population data for every star (age, metallicity, population index). Also MOCCA simulations have been published [522] hosting two generations of stars as above in tidally filling and underfilling star clusters.
They are able to reproduce the observed fractions and properties of second generation stars in the GCs of the Milky Way. For the time being there is no good way known to handle encounters and mixing of material from different populations in stellar collisions and binary mass overflow interactions in any of the codes.
We are still at the beginning of deciphering the complex picture of MSPs and their dynamical evolution using computational methods for collisional stellar systems. Nevertheless, the studies above are encouraging.

Rotation
The inclusion of initial bulk rotation in direct Nbody simulations of collisional stellar systems is still unusual (e.g. [83,48,428,365,366,367,185,182,183,184,244]), although it has been known for over a century that star clusters even today show significant imprints of rotation, which can, for example, be observed in deviations in the shapes of star clusters from sphericity [523,524,525,526,527,528,529,3,530,531,532,533,4]. Moreover, present-day detectors and data processing methods have made it possible to resolve the photometry and kinematics of individual stars all the way down to the cluster centre revealing that rotational kinematic features vary between multiple stellar populations [534,535,536,537,21,22,538,539,540,541,542,543,544,545,546]. Additionally, both observations and simulations support these results and find that star clusters show significant fractality [547,548], and internal rotation at birth in general [549,550]. Velocity anisotropy has been observed in star clusters with detected elongated structures [551,548], which might be induced by rotation. So, how do we initialise rotating collisional stellar systems? Sometimes it is assumed that there exists a kind of "Maxwell's demon" that simply switches Fig. 11: 3D scatter plots in showing the spatial distribution of the ZAMS of high mass stars (and the BHs they quickly form) from the simulations that rotate extremely quickly initially ((W 0 = 6.0, ω 0 = 1.8), i.e. rotating King models by [67]) and have stellar evolution (level C by [244]) at around 0.0 Myr, 3.68 Myr, 11.86 Myr and 100.00 Myr, respectively; The stars and compact objects are color-coded by their mass between 0.0 M and 150.0 M . The stars and BHs are also projected onto the three dimensional axes, which can be seen from the light-grey dots. We can clearly see the (rotating) triaxial structure, a bar of the BHs and their progenitor stars, at t=3.68 Myr and the spatial reconfiguration of the high mass objects (mostly BHs) to axisymmetric structures already after 11.86 Myr. At 100.00 Myr a practically spherical system of BHs remains that is much more concentrated than the system of their progenitor stars at 0.0 Myr. (Figure adapted from simulations by [46]). the direction of initial particle velocities to induce angular momentum to a N -body system and assuming the preservation of the spherical distribution function (e.g. [552,553,554]) and angular momentum in the process (e.g. [555,556]). This procedure is not physical. We instead need distribution functions that at least depend on two integrals of motion, such as the total energy E and and the total angular momentum in the z-direction L z (assuming that the system rotates around the z-axis initially and ignoring a third integral, which in some cases can be approximated by the total angular momentum of a star L 2 [73], because a third integral is generally not analytically known). Such rotating equilibrium models were developed by [557,67,74,85]. They can be considered as generalizations of standard King models [75], because their energy dependence is a lowered isothermal, and the additional term for the second independent variable is exp(−L z /L z0 ) (analogous to [558]). L z0 is a scaling constant; usually a dimensionless rotation parameter ω 0 is used (derived from L z0 , see e.g. [67]). The models are axisymmetric, with a rigid rotation of the inner parts of the cluster, a maximum of the rotation curve close to the half-mass radius, and a differentially decreasing rotation curve outside in the halo. Rotation supports only a fraction of the total kinetic energy (see Table 1 in [67], and note that the 2nd column is erroneously labeled, it contains the percentage contained in rotational energy; i.e. for ω 0 = 0.6 we have 20% of the total kinetic energy in form of ordered rotational motion). Evolved star clusters obtained from these initial models agree quite well with observed clusters [559]. Due to the 2D velocity distribution function an anisotropy is possible between the velocity dispersions in radial direction (in cylindrical coordinates) and in rotational ϕ− direction; the models are isotropic between radial and vertical direction (parallel to the rotation axis). Such models have been used as initial models in 2D Fokker-Planck (FP) modelling and direct N -body simulation [76,77,78,559,560,561,562,83,46]. Note that the models by [85] are using a different form of the distribution function based on the Jacobian of a cluster rotating around the galaxy, but they are as well generalizations of King models for rotation with similar properties. They also have been used as initial models for direct Nbody models [563,564,84,565,545,546,80,81]. Furthermore, semi-analytic models exist that [42,43,44,566] used to study the formation and evolution of rotating stellar or black hole disks in nuclear star clusters. Most of the aforementioned studies find evidence for the gravogyro catastrophe and its coupling to the gravothermal catastrophe [567,17,568,79].
It is important to note that both of these effects are entirely gravitational and disappear in the absence of gravity [569]. The gravogyro catastrophe happens on all astrophysical scales from rotating stars, to the formation of the Solar system and the dynamics of spiral galaxies. A recent study on the impact of stellar evolution on rotating star clusters, finds that the initial bulk rotation leads to a rotating triaxial structure, a bar, of black holes and their progenitor stars in early star cluster evolution that then takes the shape of an axisymmetric structure, a disk, over time [46]. For one of the simulations from [46] the results are shown in Fig. 11. This bar formation and dissolution was already found by [79] in their pioneer low particle number N -body simulations (N = 1024) of rotating star clusters. However, the simulation set-up was very different and naturally less advanced. Since rotation is fundamental to star cluster dynamics and since these environments are typically very dense [1], it will be important to understand how this process can influence the growth of intermediate mass black holes (see Sect. 6.5) and their progenitor stars going forward. Recently, a study on extremely massive and rotating star Pop-III star clusters does indicate a trend that the larger the initial bulk rotation, the more stellar and BH-BH mergers there are [82]. This might be due to the extremely efficient angular momentum transport in early star cluster evolution (0-2 Myr) and the resulting contraction of of the central region before stellar evolution begins to dominate the pre-core collapse evolution. However, more simulations are needed here to reach a conclusion. Lastly, we note that due to the assumption of spherical symmetry most of the current Monte Carlo methods for collisional dynamics are currently unable to evolve initially rotation star cluster models, such as those described above (e.g. [570,68,96,97,98,101,571,428,572,110]). A restricted Monte Carlo method for rotating, axisymmetric star clusters (usable even for general geometry, but see problem below) has been presented by [112]. It uses Spitzer's Monte Carlo method (see Sect. 4.1), which distinguishes it from currently common Monte Carlo codes. But it has some more serious approximations, because the random relaxation scatterings are applied only in v and v ⊥ obtained from a fully isotropic spherically symmetric background. In 2D Fokker-Planck models of axisymmetric rotating star clusters [67] the background distribution function used for the computation of diffusion coefficients is fully self-consistent and there are five different diffusion coefficients obtained (instead of only two).

Formation of massive objects
The formation of massive objects from collisions of massive stars and mergers involving black holes is a subject with long history and in full detail beyond the scope of this review. The first detection of a massive black hole merger from LIGO-Virgo (coalescence of two black holes with 85 and 66 M [573] has led to new interest into dynamical paths of formation of massive objects (intermediate mass black holes, IMBH) in star clusters. Numerical and theoretical works on IMBH formation in dense star clusters suggests that IMBH formation may occur through differet channels (multiple stellar mergers, accretion of stellar matter onto a stellar mass black hole, or multiple generations of relativistic black hole coalesciences (and mixtures of these processes). The idea that this may happen is not new (cf. e.g. [574,575,576,577]). On the observational side IMBHs remain elusive, a recent review about the issue can be found in [578].
Here we just want to highlight a few recent results obtained with N -body and MOCCA codes, showing the feasibility to grow an IMBH within a suitably initialized star cluster (e.g. [579,184,580,581,582,188,293,185,182,183,82,494]). Star cluster environments (particularly young massive or Pop-III clusters) can precipitate the birth and growth of IMBHs in two of the main proposed mechanisms as listed in [578]: firstly, stellar mergers in star cluster evolution and subsequent core collapse into an IMBH of the progenitor star or secondly, gravitational runaway mergers of (IM)BHs to form (even more) massive IMBHs. The mergers of (IM)BHs are associated with the emission of a GWs (see Sect. 6.2.7 and sources therein) and GR merger recoil kicks (see Sect. 6.2.5 and sources therein). The emitted gravitational radiation can be detected [583,580]. (a)LIGO [27,28,29], (a)Virgo [30,28,29] and KAGRA (e.g. [28,31,32]) currently detect stellar mass black hole inspirals and mergers during their last seconds or fractions of seconds. With improved ground based detectors (Ein-stein Telescope [33] or Cosmic Explorer [34] or space-based instruments (LISA or its Chinese counterparts, see below) inspiralling stellar mass black hole binaries could be detected months and years before their final coalescence (cf. e.g. [584]), because they are sensitive at lower frequencies. Planned space based gravitational wave detectors are the European-American LISA [585,586] project and the Chinese projects TianQin [587] and Taiji [588,589]. Also these instruments will be more suitable for detection of other black holes inspiralling and merging with IMBHs. Here we present an example from our N -body simulations, which shows that from them one can very well predict GW signal patterns and the abundances of GW events related to IMBHs (in the lower mass regime as discussed above) originating from star clusters (e.g. [590,591]). Fig. 12 shows data from direct N -body simulations of [183]; many BH-BH merger events found in the models fit in almost perfectly with the LIGO-Virgo-KAGRA detections. In these models a new parameter f c is used, describing how much mass is lost in the collision of a black hole with a main sequence star, f c = 0 means no black hole growth in the process, f c = 1 means that all the mass of the main sequence star is added to the black hole; this parameter is now present in current Nbody6++GPU, and was introduced by [183]. Also already by [446,594] it was noticed that in such a system stars would be tidally disrupted near the SMBH, leading to its growth by mass accretion. Tidally disrupted stars would be preferentially on elongated radial orbits, with low angular momentum. The accretion process would take these stars out of the stellar distribution and create a "loss cone" in phase space. In a spherically symmetric cluster the loss cone would be refilled by diffusion of angular momenta due to two-body relaxation [595,596]. In axisymmetric or triaxial nuclear star clusters angular momentum diffusion can be much faster and lead to higher star accretion rates [597,598,599,600,601]. Also orbit-averaged FP models have been used for this problem [602,445,603,604,605]. More recently a first direct N -body model using Nbody6++GPU has been published [606] (DRAGON simulation of the Galactic Center -it led to the conclusion that the central cusp areas near the SMBH are dominated by stellar mass black holes, which will be accreted to the SMBH under gravitational wave emission.

Tidal disruption events
If a star is tidally disrupted by the SMBH on a parabolic orbit only half of its mass should be accreted while the other half has enough energy to escape from the SMBH (Rees' conjecture, [607]). The characteristic light curve of such events decays with t −3/5 ; under realistic conditions, however, tidally disrupted stars will arrive on bound (eccentric) and unbound (hyperbolic) orbits, not just parabolic ones; and that will cause deviations from the standard light curve as well as change the accreted mass fraction. In [608,609] large direct N -body simulations were done to find out the distribution of accreted stars. Furthermore a partial tidal disruption of giant stars was implemented, where the envelope is removed early and a core similar to a white dwarf survives and gets tidally disrupted later [609]. Detailed stellar models were also used to determine the mass fraction of a star which is accreted to the SMBH after tidal disruption [610]. In future work this will be used in conjunction with the detailed mass function, single and binary stellar evolution as in the globular cluster simulations.

Practical tools
Here we want to highlight some software tools that have been developed and successfully applied in the context of collisional dynamics.

Multiscale and multiphysics simulation with AMUSE
The Astrophysical Multipurpose Software Environment (AMUSE) aims to provide a framework by which to simulate multiscale and multiphysics in a hierarchical fashion [611,612,613,614]. It does so by constructing new applications from the combination of known codes (solvers). For example, for gravitational dynamics the user may choose from a selection of 18 N -body codes [613]. The AMUSE framework successfully connects gravitational dynamics, radiative transport, stellar evolution and hydrodynamics [614], which results in incredibly diverse research: gas expulsion in early star cluster evolution [273,615], star cluster collisions and massive star cluster formation [616], the interaction of binary stars in gaseous filaments [617], the evolution of star clusters in a cosmological tidal field [618], evolution of triple stellar systems with Roche-lobe filling binaries [619] and many more. Particle data and relevant quantities are shared between the constituent codes through the AMUSE framework and unit conversion and other data manipulation can be done in this process [613]. AMUSE is capable of taking care of all sorts of technical problems, such as communication between codes, boundary conditions etc.. However, AMUSE is naturally limited by algorithmic complexity of the software it combines and the astrophysical research objective itself. For example, if we wanted to simulate a globular cluster of realistic size with Nbody6++ [620,118,621,137,157,134,158,622,122,204,127,623,137,145] and evolve the stars in the cluster with MESA [230,231,232,233,234,235], it would take an enormous amount of time. This is not a problem of AMUSE, but rather a statement on the plausibility of certain simulations themselves. However, we note that the operations by AMUSE and communication between codes also costs computing time.
8.2 Simulation data processing and analysis AMUSE (introduced in Sect. 8.1) provides a large number of in-built data analysis tools simulation data, such as the HOP [624] or the Kepler packages [613]. These packages can be used by AMUSE users to comfortably analyse data from their simulations without having to write their own scripts. Many N -body codes are now supported by built-in simulation data postprocessing analysis packages mostly written in python. The direct N -body code PeTar [497,149] also provides built-in data analysis tools, such as movie generators from particle data and HR diagrams. Likewise, the CMC code [99,625,341,108] uses the cmctoolkit package [107] for converting the simulation output into, e.g., velocity dispersion and surface brightness profiles. Another extremely useful tool for distributed data analysis is provided by BEANS [626], which can in principle be used for simulation output data format (e.g. csv, HDF5, or FITS [626,522]). Therefore, this tool can be used for the data analysis of output from all previously mentioned codes. It uses the industry-standard, the Apache Hadoop platform [627], for data analysis. This platform is highly optimized for processing huge amounts of (continuously generated) data and is therefore ideal not only in the "Internet of Things" and the live communication between machines, but also on-the-fly data processing and analysis from N -body simulations. BEANS is currently heavily used in data processing and analysis of MOCCA simulations (e.g. [522]). Recently, BEANS has also been configured to easily process the output data from Nbody6++GPU simulations.

Photometric mock observations from star cluster simulations
It is oftentimes desirable and useful to create photometric mock observations (CMDs or spectra of star clusters) from star cluster simulations. For this purpose various codes have been written. They mostly rely on the same principle: stellar masses, temperatures, luminosities and metallicities from the N -body output data are converted into observational magnitudes in various filter systems (e.g. HST, SDSS, 2MASS and it is typically easy to systems of new observing instruments) and into spectra in a certain wavelength range (e.g. [628]). Therefore, these mock observations are heavily influenced by the stellar evolution models that the respective simulations evolve the stars with. In the following paragraph, we summarise a number of codes for the creation of the observations outlined above. The Cluster simulatiOn Comparison with ObservAtions (COCOA) code by [629,630] has been used extensively in the creation of mock observations from N -body and Monte Carlo simulations (e.g. [428,48,631,477]). COCOA has also been extensively used in MOCCA studies concerning the evolution and detectability of IMBHs in star clusters [632,633,634]. In Fig. 13 a CMD constructed with COCOA is shown taken from [48]. They then compared this CMD with HST data of NGC 4372 (the target cluster of this study) from [635] and found that the CMDs are mostly very similar. Differences between actual and mock observations were attributed to the intrinsic observational photometry error or the presence of MSPs (see Sect. 6.4.5) in NGC 4372 that were not resolved by the simulation. Inaccurate stellar evolution modelling will have further influenced the results (in [48], the so-called level A was used, see [244]). This is a practical example of how photometric mock observations and the codes that produce these can be extremely useful in constraining both obser-vation and theory. The GAlaxy EVolutionary synthesis models (Galev) [636] was adapted for the use creating mock observations from N -body simulations (mostly intended for Nbody6++GPU and codes from that family, but it can be used for any other code that provides stellar masses, temperatures, luminosities and metallicities of the stars) by [628]. The resulting code GalevNB was applied to study the dynamical origin of MSPs (see Sect. 6.4.5) [637], the dynamical evolution of planetary systems in star clusters [638], the long-term evolution of binaries in the Dragon simulations [48,639] and PeTar simulations of open clusters and exploring their UV-excess [640]. More codes that are worth mentioning and fulfil the same underlying purpose as COCOA and GalevNB are: Flexible Stellar Population Synthesis (FSPS) [641,642,643], Simulating IFU Star Cluster Observations (SISCO) for integral field unit (IFU) observations of globular clusters [644], Make Your Own Synthetic ObservaTIonS (MYOSOTIS) for creating generic photometric observations [645] and Massive Cluster Evolution and Analysis Package (MASSCLEAN), which serves the same purposes [646].

Summary and Conclusion
The gravitational N -body problem remains to be one of the oldest and most exciting problems in astrophysics. In this review, we have focused on a branch of this incredibly diverse field: the computation of collisional stellar N -body systems. The term collisional refers to the cumulative effect of distant elastic two-body encounters here. In nature these systems are commonly represented young and massive, globular, or nuclear star clusters. These clusters are also dense and gravothermal -stellar density are high enough that direct collisions and mergers between stars occur in certain phases, and gravothermal in the sense that two-body relaxation, through distant elastic encounters, provides heat transfer and viscous angular momentum transfer.
Despite exhaustive efforts to accurately model these systems, there is still no self-consistent theory available. This is partly because N -body systems are dynamically chaotic for a significant fraction of the phase space of initial conditions if N > 2. For small N ≤ 3, there do exist some stable solutions [647,648], but these special configurations are extremely unlikely to occur in nature [649]. For all N > 3, there exist no analytical solutions and we must solve the equation of motions numerically, thus with the help of computers. For increasing N , certain configurations of such systems exhibit truly remarkable physical properties. They might exhibit negative heat capacities that result in gravothermal contraction (sometimes also coined the "gravothermal catastrophe" (e.g. [15]). If the star cluster rotates, it could exhibit a negative moment of inertia leading angular momentum transport from high mass particles to low mass particles in the system, while at the same time increasing the rotational velocity of the high mass particles, which is known as the gravogyro catastrophe (e.g. [17,78]). In a realistic star cluster simulations, both of these processes will happen at different timescales (convective angular momentum transport is more efficient than conductive heat transport in star clusters), but they are coupled and ultimately reinforce one another [568]. Importantly, these processes are entirely gravitational in nature and disappear in the absence of self-gravity. The two astrophysical methods that have been most successfully tackling collisional gravitational N -body systems are either direct N -body simulations (using for gravitational force calculations at least at some times the summation over all N − 1 other stars) or approximate methods, based on statistical mechanics, using the Fokker-Planck modelling with Monte Carlo methods. While the former method appears physically more accurate, less affected by approximations, and resolves any astrophysical object of interest (star clusters) better, the latter is much faster and computationally less expensive. And recent star-by-star Monte Carlo modelling delivers detailed data comparable to direct N -body simulations. It has been demonstrated repeatedly that for global quantities such as the time evolution of binary fractions both methods yield very similar results. Therefore, these methods ultimately complement each other and many works have employed both side by side. All fields in astrophysics have benefited from the revolutionary invention of computers and the ingenious hardware architectures that have been developed since. By the same token, new programming languages and parallelization techniques have had a significant impact on accelerating solving astrophysical problems on hardware and made certain problems possible to begin with. The collisional gravitational N -body problem is no exception here. When the astrophysics, the hardware and software domains are mastered, research in this field can culminate in ground-breaking results: the first million-body simulations of globular clusters across cosmic time, the Dragon simulations performed with Nbody6++GPU, are a product of this but more direct N -body models of Galactic and extragalactic globular and nuclear star clusters are needed to compare with both observational data and the wealth of data from Monte Carlo surveys. It follows that nowadays the gravitational N -body problem is in as much an astrophysics as it is a computer science problem. Examples here are the advent of novel simulation software such as the direct N -body codes PeTar and Frost, and new hybrid code approaches such as ETICS, which will make it possible to break into the 10 6 − 10 9 particle domain. Massive globular clusters such as 47 Tuc or ω Cen as well as most nuclear star clusters belong to this parameter range. Will a traditional workhorse such as Nbody6++GPU again keep up and compete here? The first Exaflop/s computer has been inaugurated in the US this year 5 -for the next steps we need to again adopt our codes to the new developing hardware. It is now 32 years ago that D. Sugimoto announced the million body problem as a challenge for direct N -body simulations -here we want to raise the question of when the billion body problem can be tackled and what codes and what hardware do we need? Another challenge lies in the integration of astrophysics into gravitational Nbody simulations. This involves classical single and binary stellar evolution as well as relativistic dynamics of compact objects, or proper modelling of external tidal fields. The central densities of collisional stellar systems are typically so high, that stars frequently interact with each other; they exchange mass and angular momentum or even collide and merge. In other words, the stars in the simulations do not only evolve on their own, they also interact with nearby stars (and compact objects) through diverse astrophysical processes, such as tidal interaction, dynamical mass transfer, gravitational wave in-spiral and many more. For this purpose there are two main methods that have been thoroughly tested: interpolation between tables or approximation of stellar evolution data by some interpolation (fitting) formulae as functions of mass, age and metallicity. These two approaches are also not in competition, but rather complement one another. The former has traditionally been limited by memory on hardware, since the tables are necessarily rather large. Due to advances in hardware, this problem is diminishing. Unlike the latter approach, fitting formulae, stellar parameters from the given set of detailed tracks are calculated in real time with this method. Therefore, this approach is probably the most flexible, robust and efficient today when combining detailed stellar evolution with stellar dynamics. Furthermore, fitting formulae typically take much more care and time to set up. But the effort is rewarded with rapid, robust and analytic formulae, which can be easily modified and integrated into an N -body code and that quickly give stellar luminosity, radius and core mass of the stars as functions of mass, metallicity and age for all stellar evolutionary phases. When full stellar evolution is combined with collisional dynamics and the resulting code is used to study the dynamical evolution of dense star clusters of various make-ups, exciting research can be done. From the formation and evolution of exotic stellar and compact binaries, such as blue stragglers, Cataclysmic Variables, Algols, X-ray binaries, and many others to doubledegenerate binaries, such as BH-BH and the elusive BH-NS binaries, to the abundances of compact objects and their dynamical properties an extremely diverse array of astrophysically fascinating populations can be modelled for large metallicity and mass ranges. Therefore, the simulations of collisional stellar systems with modern production codes are the ideal laboratory to study stellar evolution in dense stellar environments. From the gravitational wave in-spiral of compact objects that can be modelled nowadays up to order PN(3.5) as well as the associated general relativistic merger recoil kick with the use of fitting formulae to numerical relativity, the full compact object merger phase can be modelled. As a result, simulations can yield predictions for theoreticians and observers alike on the properties, abundances and dynamics of gravitational wave sources from star clusters (and field) across cosmic time. Additionally, with the use of photometric mock observations from star cluster simulations, simulation data can be translated into magnitudes and fluxes for most mainstream filter systems, which has been done extensively already. In summary, simulations of collisional stellar systems are important and useful tools for unravelling our cosmic history in the age of multi-messenger astronomy. They are at the crossroads of many seemingly disparate astrophysical research fields, much like the simulation target, star clusters, are a fundamental building block in a hierarchy of cosmological structure formation [1].