Monte Carlo simulations of star clusters – IV. Calibration of the Monte Carlo code and comparison with observations for the open cluster M67

We outline the steps needed in order to incorporate the evolution of single and binary stars into a particular Monte Carlo code for the dynamical evolution of a star cluster. We calibrate the results against N -body simulations, and present models for the evolution of the old open cluster M67 (which has been studied thoroughly in the literature with N -body techniques). The calibration is done by choosing appropriate free code parameters. We describe, in particular, the evolution of the binary, white dwarf and blue straggler populations, though not all channels for blue straggler formation are represented yet in our simulations. Calibrated Monte Carlo runs show good agreement with results of N -body simulations not only for global cluster parameters, but also for, for example, binary fraction, luminosity function and surface brightness. Comparison of Monte Carlo simulations with observational data for M67 shows that it is possible to get reasonably good agreement between them. Unfortunately, because of the large statistical fluctuations of the numerical data and uncertainties in the observational data the inferred conclusions about the cluster initial conditions are not firm.     

Direct N-body modelling of stellar populations: blue stragglers in M67

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 =10 000 simulations in order to quantify the rate of escape of stars from a cluster subject to the Galactic tidal field.     

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.     

Monte Carlo simulations of star clusters – II. Tidally limited, multimass systems with stellar evolution

A revision of Stodółkiewicz's Monte Carlo code is used to simulate evolution of large star clusters. The new method treats each superstar as a single star and follows the evolution and motion of all individual stellar objects. A survey of the evolution of N -body systems influenced by the tidal field of a parent galaxy and by stellar evolution is presented. 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 (Fokker–Planck, Monte Carlo and N -body). The initial rapid mass loss, resulting from stellar evolution of the most massive stars, causes expansion of the whole cluster and eventually leads to the disruption of less bound systems ( W ₀=3). Models with larger W ₀ survive this phase of evolution and then undergo core collapse and subsequent post-collapse expansion, like isolated models. The expansion phase is eventually reversed when tidal limitation becomes important. The results presented are the first major step in the direction of simulating evolution of real globular clusters by means of the Monte Carlo method.     

Single star scattering on an open cluster

We examine the scattering of single stars from an open star cluster. The probability of the capture of a star by a star cluster is dependent on the velocity and mass of the star, and the stars that are not captured experience a velocity change. For low-velocity stars there is an exponential decrease of the capture probability with the initial velocity, and the velocity change decreases almost linearly. For high-velocity stars there is a v ⁻⁶ dependence for the capture probability, and a v ⁻¹ dependence for the velocity change. Analytical estimations, Monte Carlo and full N -body simulations are all in good agreement.     

How well do starlab and nbody4 compare? I. Simple models

N -body simulations are widely used to simulate the dynamical evolution of a variety of systems, among them star clusters. Much of our understanding of their evolution rests on the results of such direct N -body simulations. They provide insight in the structural evolution of star clusters, as well as into the occurrence of stellar exotica. Although the major pure N -body codes starlab/kira and nbody4 are widely used for a range of applications, there is no thorough comparison study yet.
Here, we thoroughly compare basic quantities as derived from simulations performed either with starlab/kira or nbody4 .
We construct a large number of star cluster models for various stellar mass function settings (but without stellar/binary evolution, primordial binaries, external tidal fields, etc.), evolve them in parallel with starlab/kira and nbody4 , analyse them in a consistent way and compare the averaged results quantitatively. For this quantitative comparison, we develop a bootstrap algorithm for functional dependencies.
We find an overall excellent agreement between the codes, both for the clusters' structural and energy parameters as well as for the properties of the dynamically created binaries. However, we identify small differences, like in the energy conservation before core collapse and the energies of escaping stars, which deserve further studies.
Our results reassure the comparability and the possibility to combine results from these two major N -body codes, at least for the purely dynamical models (i.e. without stellar/binary evolution) we performed. Further detailed comparison studies for more complex systems, e.g. including stellar/binary evolution, are required.     

Cosmological evolution and hierarchical galaxy formation

We calculate the rate at which dark matter haloes merge to form higher mass systems. Two complementary derivations using Press–Schechter theory are given, both of which result in the same equation for the formation rate. First, a derivation using the properties of the Brownian random walks within the framework of Press–Schechter theory is presented. We then use Bayes' theorem to obtain the same result from the standard Press–Schechter mass function. The rate obtained is shown to be in good agreement with results from Monte Carlo and N -body simulations. We illustrate the usefulness of this formula by calculating the expected cosmological evolution in the rate of star formation that is due to short-lived, merger-induced starbursts. The calculated evolution is well-matched to the observed evolution in ultraviolet luminosity density, in contrast to the lower rates of evolution that are derived from semi-analytic models that do not include a dominant contribution from starbursts. Hence we suggest that the bulk of the observed ultraviolet starlight at z >1 arises from merger-induced starbursts. Finally, we show that a simple merging-halo model can also account for the bulk of the observed evolution in the comoving quasar space density.     

Multimass schemes for collisionless N-body simulations

We present a general scheme for constructing Monte Carlo realizations of equilibrium, collisionless galaxy models with known distribution function (DF) f ₀. Our method uses importance sampling to find the sampling DF f _s that minimizes the mean-square formal errors in a given set of projections of the DF f ₀. The result is a multimass N -body realization of the galaxy model in which 'interesting' regions of phase space are densely populated by lots of low-mass particles, increasing the effective N there, and less interesting regions by fewer, higher mass particles.
As a simple application, we consider the case of minimizing the shot noise in estimates of the acceleration field for an N -body model of a spherical Hernquist model. Models constructed using our scheme easily yield a factor of ∼100 reduction in the variance at the central acceleration field when compared to a traditional equal-mass model with the same number of particles. When evolving both models with a real N -body code, the diffusion coefficients in our model are reduced by a similar factor. Therefore, for certain types of problems, our scheme is a practical method for reducing the two-body relaxation effects, thereby bringing the N -body simulations closer to the collisionless ideal.     

Multimass spherical structure models for N-body simulations

We present a simple and efficient method to set up spherical structure models for N -body simulations with a multimass technique. This technique reduces by a substantial factor the computer run time needed in order to resolve a given scale as compared to single-mass models. It therefore allows to resolve smaller scales in N -body simulations for a given computer run time. Here, we present several models with an effective resolution of up to  1.68 × 10⁹  particles within their virial radius which are stable over cosmologically relevant time-scales. As an application, we confirm the theoretical prediction by Dehnen that in mergers of collisionless structures like dark matter haloes always the cusp of the steepest progenitor is preserved. We model each merger progenitor with an effective number of particles of approximately 10⁸ particles. We also find that in a core–core merger the central density approximately doubles whereas in the cusp–cusp case the central density only increases by approximately 50 per cent. This may suggest that the central regions of flat structures are better protected and get less energy input through the merger process.     

H i and CO observations of Arp 104: a spiral–elliptical interacting pair

We present data probing the spatial and kinematical distribution of both the atomic (H  i ) and molecular (CO) gas in NGC 5218, the late-type barred spiral galaxy in the spiral–elliptical interacting pair, Arp 104. We consider these data in conjunction with far-infrared and radio-continuum data, and N -body simulations, to study the galaxies interactions, and the star formation properties of NGC 5218. We use these data to assess the importance of the bar and tidal interaction on the evolution of NGC 5218, and the extent to which the tidal interaction may have been important in triggering the bar. The molecular gas distribution of NGC 5218 appears to have been strongly affected by the bar; the distribution is centrally condensed with a very large surface density in the central region. The N -body simulations indicate a time-scale since perigalacticon of  ∼3 × 10⁸ yr  , which is consistent with the interaction having triggered or enhanced the bar potential in NGC 5218, leading to inflow and the large central molecular gas density observed. Whilst NGC 5218 appears to be undergoing active star formation, its star formation efficiency is comparable to a 'normal' SBb galaxy. We propose that this system may be on the brink of a more active phase of star formation.     

Star clusters with primordial binaries – I. Dynamical evolution of isolated models

In order to interpret the results of complex realistic star cluster simulations, which rely on many simplifying approximations and assumptions, it is essential to study the behaviour of even more idealized models, which can highlight the essential physical effects and are amenable to more exact methods. With this aim, we present the results of N -body calculations of the evolution of equal-mass models, starting with primordial binary fractions of 0–100 per cent, with values of N ranging from 256 to 16 384. This allows us to extrapolate the main features of the evolution to systems comparable in particle number with globular clusters.
In this range, we find that the steady-state 'deuterium main sequence' is characterized by a ratio of the core radius to half-mass radius that follows qualitatively the analytical estimate by Vesperini & Chernoff, although the N dependence is steeper than expected. Interestingly, for an initial binary fraction f greater than 10 per cent, the binary heating in the core during the post-collapse phase almost saturates (becoming nearly independent of f ), and so little variation in the structural properties is observed. Thus, although we observe a significantly lower binary abundance in the core with respect to the Fokker–Planck simulations by Gao et al., this is of little dynamical consequence.
At variance with the study of Gao et al., we see no sign of gravothermal oscillations before 150 half-mass relaxation times. At later times, however, oscillations become prominent. We demonstrate the gravothermal nature of these oscillations.     

The evolution of the binary population in globular clusters: a full analytical computation

I present a simplified analytical model that simulates the evolution of the binary population in a dynamically evolving globular cluster. A number of simulations have been run spanning a wide range in initial cluster and environmental conditions by taking into account the main mechanisms of formation and destruction of binary systems. Following this approach, I investigate the evolution of the fraction, the radial distribution, the distribution of mass ratios and periods of the binary population. According to these simulations, the fraction of surviving binaries appears to be dominated by the processes of binary ionization and evaporation. In particular, the frequency of binary systems changes by a factor of 1–5 depending on the initial conditions and on the assumed initial distribution of periods. The comparison with the existing estimates of binary fractions in Galactic globular clusters suggests that significant variations in the initial binary content could exist among the analysed globular cluster. This model has been also used to explain the observed discrepancy found between the most recent N -body and Monte Carlo simulations in the literature.     

Scaling up tides in numerical models of galaxy and halo formation

The purpose of this article is to show that when dynamically cold, dissipationless self-gravitating systems collapse, their evolution is a strong function of the symmetry in the initial distribution. We explore with a set of pressureless homogeneous fluids the time evolution of ellipsoidal distributions and map the depth of potential achieved during relaxation as function of initial ellipsoid axis ratios. We then perform a series of N -body numerical simulations and contrast their evolution with the fluid solutions. We verify an analytic relation between collapse factor and particle number N in spherical symmetry, such that  ∝ N ^1/3  . We sought a similar relation for axisymmetric configurations, and found an empirical scaling relation such that  ∝ N ^1/6  in these cases. We then show that when mass distributions do not respect spherical or axial symmetry, the ensuing gravitational collapse deepens with increasing particle number N but only slowly: 86 per cent of triaxial configurations may collapse by a factor of no more than 40 as   N →∞  . For   N ≈10⁵  and larger, violent relaxation develops fully under the Lin–Mestel–Shu instability such that numerical N -body solutions now resolve the different initial morphologies adequately.     

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.     

Ratios of star cluster core and half-mass radii: a cautionary note on intermediate-mass black holes in star clusters

There is currently much interest in the possible presence of intermediate-mass black holes (IMBHs) in the cores of globular clusters (GCs). Based on theoretical arguments and simulation results it has previously been suggested that a large core radius – or particularly a large ratio of the core radius to half-mass radius – is a promising indicator for finding such a black hole (BH) in a star cluster. In this study N -body models of 100 000 stars with and without primordial binaries are used to investigate the long-term structural evolution of star clusters. Importantly, the simulation data are analysed using the same processes by which structural parameters are extracted from observed star clusters. This gives a ratio of the core and half-mass (or half-light) radii that are directly comparable to the Galactic GC sample. As a result, it is shown that the ratios observed for the bulk of this sample can be explained without the need for an IMBH. Furthermore, it is possible that clusters with large core to half-light radius ratios harbour a BH binary (comprising stellar mass BHs) rather than a single massive BH. This work does not rule out the existence of IMBHs in the cores of at least some star clusters.     

Evolution of the cluster abundance in non-Gaussian models

We carry out N -body simulations of several non-Gaussian structure formation models, including Peebles' isocurvature cold dark matter model, cosmic string models, and a model with primordial voids. We compare the evolution of the cluster mass function in these simulations with that predicted by a modified version of the Press–Schechter formalism. We find that the Press–Schechter formula can accurately fit the cluster evolution over a wide range of redshifts for all of the models considered, with typical errors in the mass function of less than 25 per cent, considerably smaller than the amount by which predictions for different models may differ. This work demonstrates that the Press–Schechter formalism can be used to place strong model-independent constraints on non-Gaussianity in the Universe.     

The binary star population of the young cluster NGC 1818 in the Large Magellanic Cloud

We determine the binary star fraction as a function of radius in NGC 1818, a young rich cluster in the Large Magellanic Cloud, using Hubble Space Telescope images in bands F336W (∼ U ) and F555W (∼ V ). Our sample includes binaries with M _primary ∼ 2–5.5 M_⊙ and M _secondary ≳ 0.7 M_primary. The binary fraction increases towards the cluster centre, from ∼ 20 ± 5 per cent in the outer parts, to ∼ 35 ± 5 per cent inside the core. This increase is consistent with dynamical mass segregation and need not be primordial. We compare our results with expectations from N -body models, and discuss the implications for the formation and early evolution of such clusters.