Scaling of N-body calculations

We report results of collisional N -body simulations aimed at studying the N dependence of the dynamical evolution of star clusters. Our clusters consist of equal-mass stars and are in virial equilibrium. Clusters moving in external tidal fields and clusters limited by a cut-off radius are simulated. Our main focus is to study the dependence of the lifetimes of the clusters on the number of cluster stars and the chosen escape condition.
We find that star clusters in external tidal fields exhibit a scaling problem in the sense that their lifetimes do not scale with the relaxation time. Isolated clusters show a similar problem if stars are removed only after their distance to the cluster centre exceeds a certain cut-off radius. If stars are removed immediately after their energy exceeds the energy necessary for escape, the scaling problem disappears.
We show that some stars that gain the energy necessary for escape are scattered to lower energies before they can leave the cluster. As the efficiency of this process decreases with increasing particle number, it causes the lifetimes not to scale with the relaxation time. Analytic formulae are derived for the scaling of the lifetimes in the different cases.     

Basic N-body modelling of the evolution of globular clusters – I. Time scaling

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; 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 that are in good agreement with N -body models of those clusters that are relaxation-dominated, at least for N -body models with several thousand particles, but that the Fokker–Planck models evolve too fast for clusters that dissipate promptly.     

Evolution of the mass function of the Galactic globular cluster system

In this work we investigate the evolution of the mass function of the Galactic globular cluster system (GCMF) taking into account the effects of stellar evolution, two-body relaxation, disc shocking and dynamical friction on the evolution of individual globular clusters. We have adopted a lognormal initial GCMF and considered a wide range of initial values for the dispersion, σ, and the mean value, 〈log  M 〉. We have studied in detail the dependence on the initial conditions of the final values of σ, 〈log  M 〉, the fraction of the initial number of clusters surviving after one Hubble time and the difference between the properties of the GCMF of clusters closer to the Galactic Centre and those of clusters located in the outer regions of the Galaxy. In most of the cases considered, evolutionary processes alter significantly the initial population of globular clusters and the disruption of a significant number of globular clusters leads to a flattening in the spatial distribution of clusters in the central regions of the Galaxy. The initial lognormal shape of the GCMF is preserved in most cases and if a power-law in M is adopted for the initial GCMF, evolutionary processes tend to modify it into a lognormal GCMF. The difference between initial and final values of σ and 〈log  M 〉 as well as the difference between the final values of these parameters for inner and outer clusters can be positive or negative depending on initial conditions. A significant effect of evolutionary processes does not necessarily give rise to a strong trend of 〈log  M 〉 with the galactocentric distance. The existence of a particular initial GCMF able to keep its initial shape and parameters unaltered during the entire evolution through a subtle balance between disruption of clusters and evolution of the masses of those which survive, suggested by Vesperini, is confirmed.     

Gravothermal catastrophe in anisotropic spherical systems

In this paper we investigate the gravothermal instability of spherical stellar systems endowed with a radially anisotropic velocity distribution. We focus our attention on the effects of anisotropy on the conditions for the onset of instability and in particular we study the dependence of the spatial structure of critical models on the amount of anisotropy present in a system. The investigation has been carried out by the method of linear series which has already been used in the past to study the gravothermal instability of isotropic systems._   We consider models described by King, Wilson and Woolley–Dickens distribution functions. In the case of King and Woolley–Dickens models, our results show that, for quite a wide range of the amount of anisotropy in the system, the critical value of the concentration of the system (defined as the ratio of the tidal to the King core radius of the system) is approximately constant and equal to the corresponding value for isotropic systems. Only for very anisotropic systems does the critical value of the concentration start to change and it decreases significantly as the anisotropy increases and penetrates the inner parts of the system. For Wilson models the decrease of the concentration of critical models is preceded by an intermediate regime in which critical concentration increases, reaches a maximum and then starts to decrease. The critical value of the central potential always decreases as the anisotropy increases.     

A stochastic Monte Carlo approach to model real star cluster evolution – II. Self-consistent models and primordial binaries

The new approach outlined in Paper I to follow the individual formation and evolution of binaries in an evolving, equal point-mass star cluster is extended for the self-consistent treatment of relaxation and close three- and four-body encounters for many binaries (typically a few per cent of the initial number of stars in the cluster mass). The distribution of single stars is treated as a conducting gas sphere with a standard anisotropic gaseous model. A Monte Carlo technique is used to model the motion of binaries, their formation and subsequent hardening by close encounters, and their relaxation (dynamical friction) with single stars and other binaries. The results are a further approach towards a realistic model of globular clusters with primordial binaries without using special hardware. We present, as our main result, the self-consistent evolution of a cluster consisting of 300 000 equal point-mass stars, plus 30 000 equal-mass binaries over several hundred half-mass relaxation times, well into the phase where most of the binaries have been dissolved and evacuated from the core. The cluster evolution is about three times slower than found by Gao et al. Other features are rather comparable. At every moment we are able to show the individual distribution of binaries in the cluster.     

Monte Carlo simulations of star clusters – III. A million-body star cluster

A revision of Stodółkiewicz's Monte Carlo code is used to simulate the evolution of million-body star clusters. The new method treats each superstar as a single star and follows the evolution and motion of all individual stellar objects. The evolution of N -body systems influenced by the tidal field of a parent galaxy and by stellar evolution is presented. All models consist of 1 000 000 stars. The process of energy generation is realized by means of appropriately modified versions of Spitzer's and Mikkola's formulae for the interaction cross-section between binaries and field stars and binaries themselves. The results presented are in good agreement with theoretical expectations and the results of other methods. During the evolution, the initial mass function (IMF) changes significantly. The local mass function around the half-mass radius closely resembles the actual global mass function. At the late stages of evolution, the mass of the evolved stars inside the core can be as high as 97 per cent of the total mass in this region. For the whole system, the evolved stars can compose up to 75 per cent of the total mass. The evolution of cluster anisotropy strongly depends on initial cluster concentration, IMF and the strength of the tidal field. The results presented are the first step in the direction of simulating the evolution of real globular clusters by means of the Monte Carlo method.     

Fluctuations in isothermal spheres

Isolated isothermal spheres of N gravitationally interacting points with equal mass are believed to be stable when density contrasts do not exceed 709. That stability limit does not, however, take into consideration fluctuations of temperature near the onset of instability. These are important when N is finite.
Here we correlate global mean quadratic temperature fluctuations with the onset of instability. We show that such fluctuations trigger instability when the density contrast reaches a value near 709×exp(−3.3 N ^−1/3). These lower values of limiting density contrasts are significantly smaller than 709 when N is not very large, and this suggests (i) that numerical calculations with small N may not reflect correctly the onset of core collapse in clusters with large N , and (ii) that a greater number of globular clusters than is normally believed may already be in an advanced stage of core collapse, because most of the observed globular clusters with parameters that fit quasi-isothermal configurations are close to marginal stability.     

On the clustering phase transition in self-gravitating N-body systems

The thermodynamic behaviour of self-gravitating N -body systems has been worked out by borrowing a standard method from molecular dynamics. The link between dynamics and thermodynamics is made in the microcanonical ensemble of statistical mechanics. Through the computation of basic thermodynamic observables and of the equation of state in the     plane, the clustering phase transition appears to be of the second-order type. The dynamical–microcanonical averages are compared with their corresponding canonical ensemble averages, obtained through standard Monte Carlo computations. The latter seem to have completely lost any information about the phase transition. Finally, our results – obtained in a 'microscopic' framework – are compared with some existing theoretical predictions – obtained in a 'macroscopic' (thermodynamic) framework: qualitative and quantitative agreement is found, with an interesting exception.     

Dynamical evolution of rotating stellar systems – II. Post-collapse, equal-mass system

We present the first post-core-collapse models of initially rotating star clusters, using the numerical solution of an orbit-averaged 2D Fokker–Planck equation. Based on the code developed by Einsel & Spurzem, we have improved the speed and the stability and included the steady three-body binary heating source. We have confirmed that rotating clusters, whether they are in a tidal field or not, evolve significantly faster than non-rotating ones. Consequences for the observed shapes, density distribution and kinematic properties of young and old star clusters are discussed. The results are compared with gaseous and 1D Fokker–Planck models in the non-rotating case.     

N-body simulations using customized potential–density pair basis sets

Potential–density pair basis sets can be used for highly efficient N -body simulation codes, but they suffer from a lack of versatility, i.e. a basis set has to be constructed for each different class of stellar system. We present numerical techniques for generating a biorthonormal potential–density pair basis set that has a general specified pair as its lowest-order member. We go on to demonstrate how the set can be used to construct N -body equilibria, which we then evolve using an N -body code that calculates forces using the basis set.     

On a time-symmetric Hermite integrator for planetary N-body simulation

We describe a P(EC) ⁿ Hermite scheme for planetary N -body simulation. The fourth-order implicit Hermite scheme is a time-symmetric integrator that has no secular energy error for the integration of periodic orbits with time-symmetric time-steps. In general N -body problems, however, this advantage is of little practical significance, since it is difficult to achieve time-symmetry with individual variable time-steps. However, we can easily enjoy the benefit of the time-symmetric Hermite integrator in planetary N -body systems, where all bodies spend most of the time on nearly circular orbits. These orbits are integrated with almost constant time-steps even if we adopt the individual time-step scheme. The P(EC) ⁿ Hermite scheme and almost constant time-steps reduce the integration error greatly. For example, the energy error of the P(EC)² Hermite scheme is two orders of magnitude smaller than that of the standard PEC Hermite scheme in the case of an N  = 100,  m  = 10²⁵ g planetesimal system with the rms eccentricity 〈 e ²〉^1/2 ≲0.03.     

Chaos in Hill's generalized problem: from the solar system to black holes

We generalize the well‐known Hill's circular restricted three‐body problem by assuming that the primary generates a Schwarzschild‐type field of the form U = A/r + B/r³. The term in B influences the particle, but not the far secondary. Many concrete astronomical situations can be modelled via this problem. For the two‐body problem primary‐particle, a homoclinic orbit is proved to exist for a continuous range of parameters (the constants of energy and angular momentum, and the field parameter B > 0). Within the restricted three‐body system, we prove that, under sufficiently small perturbations from the secondary, the homoclinic orbit persists, but its stable and unstable manifolds intersect transversely. Using a result of symbolic dynamics, this means the existence of a Smale horseshoe, hence chaotic behaviour. Moreover, we find that Hill's generalized problem (in our sense) is nonintegrable. (© 2005 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)