Table Of Content

Mon.Not.R.Astron.Soc.000,1–10(2002) Printed5February2008 (MNLATEXstylefilev2.2) Monte Carlo Simulations of Star Clusters - V. The globular cluster M4 Mirek Giersz1⋆ and Douglas C. Heggie2 8 1Nicolaus Copernicus Astronomical Centre, Polish Academy of Sciences, ul. Bartycka18, 00-716 Warsaw, Poland 0 2Universityof Edinburgh, School of Mathematics and Maxwell Institute for Mathematical Sciences, King’sBuildings, Edinburgh EH9 3JZ, UK 0 2 n Accepted ....Received...;inoriginalform... a J 5 ABSTRACT 2 We describe Monte Carlo models for the dynamical evolution of the nearby globular clusterM4.Thecodeincludestreatmentsoftwo-bodyrelaxation,three-andfour-body ] h interactionsinvolving primordialbinaries and those formed dynamically,the Galactic p 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 o a satisfactory match to the surface brightness profile, the velocity dispersion profile, r t and the luminosity function in two fields. We describe in particular the evolution s of the core, and find that M4 (which has a classic King profile) is actually a post- a [ collapsecluster,itscoreradiusbeingsustainedbybinaryburning.Wealsoconsiderthe distribution of its binaries, including those which would be observed as photometric 1 binaries and as radial-velocity binaries. We also consider the populations of white v dwarfs,neutronstars,blackholesandbluestragglers,thoughnotallchannelsforblue 8 straggler formation are represented yet in our simulations. 6 9 Key words: stellar dynamics – methods: numerical – binaries: general – globular 3 clusters: individual: M4 . 1 0 8 0 1 INTRODUCTION are concerned with issues where the dynamical history of : v the star cluster is important, where static models are un- The present paper opens up a new road in the study of i informative. Secondly, our aim is to construct a model of a X the dynamical evolution of globular clusters. We adopt the specific star cluster, rather than trying to understand the r Monte Carlo method of Giersz (Giersz 1998, 2001, 2006), general properties of the evolution of a population of star a which in recent years has been enhanced to deal quitereal- clusters. This has been done before, and a brief history is isticallywiththestellarevolutionofsingleandbinarystars, outlined in Giersz & Heggie (2008), but the present work to study the dynamical history of thenearby globular clus- takes these efforts onto a new level of realism, in terms of terM4.Anearlierversionofthecodehadalreadybeenused the description of stellar evolution, and dynamical interac- to study the dynamical history of ω Cen (Giersz & Heggie tions involvingbinary stars. 2003), but at that time the treatment of stellar evolution This problem is not easy. Not only is it necessary to wasprimitiveandtherewerenobinaries.Thenewcodehas use an elaborate technique for simulating the relevant as- been thoroughly tested on smaller systems, by comparison trophysical processes, but it is necessary also to search for with N-body simulations and observations of the old open initialconditionswhich,afterabout12Gyrofevolution,lead star cluster M67 (Giersz & Heggie 2008). There we showed toanobjectresemblingagivenstarcluster.By“resembling” that the Monte Carlo code could produce data of a similar we do not simply mean matching the overall mass, radius levelofdetailandrealismasthebestN-bodycodes.Nowfor andbinaryfractionofacluster,forexample,fortworeasons: the first time we consider much richer systems, with about halfamillionstarsinitially,whichareatpresentbeyondthe reach of N-body methods. (i) We have found that values found for data in the lit- This paper has a place within a long tradition of the erature are highly uncertain, and different sources are con- modelling of globular star clusters, but the place is a dis- tradictory. These data are usually derived, in some model- tinctiveone.First,wearenotconcernedwithastaticmodel dependentway,fromsuchdataassurface-brightnessprofiles of a star cluster at the present day, like a King model. We and velocity dispersion profiles, and we prefer to compare our models directly with this data, and not with inferred global parameters. ⋆ E-mail:[email protected](MG);[email protected](DCH) (ii) We have found that, even if one achieves a satisfac- 2 M. Giersz and D.C. Heggie tory fit to these profiles, the model may give a very poor ical star cluster is represented by its mass, energy and an- comparison with theluminosity function. gular momentum, and its stellar evolutionary state may be computed at any time using synthetic formulae for single From these considerations we conclude that a model which and binary evolution. It may be a binary or a special kind aims to fit only the mass and radius of a star cluster (say) of single star that has been created in a collision or merger may be very far from thetruth. event. Tackling this difficult problem is not just interesting, Neighbouring stars interact with each other in accor- however. We have been motivated by a number of press- dance (in a statistical sense) with the theory of two-body ing astrophysical problems. For example, the two nearby relaxation.Ifoneorbothoftheparticipantsisabinary,the star clusters M4 (the subject of this study) and NGC 6397 probability of an encounter affecting the internal dynamics (which we shall consider in our next paper in this series) iscalculatedaccordingtoanalyticcrosssections,whichalso haverathersimilarmassandradius,andyetonehasaclas- determine the outcome. This is one of the main shortcom- sicKingprofile,whiletheotherisawell-studiedexampleofa ingsofthecode,asthesecrosssectionsarenotwellknownin clusterwith a“collapsed core” (Trager et al. 1995).Among thecaseofunequalmasses,andalsothepossibilityofstellar possible explanations one may consider differences in the collisions duringlong-lived temporary captureis excluded. populationofbinaries,whichareknowntoaffectcoreprop- A star or binary may escape if its energy exceeds a erties,orintidaleffects.Indeedthepresentpaperwillshow certain value, which we choose to be lower than the energy thatthesetwo clustersmay bemuchmoresimilar thanone atthenominaltidalradius,inordertoimprovethescalingof would suppose from the surface brightness profiles alone. thelifetimewithN,asexplainedinGiersz & Heggie(2008). Asecondmotivationforourworkisourinvolvementin This is the second main shortcoming of the models, as it observational programmes aimed at characterising the bi- leadstoacutoffradiusofthemodelthatissmallerthanthe narypopulationsinglobularclusters.Whatdifferences(e.g. truetidal radius, and this lowers the surface density profile in the distributions of periods and abundances) should one in the outerparts of thesystem. expect to find between the core and the halo? These issues A difficulty in applyingtheMonte Carlo code toM4 is are important in the planning of observations, and in their that it employs a static tide, whereas the orbit of M4 ap- interpretation. pearstobeveryelliptical (Dinescu et al. 1999).Wehaveto The cluster M4 is the focus of much of this effort be- assumethataclustercanbeplacedinasteadytideofsucha cause it is nearby, making it a relatively easy target for strengththattheclusterlosesmassatthesameaveragerate deepobservationalstudy.Itwasthefirstglobularclusterto as it would on its true orbit. Some support for this proce- yield a deep sequence of white dwarfs (Richeret al. 1995). dure comes from N-body modelling. Baumgardt & Makino More recently it has been subjected to an intensive ob- (2003) show that clusters on an elliptical orbit between servational programme by the Padova group (Bedin et al. about 2.8 and 8.5kpc dissolve on a time scale intermedi- 2001, 2003; Anderson et al. 2006), which includes searches ate between that for circular orbits at these two radii, and for radial-velocity binaries in the upper main sequence thatthedissolutiontimescalesinalmostthesamewaywith (Sommariva et al 2008). It also turns out to be a cluster thesizeofthesystem.Wilkinson et al.(2003)showthatthe which (we conclude) started with only about half a million core radius of a cluster on an elliptic orbit evolves in very stars, which facilitates the modelling. Along with the open nearly the same way as in a cluster with a circular orbit at clusterM67, M4waschosen bytheinternationalMODEST thetime-averaged galactocentric distance. consortium, at a meeting in Hamilton in 2005, as the focus Allotherfreeparametersofthecode(e.g.thecoefficient for joint effort by theorists and observers, to cast light on of N in the Coulomb logarithm) take the optimal values its binary population and dynamical properties. M67 has found in theabove study. beenmodelled verysuccessfully byHurley et al. (2005),us- ing N-body techniques, and this paper represent the first theoretical step in a similar study of M4. 2.2 Initial Conditions The paper is organised as follows. First, we summarise features of the code and the models, the data we used, and TheinitialmodelsareasspecifiedinTable1.Manyofthese our approach to the problem of finding initial conditions features (e.g. theproperties of thebinaries, exceptfor their for M4. Then we present data for our best models: surface overallabundance)wereinheritedfromourmodellingofthe brightness and velocity dispersion profiles, luminosity func- openclusterM67,andthosewereinturnmainlydrawnfrom tions, theproperties ofthebinarypopulation, whitedwarfs the work of Hurley et al. (2005). Some of the parameters and other degenerate remnants, and the inferred dynami- were taken to be freely adjustable, and this freedom was cal state of the cluster. The final section summarises our exploited in the search for an acceptable fit to the current findingsanddiscussesthemin thecontextofwork onother observational data. clusters, including objects to which we will turn in future papers. 2.3 Observational data and its computation Our first task was to iterate on the initial parameters of 2 METHODS ourmodelinordertoproduceasatisfactoryfittoarangeof 2.1 The Monte Carlo Code observationaldataatanageof12Gyr.Thedataweadopted are as follows: The details of our simulation method have been amply de- scribedinpreviouspapersinthisseries.Eachstarinaspher- (i) Surface brightness profile: here we used the compila- Monte Carlo Simulations of Star Clusters - V. The globular cluster M4 3 Table 1.InitialparametersforM4 Table 2.PropertiesofM4 Distancefromsun a1.72kpc Fixedparameters DistancefromGC 5.9kpc Mass a63000M⊙ Structure Plummermodel Initialmassfunction1 Kroupa,Tout,&Gilmore(1993)inthe Coreradius 0.53pc Half-lightradius 2.3pc range[0.1,50]M⊙ Tidalradius 21pc Binarymassdistribution Kroupaetal.(1991) Binarymassratio Uniform(withcomponent Half-massrelaxationtime(Rh) 660Myr Binaryfraction a1-15% massesrestrictedasfor [Fe/H] -1.2 singlestars) Age b12Gyr Binarysemi-majoraxis Uton5if0oArmUinlog,2(R1+R2) AV a1.33 Binaryeccentricity Thermal,witheigenevolution (Kroupa1995) References:Alldataarefromthecurrentversionofthecatalogue MetallicityZ 0.002 of Harris (1996), except a Richeretal. (2004) (though this is Age 12Gyr(Hansenetal.2004) not always the original reference for the quoted number) and b Hansenetal.(2004). Freeparameters Mass M The line-of-sight velocity dispersion is computed in a Tidalradius rt similar way. Foreach particle we calculate themean square Half-massradius rh line-of-sight velocity (because the orientation of the trans- Binaryfraction fb verse component of the velocity is random), and then sum Slopeofthelowermass overallparticles.Inthissum,eachparticleisweightedbya function α(Kroupa=1.3) geometricalfactorproportionaltothesurfacedensityofthe particle’s shell along the line of sight. The result is an esti- mate of the velocity dispersion that is weighted by neither tion byTrager et al. (1995),where thesurface brightnessis mass nor brightness, but only number. To check the effects expressed in V magnitudes persquare arcsec. of this, we have sometimes calculated velocity dispersions (ii) Radial velocity profile: this came from withvariouscutoffsinthemagnitudesofthestarsincluded. (Peterson et al. 1995), and is the result of binning The effects are usually quitesmall. datatheontheradialvelocities ofnearly200stars.Strictly Computation of the luminosity function is the least we should refer to the line-of-sight velocity dispersion, as problematic. We count stars in each bin in V-magnitude, “radial” velocity has a different meaning in the Monte butlyingabovealineinthecolour-magnitudediagramjust Carlo model. below the main sequence. Again the contribution of each (iii) The V-luminosity function (Richeret al. 2004, from star is weighted by the same geometrical factor, and the V whichweconsideredtheresultsfortheinnermostandouter- magnitudeis corrected for extinction. mostoftheirfourannuli).Thesedataarelackcorrectionfor completeness,thoughthecompletenessfactorfortheouter- most field is plotted in Richeret al. (2002). For main se- 2.4 Finding initial conditions quence stars it is almost 100% down to V = 15, and drops Herewesummariseourexperienceinapproachingthisprob- steadilytolessthan50%atV =17.Fortheinnermostfield lem.Anumberofstudies(e.g.Baumgardt & Makino(2003), thecompleteness correction would belarger. Lamers et al.(2005))givesimpleformulaefortheevolution Now we consider how to compare this data with the with time of the bound mass of a rich star cluster. It is outputoftheMonteCarlocode.Thisincludesalistofeach possible to derive similar simple formulae for the evolution particleinthesimulation, alongwithdataonitsradius,ra- of thehalf-mass radiusand otherquantities. Thereare sev- dial and transverse velocities and absolute magnitude, and eral problems with inverting these formulae, however, i.e. numerousotherquantities. Toconstruct thesurface bright- using them to infer the initial parameters of a cluster from ness, we think of each particle as representing a luminous its present mass and radius. First, these present-day global shell of the same radius, and superpose the surface bright- parameters are quite uncertain, even to within a factor of ness of all shells. While this procedure involves a minimum two.Second,theseformulaedependonthegalacticorbitand of effects from binning, or randomly assigning the full po- other parameters which are equally subject to uncertainty. sition of the star, a shell has an infinite surface brightness Therefore we have adopted the more straightforward but at its projected edge. The effect of this, especially from the more laborious approach of iterating on the initial param- brightest stars, will be seen in some of the profiles to be eters of our models (Tab 1, lower half); that is, we select presented in this paper. A similar problem arises in actual values for the five stated parameters, run the model, find observations, and is often handled by simply removing the where the match with the observations is poor, adjust the brightest stars from the surface brightness data presented. parameters, and repeat cyclically. Theoutputfromthemodeliscorrectedforextinction(Table We have employed two methods to facilitate this pro- 2. The distance of the line of sight from the cluster centre cess to some extent. First, we have often carried out mini- is converted between pc(as in themodel) and arcsec (asin surveys,i.e.small,coarsegrid-searchesaroundagivenstart- theobservations) using the distance in the same table. ingmodel,tofindouthowchangesinindividualparameters 4 M. Giersz and D.C. Heggie affect the results. Second, we have used scaling to acceler- Second,thisexperienceunderscorestheimportanceofprop- atetheprocess, andwenowdescribethismethodinalittle erly normalised luminosity functions. In other words, it is detail. importanttoknowtheareaofthefieldwherethestarshave Suppose we wish to represent a star cluster which has beencounted,orsomeequivalentrepresentationofproperly mass M and radius R with a model representing a cluster normalised data. Very often, the emphasis is solely on the witha(usuallysmaller)massM∗andradiusR∗.Sincetwo- shape of the luminosity function, but we have found that body relaxation dominates muchof thedynamics,we insist theabsolute normalisation is an essential constraint. that thetwo clustershavethesamerelaxation time, and so Baumgardt & Makino(2003)showthatthelowermass function becomes flatter as the fraction of mass lost by the R∗ M 1/3 logγN∗ 2/3 cluster increases (i.e. towards the end of its life). There- = , R (cid:16)M∗(cid:17) (cid:18) logγN (cid:19) fore we could perhaps have improved the fit with the lumi- nosity function by starting with a more massive model and where N,N∗ are the corresponding particle numbers. Thus somehow ensuring a larger escape rate so as toleave a sim- themodeloflowermasshaslargerradius.Thetidalradiusis ilar mass at the present day. Instead, we elected to change scaledinthesameway.Thentheobservationalresults(sur- the slope of the low-mass IMF from the canonical value of facebrightnessprofile,etc)canbecomputedforthemodelof α = 1.3 (Kroupa 2007a) to α = 0.9. (There is some jus- lower mass andthenrescaled (byappropriate factors of the tification for a lower value for low-metallicity populations, massandradius) togivearesultforthemoremassiveclus- though it has been argued (Kroupa2007b) that there is no ter (assuming that the evolution is dominated by the pro- pristine low-metallicity population where the IMF can be cesses of relaxation, stellar evolution and tidal stripping). inferred securely.) 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)initialmass,certainlywhentheproportionof 3.2 A Monte Carlo model of M4 binaries is 10% or less. Some aspects of the evolution are Bysomeexperimentationwearrivedatamodelwhichgave not well described by this scaling technique. For example, a fair fit to all three kinds of observational data; see Ta- wedonotchangethedistributionofthesemi-majoraxesof ble 3 (where we compare with a King model developed by the binaries. In this way the internal evolution of the bina- Richeret al.(2004),andFigs1-4.Itisworthnotingthatno riesiscorrectlymodelled(providedthatthebinariesremain arbitrary normalisation has been applied in these compar- isolated dynamically), though their dynamical interactions isons between our model and the observations. The surface with the rest of the system are not. In principle one could brightness profile, for example, is computed directly from scale the semi-major axes in the same way as R, but then the V-magnitudes of the stars in the Monte Carlo simula- theinternal evolution of the binaries would be altered. tion, as described in Sec.2.3. In the construction of a King model,bycontrast,itisoftenassumedthatthemass-to-light ratio is arbitrary. 3 MODELS OF M4 3.1 Finding initial parameters 3.2.1 Surface brightness Our starting point was our work on the old open cluster M67 (Giersz & Heggie 2008), but with a larger initial mass While the overall surface brightness profile is slightly faint, andradius.(Baumgardt & Makino(2003),forexample,sug- the most noticeable feature of Fig.1 is that the model has gestedthattheinitialmassofM4(NGC6121) wasoforder a somewhat smaller limiting radius than the observational 7.5×105M⊙,thoughweusedsomewhatsmallervalues.)To data. The reason for this is explained in Giersz & Heggie begin with, our choice of initial tidal radius was inferred (2008): in short, we impose a smaller tidal radius than the from theinitialmassandthepresent-dayestimates ofmass nominal tidal radius, in an N-dependent way, which is in- andtidalradiusgiveninTable2,assumingthatr ∝M1/3. tended to ensure that the overall rate of escape from the t At first we adopted similar values for the “concentration” modelbehavesinthesamewayasinanN-bodysimulation. (r /r )and binaryfraction as in themodelling of M67, but Another point to notice is the disagreement between the t h found thatthesurface brightnessprofilefittedpoorly (with total luminosity of our model and that of the King model toolargeacore)unlessthebinaryfractionwasmuchsmaller quoted in the final column of Table 3. Trager et al. (1995) (5 to 10%) and the concentration muchhigher. givean analytic fittothesurface brightnessprofile, and we Before describing our best models, it is worth briefly havecheckedthat theintegrated value is close to ours. mentioningonewhich providedasatisfactory fittothesur- The data for the Monte Carlo model in Table 3 would face brightness profiles. The fit to the velocity dispersion be consistent with a galactocentric radius of about 1.7kpc, profile was tolerable, but indicated a model that was too in an isothermal galaxy model with a circular velocity of massivebyafactorofabout1.4.Itsmainflaw,however,was 220km/s. While this is certainly much smaller than its cur- intheluminosityfunction,whichwasgenerallytoolargeby rent galactocentric distance, a small value was also found a factor in the range 2–3. There are two reasons why this (usingasimilar argument and publishedvalues of themass is interesting. One is that, for a long time, models of star and tidal radius) by van den Bergh (1995). The orbit given clusterswereconstructedentirelyonthebasisofthesurface by Dinescu et al. (1999) has a still smaller perigalactic dis- brightness and velocity dispersion profiles. It should be re- tance,thegalactocentricdistancevaryingbetweenextremes alised that such models may be misleading in other ways. of 0.6 and 5.9 kpc. Monte Carlo Simulations of Star Clusters - V. The globular cluster M4 5 Table 3.MonteCarloandKingmodelsforM4 Quantity MCmodel MCmodel Kingmodel (t=0) (t=12Gyr) (Richeretal.2004) Mass(M⊙) 3.40×105 4.61×104 Luminosity(L⊙) 6.1×106 2.55×104 6.25×104 Binaryfractionfb 0.07 0.057 0 Low-massMFslopeα 0.9 0.03 0.1 Massofwhitedwarfs(M⊙) 0 1.81×104 3.25×104∗ Massofneutronstars(M⊙) 0 3.24×103 Tidalradiusrt (pc) 35.0 18.0 Half-massradiusrh (pc) 0.58 2.89 ∗:thisisthequotedmassof“degenerates” Surface brightness Innermost annulus 16 100 Trager et al MC Model MC Model Richer et al 18 20 10 -2ec arcs 22 N g a m V- 24 1 26 28 0.1 0.1 1 10 100 1000 10000 7 8 9 10 11 12 13 14 15 radius (arcsec) MV Figure 1. Surface brightness profile of our Monte Carlo Figure 3.LuminosityfunctionofourMonteCarlomodelat model,comparedwiththedataofTrageretal.(1995). the median radius of the innermost annulus in Richeretal. (2004),comparedwiththeirdata. Velocity dispersion Outermost annulus MC Model 100 Peterson et al MC Model 5 Richer et al 4 ec) m/s 3 V (kz N 10 S M R 2 1 0 0.01 0.1 1 10 100 1 7 8 9 10 11 12 13 14 15 radius (arcmin) MV Figure 2. Velocity dispersion profile of our Monte Carlo Figure 4.LuminosityfunctionofourMonteCarlomodelat model,comparedwiththedataofPetersonetal.(1995). the median radius of the outermost annulus in Richeretal. (2004),comparedwiththeirdata. 6 M. Giersz and D.C. Heggie 3.2.2 Velocity dispersion profile This is illustrated in Fig.2, where it is compared with the observationaldataofPeterson et al.(1995).Theshortfallat large radii is of doubtfulsignificance. 3.2.3 Luminosity Functions These are shown for our model at the median radius of theinnermostandoutermostfieldsobservedbyRicher et al. (2004). The disagreement in the inner luminosity function at faint magnitudes may be attributable to the fact that thetheoreticalresultassumes100%completeness,whilethe observationaldataareuncorrectedforcompleteness.Aplot of the completeness correction in the outer field is given by Hansen et al. (2002, Fig.3) and discussed above in Sec.2.3. 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 Figure 5.TheoreticalcoreradiusofourMonteCarlomodel. be accounted for by the estimated value of the complete- nesscorrection.Theerrorbarsinthelasttwoobservational points on this plot (not shown) almost overlap, and much 3.2.5 The colour-magnitude diagram of the mismatch may be simply sampling uncertainty. It is Figure6showsthecolour-magnitudediagramofthemodel. worthcomparingthemulti-massKingmodelconstructedby This is of interest, not so much for comparison with ob- Richeret al. (2004),which isalso problematic in theouter- servations, but for the presence of a number of interesting most field. features. The division of the lower main sequence is simply an artifact of the way binary masses were selected (a total massabove0.2M⊙andacomponentmassabove0.1M⊙.)Of particularinterestarethehighnumbersofmergerremnants 3.2.4 Core collapse on thelower whitedwarf sequence.Thereare veryfewblue stragglers.Partlythisisaresultofthelowbinaryfrequency, Fig.5 shows the evolution with time of the theoretical core but it is also important to note that some formation chan- radius. There is an early period of very rapid contraction, nelsareunrepresentedinourmodels(inparticular,collisions associated withmasssegregation,followed byaslowerreex- during triple or four-body interactions, though if a binary pansion, caused by the loss of mass from the evolving mas- emergesfromaninteractionwithappropriateparameters,it sive stars which are now concentrated within thecore. will be treated as merged.) These numbers also depend on Ourmost surprisingdiscoveryfrom ourmodelofM4is theassumedinitialdistributionofsemi-majoraxis,whichis thesubsequent behaviour. M4 is classified as a King-profile notyetwellconstrainedbyobservationsinglobularclusters. cluster (Trager et al. 1993), and such clusters are usually Note in Fig.6 that there are some white dwarf-main interpreted as being clusters whose cores have not yet col- sequence binaries (below the main sequence). Richer et al. lapsed.Buttheplotofthetheoretical coreradius(Fig.5 re- (2004)drewattentiontothepossiblepresenceofsuchbina- vealsthatthemodelexhibitedcorecollapse atabout 8Gyr. ries in theircolour magnitude diagrams of this cluster. Subsequentlyits core radius is presumably sustained by bi- naryburning.Evennon-primordialbinariesmaybeplaying a role here. To the best of our knowledge it has not pre- 3.2.6 Binaries viously been suggested that M4, is a post-collapse cluster, thoughonstatisticalgroundsDe Marchi, Paresce, & Pulone PhotometricbinariesarevisibleinFig.6,andthesearecom- (2007)havesuggestedthatsomeKing-typeclustershaveal- pared with observations in the inner field of Richer et al. readycollapsed. Thisissueisoften approachedbyreference (2004)inFig.7.Inthisfigure,themodelhistogramhasbeen tothehalf-massrelaxationtime(Table2,butforM4thisis normalised to the same total number of stars as the obser- not short enough to bedecisive. vationalone. Wemadenoattempt tosimulate photometric The presence of radial colour gradients is correlated errors,butthebinsaroundabscissa=-0.75suggestthatthe with the presence of a non-King surface brightness profile binary fractions in the model (which is under 6% globally; andhencewithcorecollapse. Wehavecomputedthecolour see Table 3) and the observations are comparable. We note profile of our model, and find that it is nearly flat except here that Richeret al. (2004), using this same data, con- for the influence of one or two very bright stars within the cluded that the binary frequency was approximately 2% in innermostfewarcsec.Mostvaluesliearound0.65inB−V, the innermost field. But they also note that the measured which is much less than the global value of 1.03 given in frequency of “approximately equal-mass binaries” is 2.2%, Harris (1996, January 2008). On the other hand roughly and we consider that the balance between this number and similarmismatchesbetweenobservationsarefoundinother ourbinaryfractioncanbemadeofbinarieswithcompanions clusters (e.g. M30, (Piotto, King, & Djorgovski 1988)). whose masses are more unequal. Monte Carlo Simulations of Star Clusters - V. The globular cluster M4 7 Binary sequence 60 CM Diagram RicMheCr eMt oadl einl,n neor ramnanliusleuds 50 0 40 N 30 Binaries 5 CollisioSnin pgrleo dsutacrtss 20 10 V 10 0-1.5 -1 -0.5 0 0.5 1 1.5 V offset from single star MS Figure 7. Histogram of V-offset from the main sequence, compared with the corresponding data from the innermost annulusstudiedbyRicheretal.(2004) Seetextfordetails. 15 Distribution of semi-major axes 1 0 Gyr 12 Gyr 0.9 0.8 20 -0.5 0 0.5 1 1.5 2 0.7 B-V 0.6 sFtiagrus;rreed6:.bTihnearcioelso;ubrl-umeapglnuisteusd:ecdolilaigsiroanmoartm12erGgyerr.rGemrenenan:tssin.gle Fraction 0.5 0.4 0.3 RadialvelocitybinariesshouldalsobedetectableinM4, 0.2 and will be part of a separate investigation. 0.1 Now we consider the distributions of the binaries, as 0 0.1 1 10 100 1000 10000 100000 predicted by the theoretical model. Fig.8 shows that bina- Semi-major axis (solar radii) ries have evolved dynamically as well as through their in- Figure8.Distributionfunctionofthesemi-majoraxesofthe ternal evolution. In particular the softest pairs been almost binariesat0Gyrand12Gyr.Units:solarradii. destroyed. By 12Gyrthebinariesexhibit segregation towards the centreof thecluster, butperhapsin more subtleways than value in the first few tens of millions of years, if valid for mightbeexpected(Figs.9,10).Whenallbinariesareconsid- M4, could have interesting consequences for the early evo- ered, there is little segregation relative to the other objects lution of the cluster and its stars. It draws attention to the inthesystem.(Mostbinariesinourmodelareoflowmass.) very high initial density of our model (Table 3), whose av- But if one restricts attention to bright binaries, which we erage value within the half-mass radius is 2×105M⊙/pc3. here take to mean those with MV < 7 (i.e. brighter than Thecentraldensityis about 106M⊙/pc3,about an orderof abouttwomagnitudesbelowturnoff),thesegregationisvery magnitude larger than in the central young cluster in the noticeable (Fig.10), with a half-mass radius smaller by al- HIIregion NGC3603 (Stolte et al. 2006).Intheabsenceof most a factor of 2 than for bright single stars. Still, bright localyoungstarclustersasmassiveasourM4progenitor,it binaries are not nearly as mass-segregated as neutron stars isdifficulttobesurewhetherourinitialmodelisimplausibly (Fig.9), which, incidentally, receive no natal kicks in our dense. model. The history of the binary fraction is, in effect, given in Fig.13. In thefirst Gyr this falls steadily from theinitial 3.2.8 Degenerate components valueof0.07toabout0.06,anditremainsclosetothisvalue untilthe present day. 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 3.2.7 Escape velocity the historical evolution of this population over the lifetime We have already referred to the fact that, in our model, of thecluster so far, according to our model. neutronstarsreceivenonatalkicks.Escapeisofcoursegov- Fig.12 shows the evolution of these populations. It has erned by the escape speed, and this is of much interest in already become well established (Vesperini& Heggie 1997; connectionwiththepossibilityofretaininggasfromthefirst Hurley & Shara 2003) that white dwarfs account for an in- generation ofrapidly evolvingstars.Forthisreason weplot creasingproportionofthemassofglobularclusters,andin- the central escape velocity in Fig. 11. The remarkably high deedit isoforder39% at thepresent day,accordingtothis 8 M. Giersz and D.C. Heggie Radial distribution 1 All objects Binaries 0.9 Neutron stars White dwarfs Collision products 0.8 0.7 0.6 Fraction 0.5 0.4 0.3 0.2 0.1 0 0.1 1 10 Radius (pc) Figure9.Radialdistributionfunctionsat12Gyr.Units:par- secs. Thekey identifies the class of object included. The dis- Figure12.Massofwhitedwarfs,andofneutronstarsandblack tributionsofwhitedwarfsandbinariesarealmostidentical. holestogether. Radial distribution of bright objects 1 All objects All bright objects (V<7) 0.9 Bright binaries 0.8 0.7 0.6 Fraction 0.5 0.4 0.3 0.2 0.1 0 0.1 1 10 Radius (pc) Figure 10. Radial distribution functions at 12Gyr, showing the extent of segregation between bright single and bright binary stars (MV < 7). Also shown for comparison is the Figure13.Numbersofallstars,whitedwarfs,neutronstarsand distributionforallobjects, asinFig.9.Units:parsecs. blackholes. model. The proportion of neutron stars is almost certainly excessive, because of ourassumption of complete retention. In order to separate two of the components in this fig- ure (black holes and neutron stars), we show in Fig. 13 the numbers of these degenerate stars, along with the numbers ofallstarsand allbinaries.Thisshows thatnostellar-mass black holes are expected to bepresent in M4 now. 4 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-bodybinaries,andtheinternalevolutionofsinglestars andbinaries.Byadjustmentoftheinitial parameters(total mass,tidalradius,initialmassfunction,binaryfraction)we havefoundamodelwhich,after12Gyrofevolution,leadsto amodelwithasurfacebrightnessprofile,velocitydispersion profile, and luminosity functions at two radii, all of which Figure 11.Centralescapevelocityasafunctionoftime. are in tolerable agreement with observational data. It also leadstoanumberofphotometricbinariesroughlyconsistent with observation. Monte Carlo Simulations of Star Clusters - V. The globular cluster M4 9 Thismodelhasacurrentmassof4.6×104M⊙,aM/LV italsomeansthatweareunabletodetermineifacollision ratio of 1.8, a binary fraction of almost 6%, and an almost occursduringalong-livedinteraction.Forthisreason,we flat lower mass function (Table3).Almost 40% of themass have said almost nothing about collision products (e.g. is in white dwarfs, and about 7% in neutron stars, though blue stragglers) in this paper. This limitation could be ourmodelassumesa100%retentionrate.Theyarestrongly overcome by direct integration of the interactions, as is mass-segregated to the centre of the cluster. There are no done by Fregeau & Rasio (2007) in their version of the stellar-massblackholesleftinthecluster.Thebinarieshave MonteCarlo scheme. experience significant dynamical evolution, almost all the (b) Long-lived triples are neglected in the model at soft pairs having been destroyed. The binary population as present. These are commonly produced in binary-binary a whole is only slightly segregated towards the centre, but encounters (Mikkola 1984), and it is desirable to include there is more evident segregation of bright binaries (such these as a third species (beyond single and binary stars). as those that would be more readily observed in a radial Theirobservableeffectsmaybesmall,butofcoursethere velocity search.) is one intriguing example in the very cluster we have fo- The most significant new result from our model is the cused on here(Thorsett et al. 1993). implication that M4 is a core collapse cluster, despite the (c) Weassumethatthetideismodelled asatidalcut- uncollapsed appearance of thesurface brightness profile. off. We have taken some care to ensure that the tidal Nowwediscussanumberofshortcomingsofourmodel radius is adjusted so that the model loses mass through and other issues related to thisstudy. escape at the same rate as an N-body model would, if immersed in a tidal field with the same tidal radius. (i) Uniqueness:Thoughwehavearrivedatabroadlysat- Thismeans,however,thattheeffectivetidalradiusofthe isfactory model, it is not at all clear how unique it is. Cer- model is somewhat too small. A better treatment of a tainlywewereunabletoconstructmodelswithmuchlarger steady tide may be possible, without thisdrawback. numbers of primordial binaries, or with significantly larger (d) Weassumethatthetideissteady,somethingwhich initialradius,assuchmodelsproducedaninsufficientlycon- is not true for M4. The effects of tidal shocks have been centrated surface brightness profile. Furthermore, the fact studied by a number of authors (e.g. Kundic& Ostriker thatthelowerslopeoftheinitialmassfunctionislesssteep 1995), and it would be possible to add the effects as thanthecanonicalvalueof1.3shouldnotbetakentoimply another process altering the energies and angular mo- that we haveinferred theinitial value. menta of the stars in the simulation. On the other hand (ii) Fluctuations:wehavenoticedthatrunsdifferingonly Aguilar, Hut,& Ostriker (1988) found, on the basis of a intheinitialseed oftherandomnumbergenerator cangive simple model, that tidal evaporation was the dominant surprisingly different surface brightness profiles. In broad mechanism in the evolution of the mass of M4 at the termsthis confirmstheimportant findingof Hurley (2007), presentday,andthatotherfactors(suchasdiskandbulge though we have not yet established (as he did) that the shocking) contributed at a level less than 1%. presence or absence of black hole binaries is the underly- (e) Rotation:ithasbeenshown(Kim et al.2004)that, ingmechanism.Weshallreturntothisissueat appropriate to the extent that rotating and non-rotating models can length in the next paper in this series, on the cluster NGC be compared, rotation somewhat accelerates the rate of 6397. corecollapse.RotationishardtoimplementinthisMonte (iii) TheInitialModel:Theinitialmodelisastonishingly dense, the central density being of order 106M⊙pc−3. It Carlo model, however. (f) The search for initial conditions is still very labo- also has to be realised that we are imposing initial con- rious, and we are constantly seeking ways of expediting ditions at a time when the residual gas from the birth this. Our most fruitful technique at present is the use of of the cluster has already dispersed. Various properties of small-scale models which relax at thesame rateas afull- a star cluster, including its binary population and mass sizedmodel.Notallaspectsofthedynamicsarefaithfully function, may change significantly during the phase of renderedbyascaled model,butthetechniqueappearsto gas expulsion, which further undermines the power of our besuccessfulaslongasthefractionofprimordialbinaries model to establish the initial conditions. Finally, we have is not too large. assumed no initial mass segregation, even though recent work suggests that this should be present already before Despitethesecaveatsandshortcomings, itisclearthat the cluster has assembled into a roughly spherical object the Monte Carlo code we have been developing is now an (McMillan, Vesperini, & Portegies Zwart 2007). This is an extraordinarilyusefultoolforassessingthedynamicsofrich aspect of themodelling that could be readily improved. starclusters.FortheclusterM4ithasledtotheconclusion (iv) Theearlyevolutionofthemodel:Itisworthdrawing thatM4canbeexplainedasapost-collapseclusters.Italso attention here to the high initial escape velocity from the provides a wealth of information on the distribution of the centreofourmodel,inthefirstfewtensofmillionsofyears. binarypopulation,whichwillbeimportantfortheplanning Clearly this is dependent on the small initial radius. The and interpretation of searches for radial velocity binaries, initial high density is also responsible for the very rapid now underway. initial evolution of thecore. N-body models will eventually supplant the Monte (v) Imperfectionsofthemodelling:Severalimportantim- Carlotechnique,butatpresentareincapableofprovidinga provementsneed to be madeto our technique. star-by-starmodelofevensuchasmallclusterasM4.They (a) At present few-body interactions are handled with can of course provide a great deal of general guidance on cross sections. The problem is not simply that these are the dynamical evolution of rich star clusters, and they un- not well known, especially in thecase of unequal masses; derpin the Monte Carlo by providing benchmarks for small 10 M. Giersz and D.C. Heggie models. Even when N-body simulations eventually become McMillan S. L. W., Vesperini E., Portegies Zwart S. F., possible, Monte Carlo models will remain as a quicker way 2007, ApJ, 655, L45 of exploring the main issues, just as King models havecon- Mikkola, S.1984, MNRAS,207, 115 tinued to dominate the field of star cluster modelling even Peterson, R. C., Rees, R. F., & Cudworth, K. M. 1995, when more advanced methods (e.g. Fokker-Planck models) ApJ,443, 124 havebecome available. Piotto G., King I. R.,Djorgovski S., 1988, AJ, 96, 1918 RicherH. B., et al., 1995, ApJ,451, L17 RicherH. B., et al., 2002, ApJ,574, L151 5 DISCUSSION Richer, H.B., et al. 2004, AJ, 127, 2771 Sommariva, V., et al, in E. Vesperini, M. Giersz, A. Sills, ACKNOWLEDGEMENTS eds, IAUS246, in press StolteA.,BrandnerW.,BrandlB.,ZinneckerH.,2006,AJ, We are indebted to Jarrod Hurley for much help with the 132, 253 BSEstellarevolutionpackage.WethankJanuszKaluznyfor Thorsett, S. E., Arzoumanian, Z., & Taylor, J. H. 1993, hiskindadviceonobservationalmatters,andHarveyRicher ApJ,412, L33 for comments on various matters, especially the observed Trager, S. C., Djorgovski, S., & King, I. R. 1993, in Djor- luminosity functions. This research was supported in part govski, S.G., Meylan G., eds, Structureand Dynamics of by the Polish National Committee for Scientific Research Globular Clusters, ASPCS 50, 347 undergrant 1 P03D 002 27. Trager, S. C., King, I. R., & Djorgovski, S. 1995, AJ, 109, 218 van den Bergh S., 1995, AJ, 110, 1171 REFERENCES Vesperini E., Heggie D.C., 1997, MNRAS,289, 898 Aguilar L., Hut P., OstrikerJ. P., 1988, ApJ, 335, 720 WilkinsonM.I.,HurleyJ.R.,MackeyA.D.,GilmoreG.F., Anderson, J., Bedin, L. R., Piotto, G., Yadav, R. S., & Tout C. A., 2003, MNRAS,343, 1025 Bellini, A. 2006, A&A,454, 1029 Baumgardt, H., & Makino, J. 2003, MNRAS,340, 227 ThispaperhasbeentypesetfromaTEX/LATEXfileprepared Bedin,L.R.,Piotto, G.,King,I.R.,&Anderson,J. 2003, bythe author. AJ, 126, 247 Bedin,L.R.,Anderson,J.,King,I.R.,&Piotto, G.2001, ApJ,560, L75 De Marchi G., Paresce F., Pulone L., 2007, ApJ,656, L65 Dinescu, D. I., Girard, T. M., & van Altena, W. F. 1999, AJ, 117, 1792 Fregeau, J. M., & Rasio, F. A. 2007, ApJ, 658, 1047 Giersz, M. 1998, MNRAS,298, 1239 Giersz, M. 2001, MNRAS,324, 218 Giersz, M. 2006, MNRAS,371, 484 Giersz, M., & Heggie, D. C. 2003, MNRAS,339, 486 Giersz, M., & Heggie, D. C. 2008, preprint Hansen, B. M. S., et al. 2002, ApJ, 574, L155 Hansen, B. M. S., et al. 2004, ApJS, 155, 551 Harris, W. E. 1996, AJ, 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, MNRAS,363, 293 Kim, E., Lee, H. M., & Spurzem, R. 2004, MNRAS, 351, 220 Kroupa P., Tout C. A., Gilmore G., 1993, MNRAS, 262, 545 Kroupa, P. 1995, MNRAS,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, MNRAS, 251, 293 Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545 Kundic,T., & Ostriker, J. P.1995, ApJ, 438, 702 Lamers,H.J.G.L.M.,Gieles,M.,Bastian,N.,Baumgardt, H.,Kharchenko,N.V.,&PortegiesZwart,S.2005,A&A, 441, 117

ebook img

Monte Carlo Simulations of Star Clusters - V. The globular cluster M4 PDF

5 MB·English
by  Mirek Giersz
#additional_collections #journals #arxiv
Save to my drive
Quick download
Download
Upgrade Premium
Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.

Preview Monte Carlo Simulations of Star Clusters - V. The globular cluster M4

Mon.Not.R.Astron.Soc.000,1–10(2002) Printed5February2008 (MNLATEXstylefilev2.2) Monte Carlo Simulations of Star Clusters - V. The globular cluster M4 Mirek Giersz1⋆ and Douglas C. Heggie2 8 1Nicolaus Copernicus Astronomical Centre, Polish Academy of Sciences, ul. Bartycka18, 00-716 Warsaw, Poland 0 2Universityof Edinburgh, School of Mathematics and Maxwell Institute for Mathematical Sciences, King’sBuildings, Edinburgh EH9 3JZ, UK 0 2 n Accepted ....Received...;inoriginalform... a J 5 ABSTRACT 2 We describe Monte Carlo models for the dynamical evolution of the nearby globular clusterM4.Thecodeincludestreatmentsoftwo-bodyrelaxation,three-andfour-body ] h interactionsinvolving primordialbinaries and those formed dynamically,the Galactic p 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 o a satisfactory match to the surface brightness profile, the velocity dispersion profile, r t and the luminosity function in two fields. We describe in particular the evolution s of the core, and find that M4 (which has a classic King profile) is actually a post- a [ collapsecluster,itscoreradiusbeingsustainedbybinaryburning.Wealsoconsiderthe distribution of its binaries, including those which would be observed as photometric 1 binaries and as radial-velocity binaries. We also consider the populations of white v dwarfs,neutronstars,blackholesandbluestragglers,thoughnotallchannelsforblue 8 straggler formation are represented yet in our simulations. 6 9 Key words: stellar dynamics – methods: numerical – binaries: general – globular 3 clusters: individual: M4 . 1 0 8 0 1 INTRODUCTION are concerned with issues where the dynamical history of : v the star cluster is important, where static models are un- The present paper opens up a new road in the study of i informative. Secondly, our aim is to construct a model of a X the dynamical evolution of globular clusters. We adopt the specific star cluster, rather than trying to understand the r Monte Carlo method of Giersz (Giersz 1998, 2001, 2006), general properties of the evolution of a population of star a which in recent years has been enhanced to deal quitereal- clusters. This has been done before, and a brief history is isticallywiththestellarevolutionofsingleandbinarystars, outlined in Giersz & Heggie (2008), but the present work to study the dynamical history of thenearby globular clus- takes these efforts onto a new level of realism, in terms of terM4.Anearlierversionofthecodehadalreadybeenused the description of stellar evolution, and dynamical interac- to study the dynamical history of ω Cen (Giersz & Heggie tions involvingbinary stars. 2003), but at that time the treatment of stellar evolution This problem is not easy. Not only is it necessary to wasprimitiveandtherewerenobinaries.Thenewcodehas use an elaborate technique for simulating the relevant as- been thoroughly tested on smaller systems, by comparison trophysical processes, but it is necessary also to search for with N-body simulations and observations of the old open initialconditionswhich,afterabout12Gyrofevolution,lead star cluster M67 (Giersz & Heggie 2008). There we showed toanobjectresemblingagivenstarcluster.By“resembling” that the Monte Carlo code could produce data of a similar we do not simply mean matching the overall mass, radius levelofdetailandrealismasthebestN-bodycodes.Nowfor andbinaryfractionofacluster,forexample,fortworeasons: the first time we consider much richer systems, with about halfamillionstarsinitially,whichareatpresentbeyondthe reach of N-body methods. (i) We have found that values found for data in the lit- This paper has a place within a long tradition of the erature are highly uncertain, and different sources are con- modelling of globular star clusters, but the place is a dis- tradictory. These data are usually derived, in some model- tinctiveone.First,wearenotconcernedwithastaticmodel dependentway,fromsuchdataassurface-brightnessprofiles of a star cluster at the present day, like a King model. We and velocity dispersion profiles, and we prefer to compare our models directly with this data, and not with inferred global parameters. ⋆ E-mail:[email protected](MG);[email protected](DCH) (ii) We have found that, even if one achieves a satisfac- 2 M. Giersz and D.C. Heggie tory fit to these profiles, the model may give a very poor ical star cluster is represented by its mass, energy and an- comparison with theluminosity function. gular momentum, and its stellar evolutionary state may be computed at any time using synthetic formulae for single From these considerations we conclude that a model which and binary evolution. It may be a binary or a special kind aims to fit only the mass and radius of a star cluster (say) of single star that has been created in a collision or merger may be very far from thetruth. event. Tackling this difficult problem is not just interesting, Neighbouring stars interact with each other in accor- however. We have been motivated by a number of press- dance (in a statistical sense) with the theory of two-body ing astrophysical problems. For example, the two nearby relaxation.Ifoneorbothoftheparticipantsisabinary,the star clusters M4 (the subject of this study) and NGC 6397 probability of an encounter affecting the internal dynamics (which we shall consider in our next paper in this series) iscalculatedaccordingtoanalyticcrosssections,whichalso haverathersimilarmassandradius,andyetonehasaclas- determine the outcome. This is one of the main shortcom- sicKingprofile,whiletheotherisawell-studiedexampleofa ingsofthecode,asthesecrosssectionsarenotwellknownin clusterwith a“collapsed core” (Trager et al. 1995).Among thecaseofunequalmasses,andalsothepossibilityofstellar possible explanations one may consider differences in the collisions duringlong-lived temporary captureis excluded. populationofbinaries,whichareknowntoaffectcoreprop- A star or binary may escape if its energy exceeds a erties,orintidaleffects.Indeedthepresentpaperwillshow certain value, which we choose to be lower than the energy thatthesetwo clustersmay bemuchmoresimilar thanone atthenominaltidalradius,inordertoimprovethescalingof would suppose from the surface brightness profiles alone. thelifetimewithN,asexplainedinGiersz & Heggie(2008). Asecondmotivationforourworkisourinvolvementin This is the second main shortcoming of the models, as it observational programmes aimed at characterising the bi- leadstoacutoffradiusofthemodelthatissmallerthanthe narypopulationsinglobularclusters.Whatdifferences(e.g. truetidal radius, and this lowers the surface density profile in the distributions of periods and abundances) should one in the outerparts of thesystem. expect to find between the core and the halo? These issues A difficulty in applyingtheMonte Carlo code toM4 is are important in the planning of observations, and in their that it employs a static tide, whereas the orbit of M4 ap- interpretation. pearstobeveryelliptical (Dinescu et al. 1999).Wehaveto The cluster M4 is the focus of much of this effort be- assumethataclustercanbeplacedinasteadytideofsucha cause it is nearby, making it a relatively easy target for strengththattheclusterlosesmassatthesameaveragerate deepobservationalstudy.Itwasthefirstglobularclusterto as it would on its true orbit. Some support for this proce- yield a deep sequence of white dwarfs (Richeret al. 1995). dure comes from N-body modelling. Baumgardt & Makino More recently it has been subjected to an intensive ob- (2003) show that clusters on an elliptical orbit between servational programme by the Padova group (Bedin et al. about 2.8 and 8.5kpc dissolve on a time scale intermedi- 2001, 2003; Anderson et al. 2006), which includes searches ate between that for circular orbits at these two radii, and for radial-velocity binaries in the upper main sequence thatthedissolutiontimescalesinalmostthesamewaywith (Sommariva et al 2008). It also turns out to be a cluster thesizeofthesystem.Wilkinson et al.(2003)showthatthe which (we conclude) started with only about half a million core radius of a cluster on an elliptic orbit evolves in very stars, which facilitates the modelling. Along with the open nearly the same way as in a cluster with a circular orbit at clusterM67, M4waschosen bytheinternationalMODEST thetime-averaged galactocentric distance. consortium, at a meeting in Hamilton in 2005, as the focus Allotherfreeparametersofthecode(e.g.thecoefficient for joint effort by theorists and observers, to cast light on of N in the Coulomb logarithm) take the optimal values its binary population and dynamical properties. M67 has found in theabove study. beenmodelled verysuccessfully byHurley et al. (2005),us- ing N-body techniques, and this paper represent the first theoretical step in a similar study of M4. 2.2 Initial Conditions The paper is organised as follows. First, we summarise features of the code and the models, the data we used, and TheinitialmodelsareasspecifiedinTable1.Manyofthese our approach to the problem of finding initial conditions features (e.g. theproperties of thebinaries, exceptfor their for M4. Then we present data for our best models: surface overallabundance)wereinheritedfromourmodellingofthe brightness and velocity dispersion profiles, luminosity func- openclusterM67,andthosewereinturnmainlydrawnfrom tions, theproperties ofthebinarypopulation, whitedwarfs the work of Hurley et al. (2005). Some of the parameters and other degenerate remnants, and the inferred dynami- were taken to be freely adjustable, and this freedom was cal state of the cluster. The final section summarises our exploited in the search for an acceptable fit to the current findingsanddiscussesthemin thecontextofwork onother observational data. clusters, including objects to which we will turn in future papers. 2.3 Observational data and its computation Our first task was to iterate on the initial parameters of 2 METHODS ourmodelinordertoproduceasatisfactoryfittoarangeof 2.1 The Monte Carlo Code observationaldataatanageof12Gyr.Thedataweadopted are as follows: The details of our simulation method have been amply de- scribedinpreviouspapersinthisseries.Eachstarinaspher- (i) Surface brightness profile: here we used the compila- Monte Carlo Simulations of Star Clusters - V. The globular cluster M4 3 Table 1.InitialparametersforM4 Table 2.PropertiesofM4 Distancefromsun a1.72kpc Fixedparameters DistancefromGC 5.9kpc Mass a63000M⊙ Structure Plummermodel Initialmassfunction1 Kroupa,Tout,&Gilmore(1993)inthe Coreradius 0.53pc Half-lightradius 2.3pc range[0.1,50]M⊙ Tidalradius 21pc Binarymassdistribution Kroupaetal.(1991) Binarymassratio Uniform(withcomponent Half-massrelaxationtime(Rh) 660Myr Binaryfraction a1-15% massesrestrictedasfor [Fe/H] -1.2 singlestars) Age b12Gyr Binarysemi-majoraxis Uton5if0oArmUinlog,2(R1+R2) AV a1.33 Binaryeccentricity Thermal,witheigenevolution (Kroupa1995) References:Alldataarefromthecurrentversionofthecatalogue MetallicityZ 0.002 of Harris (1996), except a Richeretal. (2004) (though this is Age 12Gyr(Hansenetal.2004) not always the original reference for the quoted number) and b Hansenetal.(2004). Freeparameters Mass M The line-of-sight velocity dispersion is computed in a Tidalradius rt similar way. Foreach particle we calculate themean square Half-massradius rh line-of-sight velocity (because the orientation of the trans- Binaryfraction fb verse component of the velocity is random), and then sum Slopeofthelowermass overallparticles.Inthissum,eachparticleisweightedbya function α(Kroupa=1.3) geometricalfactorproportionaltothesurfacedensityofthe particle’s shell along the line of sight. The result is an esti- mate of the velocity dispersion that is weighted by neither tion byTrager et al. (1995),where thesurface brightnessis mass nor brightness, but only number. To check the effects expressed in V magnitudes persquare arcsec. of this, we have sometimes calculated velocity dispersions (ii) Radial velocity profile: this came from withvariouscutoffsinthemagnitudesofthestarsincluded. (Peterson et al. 1995), and is the result of binning The effects are usually quitesmall. datatheontheradialvelocities ofnearly200stars.Strictly Computation of the luminosity function is the least we should refer to the line-of-sight velocity dispersion, as problematic. We count stars in each bin in V-magnitude, “radial” velocity has a different meaning in the Monte butlyingabovealineinthecolour-magnitudediagramjust Carlo model. below the main sequence. Again the contribution of each (iii) The V-luminosity function (Richeret al. 2004, from star is weighted by the same geometrical factor, and the V whichweconsideredtheresultsfortheinnermostandouter- magnitudeis corrected for extinction. mostoftheirfourannuli).Thesedataarelackcorrectionfor completeness,thoughthecompletenessfactorfortheouter- most field is plotted in Richeret al. (2002). For main se- 2.4 Finding initial conditions quence stars it is almost 100% down to V = 15, and drops Herewesummariseourexperienceinapproachingthisprob- steadilytolessthan50%atV =17.Fortheinnermostfield lem.Anumberofstudies(e.g.Baumgardt & Makino(2003), thecompleteness correction would belarger. Lamers et al.(2005))givesimpleformulaefortheevolution Now we consider how to compare this data with the with time of the bound mass of a rich star cluster. It is outputoftheMonteCarlocode.Thisincludesalistofeach possible to derive similar simple formulae for the evolution particleinthesimulation, alongwithdataonitsradius,ra- of thehalf-mass radiusand otherquantities. Thereare sev- dial and transverse velocities and absolute magnitude, and eral problems with inverting these formulae, however, i.e. numerousotherquantities. Toconstruct thesurface bright- using them to infer the initial parameters of a cluster from ness, we think of each particle as representing a luminous its present mass and radius. First, these present-day global shell of the same radius, and superpose the surface bright- parameters are quite uncertain, even to within a factor of ness of all shells. While this procedure involves a minimum two.Second,theseformulaedependonthegalacticorbitand of effects from binning, or randomly assigning the full po- other parameters which are equally subject to uncertainty. sition of the star, a shell has an infinite surface brightness Therefore we have adopted the more straightforward but at its projected edge. The effect of this, especially from the more laborious approach of iterating on the initial param- brightest stars, will be seen in some of the profiles to be eters of our models (Tab 1, lower half); that is, we select presented in this paper. A similar problem arises in actual values for the five stated parameters, run the model, find observations, and is often handled by simply removing the where the match with the observations is poor, adjust the brightest stars from the surface brightness data presented. parameters, and repeat cyclically. Theoutputfromthemodeliscorrectedforextinction(Table We have employed two methods to facilitate this pro- 2. The distance of the line of sight from the cluster centre cess to some extent. First, we have often carried out mini- is converted between pc(as in themodel) and arcsec (asin surveys,i.e.small,coarsegrid-searchesaroundagivenstart- theobservations) using the distance in the same table. ingmodel,tofindouthowchangesinindividualparameters 4 M. Giersz and D.C. Heggie affect the results. Second, we have used scaling to acceler- Second,thisexperienceunderscorestheimportanceofprop- atetheprocess, andwenowdescribethismethodinalittle erly normalised luminosity functions. In other words, it is detail. importanttoknowtheareaofthefieldwherethestarshave Suppose we wish to represent a star cluster which has beencounted,orsomeequivalentrepresentationofproperly mass M and radius R with a model representing a cluster normalised data. Very often, the emphasis is solely on the witha(usuallysmaller)massM∗andradiusR∗.Sincetwo- shape of the luminosity function, but we have found that body relaxation dominates muchof thedynamics,we insist theabsolute normalisation is an essential constraint. that thetwo clustershavethesamerelaxation time, and so Baumgardt & Makino(2003)showthatthelowermass function becomes flatter as the fraction of mass lost by the R∗ M 1/3 logγN∗ 2/3 cluster increases (i.e. towards the end of its life). There- = , R (cid:16)M∗(cid:17) (cid:18) logγN (cid:19) fore we could perhaps have improved the fit with the lumi- nosity function by starting with a more massive model and where N,N∗ are the corresponding particle numbers. Thus somehow ensuring a larger escape rate so as toleave a sim- themodeloflowermasshaslargerradius.Thetidalradiusis ilar mass at the present day. Instead, we elected to change scaledinthesameway.Thentheobservationalresults(sur- the slope of the low-mass IMF from the canonical value of facebrightnessprofile,etc)canbecomputedforthemodelof α = 1.3 (Kroupa 2007a) to α = 0.9. (There is some jus- lower mass andthenrescaled (byappropriate factors of the tification for a lower value for low-metallicity populations, massandradius) togivearesultforthemoremassiveclus- though it has been argued (Kroupa2007b) that there is no ter (assuming that the evolution is dominated by the pro- pristine low-metallicity population where the IMF can be cesses of relaxation, stellar evolution and tidal stripping). inferred securely.) 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)initialmass,certainlywhentheproportionof 3.2 A Monte Carlo model of M4 binaries is 10% or less. Some aspects of the evolution are Bysomeexperimentationwearrivedatamodelwhichgave not well described by this scaling technique. For example, a fair fit to all three kinds of observational data; see Ta- wedonotchangethedistributionofthesemi-majoraxesof ble 3 (where we compare with a King model developed by the binaries. In this way the internal evolution of the bina- Richeret al.(2004),andFigs1-4.Itisworthnotingthatno riesiscorrectlymodelled(providedthatthebinariesremain arbitrary normalisation has been applied in these compar- isolated dynamically), though their dynamical interactions isons between our model and the observations. The surface with the rest of the system are not. In principle one could brightness profile, for example, is computed directly from scale the semi-major axes in the same way as R, but then the V-magnitudes of the stars in the Monte Carlo simula- theinternal evolution of the binaries would be altered. tion, as described in Sec.2.3. In the construction of a King model,bycontrast,itisoftenassumedthatthemass-to-light ratio is arbitrary. 3 MODELS OF M4 3.1 Finding initial parameters 3.2.1 Surface brightness Our starting point was our work on the old open cluster M67 (Giersz & Heggie 2008), but with a larger initial mass While the overall surface brightness profile is slightly faint, andradius.(Baumgardt & Makino(2003),forexample,sug- the most noticeable feature of Fig.1 is that the model has gestedthattheinitialmassofM4(NGC6121) wasoforder a somewhat smaller limiting radius than the observational 7.5×105M⊙,thoughweusedsomewhatsmallervalues.)To data. The reason for this is explained in Giersz & Heggie begin with, our choice of initial tidal radius was inferred (2008): in short, we impose a smaller tidal radius than the from theinitialmassandthepresent-dayestimates ofmass nominal tidal radius, in an N-dependent way, which is in- andtidalradiusgiveninTable2,assumingthatr ∝M1/3. tended to ensure that the overall rate of escape from the t At first we adopted similar values for the “concentration” modelbehavesinthesamewayasinanN-bodysimulation. (r /r )and binaryfraction as in themodelling of M67, but Another point to notice is the disagreement between the t h found thatthesurface brightnessprofilefittedpoorly (with total luminosity of our model and that of the King model toolargeacore)unlessthebinaryfractionwasmuchsmaller quoted in the final column of Table 3. Trager et al. (1995) (5 to 10%) and the concentration muchhigher. givean analytic fittothesurface brightnessprofile, and we Before describing our best models, it is worth briefly havecheckedthat theintegrated value is close to ours. mentioningonewhich providedasatisfactory fittothesur- The data for the Monte Carlo model in Table 3 would face brightness profiles. The fit to the velocity dispersion be consistent with a galactocentric radius of about 1.7kpc, profile was tolerable, but indicated a model that was too in an isothermal galaxy model with a circular velocity of massivebyafactorofabout1.4.Itsmainflaw,however,was 220km/s. While this is certainly much smaller than its cur- intheluminosityfunction,whichwasgenerallytoolargeby rent galactocentric distance, a small value was also found a factor in the range 2–3. There are two reasons why this (usingasimilar argument and publishedvalues of themass is interesting. One is that, for a long time, models of star and tidal radius) by van den Bergh (1995). The orbit given clusterswereconstructedentirelyonthebasisofthesurface by Dinescu et al. (1999) has a still smaller perigalactic dis- brightness and velocity dispersion profiles. It should be re- tance,thegalactocentricdistancevaryingbetweenextremes alised that such models may be misleading in other ways. of 0.6 and 5.9 kpc. Monte Carlo Simulations of Star Clusters - V. The globular cluster M4 5 Table 3.MonteCarloandKingmodelsforM4 Quantity MCmodel MCmodel Kingmodel (t=0) (t=12Gyr) (Richeretal.2004) Mass(M⊙) 3.40×105 4.61×104 Luminosity(L⊙) 6.1×106 2.55×104 6.25×104 Binaryfractionfb 0.07 0.057 0 Low-massMFslopeα 0.9 0.03 0.1 Massofwhitedwarfs(M⊙) 0 1.81×104 3.25×104∗ Massofneutronstars(M⊙) 0 3.24×103 Tidalradiusrt (pc) 35.0 18.0 Half-massradiusrh (pc) 0.58 2.89 ∗:thisisthequotedmassof“degenerates” Surface brightness Innermost annulus 16 100 Trager et al MC Model MC Model Richer et al 18 20 10 -2ec arcs 22 N g a m V- 24 1 26 28 0.1 0.1 1 10 100 1000 10000 7 8 9 10 11 12 13 14 15 radius (arcsec) MV Figure 1. Surface brightness profile of our Monte Carlo Figure 3.LuminosityfunctionofourMonteCarlomodelat model,comparedwiththedataofTrageretal.(1995). the median radius of the innermost annulus in Richeretal. (2004),comparedwiththeirdata. Velocity dispersion Outermost annulus MC Model 100 Peterson et al MC Model 5 Richer et al 4 ec) m/s 3 V (kz N 10 S M R 2 1 0 0.01 0.1 1 10 100 1 7 8 9 10 11 12 13 14 15 radius (arcmin) MV Figure 2. Velocity dispersion profile of our Monte Carlo Figure 4.LuminosityfunctionofourMonteCarlomodelat model,comparedwiththedataofPetersonetal.(1995). the median radius of the outermost annulus in Richeretal. (2004),comparedwiththeirdata. 6 M. Giersz and D.C. Heggie 3.2.2 Velocity dispersion profile This is illustrated in Fig.2, where it is compared with the observationaldataofPeterson et al.(1995).Theshortfallat large radii is of doubtfulsignificance. 3.2.3 Luminosity Functions These are shown for our model at the median radius of theinnermostandoutermostfieldsobservedbyRicher et al. (2004). The disagreement in the inner luminosity function at faint magnitudes may be attributable to the fact that thetheoreticalresultassumes100%completeness,whilethe observationaldataareuncorrectedforcompleteness.Aplot of the completeness correction in the outer field is given by Hansen et al. (2002, Fig.3) and discussed above in Sec.2.3. 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 Figure 5.TheoreticalcoreradiusofourMonteCarlomodel. be accounted for by the estimated value of the complete- nesscorrection.Theerrorbarsinthelasttwoobservational points on this plot (not shown) almost overlap, and much 3.2.5 The colour-magnitude diagram of the mismatch may be simply sampling uncertainty. It is Figure6showsthecolour-magnitudediagramofthemodel. worthcomparingthemulti-massKingmodelconstructedby This is of interest, not so much for comparison with ob- Richeret al. (2004),which isalso problematic in theouter- servations, but for the presence of a number of interesting most field. features. The division of the lower main sequence is simply an artifact of the way binary masses were selected (a total massabove0.2M⊙andacomponentmassabove0.1M⊙.)Of particularinterestarethehighnumbersofmergerremnants 3.2.4 Core collapse on thelower whitedwarf sequence.Thereare veryfewblue stragglers.Partlythisisaresultofthelowbinaryfrequency, Fig.5 shows the evolution with time of the theoretical core but it is also important to note that some formation chan- radius. There is an early period of very rapid contraction, nelsareunrepresentedinourmodels(inparticular,collisions associated withmasssegregation,followed byaslowerreex- during triple or four-body interactions, though if a binary pansion, caused by the loss of mass from the evolving mas- emergesfromaninteractionwithappropriateparameters,it sive stars which are now concentrated within thecore. will be treated as merged.) These numbers also depend on Ourmost surprisingdiscoveryfrom ourmodelofM4is theassumedinitialdistributionofsemi-majoraxis,whichis thesubsequent behaviour. M4 is classified as a King-profile notyetwellconstrainedbyobservationsinglobularclusters. cluster (Trager et al. 1993), and such clusters are usually Note in Fig.6 that there are some white dwarf-main interpreted as being clusters whose cores have not yet col- sequence binaries (below the main sequence). Richer et al. lapsed.Buttheplotofthetheoretical coreradius(Fig.5 re- (2004)drewattentiontothepossiblepresenceofsuchbina- vealsthatthemodelexhibitedcorecollapse atabout 8Gyr. ries in theircolour magnitude diagrams of this cluster. Subsequentlyits core radius is presumably sustained by bi- naryburning.Evennon-primordialbinariesmaybeplaying a role here. To the best of our knowledge it has not pre- 3.2.6 Binaries viously been suggested that M4, is a post-collapse cluster, thoughonstatisticalgroundsDe Marchi, Paresce, & Pulone PhotometricbinariesarevisibleinFig.6,andthesearecom- (2007)havesuggestedthatsomeKing-typeclustershaveal- pared with observations in the inner field of Richer et al. readycollapsed. Thisissueisoften approachedbyreference (2004)inFig.7.Inthisfigure,themodelhistogramhasbeen tothehalf-massrelaxationtime(Table2,butforM4thisis normalised to the same total number of stars as the obser- not short enough to bedecisive. vationalone. Wemadenoattempt tosimulate photometric The presence of radial colour gradients is correlated errors,butthebinsaroundabscissa=-0.75suggestthatthe with the presence of a non-King surface brightness profile binary fractions in the model (which is under 6% globally; andhencewithcorecollapse. Wehavecomputedthecolour see Table 3) and the observations are comparable. We note profile of our model, and find that it is nearly flat except here that Richeret al. (2004), using this same data, con- for the influence of one or two very bright stars within the cluded that the binary frequency was approximately 2% in innermostfewarcsec.Mostvaluesliearound0.65inB−V, the innermost field. But they also note that the measured which is much less than the global value of 1.03 given in frequency of “approximately equal-mass binaries” is 2.2%, Harris (1996, January 2008). On the other hand roughly and we consider that the balance between this number and similarmismatchesbetweenobservationsarefoundinother ourbinaryfractioncanbemadeofbinarieswithcompanions clusters (e.g. M30, (Piotto, King, & Djorgovski 1988)). whose masses are more unequal. Monte Carlo Simulations of Star Clusters - V. The globular cluster M4 7 Binary sequence 60 CM Diagram RicMheCr eMt oadl einl,n neor ramnanliusleuds 50 0 40 N 30 Binaries 5 CollisioSnin pgrleo dsutacrtss 20 10 V 10 0-1.5 -1 -0.5 0 0.5 1 1.5 V offset from single star MS Figure 7. Histogram of V-offset from the main sequence, compared with the corresponding data from the innermost annulusstudiedbyRicheretal.(2004) Seetextfordetails. 15 Distribution of semi-major axes 1 0 Gyr 12 Gyr 0.9 0.8 20 -0.5 0 0.5 1 1.5 2 0.7 B-V 0.6 sFtiagrus;rreed6:.bTihnearcioelso;ubrl-umeapglnuisteusd:ecdolilaigsiroanmoartm12erGgyerr.rGemrenenan:tssin.gle Fraction 0.5 0.4 0.3 RadialvelocitybinariesshouldalsobedetectableinM4, 0.2 and will be part of a separate investigation. 0.1 Now we consider the distributions of the binaries, as 0 0.1 1 10 100 1000 10000 100000 predicted by the theoretical model. Fig.8 shows that bina- Semi-major axis (solar radii) ries have evolved dynamically as well as through their in- Figure8.Distributionfunctionofthesemi-majoraxesofthe ternal evolution. In particular the softest pairs been almost binariesat0Gyrand12Gyr.Units:solarradii. destroyed. By 12Gyrthebinariesexhibit segregation towards the centreof thecluster, butperhapsin more subtleways than value in the first few tens of millions of years, if valid for mightbeexpected(Figs.9,10).Whenallbinariesareconsid- M4, could have interesting consequences for the early evo- ered, there is little segregation relative to the other objects lution of the cluster and its stars. It draws attention to the inthesystem.(Mostbinariesinourmodelareoflowmass.) very high initial density of our model (Table 3), whose av- But if one restricts attention to bright binaries, which we erage value within the half-mass radius is 2×105M⊙/pc3. here take to mean those with MV < 7 (i.e. brighter than Thecentraldensityis about 106M⊙/pc3,about an orderof abouttwomagnitudesbelowturnoff),thesegregationisvery magnitude larger than in the central young cluster in the noticeable (Fig.10), with a half-mass radius smaller by al- HIIregion NGC3603 (Stolte et al. 2006).Intheabsenceof most a factor of 2 than for bright single stars. Still, bright localyoungstarclustersasmassiveasourM4progenitor,it binaries are not nearly as mass-segregated as neutron stars isdifficulttobesurewhetherourinitialmodelisimplausibly (Fig.9), which, incidentally, receive no natal kicks in our dense. model. The history of the binary fraction is, in effect, given in Fig.13. In thefirst Gyr this falls steadily from theinitial 3.2.8 Degenerate components valueof0.07toabout0.06,anditremainsclosetothisvalue untilthe present day. 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 3.2.7 Escape velocity the historical evolution of this population over the lifetime We have already referred to the fact that, in our model, of thecluster so far, according to our model. neutronstarsreceivenonatalkicks.Escapeisofcoursegov- Fig.12 shows the evolution of these populations. It has erned by the escape speed, and this is of much interest in already become well established (Vesperini& Heggie 1997; connectionwiththepossibilityofretaininggasfromthefirst Hurley & Shara 2003) that white dwarfs account for an in- generation ofrapidly evolvingstars.Forthisreason weplot creasingproportionofthemassofglobularclusters,andin- the central escape velocity in Fig. 11. The remarkably high deedit isoforder39% at thepresent day,accordingtothis 8 M. Giersz and D.C. Heggie Radial distribution 1 All objects Binaries 0.9 Neutron stars White dwarfs Collision products 0.8 0.7 0.6 Fraction 0.5 0.4 0.3 0.2 0.1 0 0.1 1 10 Radius (pc) Figure9.Radialdistributionfunctionsat12Gyr.Units:par- secs. Thekey identifies the class of object included. The dis- Figure12.Massofwhitedwarfs,andofneutronstarsandblack tributionsofwhitedwarfsandbinariesarealmostidentical. holestogether. Radial distribution of bright objects 1 All objects All bright objects (V<7) 0.9 Bright binaries 0.8 0.7 0.6 Fraction 0.5 0.4 0.3 0.2 0.1 0 0.1 1 10 Radius (pc) Figure 10. Radial distribution functions at 12Gyr, showing the extent of segregation between bright single and bright binary stars (MV < 7). Also shown for comparison is the Figure13.Numbersofallstars,whitedwarfs,neutronstarsand distributionforallobjects, asinFig.9.Units:parsecs. blackholes. model. The proportion of neutron stars is almost certainly excessive, because of ourassumption of complete retention. In order to separate two of the components in this fig- ure (black holes and neutron stars), we show in Fig. 13 the numbers of these degenerate stars, along with the numbers ofallstarsand allbinaries.Thisshows thatnostellar-mass black holes are expected to bepresent in M4 now. 4 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-bodybinaries,andtheinternalevolutionofsinglestars andbinaries.Byadjustmentoftheinitial parameters(total mass,tidalradius,initialmassfunction,binaryfraction)we havefoundamodelwhich,after12Gyrofevolution,leadsto amodelwithasurfacebrightnessprofile,velocitydispersion profile, and luminosity functions at two radii, all of which Figure 11.Centralescapevelocityasafunctionoftime. are in tolerable agreement with observational data. It also leadstoanumberofphotometricbinariesroughlyconsistent with observation. Monte Carlo Simulations of Star Clusters - V. The globular cluster M4 9 Thismodelhasacurrentmassof4.6×104M⊙,aM/LV italsomeansthatweareunabletodetermineifacollision ratio of 1.8, a binary fraction of almost 6%, and an almost occursduringalong-livedinteraction.Forthisreason,we flat lower mass function (Table3).Almost 40% of themass have said almost nothing about collision products (e.g. is in white dwarfs, and about 7% in neutron stars, though blue stragglers) in this paper. This limitation could be ourmodelassumesa100%retentionrate.Theyarestrongly overcome by direct integration of the interactions, as is mass-segregated to the centre of the cluster. There are no done by Fregeau & Rasio (2007) in their version of the stellar-massblackholesleftinthecluster.Thebinarieshave MonteCarlo scheme. experience significant dynamical evolution, almost all the (b) Long-lived triples are neglected in the model at soft pairs having been destroyed. The binary population as present. These are commonly produced in binary-binary a whole is only slightly segregated towards the centre, but encounters (Mikkola 1984), and it is desirable to include there is more evident segregation of bright binaries (such these as a third species (beyond single and binary stars). as those that would be more readily observed in a radial Theirobservableeffectsmaybesmall,butofcoursethere velocity search.) is one intriguing example in the very cluster we have fo- The most significant new result from our model is the cused on here(Thorsett et al. 1993). implication that M4 is a core collapse cluster, despite the (c) Weassumethatthetideismodelled asatidalcut- uncollapsed appearance of thesurface brightness profile. off. We have taken some care to ensure that the tidal Nowwediscussanumberofshortcomingsofourmodel radius is adjusted so that the model loses mass through and other issues related to thisstudy. escape at the same rate as an N-body model would, if immersed in a tidal field with the same tidal radius. (i) Uniqueness:Thoughwehavearrivedatabroadlysat- Thismeans,however,thattheeffectivetidalradiusofthe isfactory model, it is not at all clear how unique it is. Cer- model is somewhat too small. A better treatment of a tainlywewereunabletoconstructmodelswithmuchlarger steady tide may be possible, without thisdrawback. numbers of primordial binaries, or with significantly larger (d) Weassumethatthetideissteady,somethingwhich initialradius,assuchmodelsproducedaninsufficientlycon- is not true for M4. The effects of tidal shocks have been centrated surface brightness profile. Furthermore, the fact studied by a number of authors (e.g. Kundic& Ostriker thatthelowerslopeoftheinitialmassfunctionislesssteep 1995), and it would be possible to add the effects as thanthecanonicalvalueof1.3shouldnotbetakentoimply another process altering the energies and angular mo- that we haveinferred theinitial value. menta of the stars in the simulation. On the other hand (ii) Fluctuations:wehavenoticedthatrunsdifferingonly Aguilar, Hut,& Ostriker (1988) found, on the basis of a intheinitialseed oftherandomnumbergenerator cangive simple model, that tidal evaporation was the dominant surprisingly different surface brightness profiles. In broad mechanism in the evolution of the mass of M4 at the termsthis confirmstheimportant findingof Hurley (2007), presentday,andthatotherfactors(suchasdiskandbulge though we have not yet established (as he did) that the shocking) contributed at a level less than 1%. presence or absence of black hole binaries is the underly- (e) Rotation:ithasbeenshown(Kim et al.2004)that, ingmechanism.Weshallreturntothisissueat appropriate to the extent that rotating and non-rotating models can length in the next paper in this series, on the cluster NGC be compared, rotation somewhat accelerates the rate of 6397. corecollapse.RotationishardtoimplementinthisMonte (iii) TheInitialModel:Theinitialmodelisastonishingly dense, the central density being of order 106M⊙pc−3. It Carlo model, however. (f) The search for initial conditions is still very labo- also has to be realised that we are imposing initial con- rious, and we are constantly seeking ways of expediting ditions at a time when the residual gas from the birth this. Our most fruitful technique at present is the use of of the cluster has already dispersed. Various properties of small-scale models which relax at thesame rateas afull- a star cluster, including its binary population and mass sizedmodel.Notallaspectsofthedynamicsarefaithfully function, may change significantly during the phase of renderedbyascaled model,butthetechniqueappearsto gas expulsion, which further undermines the power of our besuccessfulaslongasthefractionofprimordialbinaries model to establish the initial conditions. Finally, we have is not too large. assumed no initial mass segregation, even though recent work suggests that this should be present already before Despitethesecaveatsandshortcomings, itisclearthat the cluster has assembled into a roughly spherical object the Monte Carlo code we have been developing is now an (McMillan, Vesperini, & Portegies Zwart 2007). This is an extraordinarilyusefultoolforassessingthedynamicsofrich aspect of themodelling that could be readily improved. starclusters.FortheclusterM4ithasledtotheconclusion (iv) Theearlyevolutionofthemodel:Itisworthdrawing thatM4canbeexplainedasapost-collapseclusters.Italso attention here to the high initial escape velocity from the provides a wealth of information on the distribution of the centreofourmodel,inthefirstfewtensofmillionsofyears. binarypopulation,whichwillbeimportantfortheplanning Clearly this is dependent on the small initial radius. The and interpretation of searches for radial velocity binaries, initial high density is also responsible for the very rapid now underway. initial evolution of thecore. N-body models will eventually supplant the Monte (v) Imperfectionsofthemodelling:Severalimportantim- Carlotechnique,butatpresentareincapableofprovidinga provementsneed to be madeto our technique. star-by-starmodelofevensuchasmallclusterasM4.They (a) At present few-body interactions are handled with can of course provide a great deal of general guidance on cross sections. The problem is not simply that these are the dynamical evolution of rich star clusters, and they un- not well known, especially in thecase of unequal masses; derpin the Monte Carlo by providing benchmarks for small 10 M. Giersz and D.C. Heggie models. Even when N-body simulations eventually become McMillan S. L. W., Vesperini E., Portegies Zwart S. F., possible, Monte Carlo models will remain as a quicker way 2007, ApJ, 655, L45 of exploring the main issues, just as King models havecon- Mikkola, S.1984, MNRAS,207, 115 tinued to dominate the field of star cluster modelling even Peterson, R. C., Rees, R. F., & Cudworth, K. M. 1995, when more advanced methods (e.g. Fokker-Planck models) ApJ,443, 124 havebecome available. Piotto G., King I. R.,Djorgovski S., 1988, AJ, 96, 1918 RicherH. B., et al., 1995, ApJ,451, L17 RicherH. B., et al., 2002, ApJ,574, L151 5 DISCUSSION Richer, H.B., et al. 2004, AJ, 127, 2771 Sommariva, V., et al, in E. Vesperini, M. Giersz, A. Sills, ACKNOWLEDGEMENTS eds, IAUS246, in press StolteA.,BrandnerW.,BrandlB.,ZinneckerH.,2006,AJ, We are indebted to Jarrod Hurley for much help with the 132, 253 BSEstellarevolutionpackage.WethankJanuszKaluznyfor Thorsett, S. E., Arzoumanian, Z., & Taylor, J. H. 1993, hiskindadviceonobservationalmatters,andHarveyRicher ApJ,412, L33 for comments on various matters, especially the observed Trager, S. C., Djorgovski, S., & King, I. R. 1993, in Djor- luminosity functions. This research was supported in part govski, S.G., Meylan G., eds, Structureand Dynamics of by the Polish National Committee for Scientific Research Globular Clusters, ASPCS 50, 347 undergrant 1 P03D 002 27. Trager, S. C., King, I. R., & Djorgovski, S. 1995, AJ, 109, 218 van den Bergh S., 1995, AJ, 110, 1171 REFERENCES Vesperini E., Heggie D.C., 1997, MNRAS,289, 898 Aguilar L., Hut P., OstrikerJ. P., 1988, ApJ, 335, 720 WilkinsonM.I.,HurleyJ.R.,MackeyA.D.,GilmoreG.F., Anderson, J., Bedin, L. R., Piotto, G., Yadav, R. S., & Tout C. A., 2003, MNRAS,343, 1025 Bellini, A. 2006, A&A,454, 1029 Baumgardt, H., & Makino, J. 2003, MNRAS,340, 227 ThispaperhasbeentypesetfromaTEX/LATEXfileprepared Bedin,L.R.,Piotto, G.,King,I.R.,&Anderson,J. 2003, bythe author. AJ, 126, 247 Bedin,L.R.,Anderson,J.,King,I.R.,&Piotto, G.2001, ApJ,560, L75 De Marchi G., Paresce F., Pulone L., 2007, ApJ,656, L65 Dinescu, D. I., Girard, T. M., & van Altena, W. F. 1999, AJ, 117, 1792 Fregeau, J. M., & Rasio, F. A. 2007, ApJ, 658, 1047 Giersz, M. 1998, MNRAS,298, 1239 Giersz, M. 2001, MNRAS,324, 218 Giersz, M. 2006, MNRAS,371, 484 Giersz, M., & Heggie, D. C. 2003, MNRAS,339, 486 Giersz, M., & Heggie, D. C. 2008, preprint Hansen, B. M. S., et al. 2002, ApJ, 574, L155 Hansen, B. M. S., et al. 2004, ApJS, 155, 551 Harris, W. E. 1996, AJ, 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, MNRAS,363, 293 Kim, E., Lee, H. M., & Spurzem, R. 2004, MNRAS, 351, 220 Kroupa P., Tout C. A., Gilmore G., 1993, MNRAS, 262, 545 Kroupa, P. 1995, MNRAS,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, MNRAS, 251, 293 Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545 Kundic,T., & Ostriker, J. P.1995, ApJ, 438, 702 Lamers,H.J.G.L.M.,Gieles,M.,Bastian,N.,Baumgardt, H.,Kharchenko,N.V.,&PortegiesZwart,S.2005,A&A, 441, 117

See more

The list of books you might like

Upgrade Premium
Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.
DMCA & Copyright: Dear all, most of the website is community built, users are uploading hundred of books everyday, which makes really hard for us to identify copyrighted material, please contact us if you want any material removed.
Copyright @ PDFdrive.to, 2022 | We love our users