---
abstract: 'We describe Monte Carlo models for the dynamical evolution of the nearby globular cluster M4. The code includes treatments of two-body relaxation, three- and four-body interactions involving primordial binaries and those formed dynamically, the Galactic tide, and the internal evolution of both single and binary stars. We arrive at a set of initial parameters for the cluster which, after 12Gyr of evolution, gives a model with a satisfactory match to the surface brightness profile, the velocity dispersion profile, and the luminosity function in two fields. We describe in particular the evolution of the core, and find that M4 (which has a classic King profile) is actually a [*post-collapse*]{} cluster, its core radius being sustained by binary burning. We also consider the distribution of its binaries, including those which would be observed as photometric binaries and as radial-velocity binaries. We also consider the populations of white dwarfs, neutron stars, black holes and blue stragglers, though not all channels for blue straggler formation are represented yet in our simulations.'
author:
- |
    Mirek Giersz$^{1}$[^1] and Douglas C. Heggie$^{2}$\
    $^{1}$Nicolaus Copernicus Astronomical Centre, Polish Academy of Sciences, ul. Bartycka 18, 00-716 Warsaw, Poland\
    $^{2}$University of Edinburgh, School of Mathematics and Maxwell Institute for Mathematical Sciences, King’s Buildings, Edinburgh EH9 3JZ, UK
date: 'Accepted …. Received …; in original form …'
title: 'Monte Carlo Simulations of Star Clusters - V. The globular cluster M4'
---

\[firstpage\]

stellar dynamics – methods: numerical – binaries: general – globular clusters: individual: M4

Introduction
============

The present paper opens up a new road in the study of the dynamical evolution of globular clusters. We adopt the Monte Carlo method of Giersz [@Gi1998; @Gi2001; @Gi2006], which in recent years has been enhanced to deal quite realistically with the stellar evolution of single and binary stars, to study the dynamical history of the nearby globular cluster M4. An earlier version of the code had already been used to study the dynamical history of $\omega$ Cen [@GH2003], but at that time the treatment of stellar evolution was primitive and there were no binaries. The new code has been thoroughly tested on smaller systems, by comparison with $N$-body simulations and observations of the old open star cluster M67 [@GH2008]. There we showed that the Monte Carlo code could produce data of a similar level of detail and realism as the best $N$-body codes. Now for the first time we consider much richer systems, with about half a million stars initially, which are at present beyond the reach of $N$-body methods.

This paper has a place within a long tradition of the modelling of globular star clusters, but the place is a distinctive one. First, we are not concerned with a static model of a star cluster at the present day, like a King model. We are concerned with issues where the dynamical history of the star cluster is important, where static models are uninformative. Secondly, our aim is to construct a model of a specific star cluster, rather than trying to understand the general properties of the evolution of a population of star clusters. This has been done before, and a brief history is outlined in @GH2008, but the present work takes these efforts onto a new level of realism, in terms of the description of stellar evolution, and dynamical interactions involving binary stars.

This problem is not easy. Not only is it necessary to use an elaborate technique for simulating the relevant astrophysical processes, but it is necessary also to search for initial conditions which, after about 12Gyr of evolution, lead to an object resembling a given star cluster. By “resembling” we do not simply mean matching the overall mass, radius and binary fraction of a cluster, for example, for two reasons:

1.  We have found that values found for data in the literature are highly uncertain, and different sources are contradictory. These data are usually derived, in some model-dependent way, from such data as surface-brightness profiles and velocity dispersion profiles, and we prefer to compare our models directly with this data, and not with inferred global parameters.

2.  We have found that, even if one achieves a satisfactory fit to these profiles, the model may give a very poor comparison with the luminosity function.

From these considerations we conclude that a model which aims to fit only the mass and radius of a star cluster (say) may be very far from the truth.

Tackling this difficult problem is not just interesting, however. We have been motivated by a number of pressing astrophysical problems. For example, the two nearby star clusters M4 (the subject of this study) and NGC 6397 (which we shall consider in our next paper in this series) have rather similar mass and radius, and yet one has a classic King profile, while the other is a well-studied example of a cluster with a “collapsed core” [@Tr1995]. Among possible explanations one may consider differences in the population of binaries, which are known to affect core properties, or in tidal effects. Indeed the present paper will show that these two clusters may be much more similar than one would suppose from the surface brightness profiles alone.

A second motivation for our work is our involvement in observational programmes aimed at characterising the binary populations in globular clusters. What differences (e.g. in the distributions of periods and abundances) should one expect to find between the core and the halo? These issues are important in the planning of observations, and in their interpretation.

The cluster M4 is the focus of much of this effort because it is nearby, making it a relatively easy target for deep observational study. It was the first globular cluster to yield a deep sequence of white dwarfs [@Ri1995]. More recently it has been subjected to an intensive observational programme by the Padova group [@bedinetal2001; @bedinetal2003; @andersonetal2006], which includes searches for radial-velocity binaries in the upper main sequence [@Somm]. It also turns out to be a cluster which (we conclude) started with only about half a million stars, which facilitates the modelling. Along with the open cluster M67, M4 was chosen by the international MODEST consortium, at a meeting in Hamilton in 2005, as the focus for joint effort by theorists and observers, to cast light on its binary population and dynamical properties. M67 has been modelled very successfully by @Hu2005, using $N$-body techniques, and this paper represent the first theoretical step in a similar study of M4.

The paper is organised as follows. First, we summarise features of the code and the models, the data we used, and our approach to the problem of finding initial conditions for M4. Then we present data for our best models: surface brightness and velocity dispersion profiles, luminosity functions, the properties of the binary population, white dwarfs and other degenerate remnants, and the inferred dynamical state of the cluster. The final section summarises our findings and discusses them in the context of work on other clusters, including objects to which we will turn in future papers.

Methods
=======

The Monte Carlo Code
--------------------

The details of our simulation method have been amply described in previous papers in this series. Each star in a spherical star cluster is represented by its mass, energy and angular momentum, and its stellar evolutionary state may be computed at any time using synthetic formulae for single and binary evolution. It may be a binary or a special kind of single star that has been created in a collision or merger event.

Neighbouring stars interact with each other in accordance (in a statistical sense) with the theory of two-body relaxation. If one or both of the participants is a binary, the probability of an encounter affecting the internal dynamics is calculated according to analytic cross sections, which also determine the outcome. This is one of the main shortcomings of the code, as these cross sections are not well known in the case of unequal masses, and also the possibility of stellar collisions during long-lived temporary capture is excluded.

A star or binary may escape if its energy exceeds a certain value, which we choose to be lower than the energy at the nominal tidal radius, in order to improve the scaling of the lifetime with $N$, as explained in @GH2008. This is the second main shortcoming of the models, as it leads to a cutoff radius of the model that is smaller than the true tidal radius, and this lowers the surface density profile in the outer parts of the system.

A difficulty in applying the Monte Carlo code to M4 is that it employs a static tide, whereas the orbit of M4 appears to be very elliptical [@Di1999]. We have to assume that a cluster can be placed in a steady tide of such a strength that the cluster loses mass at the same average rate as it would on its true orbit. Some support for this procedure comes from $N$-body modelling. @BM2003 show that clusters on an elliptical orbit between about 2.8 and 8.5kpc dissolve on a time scale intermediate between that for circular orbits at these two radii, and that the dissolution time scales in almost the same way with the size of the system. @wilkinsonetal2003 show that the core radius of a cluster on an elliptic orbit evolves in very nearly the same way as in a cluster with a circular orbit at the time-averaged galactocentric distance.

All other free parameters of the code (e.g. the coefficient of $N$ in the Coulomb logarithm) take the optimal values found in the above study.

Initial Conditions
------------------

The initial models are as specified in Table \[tab:ics\]. Many of these features (e.g. the properties of the binaries, except for their overall abundance) were inherited from our modelling of the open cluster M67, and those were in turn mainly drawn from the work of @Hu2005. Some of the parameters were taken to be freely adjustable, and this freedom was exploited in the search for an acceptable fit to the current observational data.

  Fixed parameters            
  --------------------------- ------------------------------
  Structure                   Plummer model
  Initial mass function$^1$   @Kr1993 in the
                              range \[0.1,50\]$\msun$
  Binary mass distribution    @Kr1991
  Binary mass ratio           Uniform (with component
                              masses restricted as for
                              single stars)
  Binary semi-major axis      Uniform in log, $2(R_1+R_2)$
                              to 50AU
  Binary eccentricity         Thermal, with eigenevolution
                              [@Kr1995]
  Metallicity $Z$             0.002
  Age                         12Gyr [@Ha2004]
  Free parameters             
  Mass                        $M$
  Tidal radius                $r_t$
  Half-mass radius            $r_h$
  Binary fraction             $f_b$
  Slope of the lower mass     
  function                    $\alpha$ (Kroupa = 1.3)

  : Initial parameters for M4

\[tab:ics\]

Observational data and its computation {#sec:data_computation}
--------------------------------------

Our first task was to iterate on the initial parameters of our model in order to produce a satisfactory fit to a range of observational data at an age of 12Gyr. The data we adopted are as follows:

1.  Surface brightness profile: here we used the compilation by @Tr1995, where the surface brightness is expressed in $V$ magnitudes per square arcsec.

2.  Radial velocity profile: this came from [@Pe1995], and is the result of binning data the on the radial velocities of nearly 200 stars. Strictly we should refer to the line-of-sight velocity dispersion, as “radial” velocity has a different meaning in the Monte Carlo model.

3.  The V-luminosity function (@Ri2004, from which we considered the results for the innermost and outermost of their four annuli). These data are lack correction for completeness, though the completeness factor for the outermost field is plotted in @richeretal2002. For main sequence stars it is almost 100% down to $V=15$, and drops steadily to less than 50% at $V=17$. For the innermost field the completeness correction would be larger.

Now we consider how to compare this data with the output of the Monte Carlo code. This includes a list of each particle in the simulation, along with data on its radius, radial and transverse velocities and absolute magnitude, and numerous other quantities. To construct the surface brightness, we think of each particle as representing a luminous shell of the same radius, and superpose the surface brightness of all shells. While this procedure involves a minimum of effects from binning, or randomly assigning the full position of the star, a shell has an infinite surface brightness at its projected edge. The effect of this, especially from the brightest stars, will be seen in some of the profiles to be presented in this paper. A similar problem arises in actual observations, and is often handled by simply removing the brightest stars from the surface brightness data presented. The output from the model is corrected for extinction (Table \[tab:m4dat\]. The distance of the line of sight from the cluster centre is converted between pc (as in the model) and arcsec (as in the observations) using the distance in the same table.

The line-of-sight velocity dispersion is computed in a similar way. For each particle we calculate the mean square line-of-sight velocity (because the orientation of the transverse component of the velocity is random), and then sum over all particles. In this sum, each particle is weighted by a geometrical factor proportional to the surface density of the particle’s shell along the line of sight. The result is an estimate of the velocity dispersion that is weighted by neither mass nor brightness, but only number. To check the effects of this, we have sometimes calculated velocity dispersions with various cutoffs in the magnitudes of the stars included. The effects are usually quite small.

Computation of the luminosity function is the least problematic. We count stars in each bin in $V$-magnitude, but lying above a line in the colour-magnitude diagram just below the main sequence. Again the contribution of each star is weighted by the same geometrical factor, and the $V$ magnitude is corrected for extinction.

  ----------------------------------- ------------------- --
  Distance from sun                   $^a$1.72kpc         
  Distance from GC                    5.9kpc              
  Mass                                $^a$63 000$\msun$   
  Core radius                         0.53 pc             
  Half-light radius                   2.3 pc              
  Tidal radius                        21 pc               
  Half-mass relaxation time ($R_h$)   660 Myr             
  Binary fraction                     $^a$1-15%           
  $\left[\right.$Fe/H$]$              -1.2                
  Age                                 $^b$12Gyr           
  $A_V$                               $^a$1.33            
  ----------------------------------- ------------------- --

  : Properties of M4

\[tab:m4dat\]

References: All data are from the current version of the catalogue of @Ha1996, except $^a$ @Ri2004 (though this is not always the original reference for the quoted number) and $^b$ @Ha2004.

Finding initial conditions
--------------------------

Here we summarise our experience in approaching this problem. A number of studies (e.g. @BM2003, @La2005) give simple formulae for the evolution with time of the bound mass of a rich star cluster. It is possible to derive similar simple formulae for the evolution of the half-mass radius and other quantities. There are several problems with inverting these formulae, however, i.e. using them to infer the initial parameters of a cluster from its present mass and radius. First, these present-day global parameters are quite uncertain, even to within a factor of two. Second, these formulae depend on the galactic orbit and other parameters which are equally subject to uncertainty. Therefore we have adopted the more straightforward but more laborious approach of iterating on the initial parameters of our models (Tab \[tab:ics\], lower half); that is, we select values for the five stated parameters, run the model, find where the match with the observations is poor, adjust the parameters, and repeat cyclically.

We have employed two methods to facilitate this process to some extent. First, we have often carried out mini-surveys, i.e. small, coarse grid-searches around a given starting model, to find out how changes in individual parameters affect the results. Second, we have used scaling to accelerate the process, and we now describe this method in a little detail.

Suppose we wish to represent a star cluster which has mass $M$ and radius $R$ with a model representing a cluster with a (usually smaller) mass $M\ast$ and radius $R\ast$. Since two-body relaxation dominates much of the dynamics, we insist that the two clusters have the same relaxation time, and so $$\frac{R\ast}{R} =
\left(\frac{M}{M\ast}\right)^{1/3}\left(\frac{\log\gamma
  N\ast}{\log\gamma N}\right)^{2/3},$$ where $N,N\ast$ are the corresponding particle numbers. Thus the model of lower mass has larger radius. The tidal radius is scaled in the same way. Then the observational results (surface brightness profile, etc) can be computed for the model of lower mass and then rescaled (by appropriate factors of the mass and radius) to give a result for the more massive cluster (assuming that the evolution is dominated by the processes of relaxation, stellar evolution and tidal stripping). In fact we have found that this is very successful, in the sense that the inferred best values for the initial conditions change little when runs are carried out with the “correct” (i.e. unscaled) initial mass, certainly when the proportion of binaries is 10% or less. Some aspects of the evolution are not well described by this scaling technique. For example, we do not change the distribution of the semi-major axes of the binaries. In this way the internal evolution of the binaries is correctly modelled (provided that the binaries remain isolated dynamically), though their dynamical interactions with the rest of the system are not. In principle one could scale the semi-major axes in the same way as $R$, but then the internal evolution of the binaries would be altered.

Models of M4
============

Finding initial parameters
--------------------------

Our starting point was our work on the old open cluster M67 [@GH2008], but with a larger initial mass and radius. (@BM2003, for example, suggested that the initial mass of M4 (NGC 6121) was of order $7.5\times10^5\msun$, though we used somewhat smaller values.) To begin with, our choice of initial tidal radius was inferred from the initial mass and the present-day estimates of mass and tidal radius given in Table \[tab:m4dat\], assuming that $r_t\propto M^{1/3}$. At first we adopted similar values for the “concentration” ($r_t/r_h$) and binary fraction as in the modelling of M67, but found that the surface brightness profile fitted poorly (with too large a core) unless the binary fraction was much smaller (5 to 10%) and the concentration much higher.

Before describing our best models, it is worth briefly mentioning one which provided a satisfactory fit to the surface brightness profiles. The fit to the velocity dispersion profile was tolerable, but indicated a model that was too massive by a factor of about 1.4. Its main flaw, however, was in the luminosity function, which was generally too large by a factor in the range 2–3. There are two reasons why this is interesting. One is that, for a long time, models of star clusters were constructed entirely on the basis of the surface brightness and velocity dispersion profiles. It should be realised that such models may be misleading in other ways. Second, this experience underscores the importance of properly normalised luminosity functions. In other words, it is important to know the area of the field where the stars have been counted, or some equivalent representation of properly normalised data. Very often, the emphasis is solely on the [*shape*]{} of the luminosity function, but we have found that the absolute normalisation is an essential constraint.

@BM2003 show that the lower mass function becomes flatter as the fraction of mass lost by the cluster increases (i.e. towards the end of its life). Therefore we could perhaps have improved the fit with the luminosity function by starting with a more massive model and somehow ensuring a larger escape rate so as to leave a similar mass at the present day. Instead, we elected to change the slope of the low-mass IMF from the canonical value of $\alpha=1.3$ [@Kr2007a] to $\alpha=0.9$. (There is some justification for a lower value for low-metallicity populations, though it has been argued [@Kr2007b] that there is no pristine low-metallicity population where the IMF can be inferred securely.)

A Monte Carlo model of M4
-------------------------

By some experimentation we arrived at a model which gave a fair fit to all three kinds of observational data; see Table 3 (where we compare with a King model developed by @Ri2004, and Figs 1-4. It is worth noting that no arbitrary normalisation has been applied in these comparisons between our model and the observations. The surface brightness profile, for example, is computed directly from the V-magnitudes of the stars in the Monte Carlo simulation, as described in Sec.\[sec:data\_computation\]. In the construction of a King model, by contrast, it is often assumed that the mass-to-light ratio is arbitrary.

  --------------------------------- ------------------ ------------------ -------------------------
  Quantity                          MC model           MC model           King model
                                    ($t = 0$)          ($t=12$Gyr)        [@Ri2004]
  Mass ($\msun$)                    $3.40\times10^5$   $4.61\times10^4$   
  Luminosity ($L_\odot$)            $6.1\times10^6$    $2.55\times10^4$   $6.25\times10^4$
  Binary fraction $f_b$             0.07               0.057              0
  Low-mass MF slope $\alpha$        0.9                0.03               0.1
  Mass of white dwarfs ($\msun$)    0                  $1.81\times10^4$   $3.25\times10^{4^\ast}$
  Mass of neutron stars ($\msun$)   0                  $3.24\times10^3$   
  Tidal radius $r_t$ (pc)           35.0               18.0               
  Half-mass radius $r_h$ (pc)       0.58               2.89               
  --------------------------------- ------------------ ------------------ -------------------------

\[tab:mc\_king\]

$\ast$: this is the quoted mass of “degenerates”

### Surface brightness

While the overall surface brightness profile is slightly faint, the most noticeable feature of Fig.\[fig:sbp\] is that the model has a somewhat smaller limiting radius than the observational data. The reason for this is explained in @GH2008: in short, we impose a smaller tidal radius than the nominal tidal radius, in an $N$-dependent way, which is intended to ensure that the overall rate of escape from the model behaves in the same way as in an $N$-body simulation. Another point to notice is the disagreement between the total luminosity of our model and that of the King model quoted in the final column of Table \[tab:mc\_king\]. @Tr1995 give an analytic fit to the surface brightness profile, and we have checked that the integrated value is close to ours.

The data for the Monte Carlo model in Table \[tab:mc\_king\] would be consistent with a galactocentric radius of about 1.7kpc, in an isothermal galaxy model with a circular velocity of 220km/s. While this is certainly much smaller than its current galactocentric distance, a small value was also found (using a similar argument and published values of the mass and tidal radius) by @vdb1995. The orbit given by @Di1999 has a still smaller perigalactic distance, the galactocentric distance varying between extremes of 0.6 and 5.9 kpc.

### Velocity dispersion profile

This is illustrated in Fig.\[fig:vdp\], where it is compared with the observational data of @Pe1995. The shortfall at large radii is of doubtful significance.

### Luminosity Functions

These are shown for our model at the median radius of the innermost and outermost fields observed by @Ri2004. The disagreement in the inner luminosity function at faint magnitudes may be attributable to the fact that the theoretical result assumes 100% completeness, while the observational data are uncorrected for completeness. A plot of the completeness correction in the outer field is given by @Ha2002 [Fig.3] and discussed above in Sec.\[sec:data\_computation\]. The mismatch between model and data is smaller in the outer field, but from the discussion above it would seem that the mismatch in the faintest bin may be too large to be accounted for by the estimated value of the completeness correction. The error bars in the last two observational points on this plot (not shown) almost overlap, and much of the mismatch may be simply sampling uncertainty. It is worth comparing the multi-mass King model constructed by @Ri2004, which is also problematic in the outermost field.

### Core collapse

Fig.\[fig:rc\] shows the evolution with time of the theoretical core radius. There is an early period of very rapid contraction, associated with mass segregation, followed by a slower reexpansion, caused by the loss of mass from the evolving massive stars which are now concentrated within the core.

Our most surprising discovery from our model of M4 is the subsequent behaviour. M4 is classified as a King-profile cluster [@Tr1993], and such clusters are usually interpreted as being clusters whose cores have not yet collapsed. But the plot of the theoretical core radius (Fig.\[fig:rc\] reveals that the model exhibited core collapse at about 8Gyr. Subsequently its core radius is presumably sustained by binary burning. Even non-primordial binaries may be playing a role here. To the best of our knowledge it has not previously been suggested that M4, is a post-collapse cluster, though on statistical grounds @dempp2007 have suggested that some King-type clusters have already collapsed. This issue is often approached by reference to the half-mass relaxation time (Table \[tab:m4dat\], but for M4 this is not short enough to be decisive.

The presence of radial colour gradients is correlated with the presence of a non-King surface brightness profile and hence with core collapse. We have computed the colour profile of our model, and find that it is nearly flat except for the influence of one or two very bright stars within the innermost few arcsec. Most values lie around 0.65 in $B-V$, which is much less than the global value of 1.03 given in @Ha1996 [January 2008]. On the other hand roughly similar mismatches between observations are found in other clusters (e.g. M30, [@Pi1988]).

### The colour-magnitude diagram

Figure \[fig:colour-magnitude\] shows the colour-magnitude diagram of the model. This is of interest, not so much for comparison with observations, but for the presence of a number of interesting features. The division of the lower main sequence is simply an artifact of the way binary masses were selected (a total mass above 0.2$\msun$ and a component mass above 0.1$\msun$.) Of particular interest are the high numbers of merger remnants on the lower white dwarf sequence. There are very few blue stragglers. Partly this is a result of the low binary frequency, but it is also important to note that some formation channels are unrepresented in our models (in particular, collisions during triple or four-body interactions, though if a binary emerges from an interaction with appropriate parameters, it will be treated as merged.) These numbers also depend on the assumed initial distribution of semi-major axis, which is not yet well constrained by observations in globular clusters.

Note in Fig.\[fig:colour-magnitude\] that there are some white dwarf-main sequence binaries (below the main sequence). @Ri2004 drew attention to the possible presence of such binaries in their colour magnitude diagrams of this cluster.

### Binaries

Photometric binaries are visible in Fig.\[fig:colour-magnitude\], and these are compared with observations in the inner field of @Ri2004 in Fig.\[fig:photometrics\]. In this figure, the model histogram has been normalised to the same total number of stars as the observational one. We made no attempt to simulate photometric errors, but the bins around abscissa = -0.75 suggest that the binary fractions in the model (which is under 6% globally; see Table \[tab:mc\_king\]) and the observations are comparable. We note here that @Ri2004, using this same data, concluded that the binary frequency was approximately 2% in the innermost field. But they also note that the measured frequency of “approximately equal-mass binaries” is 2.2%, and we consider that the balance between this number and our binary fraction can be made of binaries with companions whose masses are more unequal.

Radial velocity binaries should also be detectable in M4, and will be part of a separate investigation.

Now we consider the distributions of the binaries, as predicted by the theoretical model. Fig.\[fig:sma\] shows that binaries have evolved dynamically as well as through their internal evolution. In particular the softest pairs been almost destroyed.

By 12 Gyr the binaries exhibit segregation towards the centre of the cluster, but perhaps in more subtle ways than might be expected (Figs.\[fig:segregation\],\[fig:brights\]). When [*all*]{} binaries are considered, there is little segregation relative to the other objects in the system. (Most binaries in our model are of low mass.) But if one restricts attention to bright binaries, which we here take to mean those with $M_V < 7$ (i.e. brighter than about two magnitudes below turnoff), the segregation is very noticeable (Fig.\[fig:brights\]), with a half-mass radius smaller by almost a factor of 2 than for bright single stars. Still, bright binaries are not nearly as mass-segregated as neutron stars (Fig.\[fig:segregation\]), which, incidentally, receive no natal kicks in our model.

The history of the binary fraction is, in effect, given in Fig.\[fig:degenerate-numbers\]. In the first Gyr this falls steadily from the initial value of 0.07 to about 0.06, and it remains close to this value until the present day.

### Escape velocity

We have already referred to the fact that, in our model, neutron stars receive no natal kicks. Escape is of course governed by the escape speed, and this is of much interest in connection with the possibility of retaining gas from the first generation of rapidly evolving stars. For this reason we plot the central escape velocity in Fig. \[fig:vesc\]. The remarkably high value in the first few tens of millions of years, if valid for M4, could have interesting consequences for the early evolution of the cluster and its stars. It draws attention to the very high initial density of our model (Table \[tab:mc\_king\]), whose average value within the half-mass radius is $2\times10^5\msun$/pc$^3$. The central density is about $10^6\msun$/pc$^3$, about an order of magnitude larger than in the central young cluster in the HII region NGC 3603 [@St2006]. In the absence of local young star clusters as massive as our M4 progenitor, it is difficult to be sure whether our initial model is implausibly dense.

### Degenerate components

We have already mentioned the spatial segregation of the population of neutron stars, and the presence of a number of degenerate binaries at the present day. Now we consider the historical evolution of this population over the lifetime of the cluster so far, according to our model.

Fig.\[fig:degenerate-mass\] shows the evolution of these populations. It has already become well established [@VH1997; @HS2003] that white dwarfs account for an increasing proportion of the mass of globular clusters, and indeed it is of order 39% at the present day, according to this model. The proportion of neutron stars is almost certainly excessive, because of our assumption of complete retention.

In order to separate two of the components in this figure (black holes and neutron stars), we show in Fig. \[fig:degenerate-numbers\] the numbers of these degenerate stars, along with the numbers of all stars and all binaries. This shows that no stellar-mass black holes are expected to be present in M4 now.

Conclusions and discussion
==========================

In this paper we have presented a Monte Carlo model for the nearby globular cluster M4. This model includes the effects of two-body relaxation, evaporation across the tidal boundary, dynamical interactions involving primordial and three-body binaries, and the internal evolution of single stars and binaries. By adjustment of the initial parameters (total mass, tidal radius, initial mass function, binary fraction) we have found a model which, after 12Gyr of evolution, leads to a model with a surface brightness profile, velocity dispersion profile, and luminosity functions at two radii, all of which are in tolerable agreement with observational data. It also leads to a number of photometric binaries roughly consistent with observation.

This model has a current mass of $4.6\times10^4\msun$, a $M/L_V$ ratio of 1.8, a binary fraction of almost 6%, and an almost flat lower mass function (Table \[tab:mc\_king\]). Almost 40% of the mass is in white dwarfs, and about 7% in neutron stars, though our model assumes a 100% retention rate. They are strongly mass-segregated to the centre of the cluster. There are no stellar-mass black holes left in the cluster. The binaries have experience significant dynamical evolution, almost all the soft pairs having been destroyed. The binary population as a whole is only slightly segregated towards the centre, but there is more evident segregation of bright binaries (such as those that would be more readily observed in a radial velocity search.)

The most significant new result from our model is the implication that M4 is a core collapse cluster, despite the uncollapsed appearance of the surface brightness profile.

Now we discuss a number of shortcomings of our model and other issues related to this study.

1.  [*Uniqueness:*]{} Though we have arrived at a broadly satisfactory model, it is not at all clear how unique it is. Certainly we were unable to construct models with much larger numbers of primordial binaries, or with significantly larger initial radius, as such models produced an insufficiently concentrated surface brightness profile. Furthermore, the fact that the lower slope of the initial mass function is less steep than the canonical value of 1.3 should not be taken to imply that we have inferred the initial value.

2.  [*Fluctuations*]{}: we have noticed that runs differing only in the initial seed of the random number generator can give surprisingly different surface brightness profiles. In broad terms this confirms the important finding of @Hu2007, though we have not yet established (as he did) that the presence or absence of black hole binaries is the underlying mechanism. We shall return to this issue at appropriate length in the next paper in this series, on the cluster NGC 6397.

3.  [*The Initial Model:*]{} The initial model is astonishingly dense, the central density being of order $10^6\msun$pc$^{-3}$. It also has to be realised that we are imposing initial conditions at a time when the residual gas from the birth of the cluster has already dispersed. Various properties of a star cluster, including its binary population and mass function, may change significantly during the phase of gas expulsion, which further undermines the power of our model to establish the initial conditions. Finally, we have assumed no initial mass segregation, even though recent work suggests that this should be present already before the cluster has assembled into a roughly spherical object [@Ve2007]. This is an aspect of the modelling that could be readily improved.

4.  [*The early evolution of the model:*]{} It is worth drawing attention here to the high initial escape velocity from the centre of our model, in the first few tens of millions of years. Clearly this is dependent on the small initial radius. The initial high density is also responsible for the very rapid initial evolution of the core.

5.  [*Imperfections of the modelling:*]{} Several important improvements need to be made to our technique.

    1.  At present few-body interactions are handled with cross sections. The problem is not simply that these are not well known, especially in the case of unequal masses; it also means that we are unable to determine if a collision occurs during a long-lived interaction. For this reason, we have said almost nothing about collision products (e.g. blue stragglers) in this paper. This limitation could be overcome by direct integration of the interactions, as is done by @FR2007 in their version of the Monte Carlo scheme.

    2.  Long-lived triples are neglected in the model at present. These are commonly produced in binary-binary encounters [@Mi1984a], and it is desirable to include these as a third species (beyond single and binary stars). Their observable effects may be small, but of course there is one intriguing example in the very cluster we have focused on here [@Th1993].

    3.  We assume that the tide is modelled as a tidal cutoff. We have taken some care to ensure that the tidal radius is adjusted so that the model loses mass through escape at the same rate as an $N$-body model would, if immersed in a tidal [*field*]{} with the same tidal radius. This means, however, that the effective tidal radius of the model is somewhat too small. A better treatment of a steady tide may be possible, without this drawback.

    4.  We assume that the tide is steady, something which is not true for M4. The effects of tidal shocks have been studied by a number of authors [e.g. @Ku1995], and it would be possible to add the effects as another process altering the energies and angular momenta of the stars in the simulation. On the other hand @Ag1988 found, on the basis of a simple model, that tidal evaporation was the dominant mechanism in the evolution of the mass of M4 at the present day, and that other factors (such as disk and bulge shocking) contributed at a level less than 1%.

    5.  Rotation: it has been shown [@Ki2004] that, to the extent that rotating and non-rotating models can be compared, rotation somewhat accelerates the rate of core collapse. Rotation is hard to implement in this Monte Carlo model, however.

    6.  The search for initial conditions is still very laborious, and we are constantly seeking ways of expediting this. Our most fruitful technique at present is the use of small-scale models which relax at the same rate as a full-sized model. Not all aspects of the dynamics are faithfully rendered by a scaled model, but the technique appears to be successful as long as the fraction of primordial binaries is not too large.

Despite these caveats and shortcomings, it is clear that the Monte Carlo code we have been developing is now an extraordinarily useful tool for assessing the dynamics of rich star clusters. For the cluster M4 it has led to the conclusion that M4 can be explained as a post-collapse clusters. It also provides a wealth of information on the distribution of the binary population, which will be important for the planning and interpretation of searches for radial velocity binaries, now under way.

$N$-body models will eventually supplant the Monte Carlo technique, but at present are incapable of providing a star-by-star model of even such a small cluster as M4. They can of course provide a great deal of general guidance on the dynamical evolution of rich star clusters, and they underpin the Monte Carlo by providing benchmarks for small models. Even when $N$-body simulations eventually become possible, Monte Carlo models will remain as a quicker way of exploring the main issues, just as King models have continued to dominate the field of star cluster modelling even when more advanced methods (e.g. Fokker-Planck models) have become available.

Discussion
==========

Acknowledgements {#acknowledgements .unnumbered}
================

We are indebted to Jarrod Hurley for much help with the BSE stellar evolution package. We thank Janusz Kaluzny for his kind advice on observational matters, and Harvey Richer for comments on various matters, especially the observed luminosity functions. This research was supported in part by the Polish National Committee for Scientific Research under grant 1 P03D 002 27.

[99]{}

Aguilar L., Hut P., Ostriker J. P., 1988, ApJ, 335, 720

Anderson, J., Bedin, L. R., Piotto, G., Yadav, R. S., & Bellini, A. 2006, , 454, 1029

Baumgardt, H., & Makino, J. 2003, , 340, 227

Bedin, L. R., Piotto, G., King, I. R., & Anderson, J. 2003, , 126, 247

Bedin, L. R., Anderson, J., King, I. R., & Piotto, G. 2001, , 560, L75

De Marchi G., Paresce F., Pulone L., 2007, ApJ, 656, L65

Dinescu, D. I., Girard, T. M., & van Altena, W. F. 1999, , 117, 1792

Fregeau, J. M., & Rasio, F. A. 2007, , 658, 1047

Giersz, M. 1998, , 298, 1239

Giersz, M. 2001, , 324, 218

Giersz, M. 2006, , 371, 484

Giersz, M., & Heggie, D. C. 2003, , 339, 486

Giersz, M., & Heggie, D. C. 2008, preprint

Hansen, B. M. S., et al. 2002, , 574, L155

Hansen, B. M. S., et al. 2004, , 155, 551

Harris, W. E. 1996, , 112, 1487

Hurley J. R., 2007, MNRAS, 379, 93

Hurley J. R., Shara M. M., 2003, ApJ, 589, 179

Hurley, J. R., Pols, O. R., Aarseth, S. J., & Tout, C. A. 2005, , 363, 293

Kim, E., Lee, H. M., & Spurzem, R. 2004, , 351, 220

Kroupa P., Tout C. A., Gilmore G., 1993, MNRAS, 262, 545

Kroupa, P. 1995, , 277, 1507

Kroupa, P. 2007a, ArXiv Astrophysics e-prints, arXiv:astro-ph/0703124

Kroupa P., 2007b, arXiv, 708, arXiv:0708.1164

Kroupa, P., Gilmore, G., & Tout, C. A. 1991, , 251, 293

Kroupa, P., Tout, C. A., & Gilmore, G. 1993, , 262, 545

Kundic, T., & Ostriker, J. P. 1995, , 438, 702

Lamers, H. J. G. L. M., Gieles, M., Bastian, N., Baumgardt, H., Kharchenko, N. V., & Portegies Zwart, S. 2005, , 441, 117

McMillan S. L. W., Vesperini E., Portegies Zwart S. F., 2007, ApJ, 655, L45

Mikkola, S. 1984, , 207, 115

Peterson, R. C., Rees, R. F., & Cudworth, K. M. 1995, , 443, 124

Piotto G., King I. R., Djorgovski S., 1988, AJ, 96, 1918

Richer H. B., et al., 1995, ApJ, 451, L17

Richer H. B., et al., 2002, ApJ, 574, L151

Richer, H. B., et al.  2004, , 127, 2771

Sommariva, V., et al, in E. Vesperini, M. Giersz, A. Sills, eds, IAUS 246, in press

Stolte A., Brandner W., Brandl B., Zinnecker H., 2006, AJ, 132, 253

Thorsett, S. E., Arzoumanian, Z., & Taylor, J. H. 1993, , 412, L33

Trager, S. C., Djorgovski, S., & King, I. R. 1993, in Djorgovski, S.G., Meylan G., eds, Structure and Dynamics of Globular Clusters, ASPCS 50, 347

Trager, S. C., King, I. R., & Djorgovski, S. 1995, , 109, 218

van den Bergh S., 1995, AJ, 110, 1171

Vesperini E., Heggie D. C., 1997, MNRAS, 289, 898

Wilkinson M. I., Hurley J. R., Mackey A. D., Gilmore G. F., Tout C. A., 2003, MNRAS, 343, 1025

\[lastpage\]

[^1]: E-mail: mig@camk.edu.pl (MG); d.c.heggie@ed.ac.uk (DCH)
