Observed hyperbolic minor bodies might have an interstellar origin, but they can be natives of the Solar system as well. Fly-bys with the known planets or the Sun may result in the hyperbolic ejection of an originally bound minor body; in addition, members of the Oort cloud could be forced to follow inbound hyperbolic paths as a result of secular perturbations induced by the Galactic disc or, less frequently, due to impulsive interactions with passing stars. These four processes must leave distinctive signatures in the distribution of radiants of observed hyperbolic objects, both in terms of coordinates and velocity. Here, we perform a systematic numerical exploration of the past orbital evolution of known hyperbolic minor bodies using a full N-body approach and statistical analyses to study their radiants. Our results confirm the theoretical expectations that strong anisotropies are present in the data. We also identify a statistically significant overdensity of high-speed radiants towards the constellation of Gemini that could be due to the closest and most recent known fly-by of a star to the Solar system, that of the so-called Scholz's star. In addition to and besides 1I/2017 U1 (`Oumuamua), we single out eight candidate interstellar comets based on their radiants' velocities.

Asteroids that follow similar orbits may have a dynamical connection as their current paths could be the result of a past interaction with a massive perturber. The pair of extreme trans-Neptunian objects or ETNOs (474640) 2004 VN112-2013 RF98 exhibits peculiar relative orbital properties, including a difference in longitude of the ascending node of just 1.61 degrees and 3.99 degrees in inclination. In addition, their reflectance spectra are similar in the visible portion of the spectrum. The origin of these similarities remains unclear. Neglecting observational bias, viable scenarios that could explain this level of coincidence include fragmentation and binary dissociation. Here, we present results of extensive direct N-body simulations of close encounters between wide binary ETNOs and one trans-Plutonian planet. We find that wide binary ETNOs can dissociate during such interactions and the relative orbital properties of the resulting unbound couples match reasonably well those of several pairs of known ETNOs, including 474640-2013 RF98. The possible presence of former binaries among the known ETNOs has strong implications for the interpretation of the observed anisotropies in the distributions of the directions of their orbital poles and perihelia.

Making use of a new N-body model to describe the evolution of a moderate-size globular cluster we investigate the characteristics of the population of black holes within such a cluster. This model reaches core-collapse and achieves a peak central density typical of the dense globular clusters of the Milky Way. Within this high-density environment we see direct confirmation of the merging of two stellar remnant black-holes in a dynamically-formed binary, a gravitational wave source. We describe how the formation, evolution and ultimate ejection/destruction of binary systems containing black holes impacts the evolution of the cluster core. Also, through comparison with previous models of lower density, we show that the period distribution of black hole binaries formed through dynamical interactions in this high-density model favours the production of gravitational wave sources. We confirm that the number of black holes remaining in a star cluster at late times and the characteristics of the binary black hole population depend on the nature of the star cluster, critically on the number density of stars and by extension the relaxation timescale.

The Planet Nine hypothesis has now enough constraints to deserve further attention in the form of detailed numerical experiments. The results of such studies can help us improve our understanding of the dynamical effects of such a hypothetical object on the extreme trans-Neptunian objects or ETNOs and perhaps provide additional constraints on the orbit of Planet Nine itself. Here, we present the results of direct N-body calculations including the latest data available on the Planet Nine conjecture. The present-day orbits of the six ETNOs originally linked to the hypothesis are evolved backwards in time and into the future under some plausible incarnations of the hypothesis to investigate if the values of several orbital elements, including the argument of perihelion, remain confined to relatively narrow ranges. We find that a nominal Planet Nine can keep the orbits of (90377) Sedna and 2012 VP113 relatively well confined in orbital parameter space for hundreds of Myr, but it may make the orbits of 2004 VN112, 2007 TG422 and 2013 RF98 very unstable on time-scales of dozens of Myr, turning them retrograde and eventually triggering their ejection from the Solar system. Far more stable orbital evolution is found with slightly modified orbits for Planet Nine.

Introducing the DRAGON simulation project, we present direct $N$-body simulations of four massive globular clusters (GCs) with $10^6$ stars and 5$\%$ primordial binaries at a high level of accuracy and realism. The GC evolution is computed with NBODY6++GPU and follows the dynamical and stellar evolution of individual stars and binaries, kicks of neutron stars and black holes, and the effect of a tidal field. We investigate the evolution of the luminous (stellar) and dark (faint stars and stellar remnants) GC components and create mock observations of the simulations (i.e. photometry, color-magnitude diagrams, surface brightness and velocity dispersion profiles). By connecting internal processes to observable features we highlight the formation of a long-lived 'dark' nuclear subsystem made of black holes (BHs), which results in a two-component structure. The inner core is dominated by the BH subsystem and experiences a core collapse phase within the first Gyr. It can be detected in the stellar (luminous) line-of-sight velocity dispersion profiles. The outer extended core - commonly observed in the (luminous) surface brightness profiles - shows no collapse features and is continuously expanding. We demonstrate how a King (1966) model fit to observed clusters might help identify the presence of post core-collapse BH subsystems. For global observables like core and half-mass radii the direct simulations agree well with Monte-Carlo models. Variations in the initial mass function can result in significantly different GC properties (e.g. density distributions) driven by varying amounts of early mass loss and the number of forming BHs.

On 2013 February 15 a small asteroid rammed against the atmosphere above the region of Chelyabinsk in Russia, producing the most powerful superbolide since the Tunguska event in 1908. Lacking proper astrometric observations, the pre-impact orbit of this object has been determined using videos, satellite images, and pure geometry. Unfortunately, more than two years after the event, the published estimates vary so much that there is no clear orbital solution that could be used to investigate the origin of the impactor and the existence of dynamically, or perhaps even genetically, related asteroids. Here, we revisit this topic using a full N-body approach. A robust statistical test is applied to published solutions to discard those unable to produce a virtual impact at the observed time (03:20:20.8 s UTC). The same N-body methodology and the latest ephemerides are used to compute a new orbital solution: a=1.6247 AU, e=0.5318, i=3.9750 degrees, Omega=326.4607 degrees and omega=109.7012 degrees. This new solution --which has an impact probability > 0.99999 and uncertainties in time and space of 0.2 s and 6 km, respectively-- is utilized to explore the past orbital evolution of the impactor as well as the presence of near-Earth objects moving in similar paths. A dynamical link between asteroid 2011 EO40 and the Chelyabinsk impactor is confirmed. Alternative orbital solutions are extensively explored.

We present a hybrid code combining the OpenMP-parallel tree code VINE with an algorithmic chain regularization scheme. The new code, called "rVINE", aims to significantly improve the accuracy of close encounters of massive bodies with supermassive black holes in galaxy-scale numerical simulations. We demonstrate the capabilities of the code by studying two test problems, the sinking of a single massive black hole to the centre of a gas-free galaxy due to dynamical friction and the hardening of a supermassive black hole binary due to close stellar encounters. We show that results obtained with rVINE compare well with NBODY7 for problems with particle numbers that can be simulated with NBODY7. In particular, in both NBODY7 and rVINE we find a clear N-dependence of the binary hardening rate, a low binary eccentricity and moderate eccentricity evolution, as well as the conversion of the galaxy's inner density profile from a cusp to a a core via the ejection of stars at high velocity. The much larger number of particles that can be handled by rVINE will open up exciting opportunities to model stellar dynamics close to SMBHs much more accurately in a realistic galactic context. This will help to remedy the inherent limitations of commonly used tree solvers to follow the correct dynamical evolution of black holes in galaxy scale simulations.

Accurate direct $N$-body simulations help to obtain detailed information about the dynamical evolution of star clusters. They also enable comparisons with analytical models and Fokker-Planck or Monte-Carlo methods. NBODY6 is a well-known direct $N$-body code for star clusters, and NBODY6++ is the extended version designed for large particle number simulations by supercomputers. We present NBODY6++GPU, an optimized version of NBODY6++ with hybrid parallelization methods (MPI, GPU, OpenMP, and AVX/SSE) to accelerate large direct $N$-body simulations, and in particular to solve the million-body problem. We discuss the new features of the NBODY6++GPU code, benchmarks, as well as the first results from a simulation of a realistic globular cluster initially containing a million particles. For million-body simulations, NBODY6++GPU is $400-2000$ times faster than NBODY6 with 320 CPU cores and 32 NVIDIA K20X GPUs. With this computing cluster specification, the simulations of million-body globular clusters including $5\%$ primordial binaries require about an hour per half-mass crossing time.

Short summary of the contribution of Michel Henon to the subject of stellar dynamics.

Nearly all known extreme trans-Neptunian objects (ETNOs) have argument of perihelion close to 0 degrees. An existing observational bias strongly favours the detection of ETNOs with arguments of perihelion close to 0 degrees and 180 degrees yet no objects have been found at 180 degrees. No plausible explanation has been offered so far to account for this unusual pattern. Here, we study the dynamical evolution of comet 96P/Machholz 1, a bizarre near-Earth object (NEO) that may provide the key to explain the puzzling clustering of orbits around argument of perihelion close to 0 degrees recently found for the population of ETNOs. Comet 96P/Machholz 1 is currently locked in a Kozai resonance with Jupiter such that the value of its argument of perihelion is always close to 0 degrees at its shortest possible perihelion (highest eccentricity and lowest inclination) and about 180 degrees near its shortest aphelion (longest perihelion distance, lowest eccentricity and highest inclination). If this object is a dynamical analogue (albeit limited) of the known ETNOs, this implies that massive perturbers must keep them confined in orbital parameter space. Besides, its future dynamical evolution displays orbital flips when its eccentricity is excited to a high value and its orbit turns over by nearly 180 degrees, rolling over its major axis. This unusual behaviour, that is preserved when post-Newtonian terms are included in the numerical integrations, may also help understand the production of NEOs on retrograde orbits.

Nearly $15-20%$ of solar type stars contain one or more gas giant planet. According to the core-accretion scenario, the acquisition of their gaseous envelope must be preceded by the formation of super-critical cores with masses ten times or larger than that of the Earth. It is natural to link the formation probability of gas giant planets with the supply of gas and solid in their natal disks. However, a much richer population of super Earths suggests that 1) there is no shortage of planetary building-block material, 2) gas giants' growth barrier is probably associated with whether they can merge into super-critical cores, and 3) super Earths are probably failed cores which did not attain sufficient mass to initiate efficient accretion of gas before it is severely depleted. Here we construct a model based on the hypothesis that protoplanetary embryos migrated extensively before they were assembled into bona fide planets. We construct a Hermite-Embryo code based on a unified viscous-irradiation disk model and a prescription for the embryo-disk tidal interaction. This code is used to simulate 1) the convergent migration of embryos, and 2) their close encounters and coagulation. Around the progenitors of solar-type stars, the progenitor super-critical-mass cores of gas giant planets primarily form in protostellar disks with relatively high ($\gtrsim 10^{-7} M_\odot$ yr$^{-1}$) mass accretion rates whereas systems of super Earths (failed cores) are more likely to emerge out of natal disks with modest mass accretion rates, due to the mean motion resonance barrier and retention efficiency.

We provide a set of analytic fits to the radii of pre-mainsequence stars in the mass range 0.1 < M/Msun < 8.0. We incorporate the formulae in N-body cluster models for evolution from the beginning of pre-main sequence. In models with 1,000 stars and high initial cluster densities, pre-mainsequence evolution causes roughly twice the number of collisions between stars than in similar models with evolution begun only from the zero-age main sequence. The collisions are often all part of a runaway sequence that creates one relatively massive star.

We describe the use of Graphics Processing Units (GPUs) for speeding up the code NBODY6 which is widely used for direct $N$-body simulations. Over the years, the $N^2$ nature of the direct force calculation has proved a barrier for extending the particle number. Following an early introduction of force polynomials and individual time-steps, the calculation cost was first reduced by the introduction of a neighbour scheme. After a decade of GRAPE computers which speeded up the force calculation further, we are now in the era of GPUs where relatively small hardware systems are highly cost-effective. A significant gain in efficiency is achieved by employing the GPU to obtain the so-called regular force which typically involves some 99 percent of the particles, while the remaining local forces are evaluated on the host. However, the latter operation is performed up to 20 times more frequently and may still account for a significant cost. This effort is reduced by parallel SSE/AVX procedures where each interaction term is calculated using mainly single precision. We also discuss further strategies connected with coordinate and velocity prediction required by the integration scheme. This leaves hard binaries and multiple close encounters which are treated by several regularization methods. The present nbody6-GPU code is well balanced for simulations in the particle range $10^4-2 \times 10^5$ for a dual GPU system attached to a standard PC.

We report on results of fully consistent N-body simulations of globular cluster models with N = 100 000 members containing neutron stars and black holes. Using the improved `algorithmic regularization' method of Hellstrom and Mikkola for compact subsystems, the new code NBODY7 enables for the first time general relativistic coalescence to be achieved for post-Newtonian terms and realistic parameters. Following an early stage of mass segregation, a few black holes form a small dense core which usually leads to the formation of one dominant binary. The subsequent evolution by dynamical shrinkage involves the competing processes of ejection and mergers by radiation energy loss. Unless the binary is ejected, long-lived triple systems often exhibit Kozai cycles with extremely high inner eccentricity (e > 0.999) which may terminate in coalescence at a few Schwarzschild radii. A characteristic feature is that ordinary stars as well as black holes and even BH binaries are ejected with high velocities. On the basis of the models studied so far, the results suggest a limited growth of a few remaining stellar mass black holes in globular clusters.

We briefly report on a crater on the western slope of Mt.Ararat . It is located in an area closed to foreigners at an altitude around 2100m with geographic coordinates 39\deg 47' 30"N, 44\deg 14' 40"E. The diameter of the crater is around 60-70m, the depth is up to 15m. The origin of the crater, either of meteorite impact or volcanic, including the evaluation of its age, will need detailed studies.

We describe efforts over the last six years to implement regularization methods suitable for studying one or more interacting black holes by direct N-body simulations. Three different methods have been adapted to large-N systems: (i) Time-Transformed Leapfrog, (ii) Wheel-Spoke, and (iii) Algorithmic Regularization. These methods have been tried out with some success on GRAPE-type computers. Special emphasis has also been devoted to including post-Newtonian terms, with application to moderately massive black holes in stellar clusters. Some examples of simulations leading to coalescence by gravitational radiation will be presented to illustrate the practical usefulness of such methods.

We investigate the evolution of binary fractions in star clusters using N-body models of up to 100000 stars. Primordial binary frequencies in these models range from 5% to 50%. Simulations are performed with the NBODY4 code and include a full mass spectrum of stars, stellar evolution, binary evolution and the tidal field of the Galaxy. We find that the overall binary fraction of a cluster almost always remains close to the primordial value, except at late times when a cluster is near dissolution. A critical exception occurs in the central regions where we observe a marked increase in binary fraction with time -- a simulation starting with 100000 stars and 5% binaries reached a core binary frequency as high as 40% at the end of the core-collapse phase (occurring at 16 Gyr with ~20000 stars remaining). Binaries are destroyed in the core by a variety of processes as a cluster evolves, but the combination of mass-segregation and creation of new binaries in exchange interactions produces the observed increase in relative number. We also find that binaries are cycled into and out of cluster cores in a manner that is analogous to convection in stars. For models of 100000 stars we show that the evolution of the core-radius up to the end of the initial phase of core-collapse is not affected by the exact value of the primordial binary frequency (for frequencies of 10% or less). We discuss the ramifications of our results for the likely primordial binary content of globular clusters.

We report on the first fully consistent conventional cluster simulation which includes terms up to post^{5/2} Newtonian in the potential of the massive body. Numerical problems for treating extremely energetic binaries orbiting a single massive object are circumvented by employing the special ``wheel-spoke'' regularization method of Zare (1974) which has not been used in large-N simulations before. Idealized models containing N = 10^5 particles of mass 1 M_sun with a central black hole of 300 M_sun have been studied on GRAPE-type computers. An initial half-mass radius of r_h = 0.1 pc is sufficiently small to yield examples of relativistic coalescence. This is achieved by significant binary shrinkage within a density cusp environment, followed by the generation of extremely high eccentricities which are induced by Kozai (1962) cycles and/or resonant relaxation. More realistic models with white dwarfs and ten times larger half-mass radii also show evidence of GR effects before disruption. Experimentation with the post-Newtonian terms suggests that reducing the time-scales for activating the different orders progressively may be justified for obtaining qualitatively correct solutions without aiming for precise predictions of the final gravitational radiation wave form. The results obtained suggest that the standard loss-cone arguments underestimate the swallowing rate in globular clusters containing a central black hole.

(Abridged) In this celebratory contribution, we present some new data from standard star cluster modelling containing primordial binaries and triples, as well as results from a binary black hole simulation with two massive members. In the star cluster case, the process of mass loss from evolving stars, together with general mass segregation, promotes favourable interactions involving compact subsystems of binaries and triples in the central region. Three-body interactions often lead to energetic ejections, with one or more of the components attaining relatively large terminal velocities which can be observed in principle. The second type of stellar system to be discussed is based on the scenario of two approaching galactic cores with density cusps, each having a massive black hole. After the subsystems become well mixed, the two massive components soon form a hard binary which gains energy by ejecting other members. Such a massive binary has a large cross section and can be very effective in depleting bound stars from the core. Again high-velocity escapers are produced, with their terminal speeds related to the shrinking binary size which eventually leads to coalescence.

The formation of gas giant planets is assumed to be preceded by the emergence of solid cores in the conventional sequential-accretion paradigm. This hypothesis implies that the presence of earth-like planets can be inferred from the detection of gas giants. A similar prediction cannot be made with the gravitational instability (hereafter GI) model which assumes that gas giants (hereafter giants) formed from the collapse of gas fragments analogous to their host stars. We propose an observational test for the determination of the dominant planet-formation channel. Based on the sequential-accretion (hereafter SA) model, we identify several potential avenues which may lead to the prolific formation of a population of close-in earth-mass ($M_\oplus$) planets (hereafter close-in earths) around stars with 1) short-period or 2) solitary eccentric giants and 3) systems which contain intermediate-period resonant giants. In contrast, these close-in earths are not expected to form in systems where giants originated rapidly through GI. As a specific example, we suggest that the SA processes led to the formation of the 7.5 $M_\oplus$ planet around GJ 876 and predict that it may have an atmosphere and envelope rich in O$_2$ and liquid water. Assessments of the ubiquity of these planets will lead to 1) the detection of the first habitable terrestrial planets, 2) the verification of the dominant mode of planet formation, 3) an estimate of the fraction of earth-harboring stars, and 4) modification of bio-marker signatures.

The old open cluster M67 is an ideal testbed for current cluster evolution models because of its dynamically evolved structure and rich stellar populations that show clear signs of interaction between stellar, binary and cluster evolution. Here we present the first truly direct N-body model for M67, evolved from zero age to 4 Gyr taking full account of cluster dynamics as well as stellar and binary evolution. Our preferred model starts with 12000 single stars and 12000 binaries placed in a Galactic tidal field at 8.0 kpc from the Galactic Centre. Our choices for the initial conditions and for the primordial binary population are explained in detail. At 4 Gyr, the age of M67, the total mass has reduced by 90% as a result of mass loss and stellar escapes. The mass and half-mass radius of luminous stars in the cluster are a good match to observations although the model is more centrally concentrated than observations indicate. The stellar mass and luminosity functions are significantly flattened by preferential escape of low-mass stars. We find that M67 is dynamically old enough that information about the initial mass function is lost, both from the current luminosity function and from the current mass fraction in white dwarfs. The model contains 20 blue stragglers at 4 Gyr which is slightly less than the 28 observed in M67. Nine are in binaries. The blue stragglers were formed by a variety of means and we find formation paths for the whole variety observed in M67. Both the primordial binary population and the dynamical cluster environment play an essential role in shaping the population. A substantial population of short-period primordial binaries (with periods less than a few days) is needed to explain the observed number of blue stragglers in M67.

In Galactic open clusters, there is an apparent paucity of white dwarfs compared to the number expected assuming a reasonable initial mass function and that main-sequence stars with initial mass <= 8 M_sun become white dwarfs. We suggest that this lack of white dwarfs is due at least in part to dynamical processes. Non-spherically symmetric mass loss during the post-main-sequence evolution would lead to a few km/s isotropic recoil speed for the white dwarf remnant. This recoil speed can cause a substantial fraction of the white dwarfs formed in a cluster to leave the system. We investigate this dynamical process by carrying out high-precision N-body simulations of intermediate-mass open clusters, where we apply an isotropic recoil speed to the white dwarf remnants. Our models suggest that almost all white dwarfs would be lost from the cluster if the average recoil speed exceeds twice the velocity dispersion of the cluster.

We summarize the main results from MODEST-1, the first workshop on MOdeling DEnse STellar systems. Our goal is to go beyond traditional population synthesis models, by introducing dynamical interactions between single stars, binaries, and multiple systems. The challenge is to define and develop a software framework to enable us to combine in one simulation existing computer codes in stellar evolution, stellar dynamics, and stellar hydrodynamics. With this objective, the workshop brought together experts in these three fields, as well as other interested astrophysicists and computer scientists. We report here our main conclusions, questions and suggestions for further steps toward integrating stellar evolution and stellar (hydro)dynamics.

We discuss a new N-body simulation method for studying black hole binary dynamics. This method avoids previous numerical problems due to large mass ratios and trapped orbits with short periods. A treatment of relativistic effects is included when the associated time-scale becomes small. Preliminary results with up to N = 240 000 particles are obtained showing systematic eccentricity growth until the relativistic regime is reached, with subsequent coalescence in some cases.

We first discuss two-body and chain regularization methods for direct N-body simulations on HARP-2 and GRAPE-6. The former is used for accurate integration of perturbed binaries and hierarchies, whereas the latter deals with strong interactions involving binaries. Combined with a powerful stability criterion for hierarchical systems, this versatile treatment provides an efficient way of studying all dynamical processes in globular clusters containing a realistic population of primordial binaries. These algorithms are also ideal for modelling a variety of astrophysical effects, such as averaging over Kozai cycles, tidal circularization, spin-orbit coupling and stellar collisions, where well-defined elements are used to describe near-singular solutions. We also describe a new time-transformed leapfrog scheme which has been developed to deal with black-hole binaries in galactic centres. This formulation is valid for large mass ratios and dominant two-body motions in a compact subsystem can be treated accurately in the post-Newtonian approximation. This method has been combined with the existing regularization algorithms into a new simulation code NBODY7. Some preliminary results of a test problem illustrating possible applications are presented.

We give a full description of the code NBODY2 for direct integration of the gravitational N-body problem. The method of solution is based on the neighbour scheme of Ahmad & Cohen (1973) which speeds up the force calculation already for quite modest particle numbers. Derivations of all the relevant mathematical expressions are given, together with a detailed discussion of the algorithms. The code may be used to study a wide variety of self-consistent problems based on a small softening of the interaction potential.

(shortened) Direct N-body calculations are presented of the formation of Galactic clusters using GasEx, which is a variant of the code Nbody6. The calculations focus on the possible evolution of the Orion Nebula Cluster (ONC) by assuming that the embedded OB stars explosively drove out 2/3 of its mass in the form of gas about 0.4 Myr ago. A bound cluster forms readily and survives for 150 Myr despite additional mass loss from the large number of massive stars, and the Galactic tidal field. This is the very first time that cluster formation is obtained under such realistic conditions. The cluster contains about 1/3 of the initial 10^4 stars, and resembles the Pleiades Cluster to a remarkable degree, implying that an ONC-like cluster may have been a precursor of the Pleiades. This scenario predicts the present expansion velocity of the ONC, which will be measurable by upcoming astrometric space missions (DIVA and GAIA). These missions should also detect the original Pleiades members as an associated expanding young Galactic-field sub-population. The results arrived at here suggest that Galactic clusters form as the nuclei of expanding OB associations.

We present a state-of-the-art N-body code which includes a detailed treatment of stellar and binary evolution as well as the cluster dynamics. This code is ideal for investigating all aspects relating to the evolution of star clusters and their stellar populations. It is applicable to open and globular clusters of any age. We use the N-body code to model the blue straggler population of the old open cluster M67. Preliminary calculations with our binary population synthesis code show that binary evolution alone cannot explain the observed numbers or properties of the blue stragglers. On the other hand, our N-body model of M67 generates the required number of blue stragglers and provides formation paths for all the various types found in M67. This demonstrates the effectiveness of the cluster environment in modifying the nature of the stars it contains and highlights the importance of combining dynamics with stellar evolution. We also perform a series of N = 10000 simulations in order to quantify the rate of escape of stars from a cluster subject to the Galactic tidal field.

Multiple systems play an important role in the evolution of star clusters. First we discuss several formation mechanisms which depend on the presence of binaries, either primordial or of dynamical origin. Hierarchical configurations are often stable over long times and yet may experience evolution of the internal orbital parameters. We describe an attempt to model the eccentricity change induced by the outer component using an averaging method, together with the effects due to tidal dissipation and apsidal motion acting on the inner binary. This treatment is adopted for systems with high induced eccentricity which gives rise to some interesting outcomes of significant period shrinkage.

This paper concentrates on four key tools for performing star cluster simulations developed during the last decade which are sufficient to handle all the relevant dynamical aspects. First we discuss briefly the Hermite integration scheme which is simple to use and highly efficient for advancing the single particles. The main numerical challenge is in dealing with weakly and strongly perturbed hard binaries. A new treatment of the classical Kustaanheimo-Stiefel two-body regularization has proved to be more accurate for studying binaries than previous algorithms based on divided differences or Hermite integration. This formulation employs a Taylor series expansion combined with the Stumpff functions, still with one force evaluation per step, which gives exact solutions for unperturbed motion and is at least comparable to the polynomial methods for large perturbations. Strong interactions between hard binaries and single stars or other binaries are studied by chain regularization which ensures a non-biased outcome for chaotic motions. A new semi-analytical stability criterion for hierarchical systems has been adopted and the long-term effects on the inner binary are now treated by averaging techniques for cases of interest. These modifications describe consistent changes of the orbital variables due to large Kozai cycles and tidal dissipation. The range of astrophysical processes which can now be considered by N-body simulations include tidal capture, circularization, mass transfer by Roche-lobe overflow as well as physical collisions, where the masses and radii of individual stars are modelled by synthetic stellar evolution.

We consider the use of N-body simulations for studying the evolution of rich star clusters (i.e. globular clusters). The dynamical processes included in this study are restricted to gravitational (point-mass) interactions, the steady tidal field of a galaxy, and instantaneous mass loss resulting from stellar evolution. With evolution driven by these mechanisms, it is known that clusters fall roughly into two broad classes: those that dissipate promptly in the tidal field, as a result of mass loss, and those that survive long enough for their evolution to become dominated by two-body relaxation. The time scales of the processes we consider scale in different ways with the number of stars in the simulation, and the main aim of the paper is to suggest how the scaling of a simulation should be done so that the results are representative of the evolution of a `real' cluster. We investigate three different ways of scaling time. One of these is appropriate to the first type of cluster, i.e. those that dissipate rapidly, and similarly a second scaling is appropriate only to the second (relaxation-dominated) type. We also develop a hybrid scaling which is a satisfactory compromise for both types of cluster. Finally we present evidence that the widely used Fokker-Planck method produces models which are in good agreement with N-body models of those clusters which are relaxation dominated, at least for N-body models with several thousand particles, but that the Fokker-Planck models evolve too fast for clusters which dissipate promptly.

In this paper we study the formation, evolution and disruption of hierarchical systems in open clusters. With this purpose, N-body simulations of star clusters containing an initial population of binaries have been carried out using Aarseth's NBODY4 and NBODY5 codes. Stable triples may form from strong interactions of two binaries in which the widest pair is disrupted. The most frequent type of hierarchical systems found in the cluster models are triples in which the outer star is single, but in some cases the outer body is also a binary, giving a hierarchical quadruple. Many triple systems are non-coplanar and the presence of even a very distant and small outer companion may affect the orbital parameters of the inner binary, including a possible mechanism of significant shrinkage if the binary experiences a weak tidal dissipation. The main features of these systems are analyzed in order to derive general properties which can be checked by observations. The formation of exotic objects, such as blue stragglers and white dwarf binaries, inside hierarchical triple systems is particularly interesting.

A collisional $N$-body simulation using NBODY5 on a single CRAY YMP processor is followed well into the post-collapse regime. This is presently one of the largest particle numbers of all such models published, but some data for an even larger $N$ produced by using special-purpose computers have recently been presented. In contrast to previous ensemble-averaged $N$-body simulations the noise here is low enough to just compare this one single run with the expectations from statistical models based on the Fokker-Planck approximation. Agreement is as good as could be expected for the case of the evolution of the Lagrangian radii, radial and tangential velocity dispersions and various core quantities. We briefly discuss approximate models to understand number and energy of escapers and the question of gravothermal core oscillations; although the system exhibits post-collapse oscillations they turn out to be directly binary driven and we cannot prove the existence of a gravothermal expansion at this particle number. Finally in a detailed examination of the wandering of the density centre we find in contrast to some previous studies a clear long-time period of the order of about 14 half-mass crossing times.