FOKKER-PLANCK MODELS AND GLOBULAR CLUSTER EVOLUTIONbyGORDON ALAN DRUKIERB.Sc., The University of Toronto, 1985A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPYinTHE FACULTY OF GRADUATE STUDIES(Department of Geophysics & Astronomy)We accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIAOctober 1992© Gordon Alan Drukier, 1992B"HIn presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)Department of 6eof ki54-3^As 1-rovlol The University of British ColumbiaVancouver, CanadaDate ^Nen/ 72DE-6 (2/88)AbstractNumerical models for globular cluster evolution using the orbit-averaged Fokker-Planck equationare compared with observations of the globular clusters M71 and NGC 6397. The first set ofmodels studied includes a mass spectrum, a tidal boundary and allows for a central energy sourcefrom the formation, hardening and destruction of binaries formed in three-body interactions. Thesemodels are compared to star-count mass functions and surface density profiles for M71. It ispossible to reproduce the tidal boundary and degree of mass segregation using these models but thecentral density profiles of the fitted models are much steeper than is observed. On the basis of thedegree of mass segregation and the short relaxation time, M71 is dynamically a highly evolvedobject, but the surface density profile contradicts this. Several modifications of the model areconsidered to resolve this difficulty, and an additional energy source, which does not requireextreme densities to operate, appears to be required.New observations of the cluster NGC 6397 are presented. This cluster is presently at a similarposition in the galaxy as M71, and has a similar mass and mass function, but is much morecentrally concentrated. NGC 6397 is a low metallicity cluster with halo kinematics, in contrast tothe metal-rich disk cluster M71. The observations have been reduced to give a surface densityprofile, and mass functions at two distances from the cluster centre. The surface density profile isfit by a power-law with slope —0.8 in the central arc minute. The mass functions give someevidence for mass segregation, but the low signal-to-noise of the outer mass function limits thesignificance of this observation.A second set of models, augmented to include the effects of mass loss due to stellar evolution,are used to find matches for the observations of M71 and NGC 6397. A diagram showing theevolution with time of the amount of mass segregation and the half-density radius is shown to beuseful in identifying good candidate models. Stellar evolution appears not to solve the problem ofM71 as no satisfactorily matching model has been found. Possible solutions to this problem arediscussed. The observations of NGC 6397 can be reproduced by the models and the best suchmodel is described. This model has an initial mass function with a Salpeter-like slope and isundergoing core-collapse. This demonstrates that a numerical model based on the orbit-averagedFokker-Planck equation, and including an evolving mass spectrum, a tidal boundary and the binaryenergy source, is a useful model for the dynamical evolution of some globular clusters.llTable of ContentsAbstract^ iiList of Tables vList of Figures^ viList of Abbreviations viiiList of Symbols^ ixAcknowledgements xChapter 1. Introduction^ 1Section a. Recent Results in Fokker-Planck Modelling^ 5Section b. Comparing Fokker-Planck Models with Observations^ 11Chapter 2. The Models^ 15Section a. Basic Formalism 15Section b. Multiple Masses^ 16Section c. Binary Heating 19Section d. Tidal Boundary^ 24Section e. Stellar Evolution 29Section f. Initial Conditions^ 37Section g. Caveats^ 38Chapter 3. The Problem of M71^ 39Section a. Outline of the Problem 39Section b. Approaches to a Solution^ 48Subsection I. Gravothermal Oscillations 48hiSubsection II. Extreme Degenerates^ 49Subsection III. Modifying the Heating Rate 51Chapter 4. NGC 6397^ 56Section a. Observations and Reduction^ 58Section b. Surface Density Profile 64Section c. Mass Functions^ 75Chapter 5. Models Including Stellar Evolution ^ 91Section a. M71^ 94Section b. NGC 6397 99Chapter 6. Conclusions.^ 105References^ 110Appendix A. Initial Model Parameters^ 115Appendix B. Systematics of the S r „,--rhd Diagram^ 122ivList of Tables1. Core collapse times for two component models^ 192. Results of the sample stellar evolution grid 363. Field locations for NGC 6397 observations^ 594. Surface density profile for 1<14^ 685. Surface density profile for R15.5 726 Luminosity function for field du Pont:bf^ 767. Luminosity function for field du Pont:if 788. Luminosity function for field du Pont:n^ 799. Luminosity function for field du Pont:w 8010. Background subtracted luminosity functions for du Pont:if and du Pont:out^ 8211. Mass functions for du Pont:if and du Pont:out^ 8612. Observed mass spectral indices^ 8813. Mass spectral indices for model NGC 6397—A^ 104Al. Parameters for the model CW1^ 115A2. Parameters for the "standard" model of M71^ 116A3. Parameters for the "black hole" model of M71 117A4. Parameters for the initial set of models in Ch. 5^ 118A5. Parameters for model M71—A^ 119A6. Parameters for model M71—B 120A7. Parameters for model NGC 6397—A^ 121Bl. IMF families of Fig. B1 ^ 124B2. Results of models shown in Fig. B2^ 124List of Figures1.Projected radial profiles for the mean mass in model CW1^ 332. Projected surface density profiles for model CW1 343. Evolution of the number of stars by species in model CW1^ 354. Observations of M71^ 405. Observed segregation measure for M71^ 446. Results of the "standard" model 467. Segregation measure in the "standard" model^ 478. Results of the "black hole" model^ 509. Results of the "extra heating" model 5210.Kinematics of M71 and the "extra heating" model^ 5311.Segregation measures for the "black hole" and "extra heating" models^ 5412.Map of the fields observed in NGC 6397^ 6013.Map showing the annuli used in constructing the SDP^ 6514.Surface density profile for k 14 for the annuli in Fig. 6615. Surface density profile with merged sections^ 6916.Central surface density profiles from the four quadrants of Swope:c^ 7017. Combined /<14 and R15.5 surface density profiles^ 7318.Power-law fits to the surface density profile 7419.Comparison of the blank field star counts with those of Fahlman et al. (1989)^ 8120. Luminosity function for du Pont:if^ 8321. Mean luminosity function for du Pont:w and du Pont:n^ 8422. Composite of all deep luminosity functions for NGC 6397 8523. Mass functions for the three NGC 6397 luminosity functions^ 87vi24. Segregation measures between the three mass functions^ 9025. Comparison with observations for model M71—A 9526. Sr m rid diagram for M71 observations^ 9727. Comparison with observations for model M71—B^ 9828. S ,.„—rhd diagram for NGC 6397 observations 10029. Comparison with observations for model NGC 6397—A^ 10130. Velocity dispersion profile for model NGC 6397—A 103Bl. S r,„,—rhd diagram for M71 divided by IMF family^ 125B2. Full S r„—rhd curves for a set of models differing only in IMF^ 126viiList of AbbreviationsThere are certain phrases and references used often enough in this thesis to receive the dubioushonour of being replaced by abbreviations. Although all the abbreviations are defined where theyfirst occur, the reader might find it useful to have the definitions collected together for easyreference.JargonCCD^charge coupled deviceGTO gravothermal oscillationIMF^ initial mass functionLF luminosity functionLMC^Large Magellanic CloudMF mass functionMSI^mass spectral indexPCC post-core-collapseSDP^surface density profileReferencesCW^ Chernoff & Weinberg 1990DFR Drukier, Fahlman, & Richer 1992FRST^Fahlman, Richer, Searle, & Thompson 1989GCLM Grabhorn, Cohn, Lugger, & Murphy 1992LFR^Lee, Fahlman, & Richer 1991MCH Murphy, Cohn, & Hut 1990SOC^Statler, Ostriker, & Cohn 1985viiiList of SymbolsE^EnergyE,^Energy at tidal boundary^ §2df^Distribution functionIn A^Coulomb logarithm §2bm Stellar massm i^Mean mass for bin iInitial mass for stars in bin j^ §2eFinal mass for stars in bin j §2eM^Total massMr^Total initial mass for a model, scale mass^§2bN^Total number of starsN,^Total number of stars in bin iN,„^Number of stars per unit mass (mass function)r^Radial positionr0^Scale radius^ §2bCore radius Eq. (2.12)reff^Effective radius §4brh^Half-mass radius^ §2crhd^Half-density radius §5art^Tidal radius §2dSr^Radial segregation measure^Eq. (3.1)Srth^Mass-range normalized segregation measure^Eq. (3.2)SMm71^S r /n for M71 observations Eq. (3.3)tr,Timetth^Half-mass relaxation time^ Eq. (2.11)v^Velocity, velocity dispersionvo^Scale velocity^ Eq. (2.5)Wo^Dimensionless central density in King models^§2ex Mass spectral index Eq. (1.3)xh^Mass spectral index for high mass stars^Ch. 5x,^Mass spectral index for low mass stars Ch. 5p^Mass densityPi^Mass density for mass bin i^ Eq. (2.4)Pc^Central density^ §2ca^projected number densityan,^projected mass densityAcknowledgementsLike a globular cluster, the work of preparing a Ph. D. does not exist in isolation, and one-halfpage is permitted by university regulations for acknowledging this. Thanks are due to my supervisorDr. Greg Fahlman for many discussions, and for much advice and encouragement. Without hismany useful comments on the manuscript, there are sections of this thesis which would have beenincomprehensible. Together with Dr. Ian Thompson and the people responsible for the Observatoriesof the Carnegie Institution of Washington, he also provided the data on which Chapter 4 is based.The Fokker-Planck code that forms the foundation of this thesis was kindly given to us by Dr.Hyung Mok Lee and Dr. Peter Stetson sent me the latest copies of DAOPHOT & ALLSTAR as wellas development versions of some of his other software.My fellow grad students at UBC deserve a vote of thanks for serving as sounding boards whenI had something to talk about and for putting up with my muttering when I didn't. Special thanksgo to old-timers Claudia Mendes and Dave Woods and other sounding boards James Brewer, PhilHodder, Brad Gibson, Ted Kennelly, Steve Holland and Dave Hogg.My family were a great support through all this and my parents are especially to be thanked forletting me live 3300 km away and visit home far too infrequently. My adopted family, in the formof my friends in Vancouver, were also very supportive. I will miss them all very much.This last acknowledgement should have been first, and goes to the Creator of the Universe forgiving me the strength to do this work and for the remarkable miracle that the human mind isactually capable of making sense of the physical world he created..17 , 7 -17 inn 1 , 7 , n17.717M1 11= ^v-incn crOT`rThe heavens declare the glory of G-d; the expanse of the sky proclaims His handiwork.Psalms 19:2Chapter 1. IntroductionA spiral galaxy such as our own is composed of many parts. The spinning, flattened disk showsoff spiral arms marked by the brief, but brilliant, light of the young, massive stars. Most of themass of the disk, however, is made of their cooler, less massive, and longer-lived cousins togetherwith the gas from which the stars are built. Distinct from the disk is the spheroid; an olderpopulation of stars, less endowed with heavy elements. The spheroid, or halo, stars define aroughly spherical volume enclosing the disk. The space velocities of the halo stars are highlyrandom and show little of the orderly rotation of the disk.Along with the general field of single and binary stars, a galaxy contains distinct subsystemsarising from the preference of stars to form in groups. These groups range from small, multiplestar systems through to unbound associations and bound clusters. The disk is home to openclusters; somewhat loose groupings of hundreds of stars which, subject to the stresses of anenvironment with large density variations, are fortunate to survive together for a few billion years.Beyond the open clusters in scale are the grand dames of galactic subsystems, the globularclusters.The stars which make up a globular cluster formed at the same time from a chemically homogeneouscloud of gas. From the fitting of stellar evolution isochrones to the colour-magnitude diagrams ofglobular clusters, the ages of the globular clusters are deduced to be 15±3 Gyr (Tayler 1986;VandenBerg 1988). Much of the uncertainty reflects shortcomings in the stellar models, but it islikely that there is also a several billion year age spread amongst the clusters (VandenBerg, Bolte,& Stetson 1990). As the oldest identifiable component of the Galaxy, globular clusters potentiallyhave much to tell us about galaxy formation. In themselves, the globular clusters, which contain105 or more stars, are fascinating laboratories for studying stellar dynamics. For many decades,two-body relaxation has been identified as the force driving the dynamical evolution of globularclusters. Much effort from a number of directions has been put into understanding globular clusterevolution, but the most fruitful approach has been the numerical models based upon the Fokker-Planckequation.1There comes a time for every theoretical model when it must come up against the phenomena itis supposed to explain. In principle, the evolution of a globular cluster is simply an extension tolarge N of the N-body problem of classical mechanics. Since we are not pursuing exact predictionsfor the position and velocity of each and every star observed in a globular cluster, a statisticaldescription is sufficient. In a certain sense, even a direct N-body integration of a large N systemmay be considered statistical in that the overall behaviour of the model will depend more on theaverage properties of stars with similar orbits than on the trajectories of the individual bodies. Amore direct statistical approach is to describe the structure of the cluster in terms of a phase-spacedistribution functionf(r,v,t) for which the time evolution is governed by the Boltzmann equation,at + v • Vf — V(I) - I-a = (—df ^(1.1)^av^dtwhere (df/dt)enc is the change in the distribution function due to encounters between stars. For aself-gravitating system the potential (I)(r,t) satisfies Poisson's equationV24:0 = 47cGi fd3v.^ (1.2)If the average time between encounters is long enough that the encounter term may be neglected,then f is a solution of the collisionless Boltzmann equation,at + v • Vf —^av^Vc1)^= 0 .^(1.3)This equation gives the evolution of a stellar system under the influence of a smooth potential.Unlike the case of extremely large N, as in galaxies, the number of bodies in a globular clusteris small enough that encounters between stars are sufficiently frequent that they cannot be neglected.For a star with phase space coordinates w, let tP(w,Aw)d 3AwAt be the probability that the star isscattered by two-body encounters into a volume of phase space d3Aw around w +Aw during thetime interval At . Then the encounter term in eq. (1.1) can be written as the sum of the rates forstars to be scattered into and out of a given position in phase space, ie.(6I) =^— Aw, Aw)f (w — AW)d 3 AW — f (W)f P (W , Aw)d 3 Aw . (1.4)dt encIt can be shown (eg. Binney & Tremaine, 1987, ch. 8; Spitzer, 1987, ch. 2) that under the2conditions expected in a globular cluster, the net effect of the fluctuating forces between stars, ie.the deviation from the smooth potential, is dominated by the accumulated sum of many weakencounters rather than a few strong interactions. Goodman (1983) produced a numerical modelusing an orbit-averaged equation which included strong encounters and confirmed that, with theexception of an increase in the rate of evaporation of stars from the model, the model's evolutionwas basically unchanged. In the case of weak encounters, Awl is small and we can expand111(w – Aw,Aw)f (w – Aw) in a Taylor series in Aw :6^atp(w — Aw, Aw)f (w – Aw) =^Aw)f (w)^Aw —[t 1 f(w , Aw) f (w)11 6^o2-F— I Aw.Aw.^ [tlf(w, Aw)f (w)] + 0(Aw 3 ) . (1.5)2 awiawiThe Fokker-Planck approximation lies in ending the expansion at second order in Aw. Substitutingeq. (1.5) (to second order in Aw) in eq. (1.4) and doing the integrations over d 3Aw gives theFokker-Planck form of the encounter term^(1)FP ^{f (w)D(Aw i)1+^a2 [f (w)D(Aw iAw^(1.6)dt enc^i=i dwi^2 j_ i uwiuwiD(Awi ) and D(AwiAwi ) are referred to as diffusion coefficients and are defined asD(Awi^Aw)d3Aw^(1.7a)andD(Aw1Aw1) AwiAwjW(w, Aw)d3Aw^(1.7b)The method of Rosenbluth, MacDonald, & Judd (1957) can be used to evaluate the diffusioncoefficients given a set of coordinates. Using eq. (1.6) for the encounter term of the Boltzmannequation, eq. (1.1), gives the Fokker-Planck equation.Self-consistency requires that the Fokker-Planck equation be solved together with the Poissonequation, eq. (1.2). The coupling of these two equations makes the problem highly non-linear andcomplicated to solve numerically. Even this is an idealized situation, since a globular cluster is notan isolated system of self-gravitating point particles as has been assumed until now. To begin with,a globular cluster orbits within its parent galaxy. Its evolution will be directly affected by tidal3stresses induced by both the general gravitational field of the galaxy and by massive constituents ofthe galaxy such as giant molecular clouds. Passages through a galactic disk, or close to a massivebulge, will heighten these effects (Ostriker, Spitzer & Chevalier 1972). The individual bodies inthis N-body system are not ideal point particles of equal mass, but are stars with finite radii and arange of masses. The range in mass leads to mass segregation. Dynamical relaxation leads towardsequipartition of energy between stars with differing masses. The more massive stars will havelower velocities than the less massive stars, and so cannot travel as far from the centre of thecluster. Thus the massive stars are concentrated in the core of the cluster and leave the halo to bedominated by the low mass stars. In the presence of an external tidal field, this leads to the lowmass stars being favoured for tidal stripping.While for relaxation the strong encounters are not as important as the sum of all the weakencounters, the effects of strong encounters involving binary stars will play an important role in thedynamical evolution of a cluster. If these encounters are close enough, then the finite radius of areal star enters into the problem. The tidal forces between two passing stars will induce motions inthe stars' atmospheres. The resulting dissipation of energy may result in the two stars becomingbound, or even merging. For even closer approaches the stars will collide. In the absence ofdynamically formed binaries, there is every reason to expect that some of the original stars aremultiple. The presence of these primordial binaries must be taken into account. Finally, stellarmasses do not remain constant, but change as the star evolves. Adding in each of these complicationsto the Fokker-Planck formalism increases the complexity and time requirements of the numericalmodel. To begin with, the models were simple. Additional physical processes were added as theywere thought necessary or were numerically practical.The models in use here follow directly from those of Cohn (1980) and use his technique ofdirect integration of the Fokker-Planck equation in its orbit-averaged form (see §2a). Other approachesto solving the Fokker-Planck equation utilize a Monte Carlo approach, either with the orbit-averagedequation (116non 1971; Shapiro & Marchant 1978) or directly on the equations of motion with Avand (A v) 2 chosen to statistically satisfy the distribution expected for a Fokker-Planck process(Spitzer & Thuan 1972). StodOlkiewicz (1985), utilizing the orbit-averaged Monte Carlo method,4presented an elaborate model which included stellar evolution and several binary formationmechanisms. His model served as inspiration for much of the work that followed. A hybridN-body/Fokker-Planck approach was used by McMillan & Lightman (1984a; 1984b; McMillan1986). One problem with the Monte Carlo techniques is that the resulting models suffer from highstatistical noise, which make them difficult to compare with observations. N-body studies wouldbe an ideal modelling method, provided sufficient numerical accuracy can be retained. As yet,computing resources are insufficient to run models with realistic values of N . The reader is directedto Spitzer (1987), as well as Elson, Hut, & Inagaki (1987), for a summary of the results of thetheoretical work up to their epoch. In what follows I will discuss some of the relevant workleading up to this thesis, generally focussing on Fokker-Planck models since these are the mostdeveloped and amenable to direct comparison with observations.Section a. Recent Results in Fokker-Planck ModellingBy the mid-1980s, the mechanism of core collapse in the models had been clearly established.What was also clear was that, analogously to the situation with stars, an energy source wasrequired to reverse the core collapse and bring the cluster into an expanding phase if not anequilibrium. The favoured mechanism, binary stars, dates back to H6non's (1961) pioneeringwork. But there were other features to consider as well, including: the actual mechanisms producingthe binaries, stellar evolution, tidal effects, and gravothermal oscillations (GTOs). Many of theseprocesses were elaborated upon within the Fokker-Planck formalism in the succeeding years.Binary stars can heat a cluster via two mechanisms. Assuming the existence of a binary star,either primordial or formed during the dynamical evolution of the cluster, stars within the clusterwill interact with the binary. If the orbital velocity of the binary is higher than that of an averagefield star (ie. a "hard" binary) then the result of the interaction, on average, will be a hardening ofthe binary and an increase in the velocity of the field star. In some cases the field star will beejected from the cluster. "Soft" binaries on the other hand, tend to become even more looselybound and do not last long in the face of these perturbations. For sufficiently hard binaries, there isa good chance that the binary will receive a large enough recoil velocity in a strong encounter to be5ejected from the cluster entirely. The increased velocity acquired by the passing field stars isdistributed to other cluster members through two-body relaxation. On average the velocity of thestars in the core will increase. This is the heating effect. The ejection of binaries from the clusteralso serves to heat the cluster but through a more indirect mechanism; the removal of mass servesto decrease the binding energy of the cluster leading to expansion. For the virial theorem to besatisfied, the overall kinetic energy of the cluster must increase. The net heat input of mass AMescaping from the core is approximately (0 0 — I v02 )AM, where (1)0 and v0 are the central values ofthe potential and velocity. Mass lost from the cluster due to stellar evolution, or for any otherreason, will heat the cluster in the same way.Statler, Ostriker & Cohn (1987, hereafter SOC) used the hard binaries formed by the tidalinteractions of two passing stars to halt core collapse. The models consisted of equal mass stars,plus the binaries they form, and allowed for the heating effects of the tidal binaries, both directlyand via the ejection of the binary and the interacting field star from the cluster. Merging of the starsinto a single object was not allowed for. Much care was taken to account for all the possibleinteractions between the various kinds of stars (single and binary), resulting in a numerical schemeof some complexity. As expected, the newly formed binaries accelerate core collapse, via themass-stratification instability (Spitzer 1969) and the cooling effect of tidal binary formation, but theenergy extracted from them subsequently is able to halt the core collapse and to power a slowre-expansion. The half-mass radius and the velocity dispersion were found to scale as rec t 213 andvrh2oct -2/3 respectively for late times. These scalings are predicted on dimensional grounds (q.v.§2c). In addition, models run with short time steps undergo central oscillations. Whether these aregravothermal oscillations is unclear due to the short time span covered. SOC, as is usually the casewith groundbreaking work, raises more questions than it answers. A series of papers in subsequentyears elaborated on some of these questions.Lee (1987a) added a degenerate component to the model of SOC. Tidal binaries were allowedto form between the normal and degenerate stars and "three-body binaries" amongst the degeneratestars. In the latter mechanism, binaries are formed during the close encounter of three stars. Thethird star carries off the excess energy required to bind the other two. The population of three-body6binaries was followed explicitly in this paper and the heating coefficients were integrated from thethree-body-binary distribution function. Three-body binaries will only be important in high densityenvironments such as the collapsing core. Mass segregation ensures that the degenerates, beingmore massive than the normal stars, will dominate the core, so the environment there is ideal forthree-body binary formation to occur. Since the tidal capture process will not work with degeneratestars, the binary heating will not necessarily be dominated by tidally formed binaries as in themodels of SOC. However, in the models presented by Lee there are many more tidal binaries thanthree-body binaries at the beginning of re-expansion, and the heating is dominated by them. Thehalf-mass radius and velocity dispersion during the re-expansion follow about the same power lawdependences as in SOC. The slope of the surface density or brightness distribution for the evolvedmodels is demonstrated to be sensitive to the ratio of the masses of the degenerates and normalstars, as is expected from equilibrium theory. For a collapsed core, which is ideally an isothermalsphere, the projected surface density profile (SDP) would have a logarithmic slope of —1. In thecontext of globular clusters, a slope shallower than —1 indicates the presence of a large populationof more massive, non-luminous stars which, by mass segregation, are more centrally concentratedthan the stars producing most of the light.Lee (1987b) followed along a slightly different line from SOC by taking the extreme case thatall tidal interactions lead to mergers. The resulting stars are then allowed to evolve off the mainsequence on appropriate timescales, ejecting all of their mass from the cluster. The decrease in thecluster's binding energy associated with the mass shed by the evolving stars is sufficient to reversecore collapse. Three-body binaries were also allowed to form but, unlike Lee (1987a), they werenot followed explicitly. Instead, the heating due to three-body binaries is treated statistically byincluding a heating term in the Fokker-Planck equation. A similar procedure was employed hereand will be discussed in Chapter 2. In the case of Lee (1987b) only stars with equal masses wereallowed to form binaries. The binary subsequently interacted only with stars with the same mass asits components. Again, the post-collapse scaling of the half-mass radius and the velocity dispersionfollow the expected power-laws. In contrast to Lee (1987a), the most massive stars are also thebrightest, so the higher concentration of the more massive stars leads to the light falling off muchmore steeply with radius than does the mass.7Lee & Ostriker (1987) followed up on SOC by adding a tidal boundary. As will be discussedbelow in §2d, the tidal radius in the method of Lee & Ostriker is defined as the radius at which themean density enclosed is equal to a constant, previously specified density referred to as the "tidaldensity". The gravitational potential at this radius defines the tidal boundary in energy space. Theescape rate was taken to be zero at the tidal boundary energy and to increase for energies beyond it.As expected from dimensional arguments, the total mass declines linearly with time during there-expansion phase on a timescale proportional to the product of the total number of stars and theorbital period of the cluster about the galactic centre. The tidally-truncated Fokker-Planck modelswere shown to be similar to King (1966) models, although the Fokker-Planck models had largertidal radii for a given central surface brightness.The emphasis on dynamically formed binaries should not take away from the possible existenceof binaries in the original stellar population (Hut et al. 1992). Sufficiently hard primordial binarieswould serve as an energy source even before the cluster becomes dense enough to form tidal-captureor three-body binaries. Goodman & Hut (1989) examined the likely effects of primordial binarieson globular cluster evolution. While the primordial binaries cannot prevent core contraction (Spitzer& Mathieu 1980; McMillan, Hut & Makino 1990), they can provide the energy to reverse thecollapse. Goodman & Hut show that primordial binaries can lead to cluster with a fairly large coreradius and a core dominated by binaries. The relative importance of primordial binaries will dependon their abundance and destruction rate.Gao et al. (1991) examined this question in more detail. The models they constructed consistedof two components; one for single stars and one for the binaries. The total mass of each binary wastwice that of a single star. The distribution of internal binding energies and the exchange of energybetween the centre of mass motion and orbital motion of the binaries were treated in detail. Themodels were followed well into the post-collapse phase. The main effect of the primordial binariesis to leave the collapsed core with a fairly large radius. The ratio of core to half-mass radii isaround 0.01 to 0.04. The larger core radius suppresses the production of new binaries. Onedistinction between models with primordial binaries and those where manufactured binaries act asthe heat source is that the former have cores consisting predominantly of binaries. This should be atestable distinction between the two scenarios.8Although models with more than one mass of star have been discussed above, the additionalmass classes were not an end in themselves but served as part of the detailed treatment of thereheating mechanism. Within the direct Fokker-Planck formalism, Inagaki & Wiyanto (1984),Inagaki (1985) and Cohn (1985) have produced two species models, while Inagaki & Saslaw(1985) expanded this to 15 mass species. The latter paper emphasizes that clusters cannot ingeneral reach equipartition fully and that a range of masses accelerates core collapse.In stellar population studies, the mass function is usually described in terms of a power-lawdN oc m-(X1-1)din (1.3)When the mass spectral index (MSI), x, is defined in this way, then for the Salpeter (1955) massfunction x=1.35. The initial mass function (IMF) is the global mass function before significantevolution has taken place.Cohn, Hut & Wise (1989) took up the question posed in SOC of the post-collapse instability ofFokker-Planck models. Their single mass species models were reheated by the three-body mechanismfollowing a formalism similar to that of Lee (1987b). Gravothermal oscillations had been observedin gas sphere models, (Sugimoto & Bettwieser 1983; Bettwieser & Sugimoto 1984) but were notseen in the Fokker-Planck models until a small enough time step was chosen. Cohn et al. (1989)showed that the time steps used in previous modelling were too long to react to the gravothermalinstability. The time-step selection technique which had been used chose steps in units larger than104 central relaxation times in the epoch following core collapse. These long time steps served tosuppress the oscillations.Murphy, Cohn & Hut (1990; hereafter MCH) built multi-mass models with seven mass species.The most massive species combined heavy remnants with the upper main sequence to represent anevolved mass function. The re-expansion was powered by three-body binaries and a new expressionwas developed to describe this process statistically when multiple mass species were involved.This extended into the post-collapse phase the study of Murphy & Cohn (1988) who had confirmedthe previous findings that the presence of a range of masses hastened core collapse. MCH investigatedthe post—core-collapse (PCC) behaviour of models using small time steps, and found that multiple9mass species served to suppress the oscillations. This appears to be a result of the dependence ofGTOs on the number of stars (Goodman 1987). Since post-collapse cores will be dominated byheavy remnants, it is the number of these present, rather than the total number of stars, which willdetermine whether oscillations will take place. This conclusion is supported by the observation thatin models with power-law IMFs, those with smaller IMF MSIs (ie. those with a higher proportionof high mass stars) are more unstable to GTOs than those with a larger MSI. In general, it wasfound that the post-collapse surface brightness profiles do not differ greatly from those of core-collapseitself beyond the small radius affected by the GTOs.Lee, Fahlman & Richer (1991, hereafter LFR) produced similar models to MCH but alsoincluded a tidal limit like that of Lee & Ostriker (1987). Their focus was more on the behaviour ofthe mass function during the long-term dynamical evolution of the cluster. Within the half-massradius, significant changes are observed in the mass function. These are seen prior to core-collapse,but over longer times extending into the PCC phase, and depending on the IMF, the MF may eveninvert. Mass segregation will complicate measuring the overall MF of a cluster. It was also foundthat as a cluster approaches disruption by the Galaxy's tidal field, the half-mass relaxation timedrops dramatically. LFR commented on these points with respect to the observations of M71 ofRicher & Fahlman (1989). These considerations led to the more extensive comparisons made byDrukier, Fahlman & Richer (1992) which will be discussed in Chapter 3.Chernoff & Weinberg (1990, hereafter CW) produced pre-collapse models which includedstellar evolution along with a tidal boundary and a mass spectrum. As their models are an independentlycoded implementation of Cohn's (1980) formalism, they provide a useful check on the resultssupplied by all the other work discussed here. Their "quality control" models gave results consistentwith those Cohn (1980) and Inagaki & Wiyanto (1984). As expected from Applegate (1986), theheating effect of the mass loss due to stellar evolution served to delay core collapse in all models.For models with sufficiently low initial concentrations (the initial distribution functions were thoseof King models) or with flat enough mass functions, the coupling of stellar evolution with a tidalboundary was enough to disrupt the cluster on short order. For models which would survive to thepresent they discuss several observable quantities. Since in the advanced stages of evolution the10core is dominated by remnants more massive than the stars providing the light, the surface brightnessprofile may not give a detailed report of the cluster's dynamical state. On the other hand, masssegregation will lead to large effects on the mass function observed with radius.Section b. Comparing Fokker-Planck Models with ObservationsIn all the work described until now the main emphasis has been on the theoretical developmentof the Fokker-Planck formalism. Since the aim of the theoretical work has been to see the effects ofvarious physical processes on the models, less discussion has taken place with respect to theobserved properties of globular clusters. Murphy & Cohn (1988) built pre—core-collapse, multi-massmodels. Since the focus of Murphy & Cohn was to make realistic globular cluster models for thepurpose of comparing them to observations, projected surface densities and velocity dispersionswere presented and discussed. One of the models was compared to the U-band surface brightnessprofile of M15 by Lugger et al. (1987) and a good match was obtained. The presence of the heavyremnants was adequate to flatten the slope of the surface brightness profile from —1 expected for anisothermal sphere to the observed value of around —0.8. MCH, as mentioned above, extended thisstudy into the post—core-collapse phase and suggested that that cluster is currently undergoing acollapse rather than a re-expansion, although at that time the data did not go far enough into thecore of M15 for definitive conclusions to be drawn. LFR discussed their results with respect toM71, but no detailed comparison was undertaken.The first major comparison of the models with observations was that of Grabhorn et al.(1992a, hereafter GCLM). In that study surface brightness profiles and velocity dispersions for thecollapsed-core clusters M15 and NGC 6624 were compared with Fokker-Planck models. Thesemodels employ the methodology described in Murphy & Cohn (1988), but use the same approximation(LFR) as used here for the energy input from three-body binaries. The GCLM models have fivemass species, including heavy remnants. The important feature of these clusters is that they havevery high concentrations and so are well suited to probing the core-collapse, or post—core-collapsephase of cluster evolution. M15 is the prototypical collapsed-core cluster, its core possibly beingresolved with HST (Lauer a al. 1991, but see Grabhorn et al. 1992b). The main result of GCLM11is that both clusters can be fit by PCC models. The surface brightness profiles are best fit when themodel is in a state of maximal expansion during a gravothermal oscillation. Their preferred modelshave an IMF with x=0.9. The core in both clusters is dominated by non-luminous remnants withmasses up to 1.4 Mo . For M15, Fahlman, Richer & VandenBerg (1985) measured the massspectral index to be x =1.5 at 8.7'. More recently, Durrell & Harris (1992) measured x =1.0±0.4.These local values for the MSI are consistent with the global value in the model of GCLM. Inaddition, the gravitational potential in their M15 model is sufficient to explain the observed perioddecrease in the millisecond pulsar PSR 2127+11A (Wolszczan et al. 1989). GCLM serves as auseful sign post on the road to understanding the relationship between the models and the state ofreal globular clusters. The models are clearly capable of reproducing the surface brightness profilesand velocity dispersions of collapsed-core clusters. With some success achieved for clusters withcentral cusps, it becomes all the more necessary to test the models against clusters with a broaderrange of properties.With the theoretical models predicting the striking effects of core collapse on globular clusters,some question arose as to the whereabouts of all the clusters with collapsed cores. This questionbecame even more pressing when the accelerating effect of a mass spectrum was considered. Themodels predict that a central brightness peak would be the sign of a collapsed, or collapsing, core.The post—collapse models also predict a central power-law cusp in the stellar distribution. A surveyof globular cluster cores (Djorgovski & King 1986, Chernoff & Djorgovski 1989) has nowuncovered many cases of clusters showing this sign. In this survey the clusters are divided intothose with surface brightness profiles fitting King (1966) models, and those which do not, butinstead have power-law profiles. About 20% of the clusters they classified have power-law cusps.Deep mass functions have been observed for a number of globular clusters (Drukier et al. 1988[M13]; Fahlman et al. 1989 [NGC 6397]; Richer et al. 1990 [M13, M71]; and Richer et al. 1991[co Cen, M5, and NGC 6752]). All show increases in the MSI for stars less massive than about0.4 Mo .Velocity measurements have also been made for many clusters. Like much of the surfacebrightness and mass function data, the velocities have been used as an adjunct to fitting King12models. Examples of such velocity measurements are the pioneering study of M3 by Gunn &Griffin (1979) and more recently Peterson et al. (1989) [M15], Meylan & Mayor (1991) [NGC6397], and Pryor et al. (1989) [NGC 6624, M28, and M70].Radial colour gradients have also been observed in several central-cusp clusters (Djorgovski,Piotto, & King 1988; Piotto, King, & Djorgovski 1988; Bailyn et al. 1989; Djorgovski et al.1991). The increasing detection of X-ray binaries, millisecond pulsars and other compact objects(Bailyn 1991) as well as binaries in general (Hut et al. 1992) also have implications for globularcluster dynamics Clearly, the stellar population in a cluster is affected by the dynamical evolutionof the cluster, especially by the high densities accompanying core collapse.The motivation for the work to be described here comes from consideration of the observationaland theoretical state of globular cluster studies. The Fokker-Planck models have come a long wayfrom their initial state of idealization. The addition of various binary mechanisms, the tidal boundaryand stellar evolution have brought the modelling much closer to including all the physics presumedto be acting in real globular clusters. Still, simplifications remain in all these elaborations. What themodelling has provided are hints as to what observable properties are most sensitive to dynamicalevolution. The surface density profile (SDP) or surface brightness profile will give information onthe status of the cluster with respect to core collapse. The degree to which such profiles will reflectcore-collapse is dependent on the relative proportion of non-luminous, heavy remnants present.Mass segregation causes these to sink to the centre and reduce the slope of the profile of theluminous stars. By using star count surface density profiles rather than surface brightness profiles,questions with respect to the mass-to-light ratio of the evolved population can be avoided. Velocitydispersions give more direct information on the remnant content and the degree to which corecollapse has progressed. The radial variation in the mass function gives information on the degreeof mass segregation. Since, as will be shown below, the course of the evolution is dependent onthe IMF, observed MFs serve to constrain the range of plausible IMFs and hence the range ofpossible histories for the cluster.It is with these ideas in mind that this thesis has been undertaken. There were two phases tothe work reported on here. In the first phase, Fokker-Planck models which included (i) three-body13binaries, (ii) a mass spectrum, and (iii) a tidal limit, were constructed in an attempt to matchobservations of M71. The observations involved consisted of SDPs for two different ranges ofstellar mass, and MFs for two regions of the cluster. These models, however, did not proveadequate to match the M71 observations. The problem with the comparison lies in the contradictionbetween the observed degree of mass segregation and the short relaxation time on the one hand,and the absence of a central power-law cusp to the SDP on the other. In order to try and resolvethis contradiction several types of slightly modified models were produced. The principle result ofthe first phase comes from one of these modified models. That model suggests that an additional,continuous energy source is required in order to prevent clusters from collapsing too early. At thesame time, mass segregation must not be arrested. Two options for this energy source are primordialbinaries and the effects of mass loss during stellar evolution.Phase two followed on from this conclusion. As stellar evolution is more straightforward toimplement, and more within the statistical philosophy of Fokker-Planck models, it was added tothe code. I also felt it useful to investigate another cluster by way of comparison and so analyzedan extensive set of observations of NGC 6397 to produce a SDP and two new MFs. This nearbycluster has a central cusp in its surface brightness profile. A new series of models which includedstellar evolution was then searched for matches to the observations of M71 and NGC 6397.The rest of this thesis is arranged as follows. The details of the models used in both phases aredescribed in Chapter 2. The results of the first phase comparisons with M71 are discussed inChapter 3. Chapter 4 presents the new observations for NGC 6397. In Chapter 5, the modelsincluding stellar evolution are examined with respect to the observations of M71 and NGC 6397.In Chapter 6, I will summarize the conclusions of this thesis and discuss the implications theobservations have for future modelling. In addition, Appendix A provides lists of parameters forall the models discussed in this thesis and Appendix B discusses further the S ,.„--rhd diagramintroduced in Chapter 5.14Chapter 2. The ModelsIn this chapter I will discuss the numerical techniques used in implementing the Fokker-Planckequation for globular cluster modelling. The first two sections will describe the basic formalismintroduced by Cohn (1979; 1980) and its extension to models with more than one mass species.Subsequent sections will be concerned with issues arising from including the three-body binaryheating mechanism, the tidal boundary, and stellar evolution. Along the way, comparisons willbe made with the results of previous studies using similar techniques. Following this will be asection dealing with the initial conditions assumed for the models and a section recapitulating theassumptions made and the limitations of these models. Most of the actual computer code is notoriginal with me but was written by H. Cohn, H. M. Lee and possibly others who I am unawareof. Besides the additions I wrote, I have gone through the code to check for consistency and toimprove its efficiency. I wrote all of the routines required to analyse the output of the models .Section a. Basic FormalismThe code used in this work is based on that of Cohn (1980). The method for numericallyintegrating the isotropic Fokker-Planck equation is discussed in that paper. Three major assumptionswere made in arriving at this technique. First, it is assumed that the timescale required forsignificant changes in the energy of the stars, the relaxation time, is much longer than thedynamical time. It is then advantageous to consider separately the changes in the distributionfunction on the two time scales. Over the orbital period of a star in the cluster, two-bodyrelaxation is unimportant and the distribution function is a solution of the collisionless Boltzmannequation (eq. 1.3). By Jeans's Theorem the distribution function is a function of the isolatingintegrals of motion, /i and the long term evolution then acts as a perturbation on the distributionfunction. The appropriate equation to use in this case is the orbit-averaged Fokker-Planck equation.In essence, the diffusion constants are integrated over the region of phase space accessible tostars with integrals of motionThe second assumption made is that of spherical symmetry. In this case, f=fiE At), ie. the15distribution function is a function of the energy and angular momentum in addition to time. Cohn(1979) presented models in this formalism. Unfortunately, it proved to be impractical to run thiscode for long integrations for two reasons. One was the computational expense of the two-dimensional Fokker-Planck equation. The second was the problem that the total energy waspoorly conserved. If, as a third simplifying assumption, isotropy is imposed on the velocitydistribution, then the distribution function is reduced to a function of the energy only and bothproblems can be alleviated. One fewer dimension will obviously reduce computational expenseand in the one-dimensional form the differencing scheme of Chang & Cooper (1970), which isknown to conserve the energy well, is available.A numerical solution of these equations is obtained by first advancing the distribution functionsin time in accordance with the Fokker-Planck equation, eq. (2.6). At this point the potential isinconsistent with the density distribution. The second stage is to re-solve for the potential usingPoisson's equation and the new distribution functions. Since the changes in the potential willalso change the energies, the potential solution is constrained by requiring that the distributionfunctions remain fixed functions of q (eq. 2.2c), an adiabatic invariant. The solution of Poisson'sequation is much more expensive in computer time than the Fokker-Planck step, so four of thelatter were performed for each Poisson equation step. The details of the method, for the equal-masscase, are described in Cohn (1979; 1980).In common with the formalism as presented by Cohn (1979; 1980) the energy convention inwhat follows has the potential energy positive and falling to zero as r approaches infinity. Todenote this difference from the way the potential was defined in Chapter 1, I will use 4) instead of(13$ for the potential below.Section b. Multiple MassesOne of the first changes required in order to make detailed comparisons with observations isto provide for stars of different masses. This is easily achieved by having the distributionfunction as a function of mass as well as of energy. Although it is possible to write the distribution16function as a continuous function of the mass, from the beginning I will use a set of distributionfunctions, one for each discrete mass species, m i. For one of these distribution functions, fi(E),the number of stars per unit phase space volume with mass m i and energy we can write the1-dimensional, orbit-averaged, Fokker-Planck equation at a constant potential as (Murphy &Cohn 1988):af a---L —^rivf:DE ± DEE -adl4n 2P at - a-'[where the diffusion coefficients DE and D EE are given byDE (E) —64n4G2 in A^i rE (°) pfidE ,DEE(E)= 64n4G2 1nAEm?[q f0E fidE + L•(°) qfidE].(2.1)(2.2a)The function p is the isotropic average of the period of stars with energy E,p(E) = 4r 1(E)r2 vdr, ,^ (2.2b)where the velocity v = [24)(r) — 2412 . The integral of p, q is a similar average of the radialaction:4 r 1 (E) 2 3R'( E) = 3^r v dr .0 (2.2c)The Coulomb logarithm, In A, was taken to have a constant value of ln(0.2N), where N is thetotal number of stars in the initial model. This is to approximate the number of stars in the core.The boundary conditions on eq (2.1) require zero flux of particles past the energy limits.Since a globular cluster is a self-gravitating system we need to use Poisson's equation inorder to get 4)(r) and its inverse 4) -1 (E). Poisson's equation is used in its integral form,4 G r4)(r) = irr So p(r')r'2 dr' + 47cGirp(r')r'dr'.^(2.3)The total density p(r) is the sum over all mass species of the individual density integralsp (r) = 2 512 TCMi r (r) [00 —^fidE^(2.4)The boundary conditions on Poisson's equation are that the potential goes to zero in the limit17r--->. and that its first derivative is zero at r=0.In Cohn's method the equations are solved in a dimensionless form. The three dimensions arechosen so that a scale radius, r0, the total mass, Mt, and the gravitational constant, G, are unity.The last is done by using a scale velocity (or energy per unit mass)[v] = v =k—r°^GM,^ (2.5)For this choice of dimensions, the natural time unit is ro /vo but C,(ro /v 0 ) is used instead. Thedimensionless constant C, will be chosen below in order to eliminate certain constants from thedimensionless form of the Fokker-Planck equation. A further complication is that the masses perobject of the individual mass species, m ;, are more easily normalized by C„,M, (where thedimensionless constant is chosen such thatC„,M1 is of order 1 M o ) rather than M. Acompensating normalization 1/C,„ is applied to the distribution function to account for this. Usingthese scalings eq. (2.1) becomes,al^2^G2/W 1 a [ -= 16m In AC „,C, 2 4^+ EE^(2.6)ro v^aE^aEwhere the tildes indicate dimensionless quantities. The dimensionless diffusion integrals iE andit.E correspond to the diffusion coefficients defined in eq. (2.2a) and are given byE -==andIEE ,^ :^ tmyj 4fidEi. (2.7)For ease of calculation we may require the constant coefficient on the right side of eq. (2.6) to beunity using the constant C e. From the definition of v, G2mr2r0-2 equals vo4 , so we are left with167r2 ln ACC,. This can all be cleared away by takingC, = (16relnAC,,,) 1 .^ (2.8)After undergoing the extensive revisions described below some check runs were done tomake sure that nothing important had been changed. A single mass model starting as a Plummer)1/218model was run and compared with the results in Cohn (1980). Provided that a similar algorithmfor time step selection is employed, the resulting model was consistent with Cohn's. Two of thedimensionless quantities discussed by Cohn will be compared to the results from my model: forthe fractional rate of increase in the central density per central relaxation time 4 = 3.4x10-3 (cf.3.6x10 -3 for Cohn), and for the quantity y = d In v 2 (0)/cilnp(0)= 0.1 (Cohn's value). After stellarevolution (see §2e) was included in the models a somewhat different scheme was used for timestep selection. This leads to somewhat smaller asymptotic values of and y. I also ran twospecies models similar to those of Inagaki (1985) with stars having a mass ratio of 2:1 and arange of ratios in the total mass WM 1 . These gave similar ratios of core-collapse to initialhalf-mass relaxation times, listed in Table 1. The small differences are attributed to the differencein the time step algorithm with the somewhat larger time steps used here giving a slower corecollapse.Table 1: Core collapse times for two component models as in Inagaki (1985) Mass ratioMJMl tjtth(Inagald)15.413.610.28.59.613.60.000.010.050.111.09.015.814.810.78.89.314.8Single component modelSection c. Binary HeatingIn order to halt core contraction an energy source is required. In this study I have usedbinaries formed by three-body interactions ("three-body binaries") as the heating mechanism.While the first papers to deal with the reversal of core collapse used tidal binaries, from the timeof Lee (1987b) almost all the Fokker-Planck studies have used the statistical treatment of three-bodybinaries. There appear to be two main reasons for this change.19The first is that the overall evolution of the models is not highly dependent on the details ofthe heating mechanism. The expansion is governed by the rate of energy transfer within themodel. This can be seen qualitatively by the following dimensional argument which is based onthose of SOC and Lee (1987a).Assuming that there is an energy source in the core of the model, the rate of change inenergy, for a virialized system, is given bydv 2"Eco °c^• dt (2.9)It is assumed that the total mass is constant. The energy conduction rate outside the core isproportional to the product of the local rms velocity, v, and the kinetic temperature gradientd(mv 2 )vcond v^^ •drApproximating v locally by a power law in r, with a constant stellar mass, m, then &rid a v 3/r.Of this, r/v is a time, so for example at the half-mass radius, rh, we can take this conductiveluminosity to beV2Econd^rhtrh(2.10)where vrh is the velocity dispersion at the half-mass radius, and trh, the half-mass relaxation time(Spitzer & Hart 1971), is given bytrh = 0.1380/2rh312 (2.11)G 1/21n In AEquating Ecore and Ecohd and employing the relationship between trh and rh from eq. (2.11) and thevirial theorem, one finds that for late times v r2h a t-2/3 , rh a t2/3 and E oc t-5/3 . This result is notrestricted to the half-mass radius and will hold for any central energy source as long as thehalf-mass and central velocity dispersions are proportional.The latter condition is reasonable ifthe re-expansion is approximately self-similar as has been argued on dimensional grounds by, forexample, Inagaki & Lynden-Bell (1983).In order to check the assertion that the evolution does not depend strongly on the details ofthe heating mechanism, I have run a model with identical initial conditions to that of SOC but20with the three-body binary heating mechanism instead of the tidal binary method they employed.The behaviour at the half-mass radius is much the same in the two models. Core collapse occursat the same epoch, but it is much deeper for the model with three-body binaries. Further, thesubsequent re-expansion of the core takes place much more quickly. The latter difference can beexplained in terms of the energy argument discussed above. SOC show that for tidal binariesE cc pc2v0.8r: where p c and v2 are the central density and velocity dispersion and r, is the coreradius defined by (King 1966):3^3v2r = ^c . (2.12)47rGpcEmploying eq. (2.12) and the proportionalities given above for E, r and v, this heating rate givespc cc r°.8 . For three-body binaries Lee (1987a) shows that E cc p3cvc-71-c3 , and thus pc cc t-2 . Thedetailed structure of the cluster centre will depend on the heating mechanism, but further out thenature of the energy source is less unimportant. The actual numerical models have power-lawexponents which deviate a little from those predicted, but this is probably due to their not beingperfectly self-similar.The second reason for the change to three-body binaries, is that the statistical three-bodymethod is simpler to implement and less demanding of computer resources than the other methods.The tidal binary approach of SOC requires complex bookkeeping to handle all the possibleinteraction channels. This would become even more complex when all the possible interactionsbetween stars of different masses, including mergers as in Lee (1987b), must be considered. Lee(1987a) treats three-body binaries in detail in addition to tidal binaries. Again, the requiredformalism is quite complex.An additional support for using three-body binaries, as opposed to tidal binaries, comes fromconsideration of the stellar population in the core. Due to mass segregation, the core will bedominated by the most massive stars. At the present epoch, these will be the degenerate remnantsof stellar evolution. Such remnants do not form tidal-capture binaries, so three-body interactionswill be the favoured binary formation mechanism under certain conditions. Whether this is thecase will depend on the relative proportions of normal stars and remnants in the core and the21interaction cross sections. In the models of Lee (1987a), for example, the tidal binaries dominated.The statistical three-body approach is founded on two assumptions. First it is assumed thatthree-body binaries are formed, contribute their energy and are destroyed instantaneously. Further,it is assumed that a statistical estimate of the binary formation rate will serve as a suitablestand-in for the more realistic stochastic rate of formation and destruction. The energy input atany given time is just the product of the creation rate for hard binaries and their average energyinput.To be consistent with eq. (2.1) a heating term of the form 471 2 (k) is required, where theangle brackets indicate the orbit averaging operator, and E is a heating rate per unit mass. Basedon suggestions by P. Hut and D. Heggie, LFR useE(r) = CHG5vc2(E i (2.13)for the heating rate per unit volume due to three-body binaries. For mass class i, the mass densityis given by eq. (2.4) and Vi2 is its three-dimensional velocity dispersion. Also, v c2 is the centralvelocity dispersion. This formula may be thought of in the following manner. In order for threestars to interact strongly enough that the velocity perturbation is sufficient to bind two of themthey need to pass by each other within /v2.. The time interval between encounters out tothis distance is about (np2 v)- 1 . The probability of a third star being within distance p is pan, so thetime interval between triple encounters, at distance p or less, is (n2 p5v) - 1 . Substituting for p, theformation rate will go as G 5n2m5/v9. The energy generated per binary is CH MV: so the energygeneration rate is CHG5vc2n2m6/V9. Converting this to a rate per unit volume by multiplying by nand expressing nm as p yields CHG5ve2p 3m3/0. Since we must consider all possible combinationsof three stars from all the different mass bins it seems natural to partition this rate asCHG5v,2(p im i/O(p 2m2/v23)(p 3m3/v33) for each star involved and to sum over all possiblecombinations. This is just the extended binomial series expressed in the cubic factor of eq.(2.13).The energy release is governed by the factor CHvc2. The constant used, CH=4210 for athree-dimensional velocity dispersion, is based on the preliminary report of Hut (1985). The22detailed analysis of the gravitational scattering experiments from which this number was derived(Goodman & Hut 1992), settled on a creation rate which would give a heating rate of about 80%of that used here. The smaller heating rate implies a slightly deeper core collapse before re-expansionbegins, but otherwise should not affect the models presented in this thesis.For dimensional consistency, E must be divided by a density to get a heating rate per unitmass. If the density used is the local mean density then the energy is distributed to each massclass in proportion to its relative density at each radius. The greatest heating comes in the core atcore-collapse. Following core collapse the models remain highly concentrated and energygeneration continues.The Fokker-Planck equation, augmented by the heating term, may be writtena a^afi4n2P=W•Lmi(fiDE +167C2H^L,EE^,dEwhere the heating term is given byELL-13k 3 dr.H E(E) CHG5vc^(E) r v1.1 2^0p^vi(2.14)(2.15)In the dimensionless form used in the code, the Fokker-Planck equation including the heatingterm iswhereand the constant CL is given bya:4a [^, 2 -=^ivrt E +C IH )+Iai^at Hv^EE aill (Ell E 1 (E) r-p-213 2.17^3dr-,(2.16)(2.17)CHC2CH = 4n2 In A .^(2.18)The heating coefficients were calculated following each solution of the Poisson equation andwere used for all four of the subsequent Fokker-Planck steps. This is in contrast to the diffusioncoefficients which are calculated anew for each Fokker-Planck step. The difference is that theheating coefficient depends on the density and velocity dispersion, which is only known following23a Poisson step, while the diffusion coefficients are functions of the distribution function itself.The heating coefficients change most rapidly when the model is undergoing core collapse butthis is also the time when the time steps are the smallest, so this should not be too much of aproblem.Section d. Tidal BoundaryOne of the drawbacks of using the orbit-averaged form of the Fokker-Planck equation is thatit does not properly treat the evaporation of stars from the cluster (Henn 1960). The inclusion ofthe effects of tidal stripping will mitigate this aspect of the orbit-averaged formalism. What itwill not correct for is the fraction of the evaporation due to ejection from the core as opposed tothat due to diffusion past the edge of the cluster.The assumption of spherical symmetry prevents a thorough treatment of the galactic tide, theeffects of which on the cluster are distinctly asymmetric. Further, the assumption of isotropy maybe violated as the tidal field has more time to act on stars with low energies and high angularmomenta in the halo and preferentially removes stars from the core with high energies and lowangular momenta. As well, stars with prograde orbits are easier to strip. This leads to some netrotation, in the halo of the cluster at least. As shown by Oh and Lin (1992) the overall effect oftides is to make the velocity dispersion isotropic in the halo of the cluster. They predict that atransition zone of anisotropy in the intermediate region between the core and the halo may exist.The Oh and Lin result suggests that the assumption of isotropy will not be in extensive conflictwith tidal stripping.The method used for implementing tidal stripping follows that of Lee & Ostriker (1987) andassumes a constant tidal field with the globular cluster in a circular orbit. Rholfs & Kreitschmann(1988) point out that, to within 10%, all the proposed models of the galactic mass distribution,MG(r), obey the relationMG(r) = 2• 3 x10 -5( r j( e(r) 10 10 M 0^kpc km s -1(2.19)where 8(r) is the circular velocity at radius r in the galaxy. Innanen, Harris, & Webbink (1983)24found that for a galaxy with a flat rotation curve (C1(r) constant), the tidal radius for a clusterwith mass Mc at a distance R G from the centre of the galaxy is given byr =- 2 Mc R. (2.20)2M G3Substituting eq. (2.19) into eq. (2.20) and taking a fiducial cluster mass of 10 5 M0 givesM^e(RG ))1/3 ( 2/3 ( RG )43 (2.21)r =10.9( 105 M0^220 km s -1pckpcLee & Ostriker (1987) considered the rate of mass loss for stars with energies beyond thetidal boundary. Mass is removed from the distribution function for energies greater than E, (thegravitational potential at re when 4:$(.) = 0). The rate of mass loss, b(E), is based on the equationof motion for a particle outside the tidal boundary. In this case, for tidal strippingafi (E,t) = CTfi (E,t)b(E)t; 1 ,^ (2.22)atwhere 3121 — (f)^E Et0E> Etb(E)= (2.23)andt, —^2ic^(2.24)1Gp eis the orbital period of the cluster about the galaxy. The tidal radius is that where the enclosedmean density equals the tidal density p r. From eq. (2.20), pt is clearly constant. Cr is a dimensionlessconstant controlling the overall rate of mass loss and was always set to unity. In that case themass-loss time scale is equal to the orbital period of the cluster about the galactic centre.Equation (2.22) is applied in its dimensionless integrated form,j; (k, i + AO= jg,i)expf CT/Airi5 b(f) },16n2 ln AC„,(2.25)by assuming that the change in the tidal energy over the time interval At is small.This treatment of the tidal boundary raised several concerns. The first derivative of the tidalstripping rate is discontinuous at the tidal boundary, an unphysical situation. In Cohn's method25for solving the Fokker-Planck equation, the distribution functions are regridded in energy aftereach iteration of the Poisson equation solver. This is necessary to keep f a constant function ofthe adiabatic invariant q. The regridding is done by way of a second order Taylor expansion offin q(E). It can happen that discontinuities at the tidal boundary are reflected in the derivatives off(E) taken for the Taylor expansion and are then amplified in transferring the distribution functionsto the new energy grid. A few iterations of this can lead to a catastrophic failure of the Poissonsolver. In order to reduce the probability of this happening a cubic polynomial for the strippingrate was used over the boundary region. Let13=14,^j3^- . ^ (2.26)E,From eq. (2.23) the stripping rate goes as b(E) --= ;13 for E_ Et . For the region where 1131 S e , andfor a suitably chosen e, I find a cubic polynomial b1 03)satisfying the boundary conditions:b1 (—e) = 0,bi'(—E) = 0,bi (e) = Are-,1b1(e) = 2-1i . (2.27)The resulting cubic is3^2A re- ^a +5 a+ 3] .^k (13) = --i-[—(-vi) + (-vi)^Le^(2.28)This relation was used in eq. (2.25) instead of eq. (2.23) when Pi e . The transition regionconstant e was chosen so that two or three grid points surrounding the tidal energy would havestripping rates based on eq (2.28).One advantage of using a tidal limit is that it provides a natural outer boundary to the model.The boundary conditions require that the distribution functions go to zero at low energies. Thistranslates as the density going to zero at large radii. The expansion of the model due to relaxationdrives a population of stars out to ever larger radii. When no tidal boundary is applied, the radialgrid must be continually extended to encompass these stars. If this is not done, then eventuallythe distribution function at small energies will become large enough to contradict the boundary26conditions. With a tidal boundary, the stars at low energies will be continually removed. Oneexception to this is when the evolution of the cluster is much more rapid than the tidal strippingtimescale. In that case the grid edges may still be overwhelmed. I found that using a radial limit50 times the initial tidal radius was adequate to guard against this.On the other hand, tidal stripping will continuously reduce the mass of the cluster, so caremust be taken to adjust the limits of the energy grid appropriately. The minimum energy shouldbe taken as the current mass divided by the radial grid limit.Since the tidal stripping rate at a given energy depends on the current density profile of themodel, it can only be applied following a Poisson step despite the fact that the stripping is doneon the energy grid. One straightforward way of doing this is to simply do the tidal strippingfollowing each Poisson step before proceeding on to the next set of Fokker-Planck steps. Thiswas the approach taken by Lee & Ostriker (1987). The main problem with this is that once thetidal stripping is done the potential and densities will again be inconsistent with the distributionfunctions. If the amount of tidal stripping taking place is small, this won't present much of aworry. In the late stages of the model's evolution the mass becomes small, and since i; .. M-2/3 fora linear decrease in total mass, the rate of decrease in the tidal radius increases dramatically. Incases where the rate of decrease in I-, is significant, it becomes necessary to find a self-consistentsolution.To account for this, an iterative method of doing the tidal stripping and simultaneouslysolving Poisson's equation was used. Following the completion of the Fokker-Planck steps, thedistribution function, the functions p and q, and the previous self-consistent potential are stored. Iwill refer to this group of functions as model F and this initial model as F°. One iteration of thePoisson solver is done to estimate the tidal radius, rt°, and the tidal energy, E:D . Since the Poissonsolver changes F, F is reset to F° . Next, tidal stripping is done using E1° in eq (2.25), followed bya full solution of Poisson's equation. This gives model F1 . Based on model F1, a new tidal radius,r,1, is computed. If rf° and ri' are equal, then model F 1 is self-consistent both between thedistribution function and the potential (as a result of the Poisson solver) and with respect to thetidal boundary. If r1° and rel are not equal then the model is reset to F ° , stripped using tidal energy27E1 1 based on the tidal radius r t1 , and solved to give model F2, tidal radius r te, and tidal energy E t2 .This procedure is repeated until r7 equals and model is adopted as self-consistent. Afractional difference of 10" 3 wasdeemed "equality" in the comparison of r7 and r7 1 .The solution of Poisson's equation is the most computer intensive stage of the modellingprocedure, so it is inefficient to have to do too many iterations of the tidal stripping. If more thanone full solution of Poisson's equation was required to get r„ then the time step for the nextFokker-Planck step was arbitrarily restricted to being no more than one-half the time step justused. This procedure worked well and prevented several numerical oddities which crept intomodels run without it.CW found that when tidal stripping became rapid it was impossible to find a self-consistentsolution to the Poisson equation and their tidal boundary condition. They took this to be evidencefor an instability and confirmed the existence of such an instability using an N—body calculation.CW conclude that the nature of the instability is that the work done on the cluster by the stellarevolution mass loss is sufficient to dynamically disrupt the cluster. In my models, it was alwayspossible to find a self-consistent solution, but the mass and tidal radius go to zero in finite time.The difference in the behaviour between the two models results from the different approachesto the tidal limit. The limit of CW is much sharper with f=0 for E>E„ With E dependent on thepotential, they find it more useful to truncate f(q) at some value q0 , where q is the adiabaticinvariant defined in eq. (2.2c). In order to implement this they need to find a value of qo forwhich Alc [q0]=r,[q0] 3 in appropriate dimensionless units. The problem they run into is that asmass is lost, the implicit relationship between Mc & r as a function of qo for a given f(q) becomesparallel to Md=r,3, and the boundary condition is ill-determined. For the boundary condition ofLee & Ostriker (1987) there is no discontinuity in the density at the tidal boundary, so aself-consistent solution is easier to find. The difference here is numerical; the physical instabilityappears to be the same for both models.28Section e. Stellar EvolutionWith the exception of CW, the effects of stellar evolution have not been included in anyorbit-averaged Fokker-Planck models. StodOlkiewicz (1985) presented the results of an elaborateMonte-Carlo calculation which included a mass spectrum, tidal and three-body binaries, andstellar evolution, but stellar evolution was not an important effect in his models. For simplicity,as in LFR, it has been generally assumed that the initial state of a model was at some stage afterthe massive stars have completed their evolution and that the further evolution of lower massstars was of little importance. When this is done, degenerate remnants are included to representthe evolved massive stars, and the masses of the various mass bins is assumed to be constant. ForFokker-Planck models, which propose to model the full time evolution of a globular cluster, astatic mass spectrum is inconsistent. As will be seen below, the neglect of this fact of stellar lifecan lead to incorrect conclusions about cluster evolution.The approach used in implementing the effects of mass loss due to stellar evolution is muchthe same as that used by CW. Each mass species is assigned an initial mass based on an averageacross the stars in the bin. Considering the likely outcome of stellar evolution on stars of thatmass, a final mass is also assigned. The time range over which the stars in that species go fromtheir initial mass to their final mass is set by the main sequence lifetimes of the stars at the edgesof the bin. The mass per object for that mass species is taken to change linearly from the initialmass to the final mass over that time span. Given the degree of simplification used here it wasnot necessary to be any more exact in this scheme. Care was taken during the running of themodels that the total mass of the cluster should not change due to stellar evolution by more than1% over any time step.Let mii and mf be the initial and final mean masses for stars in bin j, and andan ti be the timerange over which stellar evolution takes place. Then, the mean mass per object for bin j at time tis given by(t) = (Am ) (t — tiAtt^t ij1 < t <'t >(2.29)29where^( Am ) rriC^(2.30)At );^tCis always negative for mass loss. The rate of change of mass with respect to time during stellarevolution was transformed into dimensionless units for consistency with the rest of the code.The effects of stellar evolution will depend on the ratio of the stellar evolution time to therelaxation time. From eq. (2.11) tth.rh'. The stellar lifetimes are independent of any of the otherscales. Hence, for otherwise identical models, models with larger scale radii (eg. half-massradius) will suffer more mass loss per relaxation time than those with smaller scale radii. Thiseffect does not occur with losses through the tidal boundary since Os also proportional to r' (cf.eq. 2.20 and 2.24).In §2b, I discussed the change of variables used to make the models dimensionless. Of thethree dimensions involved, one is set by taking of G=1. When binary heating was introduced asecond dimension, that of mass, was set by the relationship between the relaxation time and theheating time. This can be seen from eq. (2.16) where the constant CL. (eq. 2.18) depends on thetotal mass via the lnA factor. When stellar evolution is included, the final free scale, that oflength, must also be chosen. With all the scalings needing to be predetermined the usefulness ofa dimensionless model is lost, but the formalism has been retained for backward compatibility. Itshould be emphasized that in the models including stellar evolution, there are no scale parametersleft for fitting to the observations.The stellar evolution times and the final masses for stars were based on those used by CW.This was to enable me to do direct comparisons with their work. The mass bins were chosen witha logarithmic spacing with the initial mass per object per bin being the geometric mean of themasses of the bin boundary. The mass within a bin, Mt = JN(m)m dm, was integrated using apower-law mass function (eq. 1.3) between the mass bin limits and the number of stars in the bin,Ni , is defined by Ni4 1 1mi .The stellar lifetimes are based on a cubic spline of the Miller & Scalo (1979) compilation forPopulation I stars. The final masses for stars with m i<4.7 Mo , which become white dwarfs, were30from the formula of Iben & Renzini (1983) with TI= 1/3: mf=0.57+0.22(m/-1). More massivestars, to 8.0 M o , were assumed to be completely disrupted in the supernova, leaving no remnant,while those more massive than this limit were assumed to leave a 1.4 M o neutron star behind.For bins with initial masses less than 0.8 M o , which cease evolving at about 15 Gyr, noevolution was included. In other words, stellar evolution was assumed to terminate at 15 Gyr. Bythis time the amount of mass loss due to stellar evolution has dropped to a low level and is oflittle importance to the dynamical evolution of the cluster. To demonstrate this, consider a simplemodel where the mass spectral index is x, the main sequence lifetime scales as m -3 and a star withinitial mass m loses mass 6m through stellar evolution. Then the rate of mass loss is fit 0= 6m m3-".If we compare 0.9 M o stars evolving into 0.6 M o white dwarfs (&n=0.3 Mo) at 15 Gyr with 5Mo stars being completely destroyed (8m=5 M o ) at 80 Myr, then for each type of star to give thesame contribution to the mass loss rate, x=4.6 is required. Mass loss is never important for such asteep MSI. For x=1, the mass loss rate for the 0.9 Mo stars is 2 x 10 -3 that of the 5 M o stars,By way of comparison of this implementation of the Fokker-Planck equation with theindependent Fokker-Planck code of CW, I will present the results of one specific run. Theparameters for this model, designated CW1, are listed in Table Al in Appendix A. This model isanalogous to a model with W0=7 and a MSI x=1.5 in family 3 of CW. (Note that CW specifytheir power-law MFs in units where the Salpeter MF has slope 2.35. W 0 is the dimensionlesscentral potential parameter of King models. A higher value indicates a more concentrated cluster.See §2f.) Since CW have not chosen units where G is unity they retain one degree of freedom inthe scaling which they use to define their various families. In my models no such freedom existsso a specific total mass and tidal radius have been chosen.Model CW1 reaches core collapse at 46 Gyr while that of CW reaches core collapse at —36Gyr. This difference is attributable to the fact that they take In A.--..1nN while I use In A=ln(0.2 N).This corresponds to a situation where most of the relaxation takes place amongst the stars in thecore of the cluster. For this model the difference in In A makes their relaxation time about 15%shorter than mine, so a 20% earlier core collapse is not unreasonable.31Figure 1 can be compared with Fig. 23 of CW, and shows the evolution of the projectedmean mass as a function of projected radius. I have taken snapshots from my model which haveabout the same time distribution with respect to the times of greatest expansion and core collapseas in CW. The various curves are marked with their times in Gyr. Greatest expansion occurs at8.3 Gyr. The overall morphological behaviour is quite similar although the central mean massnever becomes as small as it does in CW.Figure 2 shows the projected surface densities of the various stellar populations at a time nearcore collapse. The time shown is for when model CW1 has about the same total projected centraldensity as the model shown in Fig. 25 of CW. The solid, short dashed, long dashed, and dottedlines are the projected surface densities for the unevolved stars, the neutron stars, the whitedwarfs, and all stars respectively. The overall profile is in agreement with CW but the relativeproportions of the various stellar populations is different. The reason for this will be explained inconnection with Fig. 3.Figure 3 shows the evolution of the mass function for comparison with Fig. 33 of CW. Themorphological development is much the same with the exception of the discontinuity in the MFfor masses less than 1Mo at late times. It appears that CW do not terminate stellar evolutionwhen the 0.8 Mo stars evolve, but continue it to lower masses. This difference also explains whyI have a higher proportion of luminous stars than CW do, and why the mean mass in model CW1never gets as small as does the mean mass of CW's model.I also ran a grid of models using the same stellar evolution times and final masses (for starswith m i>0.8M 0 ) as CW, with Mi=105 Mo ,W0€ 1,4,7), rrE { 15 pc, 60 pc) and XE (0.5,1.5,2.5).The results for this set of runs are shown in Table 2. All the models with W0=1 or x=0.5 weredisrupted. The remaining W0=7 models all reached core collapse. (The model with x=2.5 andre=l5pc has such a small amount of mass loss during even the first few relaxation times that theenergy input does not balance the energy loss from the core due to two-body relaxation; thecore-collapse proceeds unimpeded.) The W0=4 models present an interesting intermediate situation.The trends in this set of models is similar to that noted by CW. All in all the agreement with theresults of CW is satisfactory.321.2440A 0 . 8V0.60.427I —2 0log r (pc)1 2Fig. 1.—Projected mean mass as a function of radius for model CW1. This figure may becompared with Fig. 23 of CW. The curves are marked with the time in Gyr from the start of themodel. The minimum in the central density occurs at 8.3 Gyr.331—2^—1^0log r (pc)02C)0bb00---,Fig. 2.—Projected surface mass density profiles for the various mass components in modelCW1. The curves show the projected surface density for the neutron stars (short dash), whitedwarfs (long dash), main seqence stars (solid) and for all stars (dotted). The time for this figurewas chosen so that the total central density is about the same as that of Fig. 25 of CW. Therelative proportions of white dwarfs and luminous stars is different than in CW since no starswith masses less than 0.8 M o were allowed to evolve in model CW1.3410.8IIz0.6z0.4020^0.5^1log mi (M0)I^I^I^I^I^I^I^I^I^I^I^I^l^I^I^I1^1^1^I^1^1^1^1^1^1^I^I^I^1^I^I —0.4 —0.2 0^0.2 0 4log rni (M0)zrap0—1.5Fig. 3.—(a) The number of stars as a function of mass normalized by their initial values formodel CW1. The times were selected to have the same normalized number of stars for the leastmassive species as in Fig. 33a of CW. The break in the lines connecting the points differentiatebetween the evolved and unevolved mass classes. The three initially heaviest mass classes havebeen merged into one 1.4 Mc, mass class, and the stars in three next mass classes were entirelydestroyed. (b) The number of stars as a function of mass normalized to the total initial number ofstars. The counts are shown for only the unevolved mass bins. The line without the circlesmarking the masses of the various mass species is the initial number of stars. This figure iscomparable with Fig. 33 of CW.35Table 2: Results of sample grid of models including stellar evolutionInitial mass(M)Wc,^r,(pc)MSI Fate a t rr b(Gyr)i.e.:(Gyr)105 1 15 0.5 D 0.04105 1 15 1.5 D 0.21105 1 15 2.5 D 2.3105 1 60 0.5 D 0.07105 1 60 1.5 D 0.70105 1 60 2.5 D 7.1105 4 15 0.5 D 0.07105 4 15 1.5 D PCC 2.6 2.4105 4 15 2.5 D C 6.9105 4 60 0.5 D 0.20105 4 60 1.5 D 7.0105 4 60 2.5 C >45.0 >45.0105 7 15 0.5 D C 0.54105 7 15 1.5 D C 3.9105 7 15 2.5 D C 6.3105 7 60 0.5 D PCC 6.3 6.0105 7 60 1.5 PCC >29.9 17.0105 7 60 2.5 PCC >23.1 19.84x105 1 23.8 2.5 D 2.84x105 1 95.2 2.5 D 7.84x105 4 23.8 1.5 D 3.04x105 4 23.8 2.5 D C 23.54x105 4 95.2 1.5 D 6.44x105 4 95.2 2.5 C >30.2 >30.24x105 7 23.8 0.5 D C 1.54x105 7 23.8 1.5 D C 12.84x105 7 23.8 2.5 D C 21.64x105 7 95.2 0.5 D C 12.84x105 7 95.2 1.5 C >30.2 >30.24x105 7 95.2 2.5 C >30.4 >30.4aD: Disrupts with expanding core, C: Core is contracting,D C: disruption with contracting core, PCC: Post-core-collapse, willeventually disrupt, D PCC: Disruption in post-core-collapse phase.b^•Time of disruptionTime of maximum central densitySection f Initial ConditionsThe formation mechanism for globular clusters is still a matter of great debate. It can probablybe assumed that whatever the means by which the initial gas cloud fragments into stars theresulting system will not be in dynamical equilibrium. Such a system is expected to undergoviolent relaxation (Lynden-Bell 1967). The velocity distribution for the stars following violentrelaxation is independent of the mass of the stars. A corollary of this is that initially there will beno mass segregation.One common choice for the initial model is an n=5 polytrope (Plummer's model) which hasthe distribution functionf (E)= 24V-27 E3.5 .77E 3 (2.31)The Plummer model has the disadvantage of being unbounded, though of finite total mass. Analternative choice are the models of King (1966). This family of models is based on an approximatesolution of the Fokker-Planck equation in a square-well potential. The distribution function takesthe form of a lowered-Maxwellianf (E) = C(e -BE —1). (2.32)If, using Jeans's Theorem, it is assumed that the same form for the distribution function applieseverywhere in the cluster, this distribution function can be used in a solution of Poisson'sequation to give a full cluster model. The King models form a one parameter family of dimensionlessmodels. The parameter usually used is the dimensionless central potential, Wo . These models arenaturally truncated and cover a range of concentrations depending on the value of W 0, with theconcentration increasing with Wo . A more detailed description can be found in King (1966). AKing model was used for the initial distribution function in the models to be described here.The mass spectrum was divided up into 15 to 25 mass bins, generally using a logarithmicspacing. The IMF was then determined by the number of stars placed in each bin. While a singlepower-law is the most common parametrization of the IMF, models were run with other ways ofsetting the IMF. In the models described in Chapter 3, where stellar evolution was not included,the IMF was taken to represent an already evolved population of stars containing dark remnants37along with unevolved, luminous stars. For the runs involving stellar evolution, specifying theIMF is much simpler, but an mf(m i) relationship and a set of stellar evolution times must beadopted. I used the relations used by CW and described above.One problem with using these models for comparison with observations is the large parameterspace involved. The structural parameters of the initial model are: the dimensionless centralpotential, Wo, of the King model, which gives the initial distribution function and potential; theinitial total mass, which sets the relative rates for relaxation and binary heating; and the tidalstripping timescale, Cr The stellar populations are defined by the initial grid of bins and theIMF. The inclusion of stellar evolution also requires the specification of the final mean mass foreach bin and the time interval over which stellar evolution takes place.Section g. CaveatsBefore going on to discuss the applications of these models to actual observations it is usefulto summarize the potential problems with models of this type. These Fokker-Planck models areby their very nature numerical approximations to statistical approximations to the dynamicalevolution of a globular cluster. They assume that the number of stars is sufficiently large thatthere are a reasonable number in each bin of the distribution function. The inclusion of the tidalboundary ensures that if the simulation is run long enough this assumption will be violated. Thevery high central densities sometimes reached may also indicate a breakdown in the approximationsused. Since, as described by Makino & Hut (1991), the core mass decreases with decreasing coreradius in self-gravitating systems, we see the curious situation of the number of stars in the coredropping to zero as the central density goes to infinity. Since it is the density and velocitydispersion as a function of radius that go into the heating rate, there may be too few stars withinthe core to even form a three-body binary. In a certain sense, three-body binaries may not besufficient to reverse core collapse in these cases. Heating due to tidal and primordial binaries,which have not been included in these models, still leaves a way of reversing core collapse.If Fokker-Planck models can be considered as averages over possible clusters and the discreteevents in their histories, as the distribution function is often considered a probability density,then it is reasonable to use them in comparisons with observations despite these concerns.38Chapter 3. The Problem of M71Section a. Outline of the problemThe work described in this chapter has previously been published in Drukier, Fahlman, andRicher (1992, hereafter DFR). While the observations of M71 described in §2 of DFR are theresponsibility of my co-authors, the models and their analysis, which are described in the remainderof the paper, is my work. This chapter will provide a summary of DFR in order to lay thegroundwork for the remaining chapters of this thesis. All of the models which are presented in thefigures of this chapter have been rerun using the final version of the Fokker-Planck code. Somesmall differences will be seen compared to the diagrams in DFR, mostly due to the self-consistenttreatment of the tidal boundary and a somewhat differenct collection of saved times.The observational data is based on a series of overlapping UBV CCD (charge coupled device)frames covering the centre of M71 (denoted the "core field", Richer & Fahlman 1988; Richer &Fahlman 1989), and a single, deeper field situated 3' from the cluster centre (denoted the "3' field",Richer et al. 1990). Each of the two data sets provided a mass function, and the core field alsoprovided the information required to produce surface density profiles for the brighter stars. In DFRthe star counts were reanalysed on a finer spatial grid than was used in Richer & Fahlman (1989).The observations are shown in Fig. 4. There are two features to note in these data. First, thesurface density profiles are basically King-model like, but the one for the brighter stars does showa slight central enhancement over that expected from a King model. Second, a significant degree ofmass segregation is seen between the two fields considered.All of the models described in this chapter were run before the effects of stellar evolution wereadded to the Fokker-Planck code. The initial models are as described in §2f. In most of the modelsdiscussed in this chapter the main sequence was represented by ten mass classes covering starswith masses between 0.16 and 0.89 M o . The low mass cutoff was somewhat arbitrary in that itwas the low mass limit to the VandenBerg isochrone used, but the initial model MFs do extend intothe range of lower limits for observed globular cluster mass functions (Richer et al. 1991), andpast that observed in the 3' field. The model MF could be extended to even lower masses, but masssegregation makes them unimportant to the dynamics of the core of the model. One model, to be39314.5—30.8 —0.6 —0.4 —0.2 0log m (M0)—0.5^0^0.5log r (pc)I iiii i mi l i,11111.-^ (a) -_^ixI1^1^1^1^1^1^1^1^1^1^l^1^1^1^1^1^1^1^1^1Fig. 4.—(a) The surface density profiles for the two brightest magnitude ranges observed inM71. These correspond to the mass ranges 0.85 <m/M o <0.89 (filled circles) and 0.71< m/Mo<0.85(open circles). The crosses show the total density over the entire luminosity range. (b) The massfunctions observed for M71 from Richer & Fahlman (1989) (filled circles) and from Richer et al.(1990) (open circles). The former covers the inner 3.4' and the latter a rectangular field with amean distance of 3' from the cluster centre.40mentioned below, which was calculated with an additional low mass bin that extended the lowmass cutoff to 0.1 M o , supports this conclusion. The mass bins for the higher mass stars werechosen to be the same as the observational bins for the core ME In addition, three bins were usedto represent a distribution of white dwarf stars with masses from 0.9 to 1.14 M o . The three whitedwarf bins were of equal width and an equal number of stars was put in each, making up apredetermined number fraction of the stars in the cluster. The main sequence stars were distributedas a sum of two overlapping power-laws. The deep MFs compiled by Richer et al. (1991),including those for M71, show that the observed MFs cannot be described by a single power law.In view of these observations, a combination of two power laws, while still an approximation,may be closer to the truth than a single power law. Undoubtedly, the true shape of the IMF is yetmore complicated. The mass spectral indices of the two power laws, their relative proportions andthe number fraction of white dwarfs determined the IMF used for each model.The philosophy followed for all the comparisons between models and observations in thisthesis is to transfer the results of the Fokker-Planck models onto the observational plane. Themodelled data were then compared directly with the observations. The selection of the time when amodel best matched the observations was done subjectively. As will be seen, the models did notmatch the M71 observations particularly well, so more objective matching procedures wereinappropriate.Most of the comparisons involve surface densities, either directly, in the surface densityprofile, or via a mass function. Periodically when it is run, the Fokker-Planck code saves a fulldescription of the current state of the model. For each of these save times the density profiles foreach species were extracted and the projected densities calculated. In order to compute a modelSDP, the projected surface density for the appropriate species is integrated across the same regionsas were the observations. For example, if the observed SDP was calculated for 10 logarithmicallyspaced annuli with width 0.2 dex and the innermost boundary lying at a projected distance of 1 pc,then the projected model surface-densities would be integrated over the same set of annuli. Formass functions, where the observations were made in a region of given dimensions and positionwith respect to the cluster centre, the projected model surface -density for each mass species would41be integrated over the same region. This procedure ensures that the observations and the modelledobservations are subject to the same degree of averaging. In the case of velocity dispersions, theobervations were just compared with the line-of-sight velocity dispersion profile. A more detailedcomparison could be done if the radial distribution of the observed stars were known, but thisinformation was not available for the observations used here.In the case of M71, the model SDP was calculated using the same set of annuli used in DFR.The 0.869 Mo mass bin (see Table A2) was used for the 0.85<m/Mo <0.89 SDP, and the sum ofthe 0.816 Mo and 0.743 M o mass bins was used for the other SDP. The model core MF wascalculated by integrating the projected density profile for each species over a circular area of radius3.66 pc (This assumes a distance of 3.68 kpc [Richer & Fahlman 1988]), and dividing by widthsof the mass bins to get numbers per unit mass. The model 3' MF was calculated the same way,except that the integration region was a rectangular area with the same dimensions and positionwith respect to the cluster centre as that observed by Richer et al. (1990). The areas covered by thetwo MFs overlap to a certain extent but the core MF is dominated by the higher density of stars inthe core of the cluster. The stellar density in the 3' field is in good agreement with the SDPs fromthe core field.In order to discuss the degree of mass segregation more quantitatively I define a segregationmeasure, S, as the logarithm of the ratio of two mass functions. In this way, the effects of masssegregation can be examined independently of the mass function itself. Segregation measures couldbe defined for changes in the mass function with time (eg. with respect to an IMF) or with space(eg. with respect to a global mass function or a central mass function). In order to remove thedifference in absolute numbers due to changes in density the segregation measure may also benormalized in some manner. Being the ratio of two quantities, the errors in observed segregationmeasures tend to be large, but it still serves as a useful tool. In this thesis I will discuss radialsegregation measures of the form(m, a; b) = log{Nm (m)[field a]/Nm(m)[field , (3.1)where Nm(m)[field a] is the number of stars per unit mass of stars with mass m in the fielddesignated 'a'. A larger value of S r indicates a relative surplus of stars of that mass in field 'a' overfield 'b'.42For M71 I consider the segregation measure S r(m, 3'; core) which is shown in Fig. 5. That thesegregation measure is negative for all masses reflects the lower stellar density in the more distant3' field with respect to the core field. Sr(m, 3'; core) becomes less negative with decreasing mass.This indicates a relative surplus of low mass stars over higher mass stars in the more distant field.This is the type of behaviour expected as a result of mass segregation where the more massive starsdominate the core, leaving relatively more low mass stars further out. On the other hand, thedeficiency of stars with m=0.55 M o in the 3' field is difficult to understand. In is hard to see whysuch a short mass range should go against the general trend determined by mass segregation. Theexistence of this dip in the data and the large errors associated with S,. should serve as a warningagainst over-interpreting the observed segregation measure.To reduce the differences in the MFs to one number we can consider a segregation measurenormalized by mass over the mass range observed. This is defined asS,.,.(mp a;m2 ,b) S,.(m„a;b)— Sr (m2 ,a;b).^(3.2)For the M71 observations the normalized segregation measure isSrM: Sr,m (0.35 Mo , 3'; 0.85 Mo , core) = 0.4.^(3.3)If the mass function is a power-law form thenxa — xb =^(ml , a; m2 , b)/log(mi /m2 ) .^(3.4)For M71, x3 , — xco,. — 1In order to try and find a match to the observations of M71 I ran a large number of models withvarying MFs. The MFs differed in the MSIs and relative weights of the two power laws, in thenumber fraction of white dwarfs, and, for some special cases, in other ways as well. The MSIsranged from -1 to 5 in various combinations. It was found that in order to match the steepness ofthe mass function for the 3' field a very steep IMF, with MSI around 5, was required for the leastmassive stars. A more moderate MSI, around 1, was helpful in matching the brighter stars.I also ran some models that had one bin of 1.4 M o to represent the degenerates as opposed tothe three, lower-mass, white dwarf bins discussed above. This change did not greatly affect theresults. I also ran a model, identical to the standard model to be discussed in some detail below,43—1.4—1.2-0.8I^I^I^I^I^I 0.4 0.5 0.6 0.7^0.8^0.9m (M0)Fig. 5.—The segregation measure defined in eq (3.1) for the 3' field with respect to the corefield.44but which had an extra low-mass bin extending the mass function down to 0.1 M o . This bincontained 58% of the cluster's mass and 80% of its stars; the evolution of the core was not verydifferent from that in the standard run. The main difference was that the model evolved somewhatfaster, reaching core collapse in only 80% of the time required for the model without the extralow-mass stars.The best-fitting (in an eyeball sense) of the models run was the one with the MF and initialconditions listed in Table A2. I shall refer to this model as the "standard run". In Figures 6 and 7 Ishow a comparison of the standard run with the observations. The models are shown at four times.The dot-dash line shows the initial model, the thick line the model at the reversal of core collapse,and during the post collapse phase the thin solid and dashed lines are for when the model has 9%and 2%, respectively, of its initial mass. The standard model best fits the outer SDPs (beyond 1pc) when the mass is 0.09 Af t (thin solid line). The main problem in comparing this model withM71 is obvious. The model SDPs at small radii do not flatten off as do the observed SDPs butcontinue in a power-law fashion to the inner limit of the observations. The slope of the power lawsection of the SDP does not flatten significantly during the entire PCC evolution. This behaviour iscommon to all the models with initial conditions discussed thus far. Adjusting the scale length doesnot noticeably improve the fit to the outer points of the SDPs and certainly does not do anything forthe core.While the model does not reproduce the detailed structure of the MFs the segregationmeasure is similar to that observed.These Fokker-Planck models show the following general behaviour. First, the models reachedcore collapse in of order ten or less of their current half-mass relaxation times. Second, only verymild amounts of mass segregation are seen before core collapse, certainly an insufficient amount incomparison to that seen in M71. Third, the PCC models all have very steep central SDPs. Thissteep profile has been been considered to be the identifying feature indicating that core-collapse hastaken, or is taking, place (Djorgovski & King 1986).In view of these models we may expect to see two classes of globular clusters. In the first classwould be those relatively unevolved clusters with long relaxation times, little mass segregation andKing-model-like surface brightness profiles. I will refer to clusters with these signs as "unevolved"450—0.6 —0.4 —0.2log m (M0 )0.85<m/Mo<0.89—0.5^0^0.5^1log r (pc)No 3Q) 2b0 01—0.5^0^0.5log r (pc)log m (M0 )4Fig. 6.—The "standard" Fokker-Planck model of M71 compared with the observations of Fig.4. The four panels depict the two SDPs and the two MFs. The four lines shown are the initial Kingmodel (dot-dash line), at maximum core collapse (thick line), at the best fitting time (when themass is 9% of its initial value) (solid line), and when the mass is 2% of its initial value (dashedline). The detached line segments in the MF panels represent the white dwarfs.46-0.8—1.2—1.40.2^0.3^0.4^0.5^0.6 0.7 0.8 0.9 1m (M0 )Fig. 7.—A comparison of the segregation measures for the "standard" Fokker-Planck modelwith the observed segregation measure of M71. The lines have the same meanings as in Fig. 6.47although some dynamical evolution will have taken place. On the other hand, a highly evolvedcluster would be expected to show strong signs of mass segregation, a short relaxation time and asteep power-law cusp in the central distribution of stars. Such clusters I will refer to as "core-collapse"clusters. They may be either currently undergoing core collapse or be in the PCC re-expansionphase.Into which class of cluster should M71 be assigned? As discussed in §4 of DFR, the presenthalf-mass relaxation time of M71 is of order 10 8 years. Assuming an age of order 15 Gyr, thecluster has been around for about 100 of its present half-mass relaxation times. As well, we seeevidence for a large amount of mass segregation. These two signs would suggest that M71 is acore-collapse cluster and seem to preclude M71 from being bundled in with the unevolved clusters.On the other hand, the lack of a central cusp to the stellar distribution keeps it out of the core-collapsegroup. Identification with either class presents clear problems and all that can be concluded is thatthe observed properties of M71 are inconsistent with the models as they stand. If we are to find amodel that matches M71 then additional physical processes will need to be included. Even then,there is no guarantee that a match will be found. If such is the case, then an approach beyond thatof the orbit-averaged Fokker-Planck equation will be needed.In DFR we attempted, working within the Fokker-Planck formalism discussed in Chapter 2(but without stellar evolution), to find modified models which might give a hint as to what newphysics is required to resolve the question posed by M71. These attempts are reviewed in the nextsection.Section b. Approaches to a solutionSubsection I. Gravothermal Oscillations. As discussed in the introduction, PCC globularcluster models are susceptible to a dynamic instability known as gravothermal oscillations (GTOs).These oscillations are characterized by extreme fluctuations in the central density. They affect onlythe very central regions of the model, and it is as yet unclear whether they occur in real systems.GCLM found that the best fit to M15 was when the models were in a state of maximal expansion.Since the models spend most of their time in the expanded state, this is not so surprising. Perhaps48we are also observing M71 at the low density extreme of a GTO, and it is for this reason that thecentral cusp in the SDP is absent. However, if the observations are sufficiently removed from thecentre of the cluster, or have insufficient spatial resolution to isolate the core, then the data shouldnot depend significantly on the existence and amplitude of GTOs. The models show that atmaximum expansion only the inner 0.02 pc of the SDP are flattened to any extent (DFR, MCH).At the distance to M71 this is only 1" and will not solve the problem.Subsection II. Extreme degenerates. One of the effects of mass segregation during corecollapse is to reduce the power-law slope in the density profile of species less massive than the onethat dominates the core. In an attempt to see if this effect would be of aid in resolving the problempresented by M71, I ran a series of models which contained an extra class of massive degeneratestars. In particular, one model had 3% of its mass in 2.5 M o , non-luminous objects (Table A3).These might be thought of as black holes except that no relativistic effects or tidal effects wereincluded in the models. Fig. 8 shows this "black hole" model at a late time when only 7% of theinitial mass of 6x10 5 M0 remained, and some 47% of the current cluster mass was in the form ofblack holes. For the more massive stars, the fit to the surface density profile is very good. Theproblem with this model lies in the degree of mass segregation between species, and with radius.The model MFs are almost the same, allowing for the difference in absolute density, while themorphology of the observed MFs are quite different from each other. This case may seem somewhatextreme, but when I reduced the fraction of black holes modestly, the slope of the inner SDPbecame too steep, while the degree of mass segregation was still insufficient to match the observations.Reducing the mass of the heavy component allows an increased degree of mass segregation but themass ratio between the heavy component and the visible stars is then too small to provide sufficientflattening.A further problem with the heavy remnant hypothesis arises when the kinematics of the clusterstars is considered. Pryor (1990) has measured the velocities of 79 members of M71 and gets avelocity dispersion of 2.15±0.17 km/s at a mean radius of 1.5 pc. In the "black hole" model, at thetime which best matched the SDPs, the central velocity dispersion for this mass class is 8.2 km/sdropping to 2 km/s only at 4.5 pc. This observation appears to rule out the presence of a substantial,massive, dark component in the cluster.493 ^1^1^1^1^1^1 —0.8 —0.6 —0.4 —0.2 0log m (M0)5.553.5320—0.5^0^0.5^1log r (pc)Fig. 8.—Fit of the "black hole" model to the observations. (a) Surface density profiles. (b)Mass functions. The detached line segments are the white dwarf bins. The time was selected as thatbest matching the more massive stars' SDP. The symbols for the observations are as in Fig. 4.50Subsection III. Modifying the heating rate. The reheating source in these models is due tothe statistical formation and destruction of three-body binaries. Following Lee (1987b), whoadopted a heating rate based on the work of Hut (1985), I assumed that the amount of energygenerated per binary formed was some 100 times the binding energy of the binary. If a morevigorous energy source is available to reverse core collapse then shallower-sloped SDPs arepossible. I therefore ran a series of models identical with the standard model except that the heatingconstant was set to various multiples of the standard one ranging from 10 to 106 times the energyinput per binary.As the heating rate was increased the depth of core collapse decreased. Each model waschecked against the observations and only in the most extreme case, with a heating rate 10 6 timesthe standard value, was a good match to the SDPs achieved. Figure 9 shows the results of this"extra heating" model. Both the SDPs are matched well—if anything, the model SDPs are nowinsufficiently steep—although the model SDP for the less massive stars has somewhat too manystars at all radii. This can be accounted for with a slightly different IMF. The kinematics of thismodel are also consistent with observations, as is shown in Figure 10.On the other hand, the amount of mass segregation between the MFs is not reproduced verywell. The problem with regard to mass segregation, in both the "extra heating" and "black hole"models, can be seen in Fig. 11 where the observed segregation measure is plotted together withthose for the best fitting times for the two models. The amplitude of the mass segregation in themodels is only one-quarter of that observed.From the results of this model, it appears that the problem in fitting M71 may lie in the natureof the heating source. As discussed in DFR, the total amount of energy released in the enhancedheating model is only slightly less than the total amount of energy released in the standard run. Thedifference is that in the standard run the heating rate is effectively zero until the cluster corebecomes sufficiently dense for binary formation to occur, while in the enhanced heating modelenergy is released at a more gradual rate, but without requiring such extreme densities. The highconcentration reached in the standard models then remains imprinted on the density distribution forall later times. Models using tidal binaries as an energy source do not need to reach such high51II I III I III I III I II —0.8 —0.6 —0.4 —0.2 0log m (M0)1—0.5 0^0.5log r (pc)32011 1 111 1 11 11 11^1^1 1(b)34.5 /-Fig. 9.—Comparison with the observations of the "extra heating" model with a heating constant106 times that of the standard model. (a) Surface density profiles. (b) Mass functions. Thedetached line segments are the white dwarf bins. At the time shown the model has shrunk to 3% ofits initial mass. The symbols for the observations are as in Fig. 4.5210^0.5log r (pc)—0.52.502U)00Q)1.5Fig. 10.—Line of sight velocity dispersion from Pryor (1990) for the giants of M71 comparedwith the predictions from the "extra heating" model for the same time as Fig. 953-0.8—1.2—1.40.2^0.3^0.4^0.5^0.6 0.7 0.8 0.9 1m (M0)Fig. 11.—Comparison of the segregation measures for the best fitting times of the "black hole"(short dash) and "extra heating" (long dash) models with the observed segregtion measure. Neithermodel shows the observed amplitude of mass segregation.54concentrations before core collapse is reversed (SOC, Lee 1987a). Even so, for the model presentedby SOC, a core radius greater than 0.1 pc is not reached until hundred of billions of years havepassed. These models are also too concentrated to represent clusters like M71.The conclusion that can be drawn from this is that a continuous source of energy in a globularcluster will prevent a deep core-collapse from taking place while still allowing for dynamicalevolution. Whether the amount of dynamic evolution would be sufficient to produce the observedmass segregation is a question which requires an answer.In DFR I discussed two possible sources of this extra energy, both of which must occur tosome extent in globular clusters. The first is due to interactions of the cluster stars with primordialbinaries in the cluster. Gao et al. (1990) have run Fokker-Planck simulations which explicitlyinclude a population of initial binaries along with one mass class of single stars. In these modelsthe time to initial core collapse increases as the fraction of binaries increases. When 20% of theinitial model is in binary stars, core collapse is delayed to 100 initial half-mass relaxation times.The effects of a population of primordial binaries on mass segregation cannot be estimated from theGao et al. models. Much would depend on the mass distribution of the initial binary population.Goodman & Hut (1989), on theoretical grounds, suggest that, depending on the initial abundanceof binaries and their rate of destruction, the presence of primordial binaries would result in arelatively large core and hence a SDP which does not show a cusp. On the other hand, such a coreis likely to consist mostly of binaries. Core collapse in the models of Gao et al. terminates whenthe mass segregation between the binaries and singles is complete, again with a core dominated bybinaries. The high binary fraction expected in the core is thus a testable prediction of this scenario.A second heating source is that due to mass loss during stellar evolution. It was this approach that Ihave chosen to pursue here, and a discussion of the models which include stellar evolution followsin Chapter 5. To further understand the problem posed by M71 it would be helpful to have datafrom another cluster for comparison. In the next chapter I present observations of NGC 6397 withthis objective.55Chapter 4. NGC 6397One question raised by the contradiction between the models and observations discussed in thecontext of M71 is the prevalence of this discrepancy. In the interests of getting a slightly widerview, I took advantage of an extensive set of observations of the cluster NGC 6397 acquired byG. G. Fahlman and I. B. Thompson in 1989.NGC 6397 is in many ways a useful foil to M71. Both lie at about the same galactocentricdistance (7.4 kpc for M71, 6.9 kpc for NGC 6397 [Webbink 1985]) but M71 has the kinematicsof a disk cluster while NGC 6397 belongs to the galactic halo (Cudworth 1992). NGC 6397 ismore massive and more concentrated than M71 and is classified as a PCC cluster by Djorgovski &King (1986). We would expect that it would be easier to find matching models in the case of NGC6397 than in that of M71, in that the former has a surface density profile with a power-law core.Based on the models discussed with respect to M71, this is an indication that NGC 6397 isdynamically evolved. From the multi-mass King models of Meylan & Mayor (1991) a half-massrelaxation time may be estimated using eq. (2.11). With their mean model values for the total mass,1.0±0.1x105 Mo , and for the half-mass radius, 5.2±0.3 pc, and taking the mean stellar mass to be0.4 Mo , trk = 7x10 8 yr. This lies in the middle of the distribution times of Fig 8-2 of Binney &Tremaine (1987). Djorgovski et al. (1991) examined colour gradients in several highly concentratedglobular clusters including NGC 6397. They found that the central cusps become bluer towards thecluster centres and attribute this to the destruction of red giants and subgiants, presumably bydynamical processes.NGC 6397 (a 1950=17h 36m 38 8 , 81950= –53° 38.9') lies at 1=338.165°, b=-11.959°. Reddeningmeasurements give consistent values—Cannon (1974) found E(B–V)= 0.18±0.01, van den Bergh(1988) E(B–V)=0.19±0.02, and VandenBerg, Bolte & Stetson (1990) E(B–V)=0.19. VandenBerget al. also find evidence for a small amount of variable reddening. From the data in their paper thisappears to be at the 0.013 magnitude level, and should not affect the results below where aconstant reddening of E(B–V)=0.19 is assumed.56There are two sources for the distance modulus to NGC 6397. The first is a venerablemeasurement of the magnitude of the horizontal branch by Cannon (1974). Depending on thechoice of horizontal branch calibration, a distance modulus of (m —M)v = 12.3 ± 0.3 is obtained(q.v. Harris 1980, Zinn 1985). More recently, Anthony-Twarog, Twarog, & Suntzeff (1992)derived an independent distance using StrOmgren photometry. Using main sequence fitting and theMv, (b —y) relationship method of Laird, Carney, & Latham (1988) they got (m —M)v = 12.1±0.3.With this distance modulus, the colour-magnitude diagram does not give a good fit to the theoreticalisochrones of VandenBerg & Bell (1985) as modified for enhanced oxygen abundance by McClureet al. (1987). A better fit is achieved if (m —M)v = 12.4 is used. Anthony-Twarog et al. attributesome of the discrepancy to the bolometric corrections used to transfer the luminosities onto theobservable plane. In any event, their distance moduli are consistent with (m —M)v = 12.3 adoptedhere. These values for the distance modulus and reddening gives a heliocentric distance of 2.2 kpcto NGC 6397. At this distance 1 pc subtends 1.6'±0.2, given the uncertainty in the distancemodulus. The relatively short distance to the cluster allows for very faint stars to be observed withonly moderate investments of telescope time. Fahlman et al. (1989, hereafter FRST) have previouslydone deep starcounts in NGC 6397 and observed the luminosity function to MI = 11.5. Thiscorresponds to a mass of 0.12 M o .Since the present observations are in the /—band, an apparent distance modulus in that colour isrequired. The formula in Dean, Warren, & Cousins (1978), extrapolated from 0 and B stars,implies E(V—I)IE(B—V)=1.35 for the (B—V)-1.3 typical for the lower main sequence of NGC 6397(Alcaino et al. 1987). Consideration of the entry for K3 III stars in Table 3 of Taylor (1986) givesa value of 1.38 for this ratio. On the other hand, Grieve's (1983) calibration of the colour excessratio for F and G supergiants gives a higher value, more like 1 6. This difference will only affect(m—M), by 0.05 mag, well within the uncertainties in the distance modulus itself. Given thesevalues for the ratio of the reddenings, I adopt an apparent I distance modulus (m—M)=12.0. Thisis the same as the distance modulus used by FRST.The observations to be discussed in §4a have been used to produce surface density profiles and57mass functions. One SDP is for stars with R14 at radii from the centre out to 12'. A second SDP,reaching to 1=15.5, approximately the main sequence turn-off, only includes stars more than 20"from the cluster centre, since, for the area interior to this, the counts of the fainter stars are veryincomplete because of crowding. The SDPs will be discussed in §4b. The MFs, to be discussed in§4c, are for fields at radii bracketing that of FRST. The broad range of radii and masses observedgive useful leverage on the question of mass segregation in this cluster.Section a. Observations & reductionThe observations of NGC 6397 were obtained during an observing run at the Las CampanasObservatory during May and June 1989. There were two sets of observations, both in the Cousins/---band and using the TI #1 800x800 pixel CCD. The first data set consists of deep observations ofseveral fields using the 2.5m du Pont Telescope. Observations of four fields will be discussedhere. The second set of observations were obtained with the 1 m Swope Telescope and covers bothfour overlapping fields extending from the cluster centre, and two regions further out whichoverlap 2.5 m fields. The positions and designations of the fields are listed in Table 3. Figure 12shows the locations of the program fields with respect to the cluster centre. The fields are labelledwith their names from Table 3 except that "du Pont" has been shortened to "duP". All stars in theprogram frames observed with R13 have been marked by points. The relative sizes of the pointsgives an indication of their magnitudes and most of them can also be seen in Fig. 4 of Cannon(1974). In addition to the program fields, a field about 1° away from the cluster, and at the samegalactic latitude, was observed with both telescopes. The star counts from this "background field"were used to measure the contribution of non-cluster objects to the counts. Observations of fieldswith standard stars were also take on nights with photometric conditions. The 2.5 m observationswere used to construct luminosity and mass functions, while the purpose of the 1 m data was toobserve the surface density profile. On the 2.5 m, the plate scale was 0.162"/pixel. As the seeingwas usually around 1" FWHM, this resulted in heavily over-sampled star images. In order toimprove the signal-to-noise ratio, the frames were boxcar averaged 2x2 following debiasing and58flat fielding, giving an effective plate scale of 0.324"/pixel. For the 1 m frames the plate scale was0.435"/pixel and the seeing was typically around 1.7" FWHM.Table 3: Field locations for NGC 6397 observationsField Name Xa^Ya^dist.a^width^length(all positions and dimensions in arc minutes)ExposurebtimesSwope:c -0.12 0.03 0.12 5.775 5.745 1x10s, lx1OsSwope:sa 1.29 -4.31 4.50 5.771 5.771 2x5sSwope:sb 3.64 -6.66 7.60 5.786 5.757 2x5sSwope:na -0.25 3.32 3.33 5.778 5.778 2x30sSwope:w 11.09 -0.06 11.1 5.778 5.778 2x30sSwope:nb -1.98 11.73 11.9 5.778 5.778 2x30sSwope:bf 60 5.778 5.778 2x30sdu Pont:if 0.67 3.79 3.86 2.138 2.111 14x200s, 13x200sdu Pont:n -2.16 10.48 10.7 2.155 2.130 8x450s, 8x450sdu Pont:w 11.15 -0.25 11.2 2.144 2.122 4x600s, 5x600sdu Pont:bf 60 2.140 2.140 7x600s, 6x600sFRST 3.24 -5.59 6.46 2.160 2.160 12x450saPosition of field centre with respect to cluster centre.bNumber and length of exposures going into each reduced image are listed.Debiasing was done using routines at the Carnegie Institution. Flat fielding; image registrationand averaging, where required; and image trimming were done using procedures within IRAF. Theflat fields for most of the 2.5 m data were exposures of clouds from the first night of the run. Thisflat field did not work very well for the du Pont:bf field so a combination of this flat and dome flatswas used. Since the 2.5 m data were acquired under non-photometric conditions, the exposureswere weighted by their mean flux levels during the averaging. The sigma clipping algorithm of theIRAF task imcombine was used in order to remove cosmic rays. The deepest set of exposures ineach field were combined into two independent images for that field. While somewhat fainter starscould have been recovered if all the data had been combined into a single frame, the advantage ofworking with two images is that one can be more certain of the reality of the recovered stars(Stetson 1991). In observing a mass function the additional confidence conferred by this procedurefar outweighs the loss in limiting magnitude.59I^-i ^I— duP:nduP:w1510Swope:nbESwope:naSwope•c•_______e_duP:if.^.. *...•^..• .. .• II, ..••,,- V.' • ,.... i ..!ti, ; • ••^••Swope•w••• ,FRSTSwope:saSwoper^IIIIII^ 1---1 ^I^1^1^L —5^0 5 10 ^15X wrt centre (arc minutes)— 10 10— Fig. 12.—Positions of the program fields observed in this study. The coordinates are withrespect to the adopted cluster centre marked by a `+'. The dots are all stars with R13 observed inthe program fields. Most of these also appear in the plate of Cannon (1974) and can be used tolocate these fields. The abbreviation "duP" indicates the du Pont fields which were selected toavoid any bright stars.They are rotated slightly with respect to the Swope fields. "FRST" is thefield observed by Fahlman et al. (1989).60The data in fields Swope:nb, Swope:w and Swope:bf were taken under photometric conditionsand calibrated using standard stars in fields SA 110, SA 112, SA 113 and in the field of PG 1323-085.Good results were obtained using the growth curve program DAOG ROW (Stetson 1990) andseveral other photometry calibration programs kindly provided by P. B. Stetson. The transformationto the standard system is good to within about ±0.03 mag. The calibration for the corresponding2.5 m frames was transferred from the 1 m observations. While the remaining four 1 m fields wereapparently observed under photometric conditions, it did not prove possible to use the standardsobserved that night. The fluxes measured for these stars were smaller than expected by a factor ofabout three for unknown reasons. Fortunately, the Swope:sa and Swope:sb fields overlapped thatof FRST and a calibration, good to within the ±0.05 quoted by FRST was obtained from magnitudesof stars in that field (Fahlman 1991). This calibration was also transferred onto field du Pont:if.Auriere, Ortolani & Lauzeral (1990) found the centre of NGC 6397 to be within 1" of theirstar 1. This star is identified with star 335 of Woolley et al. (1961) which was taken as the centreby Auriere (1982). Having identified this star in field Swope:c, I took as the centre a position justoff of this star. The relative positions of all the other frames within the inner region (ie. Swope:na,Swope:sa, Swope:sb, du Pont:if, and FRST) were found with respect to Swope:c. The 2.5 mimages were found to be rotated slightly with respect to the 1 m frames. The positions of theSwope:nb and Swope:w fields with respect to the central four 1 m fields were found usingphotographic images of NGC 6397. The positions of several stars in the two pairs of regions weremeasured on the photograph with a measuring engine and identified on the program frames. Therequired offset between the origins of the Swope:sb and Swope:w fields, and between the originsof the Swope:na and Swope:nb fields, were solved for in a least squares fashion together with aconstant scale factor. No rotation was apparent in the transformation relationship and in the finalcalculation was assumed to be zero. The resulting positions for these fields with respect to thecluster accord well with a sketch of the field positions made at the telescope. The positions of thedu Pont:n and du Pont:w fields, which are contained within the Swope:nb and Swope:w fieldsrespectively, were calculated with respect to the positions of these 1 m fields.61The data was reduced using versions of DAOPHOT and ALLSTAR (Stetson 1987) in the usualmanner. Two finding passes were made on the images with a find threshold of 46. Generally, aMoffat function (Moffat 1969) with exponent 1.5 was used for the analytic portion of the PSF,plus a constant residual look up table. Occasionally, a linear variation in the look up table was usedbut it did not materially affect the photometry. Each of the two independent images for each fieldwas reduced separately. Stars on the two lists were matched and the mean offsets in magnitude andposition were calculated. These offsets were then applied and a final star list for each field wasproduced containing those pairs of stars with centroids closer than 1 pixel.In order to correct for incompleteness in the star counts, tests were done using the usualprocedure of adding artificial star images to the program frames and then re-reducing them. Theartificial stars were added with a magnitude distribution based on that estimated from the initialreduction of the frame. Ten percent of the number of stars observed in the magnitude range ofinterest were added at random to each artificial star frame. The same artificial stars, with appropriatemagnitude and position offsets, were added to both images in a pair. Ten to thirty such pairs offrames were reduced for each field.The final star lists were prepared in the following manner. For each field, regions surroundingsaturated stars, diffraction spikes, charge overflow columns and other defects, were identified.These areas have a significantly higher rate of spurious detections. All artificial stars added, andobjects detected, in these regions were ignored during further analysis. Also ignored were starsfound in the region lying outside that where artificial stars were added. The effects of theserestrictions are admittedly small. The detection lists for the two images in a pair were then shiftedonto common position and magnitude zero points. The stars on the lists were matched and allobjects found on both frames with position centroids within 1 pixel were deemed valid detections.The list of valid detections was then compared with the list of positions and magnitudes of theartificial stars added to that pair of frames. All valid detections with positions within 1 pixel andmagnitude differences less than 0.7 mag with respect to the artificial star list were counted asrecoveries of artificial stars. The remaining, presumably real, stars were put into a separate list.62Star counts were performed on the list of real stars from each artificial star test separately and theresults averaged. Small variations in the number of stars were seen from test to test due tovariations in the local crowding conditions. Errors in the star counts are a combination of thePoisson error and the standard deviation in the number of stars recovered.The observed luminosity function is a convolution of the true luminosity function with theprobability that a star of a given magnitude will be found, and if so, with what observed magnitude.The artificial star test give an approximation to the probability function and allow for the inversionof the observed luminosity function to give the true one. The luminosity functions were correctedfor incompleteness using two techniques. In the simplest approach, all stars added in a givenmagnitude bin and subsequently recovered were counted as recoveries of stars with that magnitude.The ratio of the number of recovered stars to the number added in that bin gives the recovery rate,and the incompleteness correction factor is the inverse of this. This method tends to underestimatethe uncertainties in the corrected counts and does not take into account the effects of stars with truemagnitudes lying in one bin being recovered in a neighbouring bin. This bin jumping can becorrected for in an approximate fashion using techniques such as those described by Drukier et al.(1987) and Stetson & Harris (1988).The second method I used is a simplified version of the matrix inversion method of Drukier etal. (1987). As an intermediate method between using the full matrix, as in Drukier et al., andcompressing all the elements onto the diagonal, as in the approach discussed above, a tri-diagonalrecovery matrix was constructed. The advantage of a tri-diagonal matrix is that it is straightforwardto invert and calculate the errors for. The inverted matrix, together with some assumptions aboutthe recovery rates and magnitudes for stars in the faintest bins, was used to calculate the completenesscorrected counts. The results of this method are for the most part consistent within the uncertaintieswith those of the simpler approach discussed above, except that the error estimates and magnitudelimit to the star counts are more conservative using the matrix approach.63Section b. Surface Density ProfileThe surface density profiles were produced from the 1 m data. Using the artificial star tests forthe Swope:c field, the completeness was calculated as a function of magnitude for annuli about thecentre . In the central region, within a radius of 29 pixels (0.21% the recovery rate was greaterthan 90% to This magnitude was adopted as the cutoff for the SDP. A similar result wasachieved for a pair of images of the same field with exposure times of 5s. On the longer exposurestaken for fields Swope:na, Swope:sa and Swope:sb the brighter stars are saturated and could notbe counted in the luminosity function. Therefore, the single short exposure images (averages of thetwo 5s exposures obtained for each field) were used to produce the surface density profiles. In theSwope:nb, Swope:w and Swope:bf fields the shallowest images, averages of two 30s exposures,were used. These fields contained no stars bright enough to have poor photometry over this lengthexposure. For the Swope:c data the detected stars were counted from the output of all 10 artificialstar runs and an average taken. The variation between star lists was quite small, confirming theresults of the artificial star tests.The stars were counted in two sets of overlapping annuli to reduce the effects of binning. Theouter limits of the annuli were spaced logarithmically with each annulus having a width of 0.2 dex.The radii of the inner disks of the two sets of annuli were 0.21' and 0.17'. The positions of theseannuli with respect to the fields are shown in Figure 13. The dashed and solid circles are theboundaries of the two sets of annuli. The areas of the intersections of the annuli with the rectangularfields (I will refer to these as "sections") were calculated and the star counts were converted tosurface densities. The adopted effective radius for each section is the radius which divides the areaof that section in half. This radius will vary between sections on the same annulus, but fromdifferent fields, depending on the geometry of the section. The field star density was measured onSwope:bf. Twenty-three stars brighter than 1=14 were found on that frame giving a surface densityof 0.69±0.14 stars/square arc minute. The resulting surface density profile, with the backgroundsubtracted, is shown in Figure 14. The different symbols designate which series of annuli aparticular density comes from. The errors are Poisson errors except in the case of field Swope:c64I —5^0^5X wrt centre (arc minutes)Fig. 13.—Map showing the annuli used in constructing the surface density profile. The rectanglesare the Swope fields shown in Fig. 12. The solid circles delimit the set of annuli starting at 0.21'(group B), with a logarithmic width of 0.2 dex and the broken circles the set of annuli starting at0.17' (group A) and having the same width.650background—0.5^0^0.5log r (arc minutes)iFig. 14.—Surface densities for 1<14 for the annuli in Fig. 13. The squares are for the 0.21'(group B) annuli and the x's for the 0.17' (group A) annuli. The densities are shown at the radiuswhich bisects the area of the section. The background density has been subtracted and is shown byhorizontal line. The points with large error are those for small sections at the edge of a field withvery few stars.66where a small contribution due to variations in the number recovered on the artificial star frameshas been included.In order to account for areas in which no stars were seen and eliminate the apparent anomaliesassociated with small number statistics, any section of a field which contained fewer than 10 starswas merged with its neighbouring section. The area and effective radius for the new, combinedsection were calculated, and the density and error re-evaluated. The SDP after this merging isshown in Fig. 15 and listed in Table 4. The various symbols indicate the frame of origin for eachdata point.A set of surface density profiles were produced in the same way for a set of centres offset, byvarying amounts, from that adopted. Little difference was seen for centres offset by up to ±20pixels (8.7") in either or both directions. The outer SDP was little affected by all these variations incentre position and I conclude that the errors in the determination of the relative positions of theouter frames and the inner ones will have a negligible effect on the surface densities.In the days when star counting was done on photographic plates, a reseau consisting ofsectored annuli was placed on the plate and the stars were counted by sector. The density anduncertainty for a given annulus was calculated from the counts for all the sectors in the annulus.Here, with the exception of the central region, only a limited area of each annulus was surveyed. Inorder to check the consistency of the star counts and the assigned errors I split the lists ofrecovered stars for Swope:c into quadrants about the cluster centre and repeated the countingprocedure for each quadrant separately. Figure 16 shows the results of this test. The symbolswithout error bars are the densities for the four quarters of field Swope:c. The line connects themean of the four measurements in each annulus and the error bars are the standard deviations of thesamples. The solid circles offset +0.025 dex in radius are the overall densities from the full fieldtogether with the uncertainties computed from Poisson statistics and frame-to-frame variations inthe counts. With the exception of the innermost point, the quarter-to-quarter errors are within afactor of two of the Poisson errors. There is excellent agreement in the errors for mean pointsbetween 0.25' and but a larger quadrant-to-quadrant variation is seen for the annuli beyond 1'.67The large scatter in the innermost point is an artifact of the linear arrangement of the bright stars inthe core of NGC 6397 (q.v. Auriere et al. 1990).Table 4: Surface Density Profile for/ <14reff N (TN Field^Annulusa(arc min.) (arc min.-) Group0.149 99. 28. Swope:c^B0.187 89. 21. Swope:c^A0.279 53. 17. Swope:c^B0.350 40. 11. Swope:c^A0.442 35.3 8.3 Swope:c^B0.555 32.7 6.3 Swope:c^A0.700 23.7 7.7 Swope:na^B0.700 27.9 4.7 Swope:c^B0.847 26.3 5.6 Swope:na^A0.880 19.0 3.1 Swope:c^A1.109 15.6 2.2 Swope:c^B1.119 23.8 4.5 Swope:na^B1.394 11.4 1.5 Swope:c^A1.403 12.7 2.5 Swope:na^A1.758 6.8 1.0 Swope:c^B1.766 6.2 1.4 Swope:na^B1.866 6.3 2.2 Swope:sa^B2.198 6.0 1.4 Swope:sa^A2.210 4.76 0.66 Swope:c^A2.218 4.70 0.98 Swope:na^A2.737 2.97 0.46 Swope:c^B2.749 3.45 0.70 Swope:na^B2.782 2.77 0.81 Swope:sa^B2.996 2.13 0.53 Swope:c^A3.373 2.26 0.54 Swope:na^A3.498 2.34 0.61 Swope:sa^A4.275 1.55 0.45 Swope:na^B4.345 2.49 0.52 Swope:sa^B5.297 1.41 0.42 Swope:na^A5.737 0.78 0.43 Swope:sb^A5.786 0.92 0.31 Swope:sa^A5.857 1.45 0.57 Swope:na^B6.379 0.40 0.32 Swope:sa^B6.867 0.83 0.31 Swope:sb^B8.475 0.46 0.26 Swope:sb^A9.350 0.09 0.28 Swope:sb^B9.442 0.18 0.30 Swope:w^A11.212 0.03 0.21 Swope:w^B12.020 0.15 0.21 Swope:nb A & B12.296 -0.06 0.23 Swope:w^AaGroup "A" is the set of annuli starting at 0.17' with width 0.2dex. Group "B" annuli start at 0.21' and have the same logarithmicwidth. Section merging as described in the text has been applied.68background• Swope:c• Swope:na■Swope:saOSwope:sbA Swope:nb* Swope:wI^IiIii^i^I I^I —1 —0.5 0 0.5log r (arc minutes)Fig. 15.—The surface densities shown here are for the same sections as in Fig. 14 except thatall sections which contained fewer than ten stars (including those with no stars) have been mergedwith their neighbours on the same field. As the relatively empty sections are always at the edge of afield there is no ambiguity in the direction of merging. The effective radii are those that bisect theareas of the (possibly composite) sections. The various symbols indicate the field of origin for eachsurface density measurement and the magnitude limit is 1=14.69AI —0.5^0log r (arc minutes)2.5'-"a;5z0 1 5cra)bCDO0.5Fig. 16.—A comparison of the central surface density profile for R14 derived from the fourquadrants of field Swope:c with that derived from the entire field. The various symbols withouterror bars are the surface densities derived from the star counts in the four quadrants. The lineconnects the mean densities at each radius and the solid circles are the surface densities measuredfrom the entire field and are offset outward by 0.025 dex for clarity. The error bars on the meancounts are the standard deviation of the sample while the error bars on the counts from the full fieldare the sum of Poisson and frame to frame errors as described in the text.70A second surface density profile, containing all stars brighter than /=15.5 and more than 20"from the centre was produced in the same manner as the 1<14 SDP. The inner radial limit is due toincompleteness. Beyond this point the counts from Swope:c are more than 90% complete between1=15.0 and /=15.5, and greater than 95% complete for all stars brighter that /=15.5. Theseincompletenesses would change the final surface density profile by less than 0.05 dex, well withinthe estimated errors. The magnitude cutoff of 1=15.5 was chosen since this is the magnitude of theturnoff on a roughly calibrated I —(V—/) colour-magnitude diagram produced with data taken at thesame time as that of FRST. A field density of 2.06 ± 0.25 stars/sq. arc minute was derived fromthe Swope:bf data. Since all the stars on both the 1<14 and R15.5 SDPs lie at or above the turnoffthey should have the same mass (or did have until recently) and so should have the same radialdensity distribution. This SDP is listed in Table 5. A comparison between the two SDPs is shownin Fig. 17 where an offset of 0.49 dex has been applied to the 1<14 SDP. The background level forthe 1<14 SDP has been shifted by the same amount. This combined SDP has been used for allsubsequent comparisons.The SDP of NGC 6397 seems to be best characterized by two power laws with a break at 1'.For the 1<14 SDP, the inner section (log r 50.1) has a slope of —0.8±0.1 and the intermediatesection (0..5 log r 50.6) has a slope of —1.7±0.1. The slopes are the mean of those calculated by aweighted linear regression for the two sets of annuli independently. Over the region 0.5 log r 50.6the R15.5 SDP has a slope of —1.60±0.06 and the mean slope in this region for both SDPs is—1.67±0.07. The intercepts are better determined and give the offset of 0.49±0.04 dex used inproducing Fig. 17. The fits are shown in Fig. 18. The dot-dash line is a fit to the entire regionwith r 1' and has slope —1.32±0.03 (mean of least-squares fits to all four SDPs).Given the low galactic latitude of NGC 6397 determination of a limiting radius for this clusteris a difficult proposition. For example, Da Costa (1979), from a single-mass King model fit, findsa limiting radius of 38' and Meylan & Mayor (1991) using multi-mass King models find a meanlimiting radius of 95'. These values are much larger than the limit of the new SDP (-13'). The SDPsteepens at the limit of the observations, but the densities are quite uncertain.71Table 5: Surface Density Profile for / <15.5 rdr N oN Field^Annulus a(arc min.) (arc min.-2) Group0.442 119. 15. Swope:c^B0.555 108. 12. Swope:c^A0.583 128. 28. Swope:na^A0.700 113. 17. Swope:na^B0.700 92.0 8.5 Swope:c^B0.880 70.1 6.0 Swope:c^A0.891 96. 12. Swope:na^A1.109 48.8 3.9 Swope:c^B1.119 70.5 7.7 Swope:na^B1.394 36.1 2.7 Swope:c^A1.403 43.5 4.7 Swope:na^A1.581 34. 11. Swope:sa^A1.758 23.6 1.8 Swope:c^B1.766 25.4 2.8 Swope:na^B1.866 22.9 4.2 Swope:sa^B2.210 15.7 1.2 Swope:c^A2.218 17.1 1.8 Swope:na^A2.236 15.9 2.3 Swope:sa^A2.665 11.6 0.9 Swope:c^B2.749 13.3 1.4 Swope:na^B2.782 9.0 1.4 Swope:sa^B2.996 8.7 1.0 Swope:c^A3.373 8.4 1.0 Swope:na^A3.498 7.6 1.1 Swope:sa^A3.544 6.1 1.9 Swope:c^B4.275 4.16 0.75 Swope:na^B4.345 6.29 0.85 Swope:sa^B4.815 4.9 1.6 Swope:sb^B5.297 3.36 0.67 Swope:na^A5.368 3.99 0.67 Swope:sa^A5.737 3.82 0.85 Swope:sb^A5.857 3.50 0.92 Swope:na^B6.379 2.15 0.61 Swope:sa^B7.072 0.21 0.70 Swope:sa^A7.085 2.36 0.56 Swope:sb^B8.475 1.33 0.45 Swope:sb^A9.350 0.88 0.54 Swope:sb^B9.442 0.79 0.54 Swope:w^A9.860 1.15 0.69 Swope:nb^A10.809 0.56 0.39 Swope:w^B11.234 0.48 0.41 Swope:nb^B12.296 0.35 0.42 Swope:w^A12.678 0.09 0.38 Swope:nb^A13.670 0.28 0.75 Swope:w^B14.058 -0.02 0.53 Swope:nb^BaGroup "A" is the set of annuli starting at 0.17' with width 0.2dex. Group "B" annuli start at 0.21' and have the same logarithmicwidth. Section merging as described in the text has been applied.72Fig. 17.—A comparison of the surface density profiles for all stars brighter than 1=14 (opencircles) and for all stars brighter than1=15.5 (filled circles). The former has been increased by a+0.49 dex on the basis of least-square fits to the region 0.01og r The relative backgroundlevels are: (dashed) 1<15.5, (dotted) k14. The background level for the 1<14 SDP has also beenshifted by +0.49 dex.73—0cda)0bOD01^1^1^1^1^1^I^1^1^1^1^I^1—0.5 0 0.5log r (arc minutes)Fig. 18.—The various lines show weighted least-squares, power-law fits to various radialranges of the combined SDP: dashed line, log r _0.1; solid line; 0.01og r 0.6;dot-dashed line,log r The lines are the means of the fits to each of the two sets of annuli. The means for thetwo radial ranges with the lower limit at log r=0 include points from both SDPs.74Section c. Mass functionsThe luminosity functions for fields du Pont:if, du Pont:n, du Pont:w, and du Pont:bf wereproduced in accordance with the procedures outlined in §4a. The stars were counted in two sets ofoverlapping 0.5 mag bins offset by 0.25 magnitudes in order to reduce the distortions due tobinning. For each field the true observed area, exclusive of the ignored regions, was computed andthe star counts adjusted to be the number of stars which would be observed in a 4.67 square arcminute field (800x800 pixels at a plate scale of 0.1627pixel).At first glance the luminosity functions in Tables 6-9 may appear to present some anomalies.When the incompleteness corrections do not take into account bin jumping, the corrected countscannot be less than the raw counts. The inclusion of the effects of bin jumping can lead to thecorrected number of stars in a bin being smaller than the original number. For a given bin, bin i,the next fainter bin, bin i+1, generally has more stars. The errors in the magnitudes for stars whichbelong in bin i+1 are larger than those in bin i. Even if the scattering is symmetric about the truemagnitude, bin i will receive more stars from bin i+1 than it loses to bins i-1 and i+1. Further,when bin jumping is important, even if all the artificial stars are recovered with some magnitude,the net effect of the corrections on the counts in bin i could very well be to decrease the number ofstars in that bin. This result is even more likely when it is noted that stars have a tendency to befound brighter than they really are. The point is that of the stars observed with magnitudes in bin i,it is possible that more of them belong in bin i+1 than have been lost from bin i. The correction forincompleteness in bin i may not be large enough to overcome this.The background counts, converted to the number expected on a 4.67 sq. arc minute field, aretabulated in Table 6. In the bin centred at 1=22.5, 50% of the artificial stars were recovered withmagnitudes within 0 7 mag of their true values. For the next independent 0.5 mag bin, centred at1=23.0, only 8% of the artificial stars were recovered. Since it is impossible, because of the poorstatistics, to estimate the number of stars in the 1=22.5 bin which originated in the 1=23.0 bin, the1=22.5 bin could not be used in the luminosity function. Similarly, for the overlapping bins centredon the quarter magnitude, the last useful bin is at 1=21.75. Although the next bin has a recoveryrate of 78%, the one following it, centred at 1=22.75, has a recovery rate of only 26%, and of75these, 44% were found in the next higher bin, 12% in the the next fainter bin, and the remaining44% were found in the same magnitude bin they were added in. The recovery rate drops quitequickly, so the last usable bin usually has a rather high recovery fraction. As a result of theseconsiderations, the cutoff magnitude for the du Pont:bf luminosity function is 1=22.25; the last binused being the one centred at I=22.0.Table 6: Raw and completeness corrected luminosity function for thebackground field du Pont:bf.I aNcounte d error aNcomplete error15.00 1.1 1.1 1.1 1.115.25 4.3 2.1 4.3 2.115.50 4.3 2.1 4.3 2.115.75 3.2 1.8 3.2 1.816.00 5.3 2.5 5.3 2.516.25 6.4 2.6 6.4 2.616.50 9.6 3.2 9.6 3.216.75 11.7 3.6 11.7 3.617.00 12.8 3.7 12.8 3.717.25 10.6 3.4 10.6 3.417.50 15.9 4.1 15.9 4.117.75 18.1 4.4 18.1 4.418.00 13.8 4.0 13.8 4.018.25 21.3 4.8 21.3 4.818.50 22.3 4.9 22.3 4.918.75 31.0 5.9 32.6 6.419.00 33.7 6.1 35.0 6.419.25 32.5 6.0 33.3 6.219.50 37.3 6.4 38.1 6.619.75 34.5 6.1 34.5 6.120.00 30.0 5.8 34.9 7.020.25 48.4 7.3 46.7 8.020.50 65.2 8.5 63. 10.20.75 61.8 8.3 69. 11.21.00 73.9 9.1 76. 12.21.25 83.3 9.6 85. 14.21.50 75.8 9.2 75. 14.21.75 78.0 9.3 72. 16.22.00 84.3 9.8 103. 21.'Number per 0.5 magnitudes per 4.67 sq. arc minute field.76The du Pont:if starcounts were found to have a 50% recovery rate at 1=21.7 and when the binjumping effects are taken into account the limiting magnitude is /.21.5. The 50% recovery ratemagnitude for the du Pont:n data is 1=22.8, and for the du Pont:w data it is /=22.1. The magnitudelimits of the incompleteness corrected star counts are 22.5 and 21.75 respectively. The du Pont:nLF goes 0.75 mag deeper than the du Pont:w LF because of better seeing and a longer totalexposure. Since the du Pont:w and du Pont:n are at about the same distance from the cluster centre,and have low signal-to-noise, the LFs for these two fields have been averaged where they overlapand a final, outer LF produced. The bins with centres with 121.75 are from the du Pont:n dataonly. The magnitude limit of the background LF cuts off this combined LF, which will be referredto as the du Pont:out LF, at /=22.25. Tables 7-9 contain the raw and incompleteness correctedluminosity functions for the three program fields.Figure 19 shows a comparison of the luminosity functions for the background fields observedin this work and in FRST. The agreement is good despite the two fields being on opposite sides ofthe cluster. This confirms the observation by Da Costa (1982) that the field stars are evenlydistributed in the vicinity of the cluster. The final, background subtracted, luminosity functions atthe two distances are listed in Table 10 and are plotted in Figs. 20 and 21. These two figures showthe two new LFs and compares them with the background star counts. At the distance of the twoouter fields the number of field stars is equal to the number of cluster stars, so deriving a moreprecise luminosity function will require observing a much larger area.Combined with the inner LF from field du Pont:if and the LF of FRST, the du Pont:out LFgives three deep LFs for NGC 6397. Figure 22 brings together all three of the observed deep LFsfor NGC 6397. The two types of filled symbols are the new LFs from this study, while the opensymbols are the deeper LF of FRST. No further normalizations have been applied; the offsetsbetween the LFs are a reflection of the overall decrease in density with radius. For the magnituderange shown, all the stars should be on the main sequence, since the turn-off lies at about 1=15.5.The FRST LF goes deeper because of the better seeing for that data, and because of the moreconservative approach to incompleteness corrections used here. The origin of the dip around /=1977in the du Pont:out LF is unclear, but is seen, at about the same magnitude, in both fields. Thegeneral depression in the counts between 1-18.5 and 1-20 suggests that whatever is causing the dipis not restricted to just the two low bins.Table 7: Raw and completeness corrected luminosity function for innerfield du Pont:if.Ncounteda error N^, acomp.ete error13.75 1.1 1.1 1.1 1.114.00 1.1 1.1 1.1 1.114.25 5.2 2.5 5.2 2.514.50 7.3 2.9 7.3 2.914.75 5.2 2.5 5.2 2.515.00 7.5 3.0 7.5 3.015.25 9.4 3.1 9.4 3.115.50 19.0 4.5 19.0 4.515.75 30.5 5.8 30.5 5.816.00 40.0 6.5 40.0 6.516.25 53.4 7.6 53.4 7.616.50 67.9 8.6 67.9 8.616.75 79.1 9.5 79.1 9.517.00 103. 11. 103. 11.17.25 118. 12. 118. 12.17.50 110. 11. 113. 12.17.75 110. 11. 114. 12.18.00 128. 12. 115. 13.18.25 133. 12. 122. 14.18.50 142. 13. 153. 17.18.75 153. 14. 163. 18.19.00 168. 14. 183. 19.19.25 171. 14. 193. 20.19.50 181. 15. 197. 19.19.75 214. 17. 217. 25.20.00 232. 17. 237. 29.20.25 241. 18. 293. 33.20.50 241. 17. 307. 42.20.75 234. 17. 314. 40.21.00 222. 17. 341. 55.21.25 192. 16. 322. 53.aNumber per 0.5 magnitudes per 4.67 sq. arc minute field.78Table 8: Raw and completeness corrected luminosity function for outerfield du Pont:n.I aNcounted erroraNcomplete^error15.00 4.1 2.0 4.1 2.015.25 8.2 3.0 8.2 3.015.50 4.1 2.0 4.1 2.015.75 2.0 1.4 2.0 1.416.00 5.1 2.4 5.1 2.416.25 11.1 3.5 11.1 3.516.50 18.1 4.3 18.1 4.316.75 19.7 4.5 19.7 4.517.00 22.9 4.9 22.9 4.917.25 22.6 4.8 22.6 4.817.50 28.7 5.4 28.7 5.417.75 30.3 5.7 30.3 5.718.00 34.3 6.0 34.3 6.018.25 48.1 7.2 48.1 7.218.50 37.5 6.3 37.5 6.318.75 36.5 6.2 36.8 6.319.00 49.9 7.2 50.3 7.319.25 66.8 8.3 62.5 8.619.50 66.5 8.3 69.0 8.719.75 77.2 9.0 85. 10.20.00 104. 10. 100. 11.20.25 103. 11. 99. 11.20.50 99. 10. 106. 13.20.75 108. 11. 119. 13.21.00 110. 12. 112. 15.21.25 106. 11. 104. 14.21.50 115. 11. 125. 17.21.75 118. 12. 125. 21.22.00 129. 12. 133. 24.22.25 132. 12. 166. 28.'Number per 0.5 magnitudes per 4.67 sq. arc minute field.Table 9: Raw and completeness corrected luminosity function for theouter field du Pont:w.N-counteda error Ncompletea error15.25 1.1 1.1 1.1 1.115.50 4.2 2.1 4.2 2.115.75 5.3 2.5 5.3 2.516.00 5.3 2.5 5.3 2.516.25 11.3 3.5 11.3 3.516.50 16.8 4.2 16.8 4.216.75 16.0 4.1 16.0 4.117.00 20.1 4.6 20.1 4.617.25 20.0 4.6 20.0 4.617.50 25.2 5.2 25.2 5.217.75 40.0 6.5 40.0 6.518.00 40.4 6.6 41.2 6.718.25 30.8 5.8 31.8 6.118.50 34.5 6.1 35.5 6.318.75 42.8 6.8 43.5 6.919.00 39.2 6.5 35.5 7.119.25 44.8 7.0 45.9 7.219.50 58.5 8.0 60.6 9.419.75 72.9 9.0 72. 11.20.00 78.6 9.2 82. 12.20.25 80.5 9.3 81. 13.20.50 90.4 9.9 88. 16.20.75 105. 11. 113. 15.21.00 120. 12. 153. 22.21.25 129. 12. 152. 21.21.50 118. 12. 125. 26.aNumber per 0.5 magnitudes per 4.67 sq. arc minute field.16 18 20I2220.5Fig. 19.—Comparison of the background field star counts for FRST (histogram) and duPont:bf (points). The two fields are at a distance of 1° on either side of the cluster and at the samegalactic latititude. The uncertainties in the FRST counts should be similar to those shown for the duPont:bf counts. All of the new luminosity functions are shown for two sets of overlappinghalf-magnitude bins.81Table 10: Background subtracted luminosity functions for du Pont:if(inner) and average of fields du Pont:n and du Pont:w (outer).I Ninne: error Noutera error15.25 5.2 3.8 2.0 2.215.50 14.7 5.0 0.0 2.115.75 27.3 6.1 1.0 1.916.00 34.6 7.0 0.0 2.516.25 47.0 8.0 4.8 3.116.50 58.3 9.2 7.9 3.816.75 67. 10. 6.2 4.017.00 90. 11. 8.7 4.317.25 108. 12. 10.7 4.117.50 97. 13. 11.0 4.817.75 96. 12. 17.1 5.318.00 102. 14. 23.9 5.318.25 101. 15. 18.7 5.818.50 130. 18. 14.2 5.618.75 131. 19. 7.60 6.519.00 148. 20. 7.90 6.819.25 160. 21. 20.9 7.119.50 159. 21. 26.7 7.919.75 183. 26. 44.2 8.620.00 202. 30. 56.5 9.420.25 246. 34. 43. 10.20.50 244. 43. 35. 12.20.75 246. 41. 47. 13.21.00 265. 57. 56. 16.21.25 236. 55. 43. 16.21.50 50. 19.21.75 53. 26.22.00 30. 32.'Number per 0.5 magnitudes per 4.67 sq. arc minute field. I I 1^1^1 1 i i Iii^1^if f i 1 f ff f f f f f^.--,^.f .^background field-..___..__:---- 2^------ .0I^I16 18^20^22II ,^IFig. 20.—Luminosity function for the inner field du Pont:if. The background counts areshown for comparison.83Fig. 21.—Luminosity function for the outer fields du Pont:n and du Pont:w. The counts havebeen averaged to improve the signal-to-noise. The background counts are also shown. These fieldsare about as far from the centre of the cluster as it is practical to go.84If i f ififI^I^II^I16 18 20^22^24IFig. 22.—All the NGC 6397 LFs are compared. The three LFs are for the inner field (filledcircles), FRST (open circles) and the averaged outer fields (squares). For the two LFs from thisstudy, two sets of overlapping half-magnitude bins are shown.85In order to convert the luminosity functions into mass functions a mass-luminosity law isrequired. The I—band mass-luminosity law of FRST, with an extension to more massive starstaken from VandenBerg & Bell (1985), was used. This mass-luminosity relationship is for Y=0.2,Z=10-4, and a time of 16 Gyr. Table 11 lists the two new mass functions. These are shown,together with the MF of FRST, in Fig. 23.Table 11: Mass functions for du Pont:if (inner) and average of fields du Pont:n anddu Pont:w (outer). Mr m Nma error Nma error(M0) (inner) (outer)3.50 0.804 860 2903.75 0.794 1240 280 47 884.00 0.782 1310 270 0 944.25 0.768 1500 250 154 984.50 0.751 1560 250 210 1004.75 0.731 1570 240 143 935.00 0.708 1890 230 183 895.25 0.684 2160 240 214 825.50 0.657 1830 240 207 895.75 0.631 1740 230 309 966.00 0.605 1920 260 450 1006.25 0.576 1820 270 340 1006.50 0.549 2260 310 246 986.75 0.517 2070 310 120 1007.00 0.484 2040 280 110 947.25 0.449 2170 290 285 977.50 0.410 2100 270 350 1007.75 0.371 2250 320 540 1108.00 0.332 2510 370 700 1208.25 0.293 3310 460 580 1408.50 0.260 3950 690 560 2008.75 0.232 5150 870 990 2709.00 0.209 6600 1400 1400 3909.25 0.192 7000 1600 1260 4709.50 0.176 1890 7009.75 0.164 2400 120010.00 0.154 1600 1700aNumber per unit mass per 4.67 sq arc minute field.86 J^I —0.2 —0.4^—0.6log m (M0)—0.8-0cdzt!. 3cCfC/)N1.0 2Fig. 23.—All the NGC 6397 MFs are compared. The symbols and bins are as in Fig. 22.87From the more massive, upper main-sequence, stars it is difficult to see any sign of masssegregation, although there is some indication that the outer field is more deficient in higher massstars than the two fields closer to the centre. For the less massive stars, there is certainly animpression that the outer field contains relatively higher proportions of these than do the two innerfields. For these stars there is a somewhat clearer signal for mass segregation. Table 12 gives theslopes of weighted, least-square, power-law fits to the mass bins with masses lower than 0.4,0.32, and 0.28 M 0 , together with the numbers of points fit and error estimates. The final point inthe du Pont:out MF has been excluded from this fitting due to its large uncertainty. For the duPont:if and FRST MFs the three slopes are quite consistent with one another, and are alsoconsistent with there being little mass segregation between the two fields. As the high-mass cutoffis moved to lower masses the MF from the outer fields gets steeper. The evidence for masssegregation is strongest for the stars with M<0.32 M . but is not highly significant. There is alsoan indication in the three MFs that the mass of the break in the MSI decreases with distance fromthe cluster centre. For the du Pont:if MF the upturn in the MF takes place at around 0.4M 0 , for theFRST MF this occurs between 0.4M 0 and 0.3M0 , and for the du Pont:out MF the data suggestsan upturn nearer to 0.3M 0 .To some extent this is borne out by the power laws fit to the MFs.Table 12: Observed mass spectral indicesm<0.4M0^m<0.32M0 m<0.28M0Mass Function limitingMSI #^MSI #^MSI # mass (M0 )du Pont:if 0.8±0.3 7 0.9±0.5 5 1.0±0.9 4 0.19FRST 1.0±0.2 8 1.1±0.3 7 0.9±0.4 6 0.12du Pont:out 0.5±0.3 9 1.4±0.6 7 1.8±1.0 6 0.16Constructing the segregation measures S„ as defined in eq. (3.1), might aid in quantifying thedegree of mass segregation by focusing on the relative differences in the the mass functions. Thesegregation measures between the three pairs of fields are shown in Fig. 24. One difficulty withthese diagrams are the large errors associated with taking the ratio of two, already somewhat88poorly determined, quantities. To a certain extent the conclusions about mass segregation madefrom the MFs themselves are supported by the segregation measures, but most of the scatter in thesegregation measures is probably due to observational uncertainties and only weak conclusions canbe drawn from them.8911,,,I,III,,,,,III lllll 'i l l^.. 1,11),,,,i1 ^ 1^1S,(m,du Pont:out; du Pont:if)_. S r(m. du Pont:out; FRST)^Sr(rn. FRST; du Pont:in)11^111 1 111111V111111111111—0.5—1.5^ ItititItIttittlil^[till ^ Ilititill ^ Id—0.2 —0.4 —0.6 —0.8 —0.2 —0.4 —0.6 —0.8 —0.2 —0.4 —0.6 —0.8log m (Me)^log m (MG)^log m (M0)Fig. 24.—The segregation measures between the three pairs of fields are shown. The segregationmeasure, Sr, is defined by eq. (3.1).90Chapter 5. Models Including Stellar EvolutionThe focus of this chapter will be on attempting to find Fokker-Planck models, employing allthe features discussed in ch. 2, which can reproduce the observations of M71 and NGC 6397. Noattempt will be made at surveying the set of all possible parameters. Rather, a set of models withplausible initial conditions, and which survive to an age of approximately 15 Gyr, are to becompared with the observations. The properties listed in Table A4 for the first grid of modelsseem quite reasonable, as will be discussed presently. When actually compared with the observationsof M71 and NGC 6397, they quickly prove to be wide of the mark. The lack of remotelysatisfactory matches to either cluster necessitated a wider search through the parameter space,especially in terms of the IMF. Before turning to the results of the modelling, I will discuss someof the motivation behind the initial choice of parameters in order to give a sense of the issuesinvolved.An initial mass of order 106 Mo is suggested by the Jeans mass of primordial gas at 10 4K, thestarting point of a proto—globular-cluster according to some theories of cluster formation (Fall &Rees 1985). The typical present-day mass of a galactic globular cluster is few times 105 M 0 . Theyoung populous clusters in the LMC have masses estimated to be 10 4 to 106 M 0 (Elson, Fall, &Freeman 1987). Holtzman et al. (1992), using the Hubble Space Telescope, recently detected apopulation of luminous blue objects in the galaxy NGC 1275. A new set of observations by Richeret al. (1992) confirms the existence of these objects. The most probable interpretation of theseobservations is that the objects are young globular clusters. Richer et al. estimate the mean clustermass to be around 10 6 Mo . Allowing for inefficiencies in star formation and the amount of masslost via stellar evolution and the tidal boundary, 10 6 Mo is a reasonable initial mass for clusterssuch as M71 and NGC 6397.Observations of the young populous clusters in the LMC are also providing information on theIMF for stars more massive than those remaining in the galactic globular clusters. Mateo (1992)finds that, on the basis of the available evidence, all of the MFs observed for LMC clusters areconsistent with a mean mass spectral index <x>=1.7±0.2 for stars more massive than about1.5 Mo . Furthermore, the same mean slope, with a somewhat larger spread, appears to fit the91MFs in many diverse environments. Based on this observation, the initial high mass IMF had aMSI of xh=1.7. A flatter IMF, which ultimately was required for M71 and NGC 6397, is supportedby the observations of Elson, Fall, & Freeman (1989). In six young LMC clusters, they found thatbetween 1.5 and 6 Mo the MSI lay between —0.2 and 0.8.The clusters with observed deep MFs all show a break in their MFs, with an increase in theMSI for masses below about 0.4 M o (Richer et al. 1991). These show a variety of slopes whichhave been affected by mass segregation and tidal stripping. The steepest of these have a MSI ofx=2.7. A recent observation of the MF in the galactic halo shows that the MSI there is 3.5 (Richer& Fahlman 1992). On the basis of these observations, and assuming such steep slopes areappropriate to all clusters, the initial MSI for the stars with M<0.4M0 was taken to be x 1=3.5. Theobserved IMFs for stars from the break up to the turnoff have smaller mass spectral indices. Forconvenience the high mass IMF is continued to 0.4 M o . Assuming that the transitions in the MFscollected in Richer et al. (1991) are not due to a problem with the mass-luminosity relation, theMFs do not have the appearance of being the sum of two overlapping power-laws, but that a muchsharper transition takes place. Hence the IMF has been constructed to be continuous but notsmooth at the mass bin with initial mass 0.44 M o (q.v. Table A5), with one power-law on eitherside of the break.The mass range has been taken to be between 20 Mo and 0.1 M 0 , and 25 mass bins wereused. The mass limits are somewhat arbitrary, but are governed by the observed range of masses.Constructing the IMF with an upper limit of 100 Mo and xn=1.7 gives an additional 0.4% of massloss in the first 3 to 6 Myr, and 10% more neutron stars. Extending the low mass cutoff, but withthe same total mass, will reduce the fraction of mass lost through stellar evolution and will lengthenthe relaxation time, but these are not large effects. The final masses and stellar lifetimes used arethe same as those used by CW. It is assumed that all the neutron stars are retained by the cluster.The models start as King models with various concentrations and at several galactocentric radii.They are assumed to already fill their Roche volumes when the models are started. Questions can,and will, be raised about all of these choices of parameters, but the choices are not implausible.92The analysis of the models were done in the following manner. In view of the uncertainty inthe relaxation timescale, any model which survived past 10 Gyr was examined for comparisonwith the two clusters. All the saved steps between 10 Gyr and 20 Gyr were used for this. Themodels were converted to the observable plane for M71 in the same manner discussed in Ch. 3,and examined visually to ascertain the quality of the match. The mass bin with mass 0.87 M o wasused for the SDP for stars between 0.85 and 0.89 M 0 , and the next bin, with mass 0.78 M 0 wasused for the SDP from the stars between 0.71 and 0.85 M o . I will refer to the bin with mass0.87 Mo as the "giant bin" as it contains the stars currently evolving off the main sequence andwhich dominate the light of the cluster. The segregation measure S,,(m, 3'; core) was also calculated.A similar procedure was followed for NGC 6397. Mass functions were calculated for rectangularareas at the locations of the du Pont:if, du Pont:n, du Pont:w and FRST fields (Table 3) and withareas of 1.35 pc x 1.35 pc. The model mass functions for the du Pont:n and du Pont:w 'fields'were averaged as was done with the observed MFs. One problem with using these models forNGC 6397 is that low metallicity stars (for NGC 6397 [Fe/I-1]=-2.0 vs. [Fe/1-1]=-0.5 for M71[Webbink 1985]) evolve faster than those with higher metallicities. The mass at the turnoff for themass-luminosity relationship used for NGC 6397 is only 0.83 M o as opposed to 0.87 M 0 used inthe models. For the purposes of comparison the giant bin has been used to model the shape of theSDP of NGC 6397, but an offset between the observed and model densities should be expected.The segregation measure S,.(m, du Pont:out; du Pont:if) has also been computed for each model.A diagram comparing the simulated observations and real data was constructed for each modeland cluster and included all the saved time steps between 10 and 20 Gyr. Simplified versions ofthis diagram are shown below. The full diagram was examined by eye and the time which appearedto best fit the observations was selected. With, the exception of the NGC 6397 SDP, where adensity rescaling was permitted because of the mismatch in the bin mass and width and theuncertainty in the fraction of the turnoff and giant stars included in the SDP, there are no freescaling parameters in the models to actually "fit". When the term "fit" is used in what follows, onlythis sort of visual matching is implied; no other adjustments are possible.93Section a. M71Figure 25 shows the most reasonable of the models for M71 from the series of Table A4. Thespecific parameters for this model, designated model M71—A, are listed in Table A5. The degree ofmass segregation in the model is in satisfactory agreement with the observations. The SDPscontain too many stars, indicating too high an initial mass, and have a peculiar shape, but the innerregion has a shape similar to that of the observed SDPs. The model MFs are far too steep to matchthe observations. For the models from Table A4 the effects of stellar evolution are quite irrelevant;only 3% of the initial mass is lost through stellar evolution mass loss. A model run with the MFpre-evolved to its value at 15 Gyr gave very similar results. If M71 is to be fit with a stellarevolution model it appears to require a much flatter IMF and a smaller initial mass.With this in mind several dozen further models were generated in an attempt to find a better fitto M71. Some of the choices in parameters to try were driven by the desire to fit NGC 6397, but asthe observed MFs and galactic positions of the two clusters are similar, such choices will serveequally well for M71. In addition to the initial concentration and radial size, the mass spectralindices and the initial mass were also varied. Besides 10 6 M 0 , initial masses of 7.5x105 M, and5x105 Mo were used. For the high-mass stars, mass spectral indices, xh , down to 0.2 were tried.For the low-mass stars the MSI, x,, ranged from 1 to 3.5. The strongest sensitivity shown by themodels was with respect to the IMF. Table B1 gives details of the IMFs used and the sets of initialparameters which gave models surviving to 10 Gyr or more. The initial high-mass MSI for M71was probably not as flat as xh =0.2. All of the models with this value of xh were destroyed fairlyquickly despite having an x i=3.0. For this IMF, 59% of the total initial mass is lost through stellarevolution. The models with x h=0.5 survived only if they were initially sufficiently concentrated:models with W05_4 were destroyed, models with Wo=6 or higher survived. Tidal radii (as definedfor a cluster with mass 10 5 Mo [q.v. eq. 2.21]), near 20 pc appear to give the best results both forM71 and NGC 6397.As an aid to identifying models with properties similar to those observed in M71, I calculatedfor each model the normalized segregation measure S *f7 1 (eq. 3.3) and the radius where the density94—0.5^0^0.5log r (pc)1—0.5^0^0.5log r (pc)!mil), 1 1^111 1 _0.71<m/Mo<0.851^iirr^11111111^1110.85<M/Mo<0.89 -I^1^I^I^I^1^I^I^I^I^I^I^I^I^I^I^ICaE 6• 5a)E• 4Z 3—1 —0.8-0.6-0.4-0.2 0log m (M0)—0.8—1.4—1 —0.8-0.6-0.4-0.2 0log m (M0 )Fig. 25.—Comparison of model M71—A with the observations of M71. The top panels showthe two observed SDPs together with the model SDPs. The lower-left panel shows the observedand model MFs for the two fields (filled circles, core; open circles, 3'). The lower-right panelshows the observed and model segregation measures.95of stars in the giant bin reaches half its central value. This is being used as a stand-in for theclassical 'core radius', the radius at which the surface brightness reaches half its central value(King 1962). This definition of core radius has since been overwhelmed by the use of that namefor King's scaling radius (eq. 2.12), so I will refer to the radius used here as the half-densityradius, rm. From the SDP in Fig. 4, the observed value for M71 is about 1 pc.Figure 26 shows the evolution of .SrMm7 1 and log rhd between 10 and 20 Gyr for all of the modelscomputed. Also indicated is the observed values of Srm,h71 =0.4 and rhd=1. It appears from thisfigure that none of the models can be expected to provide a good match with M71, and this isborne out by more detailed comparisons. The systematics of S,,,h—rhd diagram are discussed morefully in Appendix B. Two conclusions can be draw from Fig. 26 and the diagrams in Appendix B.The first is that IMFs with MSIs that can evolve into the MFs observed in M71 do not produce theobserved amount of mass segregation between the two fields. The breakdown of Fig. 26 by IMFshown in Fig. B 1, shows that the only IMFs producing a large value of .5?. "in in are those withxh =1.7. It is one of these models which are shown in Fig. 25, and the required reduction in theMSI for the upper main sequence stars does not occur.The second conclusion can be seen in Fig. B2 which shows the entire SM71 vs. rhd curves for ar,series of models. During the epoch when stellar evolution is important the models evolve to largerrhd at almost constant SrMn . Stellar evolution does not provide a way for mass segregation to takeplace before extensive core contraction has occurred.Model M71—B (Table A6), shown in Fig. 27, is typical of the models with the observed rhdand a more appropriate MF. This model has too high a mass to fit M71 and the MFs are still toosteep. The segregation measure is only half that observed. It should be possible to find a modelwith lower mass and flatter IMF to fit the SDPs, but both of these changes tend to decrease therelaxation time and bring core collapse sooner. Having a higher proportion of high mass starscauses a greater initial expansion but also accelerates the loss of mass through the tidal boundary.One way around this problem is to begin the evolution of the model with a limiting radius smallerthan the tidal radius at the desired galactocentric distance. This will reduce the rate of mass lossthrough the tidal boundary and slow the evolution of the model. A model run with the type of960.60.4U)0.20—3^—2^—1^0^1log rhdFig. 26^S,..—rhd diagram for M71. The model segregation measure has been calculated for thesame fields as were observed. The diagram depicts the evolution of the segregation measure andhalf-density radii with time. The evolutionary tracks are plotted for all models which have notdisrupted before 10 Gyr, and extend to 20 Gyr. Evolution begins at zero segregation measure andgoes in the direction of increasing S r and decreasing rm. Evoldtion at constant S r . occurs duringcore collapse, and the post—core-collapse evolution is generally in the direction of decreasing S,and increasing rhd. The observed value of S r,. and rhd for M71 is indicated by the snowflake. Thetrack for model M71-A lies to the left of the M71 point.971—0.5^0^0.5log r (pc)1—0.5^0^0.5log r (pc)_IIII0.71 <m/Me<0.850.85<m/M0<0.89 -I^I^t^I^1^I^I^l^I^I^I^I^I^I^I^1^I^IE 64)5a)43—1 —0.8-0.6-0.4-0.2 0log m (Me)—0.8—1.4—1 —0.8-0.6-0.4-0.2 0log m (M0)Fig. 27—As Fig. 25 for model M71—B. The shapes of the surface density profiles are reproduced,although at higher density, but the mass segregation is too small. The mass function slopes for the3' field are not too different, but the core MF slopes are poorly matched. The S r m vs. rhd curve forthis model is among those lying below the point for M71 in Fig. 26.98initial parameters just outlined will not have any higher a segregation measure so the problem ofM71 still stands.There are still possible solutions which could allow M71 to be matched by a Fokker-Planckmodel. It is yet to be seen how primordial binaries affect a multi-mass model. Tidal shocking, bythe disk and the bulge, has been proposed as method of heating globular clusters (Ostriker, Spitzer& Chevalier 1972) and could have played an important role in M71' s evolution. Alternatively, theproblem may not lie so much in the Fokker-Planck modelling as in the initial conditions. Ifglobular clusters start mass segregated, then the conclusion derived here for models which startwithout mass segregation will not apply. In the face of continuing difficulty in finding a model forM71, additional observations of the mass function are required in order to confirm the large degreeof mass segregation in the present data. The file is not closed on M71, and the possibility ofbuilding a Fokker-Planck model, perhaps with additional physics, still remains.Section b. NGC 6397One might reasonably expect that it would be easier to find a model matching NGC 6397 thanone for M71, and this is the case. As a first step, a S r,,h—rhd diagram from the same set of modelshas been constructed for NGC 6397 using the normalized segregation measure S,...(0.22M o , duPont:out; 0.77M 0 , du Pont:if). The observed segregation measure is, from Fig. 24, between 0.3and 0.4. An upper limit on rhd of 0.2 is based on the mean core radius of the multi-mass Kingmodels of Meylan & Mayor (1991). The S ,„,--rhd diagram for the NGC 6397 observational windowis shown in Fig. 28 with the estimated observed range delimited by a dashed box. There are plentyof models with properties appropriate to the observations. If the density of models within the boxseems high, this is due to the many models run in refining the model IMF for the NGC 6397.The most satisfactory looking model is model NGC 6397—A, with initial parameters listed inTable A7. The best matching of the times stored for this model is shown in Fig. 29. It should benoted that this is just a representative comparison and that a better match is possible in the intervalbetween this saving of the model and the next. The model accords well with the observed MFs.The model SDP has been normalized by eye to the observed profile out to 1 pc and has been99Az0 . 84.;0P-4 0 . 60:is 0.40a.,0 . 2OC\1C\20—3 —2^—1^0log rhd1Fig. 28.—As Fig. 26, but for the NGC 6397 normalized segregation measure S„(0.22M 0 , duPont:out; 0.77M 0 , du Pont:if). The difference between the tracks followed by the models in thisfigure and in Fig. 26 lies in the different windows through which the models are observed. Theoutlined box is the range of estimated values of S„ and TM for NGC 6397. Many good candidatesfor a NGC 6397 model lie within this region.1001^1^11^1^1^1^1^1^1^1^1^1^1^1^1^1—1 —0.5 0^0.5log r (pc)- 1^1 1^1^1^1^1^1^1^1^1^1^1^1^1^1FRST MF1^1^1^I^1^1^1^1^1^1^1^1^1^1^1^1^1 —1 —0.8 —0.6 —0.4 —0.2 0log m (M0)Ez0 — 0.4*r;r: —0.6-o —0.8-1— 1.2I.—1.4U)I^1^1^1^1^1^I^I,du Pont MFs -1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1__1—1 —0.8 —0.6 —0.4 —0.2 0log m (M0)—1 —0.8 —0.6 —0.4 —0.2 0log m (M0)Fig. 29.—Comparison of model NGC 6397—A with the observations described in Ch. 4.Top-left: Composite surface density profile including both the R14 and R15.5 data. The modelSDP has been decreased by 0.6 dex to match the inner profile. Top-right: MF of FRST. Bottom-left:New MFs. The upper one is for the du Pont:if field and the lower one (open circles) is the mean,du Pont:out MF. No normalizations have been applied to the model MFs and the fits are quitegood. Bottom-right: Segregation measure between the du Pont:if and du Pont:out MFs. Again, themodel presents a good representation of the observations.101decreased by a factor of four. With the exception of this normalization, none of the model resultshave been adjusted in any way to improve the fit. The absolute normalizations of three massfunctions are very well reproduced, as are their overall shapes. The good match for the segregationmeasure confirms the more detailed MF comparisons.Meylan & Mayor (1991) presented a line-of-sight velocity dispersion profile for NGC 6397. InFig. 30 the velocity dispersion profile for model NGC 6397—A is compared with these observations.A satisfactory match is found, although it appears that the model velocity dispersion is increasingmore rapidly towards the centre than the observed velocity dispersion. If this flattening is supportedby further observations, the discrepancy in the model may be explained as being due to theassumption that all neutron stars were retained by the cluster. The result of this assumption is togive the model a higher central velocity dispersion than would be the case if a sizable fraction of theneutron stars were ejected. That the model velocity dispersions lie systematically higher than theobservations is not a concern since continuing mass loss in the model will bring them into betteragreement.At the preferred time, 16.6 Gyr, model NGC 6397—A is undergoing core collapse and is 40Myr before the density maximum. The current mass is 10 5 Mo , of which 9% is in neutron starsand 36% in white dwarfs which have a mean mass of 0.8 M o . There are 7,000 neutron stars and48,000 white dwarfs present out of 2.6x105 stars in total. The global mean mass is 0.4 M o and inthe centre, dominated by degenerates, the mean mass is 1.3 M 0 . The half-mass radius is 4 pc andthe tidal limit is at 20 pc. The half-mass relaxation time is 1 Gyr.It is of interest to look at the changes in the MSIs and these are listed in Table 13. The globalMF has flattened considerably and this is reflected in the MFs for the three fields individually. Themodel MF for du Pont:if is flatter than the global MF but the model MFs for FRST and du Pont:outhave steeper values of xh than the does the global MF. The values of xh at these two positions aresimilar to that of the global MF. Mass segregation is seen most strongly for the more massivestars. This is a difficulty with respect to the observations, as the MFs for the high mass stars aremore uncertain than that for the lower mass stars due to the lower numbers observed.102.^I,^.^.^I^i^i^i.^I^i I^1—0.5 0 0.5 1log r (pc)Fig. 30.—Comparison of the velocity dispersion profile of model NGC 6397—A with theobserved velocity dispersions of Meylan and Mayor (1991). The model reproduces the observationsadequately although there is some indication that it is rising more quickly towards the centre than isobserved. Differences in the handling of the heavy remnants, for example ejecting some fraction ofthe neutron stars, will reduce this difference.103Table 13. Mass spectral indices for model NGC 6397—AMass Function x, xhIMF 1.5 0.9Current global 1.0 —0.1du Pont:if field 0.8 —0.2FRST 0.9 0.3du Pont:out field 1.1 0.7The apparent shift in the mass of the break in the observed MFs is not seen in the model. Theevolved MF remains fairly faithful to a power-law form for masses both above and below thebreak and does not show any movement in the position of the break. It is difficult to see how theshape of the MF can be changed by mass segregation to the extent that a feature like the break canmove about in magnitude.There were several other models that gave similarly satisfactory fits to the SDP and thesegregation measure. Where these models failed was in the reproduction of the details of MFs.Differences of a few tenths in the mass spectral indices were enough to spoil the match. Someleeway is also permitted in the distance of the initial tidal boundary.This model is not the final word on NGC 6397 by any means. Many details of the comparisonwith the observations are still unsatisfactory. For r >lpc the observed SDP is much straighter thanthe model SDP. The detailed structure of the MFs are not reproduced using such a simple form forthe IMF and it will be a delicate fine-tuning problem to get a better fit. With respect to the input datain this model, the stellar lifetimes are inappropriate for low metallicity stars, lending some uncertaintyto the result. Adjusting the stellar lifetimes will have small quantitative effects but the qualitativebehaviour of the models should be the same. In the face of the inability to find a model matchingthe M71 observations, it is an encouraging sign that such a good match can be found for NGC6397. The wide range of observations the model is able to reproduce is certainly an indication thatsuch Fokker-Planck modelling does represent an aspect of globular cluster evolution.104Chapter 6. ConclusionsThe aim of this thesis was to test a class of Fokker-Planck model with observations of globularclusters. The two clusters used, M71 and NGC 6397, have different density profiles but havesimilar masses and mass functions. From the perspective of the Fokker-Planck models, the differencesare more important than the similarities. I was unable to find a model which matched the observationsof M71, but the results of this comparison provide some illuminating insights into how to proceed.The comparison with NGC 6397 was a great success and confirms that these Fokker-Planckmodels are appropriate for some classes of globular clusters.The original problem encountered for M71 in the models with pre-evolved MFs (ch. 3) remainswhen models including stellar evolution are used. The addition of stellar evolution serves to delaycore collapse beyond the time expected from relaxation time arguments alone. On the other hand,as discussed in Appendix B, the models which suffer the largest degree of expansion due to massloss from stellar evolution are also the models with the lowest segregation measures for the fieldsobserved in M71. It appears that stellar evolution is not the solution to the puzzle posed by M71.There are several parameters in these Fokker-Planck models which can be traded off againstone another in attempting to improve the match with a set of globular cluster observations. One ofthe key problems lies in the realm of time scales. In order to even be tested with the observations Irequired a model to survive until at least 10 Gyr. Consider a model which is assigned an initialmass, an initial concentration, a scale radius, an IMF, and a position in the galaxy. The mass,concentration, and scale radius set the relaxation time. The IMF contributes to two effects. The firstis the degree of expansion and subsequent delay in the time of core collapse due to mass lossduring stellar evolution. The second is the degree to which mass segregation accelerates corecontraction. The position in the galaxy governs the size of the tidal boundary and the rate of tidalmass loss. When the initial limiting radius of the model is smaller than the initial radius of the tidalboundary, the tidal mass loss rate will be reduced, leaving the model with more mass at a giventime than if the model initially fills its Roche volume.Changing any of these parameters has side effects which can be offset somewhat by varyingother parameters. In model M71 —B (Fig. 27) we saw that the model SDPs were much higher than105the observed SDPs. It would be desirable, in order to improve the match, to reduce the total massof the model, say by a factor of two. If the initial limiting radius is constrained to be the same asthe initial tidal boundary, then the relaxation time also is reduced by about a factor of two. (Thisfollows from the proportionality t recm1/2rh3/2 [eq. 2.11]. The constraint on the limiting radius andthe definition of the tidal radius [eq 2.20] requires recM1/3, hence trhOCM. with this constraint.) Fora model like M71-B, but with half the mass, core collapse occurs at 12 Gyr; at 10 Gyr the model isalready too concentrated to give a good match. Instead of reducing the total mass, we could justredistribute it to stars with masses below the minimum mass observed. This adjustment in the IMFwould also cause there to be less expansion due to stellar evolution, and possibly more mass lossthrough the titdal boundary. These two effects would shorten the time to core collapse somewhatbut the full effects of such adjustments to the IMF have not been explored. Having a larger tidalboundary reduces the rate of tidal mass loss, lengthening the time to core collapse.These are relatively minor adjustments compared to the effects of IMF slopes. The IMF can bepartially constrained by the observed MFs, but large ranges, especially the high mass end, aredifficult to observe. Some constraint is possible through counts of degenerate stars, eg. millisecondpulsars, but these will be incomplete. Correction of the observed counts for things such asselection effects, rate of formation from the parent population, and the fraction of the remnantpopulation retained by the cluster, will make interpreting these counts difficult. Given the ability tomake trade-offs between all these parameters, it should not be expected that the model matchingNGC 6397 is unique. More work is required to assess the size of the allowable adjustments.For NGC 6397 a well fitting model was indeed found. Many of the models matched the SDPand the segregation measure adequately, but some fine tuning was required to reproduce the MFs.The model NGC 6397—A also gave a good match to the velocity dispersion profile. With theirmodels of M15 and NGC 6624, GCLM showed that Fokker-Planck models were capable of fittingthe surface brightness profile and velocity dispersion profile of a core-collapsed cluster. The modelpresented here for NGC 6397 extends the range of data which can be fit simultaneously to includeseveral local mass functions and the SDP for the brightest stars. In the face of the problem of M71,this is an important achievement. The success in modelling NGC 6397 is all the more remarkable106in that, with the exception of the SDP offset, no adjustments were allowed in the model output inorder to fit (literally, this time) the observations. The model results were transferred onto theobservational plane and plotted together with the observations. The mass function results areparticularly satisfying in view of the sensitivity to the IMF shown by the models.While including stellar evolution in the models has not helped the situation with respect toM71, it is a more realistic treatment of one of the physical processes occurring in globular clusters.The effect on the models of varying the stellar lifetimes and remnant masses has not been explored.These will have an impact on the initial evolution of the model, but, by the present epoch, the maineffects will be seen through the relative masses of the turnoff stars and the remnants. It is ashortcoming in the model of NGC 6397 that the mass bin used to model the SDP has a slightlyhigher mass than the currently evolving stars in the cluster, but this can be rectified with moreappropriate data for the stellar evolution of low metallicity stars.In ch. 4 I noted several similarities and differences for the two clusters M71 and NGC 6397.The answer to the question of why the Fokker-Planck models can easily fit the latter, but not theformer, will lie in the differences in the initial conditions and in the differing environments. Thepresent mass functions in the two clusters are similar and suggest that the IMFs were similar. Itmight be that M71 started with a higher proportion of primordial binaries than did NGC 6397. Theadvantage of primordial binaries over stellar evolution as a heating source is that they promotemass segregation even while preventing an extreme core collapse. The models with primordialbinaries of Gao et al. (1992) support this idea. Mass loss due to stellar evolution, acting as it doesthrough the potential, causes a more indiscriminate expansion. The effect of primordial binaries iscertainly worth pursuing in multi-mass models.Another possible intrinsic solution would be if the initial conditions used here are not appropriateto globular clusters. For example, a globular cluster may start mass segregated. In such a case therelationship between mass segregation and concentration is not as straightforward as the Sr,„—rhddiagram would suggest. As well, if the massive stars start centrally concentrated, then the effectsof their losing mass will be amplified. More observations of candidates for young globular clusters,107such as the ones in the LMC or NGC 1275, and a theory for globular cluster formation, wouldhelp to resolve these issues.Since M71 is a disk cluster, as opposed to NGC 6397 with its halo orbit, it will be subject tothe tidal effects of the disk and bulge to a much greater extent. For example, it is more likely toencounter a giant molecular cloud than NGC 6397 is. One might speculate that M71 has sufferedfrom some catastrophic interaction which has led to the rapid expansion of the core. If the expansionoccurs over a time much shorter than the relaxation time, then the previous degree of masssegregation will be preserved even while the core radius increases.Alternatively, errors in the observed mass functions could be giving a much higher measurementof the mass segregation than is really present in M71. Additional observations of the cluster,especially additional mass functions, would go a long way to helping to understand its physicalstate.Of the clusters which have been compared with the type of detailed Fokker-Planck model usedhere, three, M15, NGC 6624, and NGC 6397, have central cusps in their surface densities, whilethe other, M71, is reasonably well fit by a King model. Certainly M71 does not show the strongcentral cusp predicted by the models for a core-collapse or post—core-collapse cluster. It is the threecusp clusters which have successfully been matched by Fokker-Planck models but which are notwell fit by multi-mass King models (e.g. Meylan & Mayor 1991 for NGC 6397). It is premature todraw any conclusions from such a small sample, but the contrast is striking The King modeldistribution function (eq. 2.32) is based on an approximate solution of the Fokker-Planck equation(King 1965). The extension to multiple masses introduced by Gunn & Griffin (1979) is based onan assumption of equipartition of energy between the mass classes at the cluster centre. As pointedout by Pryor et al. (1986), the models don't satisfy this assumption. The low mass stars are out ofequipartition with the high mass stars, but the discrepancy is similar to that seen in the Fokker-Planckmodels. The two types of model are based on the same equation but under different sets ofassumptions and each appears to work best with different types of clusters. The static snapshotgiven by multi-mass King models accords better with the apparently unevolved clusters, but cannotproduce models with the high concentrations that the dynamic Fokker-Planck models can. However,108it has yet to be shown that King models can also reproduce a set of mass functions for a cluster ashas been done here for NGC 6397. The entire question of the relationship between the type ofFokker-Planck model presented here and King models has yet to be explored.With the mixed results obtained here, there remains plenty of scope for working with otherclusters. The use of an appropriate S r,m—rh d diagram from a grid of models will be a useful tool infinding candidate models. Fine tuning, especially of the IMF, will still be required in order toobtain an optimal match. The results for NGC 6397 are very encouraging but also indicate thatmore work is required on the effects of varying the model parameters. The sensitivity of themodels to different stellar lifetimes and remnant masses needs to be investigated. Also required is amore thorough understanding of the effects that variations in the IMF such as changing the upperand lower mass limits and the shape of the mass function have on the models.The effects of primordial binaries will need to be added to the models. Gao et al. (1991) havemade a start on this, but it will need to be extended to the multi-mass case. For a full treatment, itwill be necessary to put the finite radius effects of tidal-capture binaries and mergers back into thecode, but this time for the multi-mass case. The effects of tidal shocks should also be investigated.The success here in fitting NGC 6397, and that of GCLM in fitting M15 and NGC 6624, givessupport to the theory that orbit-averaged Fokker-Planck models have something to tell us aboutglobular cluster evolution. This gives us confidence in using such models in interpreting observationsof individual globular clusters and of the globular cluster system. With a clearer understanding ofthe effects of IMF variations, Fokker-Planck models can be used in relating local to global MFs.The models will also be of aid in relating the present properties of globular clusters both to therange of initial conditions which can lead to them, and to their ultimate fates. Studies such asRicher et al. (1991) and Capaccioli, Ortolani, & Piotto (1991) have suggested correlations betweencluster properties and the galactic environment. Fokker-Planck models may prove useful in explainingthese correlations. However, the utility of the present Fokker-Planck models in this regard may besomewhat limited due to the approximations of isotropy and spherical symmetry and the simplifiedtreatment of tidal effects. In pursuing such enterprises the difficulty in modelling M71 stands as acaution; the success in modelling NGC 6397 as an encouragement.109ReferencesAlcaino, G., Buonanno, R., Caloi, V., Castellani, V., Corsi, C. E., Iannicola, G., & Liller, W.1987, AJ, 94, 917Applegate, J. 1986, ApJ, 301, 132Anthony-Twarog, B.J., Twarog, B.A., & Suntzeff, N.B. 1992, AJ, 103, 1264Auriere, M. 1982, A &A, 109, 301Auriere, M., Ortolani, S., & Lauzeral, C. 1990, Nature, 344, 638Bailyn, C. D. 1991 in The Formation and Evolution of Star Clusters, ed. K. Janes (ASP Conf.Ser., 13), 307Bailyn, C., Grindlay, J., Cohn, H., Lugger, J., Stetson, P., & Hesser, J. 1989, AJ, 98, 882Bettwieser, E. & Sugimoto, D. 1984, MNRAS, 208, 439Binney, J. & Tremaine, S. 1987, Galactic Dynamics (Princeton: Princeton Univ. Press)Cannon, R. D. 1974, MNRAS, 167, 551Capaccioli, M., Ortolani, S., & Piotto, G. 1991, A &A, 244, 298Chang, J. S., & Cooper, G., 1970, J. Comp. Phys., 6, 1Chernoff D. F. & Djorgovski, S. 1989 ApJ, 339, 904Chernoff, D. F. & Weinberg, M. D. 1990, ApJ, 351, 121 (CW)Cohn, H. 1979, ApJ, 234, 1036. 1980, ApJ, 242, 765. 1985, in IAU Symp. 113, The Dynamics of Star Clusters, ed. J. Goodman & P. Hut(Dordrecht: Reidel), 161Cohn, H., Hut, P., & Wise, M. 1989, ApJ, 342, 814Cudworth, K. M., 1992, to appear in of The Globular Cluster-Galaxy Connection, ed. J. P.Brodie & G. Smith (ASP Conf. Ser.)Da Costa, G. S. 1979, AJ, 84, 505. 1982, AJ, 87, 990Dean, J. F., Warren, P. R., & Cousins, A. W., 1978, MNRAS, 183, 569Djorgovski, S. & King, I.R. 1986, ApJ, 305, L61Djorgovski, S., Piotto, G., & King, I. R. 1988, in Dynamics of Dense Stellar Systems, ed. D.Merritt (Cambridge: Cambridge Universty Press), 333110Djorgovski, S., Piotto, G., Phinney, E. S., & Chernoff, D.F 1991 ApJ, 372, L41Drukier, G. A., Fahlman, G. G., & Richer, H. B. 1992, ApJ 386, 106 (DFR)Drukier, G. A., Fahlman, G. G., Richer, H. B., & VandenBerg, D. A. 1988, AJ, 95, 1415Durrell, P. & Harris, W. 1992, to appear in proceedings of The Globular Cluster-GalaxyConnection, ed. J. P. Brodie & G. Smith (ASP Conf. Ser.)Elson, R. A. W., Fall, M. S., & Freeman, K. C. 1987, ApJ, 323, 54^. 1989, ApJ, 336, 734Elson, R., Hut, P., and Inagaki, S. 1987, ARA&A, 25, 565Fahlman, G. G. 1991, private communicationFahlman, G. G., Richer, H. B., & VandenBerg, D. A. 1985, ApJS, 58, 225Fahlman, G. G., Richer, H. B., Searle, L., & Thompson, I. B. 1989, ApJ, 343, L49 (FRST)Fall, M. S. & Rees, M. J. 1985, ApJ, 298, 18Gao, B., Goodman, J., Cohn, H. & Murphy, B. 1991, ApJ, 370, 567Goodman, J. 1983, ApJ, 270, 700. 1987, ApJ, 313, 576Goodman, J. & Hut, P. 1989, Nature, 339, 40. 1992, IAS preprint Ast 92/19, ApJ, submittedGrabhorn, R. P., Cohn, H. N., Lugger, P. M., & Murphy, B. W. 1992a, ApJ, 392, 86(GCLM)Grabhorn, R. P. et al., 1992b, to appear in Dynamics of Globular Clusters, ed. S. Djorgovski &G. Meylan (ASP Conf. Ser.)Grieve, G. 1983, Ph.D. thesis, University of TorontoGunn, J. E. & Griffin, R. F. 1979, AJ, 84, 752Harris, W. E. 1980, in IAU Symp. 85, Star Clusters, ed. J. E. Hesser (Dordrecht: Reidel), 81Henon, M. 1960, Ann. Astr., 23, 668. 1961, Ann. Astr., 24, 369. 1971, Ap&SS,14, 151Holtzman, J. A., Faber, S. M., Lauer, T. D., Groth, E. J., Hunter, D. A., Baum, W. A., Ewald,S. P., Hester, J. J., Light, R. M, & Lynds, C. R. 1992, AJ, 103, 691111Hut, P. 1985, in IAU Symp. 113, The Dynamics of Star Clusters, ed. J. Goodman & P. Hut(Dordrecht: Reidel), 231Hut, P., McMillan, S., Goodman, J., Mateo, M., Phinney, E. S., Pryor, C., Richer, H. B.,Verbunt, F., & Weinberg, M. 1992, PASP, in pressIben, I., Jr. & Renzini, A. 1983, ARA &A, 21, 271Inagaki, S. 1985, in IAU Symp. 113, The Dynamics of Star Clusters, ed. J. Goodman & P. Hut(Dordrecht: Reidel), 189Inagaki, S. & Lynden-Bell, D. 1983, MNRAS, 205, 913Inagaki, S. & Saslaw, W. C. 1985, ApJ, 292, 339Inagaki, S. & Wiyanto, P. 1984, PASJ, 36, 391Innanen, K. A., Harris, W. E., & Webbink, R. F. 1983, AJ, 88, 338King, I. 1962, AJ, 67, 461^. 1965, AJ, 70, 376^. 1966, AJ, 71, 64Laird, J. B., Carney, B. W., Latham, D. W. 1988, AJ, 95, 1843Lauer, T. R., Holtzman, J. A., Faber, S. M., Baum, W. A., Currie, D. G., Ewald, S. P., Groth,E. J., Hester, J. J., Kelsall, T., Light, R. M., Lynds, C. R., O'Neil, E. J., Jr., Schneider, D.P., Shaya, E.J., and Westphal, J. A. 1991, 1991, ApJ, 369, L45Lee, H.M. 1987a, ApJ, 319, 7721987b, ApJ, 319, 801Lee, H. M. & Ostriker, J. P. 1987, ApJ, 322, 123Lee, H. M., Fahlman, G. G., & Richer, H. B., 1991, ApJ, 366, 455 (LFR)Lynden-Bell, D. 1967, MNRAS, 136, 101Lugger, P. M., Cohn, H. N., Grindlay, J. E., Bailyn, C. D., & Hertz, P. 1987, ApJ, 320, 482Makino, J. & Hut, P. 1991, ApJ, 383, 181Mateo, M. 1992, to appear in The Globular Cluster-Galaxy Connection, ed. J. P. Brodie & G.Smith (ASP Conf. Ser.)McClure, R. D., VandenBerg, D. A., Smith, G. H., Fahlman, G. G., Richer, H. B., Hesser, J.E., Harris, W. E., Stetson, P. B., & Bell, R. A. 1986, ApJ, 307, L49McMillan, S. L. W. 1986, ApJ, 307, 126McMillan, S. L. W. & Lightman, A. P. 1984a, ApJ, 283, 801112^. 1984b, ApJ, 283, 813McMillan, S. L. W., Hut, P., & Makino, J. 1990, ApJ, 362, 522Meylan, G., & Mayor, M. 1991, A&A, 250, 113Moffat, A. F. J. 1969, A&A, 3, 455Murphy, B. W. & Cohn, H. N. 1988, MNRAS, 232, 835Murphy, B. W., Cohn, H. N., & Hut, P. 1990, MNRAS, 245, 335 (MCH)Oh, K. S., & Lin, D. N. C., 1992, ApJ, 386, 519Ostriker, J. P., Spitzer, L. & Chevalier, R. A. 1972, ApJ, 176, L51Peterson, R. C., Seitzer, P., & Cudworth, K. M. 1989, ApJ, 347, 251Piotto, G., King, I. R., & Djorgovski, S. 1988, AJ, 96, 1918Pryor, C. 1990, private communicationPryor, C., McClure, R. D., Fletcher, J. M., Hartwick, F. D. A., & Kormendy, J. 1986, AJ, 91,546Pryor, C., McClure, R. D., Fletcher, J. M., & Hesser, J. E. 1989, AJ, 98, 596Richer, H. B. & Fahlman, G. G. 1988, ApJ, 325, 218. 1989, ApJ, 339, 178. 1992, Nature, 358, 383Richer, H. B., Crabtree, D. R., Fabian, A. C., & Lin, D. N. C. 1992, AJ, submittedRicher, H. B., Fahlman, G. G., Buonanno, R., & Fusi Pecci, F. 1990, ApJ, 359, L11Richer, H. B., Fahlman, G. G., Buonanno, R., Fusi Pecci, F., Searle, L., & Thompson, I. B.,1991, ApJ, 381, 147Rohlfs, K. & Kreitschmann, J. 1988, A &A, 201, 51Rosenbluth, M. N., MacDonald, W. M., & Judd, D. L. 1957, Phy. Rev., 107, 1Salpeter, E. 1955, ApJ, 121, 161Shapiro, S. L. & Marchant, A. B. 1978, ApJ, 225, 603Spitzer, L. 1969, ApJ, 158, L139^ 1987, Dynamical Evolution of Globular Clusters (Princeton: Princeton Univ. Press)Spitzer, L. & Hart, M. H. 1971, ApJ, 164, 399Spitzer, L. & Mathieu, R. D. 1980, ApJ, 241, 618113Spitzer, L. & Thuan, T.X. 1972 ApJ, 175, 31Statler, T. S., Ostriker, J. P., and Cohn, H. N. 1987, ApJ, 316, 626 (SOC)Stetson, P. B. 1987, PASP, 99, 191^. 1990, PASP, 102, 932. 1991, in The Formation and Evolution of Star Clusters, ed. K. Janes (ASP Conf. Ser.,13), 88Stetson, P. B. & Harris, 1988, AJ, 96, 909StodOlkiewicz, J. S. 1985, in IAU Symp. 113, The Dynamics of Star Clusters, ed. J. Goodman &P. Hut (Dordrecht: Reidel), 361Sugimoto, D. & Bettwieser, E. 1983, MNRAS, 204, 19PTayler, R. J. 1986, QJRAS, 27, 367Taylor, B. J. 1986, ApJS, 60, 577van den Bergh, S. 1988, AJ, 95, 106VandenBerg, D. A. 1988 in The Extragalactic Distance Scale, ed. S. van den Bergh & C. J.Pritchet (ASP Conf. Ser., 4), 187VandenBerg, D. A. & Bell, R. A. 1985, ApJS, 58, 561VandenBerg, D. A., Bolte, M., & Stetson, P. B. 1990, AJ, 100, 445Webbink, R. F. 1985, in IAU Symp. 113, The Dynamics of Star Clusters, ed. J. Goodman & P.Hut (Dordrecht: Reidel), 541Woolley, R., Alexander, J. B., Mather, L., & Epps, E. 1961, Roy. Obs. Bull., 43, E303Wolszczan, A., Kulkarni, S. R., Middleditch, J., Backer, D. C., Fruchter, A. S., & Dewey,R. J. 1989, Nature, 337, 531Zinn, R. 1985, ApJ, 293, 424114Appendix A. Model ParametersTable Al. Parameters for model CW1Structure parameters:King model central potential, Wo^ 7 .Tidal radius^ 94.2 pcTotal mass 2.82x105M0Total number of stars^ 2 79x105Half-mass radius 10.97 pcMean stellar mass 1 01 MoHalf-mass relaxation time 3.59 GyrIMF mass spectral index^ 1 5IMF:Initialmass(M0)NumberFractionMassFractionMass binwidth(M0)0.438 0.239 0.103 0.0790.525 0.182 0.0946 0.0960.629 0.139 0.0864 0.1140.754 0.106 0.0789 0.1370.904 0.0806 0.0721 0.1641.08 0.0614 0.0658 0.201.30 0.0468 0.0601 0.241.56 0.0357 0.0549 0.281.87 0.0272 0.0502 0.342.24 0.0207 0.0458 0.412.68 0.0158 0.0418 0.493.21 0.0120 0.0382 0.583.85 0.00916 0.0349 0.704.62 0.00698 0.0319 0.845.54 0.00532 0.0291 1.006.64 0.00405 0.0266 1.207.96 0.00309 0.0243 1.449.54 0.00235 0.0222 1.7311.4 0.00179 0.0203 2.0713.7 0.00137 0.0185 2.49Remnant Stellar evolutionmass^epoch(M0) (Gyr)not evolvednot evolvednot evolvednot evolved0.559 8.13 15.50.598 4.07 8.130.646 2.04 4.070.703 1.07 2.040.771 0.589 1.070.852 0.339 0.5890.950 0.204 0.3391.07 0.123 0.2041.21 0.0794 0.1231.38 0.0513 0.07940.00 0.0363 0.05130.00 0.0257 0.03630.00 0.0191 0.02571.40 0.0141 0.01911.40 0.0110 0.01411.40 0.00851 0.0110Table A2. Parameters for the "standard" model of M71Structure parameters:King model central potential, Wo^ 4.Tidal radius^ 21.6 pcTotal mass 3 0x105 M0Total number of stars^ 1 27x106Half-mass radius 4 99 pcMean stellar mass 0 237 MoHalf-mass relaxation time 4.25 GyrIMF:Sum of two power-laws with MSIs 0.5 and 5.0 weighted 1:5 by number, and 1%by number in white dwarfs.Mass per star Number Mass Width of bin(M0) Fraction Fraction (M0)0.181 0.712 0.545 0.0630.247 0.140 0.146 0.0630.312 0.0488 0.0643 0.0630.398 0.0390 0.0657 0.1170.506 0.0162 0.0347 0.0880.589 0.0102 0.0255 0.0750.665 0.00846 0.0238 0.0770.743 0.00727 0.0228 0.0790.816 0.00523 0.0181 0.0660.869 0.00272 0.00999 0.0380.94 0.00333 0.0132 0.0801.02 0.00333 0.0144 0.0801.10 0.00333 0.0155 0.080Table A3. Parameters for the "black hole" model of M71Structure parameters:King model central potential, Wo^ 4.Tidal radius^ 13.1 pcTotal mass 6.0x105MoTotal number of stars^ 2 45x106Half-mass radius 3 01 pcMean stellar mass 0 245 MoHalf-mass relaxation time 2.59 GyrIMF:Sum of two power-laws with MSIs 0.5 and 5.0 weighted 1:5 by number, 1% bynumber in white dwarfs and 0.37% by number in "black holes".Mass per star Number Mass Width of bin(M0) Fraction Fraction (M0)0.181 0.709 0.529 0.0630.247 0.140 0.142 0.0630.312 0.0486 0.0624 0.0630.398 0.0389 0.0637 0.1170.506 0.0162 0.0337 0.0880.589 0.0102 0.0247 0.0750.665 0.0084 0.0231 0.0770.743 0.0072 0.0222 0.0790.816 0.0052 0.0175 0.0660.869 0.0027 0.0097 0.03800.94 0.0033 0.0128 0.08001.02 0.0033 0.0139 0.08001.10 0.0033 0.0150 0.08002.50 0.0037 0.0300Table A4. Parameters for initial grid of models in Ch. 5.Structure parameters:King model central potential, W o^ 2.,4.,6.Tidal radius (pc)^ 43.1,75.4,86.2Total mass 1.x106M0Total number of stars 6 71x106Mean stellar mass 0 15 MeIMF and stellar evolution:see Table A5.Table A5. Parameters for model M71-AStructure parameters:King model central potential, Wo^ 6.Tidal radius^ 43.1 pcTotal mass 1.x106MoTotal number of stars^ 6 71x106Half-mass radius 6 39 pcMean stellar mass 0 15 MoHalf-mass relaxation time^ 15.8 GyrIMF family (Table B1) AIMF and stellar evolution:Initialmass(M0)NumberFractionMassFractionMass binwidth(M0)0.112 0.551 0.415 0.0260.141 0.245 0.233 0.0330.178 0.110 0.131 0.0410.225 0.0487 0.0734 0.0520.283 0.0217 0.0412 0.0660.356 0.00966 0.0231 0.0830.440 0.00380 0.0112 0.0840.533 0.00275 0.00983 0.1020.645 0.00199 0.00860 0.1240.777 0.00136 0.00710 0.1400.870 0.000287 0.00168 0.0401.00 0.00117 0.00785 0.241.27 0.000779 0.00665 0.301.61 0.000520 0.00563 0.382.04 0.000347 0.00476 0.492.59 0.000232 0.00404 0.623.29 0.000155 0.00342 0.784.17 0.000103 0.00289 0.995.14 5.42e-05 0.00186 0.916.13 4.00e-05 0.00165 1.097.32 2.96e-05 0.00145 1.308.97 2.71e-05 0.00163 2.0611.3 1.84e-05 0.00139 2.5914.2 1.24e-05 0.00119 3.2617.8 8.42e-06 0.00101 4.09Remnant Stellar evolutionmass^epoch(M0) (Gyr)not evolvednot evolvednot evolvednot evolvednot evolvednot evolvednot evolvednot evolvednot evolvednot evolvednot evolved0.581^4.97^11.90.640^1.97^4.970.715^0.862^1.970.810^0.409^0.8620.931^0.205^0.4091.08^0.108^0.2051.28^0.0609 0.1080.00^0.0418 0.06090.00^0.0299 0.04180.00^0.0220 0.02991.40^0.0151 0.02201.40^0.0107 0.01511.40^0.00800 0.01071.40^0.00618 0.00800Table A6. Parameters for model M71-BStructure parameters:King model central potential, Wo^ 6.Tidal radius^ 43.1 pcTotal mass 1.0x106 MoTotal number of stars^ 3 18x106Half-mass radius 6 39 pcMean stellar mass 0 31 MoHalf-mass relaxation time^ 7 91 GyrIMF family (Table B1) IIMF and stellar evolution:Initialmass(M0)NumberFractionMassFractionMass binwidth(M0)0.112 0.353 0.126 0.0260.141 0.222 0.0997 0.0330.178 0.140 0.0793 0.0410.225 0.0882 0.0629 0.0520.283 0.0555 0.0499 0.0660.356 0.0350 0.0396 0.0830.440 0.0190 0.0265 0.0840.533 0.0157 0.0265 0.1020.645 0.0129 0.0265 0.1240.777 0.0101 0.0250 0.1400.870 0.00231 0.00638 0.0401.00 0.0104 0.0330 0.241.27 0.00816 0.0330 0.301.61 0.00643 0.0330 0.382.04 0.00507 0.0330 0.492.59 0.00400 0.0330 0.623.29 0.00315 0.0330 0.784.17 0.00248 0.0330 0.995.14 0.00151 0.0246 0.916.13 0.00126 0.0246 1.097.32 0.00106 0.0246 1.308.97 0.00112 0.0318 2.0611.3 0.000886 0.0318 2.5914.2 0.000705 0.0318 3.2617.8 0.000560 0.0318 4.09Remnant Stellar evolutionmass^epoch(M0) (Gyr)not evolvednot evolvednot evolvednot evolvednot evolvednot evolvednot evolvednot evolvednot evolvednot evolvednot evolved0.581^4.97^11.90.640^1.97^4.970.715^0.862^1.970.810^0.409^0.8620.931^0.205^0.4091.08^0.108^0.2051.28^0.0609 0.1080.00^0.0418 0.06090.00^0.0299 0.04180.00^0.0220 0.02991.40^0.0151 0.02201.40^0.0107 0.01511.40^0.00800 0.01071.40^0.00618 0.00800Table A7. Parameters for model NGC 6397-AStructure parameters:King model central potential, W o^ 6.Tidal radius^ 39.1 pcTotal mass 7.5x105MoTotal number of stars^ 1 63x 106Half-mass radius 5 81 pcMean stellar mass 0 46 MoHalf-mass relaxation time^ 5.81 GyrIMF family (Table B1) JIMF and stellar evolution:Initialmass(M0)NumberFractionMassFractionMass binwidth(M0)0.112 0.275 0.0669 0.0260.141 0.194 0.0595 0.0330.178 0.138 0.0532 0.0410.225 0.0971 0.0473 0.0520.283 0.0687 0.0422 0.0660.356 0.0485 0.0376 0.0830.440 0.0293 0.0280 0.0840.533 0.0246 0.0285 0.1020.645 0.0208 0.0291 0.1240.777 0.0165 0.0279 0.1400.870 0.00381 0.00720 0.0401.00 0.0174 0.0377 0.241.27 0.0140 0.0387 0.301.61 0.0113 0.0396 0.382.04 0.00912 0.0405 0.492.59 0.00738 0.0416 0.623.29 0.00595 0.0425 0.784.17 0.00480 0.0435 0.995.14 0.00297 0.0331 0.916.13 0.00253 0.0337 1.097.32 0.00216 0.0343 1.308.97 0.00233 0.0453 2.0611.3 0.00189 0.0463 2.5914.2 0.00154 0.0475 3.2617.8 0.00125 0.0484 4.09Remnant Stellar evolutionmass^epoch(M0) (Gyr)not evolvednot evolvednot evolvednot evolvednot evolvednot evolvednot evolvednot evolvednot evolvednot evolvednot evolved0.581^4.97^11.90.640^1.97^4.970.715^0.862^1.970.810^0.409^0.8620.931^0.205^0.4091.08^0.108^0.2051.28^0.0609 0.1080.00^0.0418 0.06090.00^0.0299 0.04180.00^0.0220 0.02991.40^0.0151 0.02201.40^0.0107 0.01511.40^0.00800 0.01071.40^0.00618 0.00800Appendix B. Systematics of the Srm—rhd DiagramThe S r, m—rha diagram was introduced in Fig 26 in Ch. 5. It plots the evolution of a normalizedsegregation measure (eq. 3.2) against the half-density radius for the stars dominating the clusterlight. The advantage of such a diagram is that it reduces a set of models into observable parameterswith which the observed values for a particular cluster can be compared. The segregation measureused is dependent on the locations of the fields observed to produce the MFs (the "observationwindow") and so a new diagram is required for each set of observations. To further constrain themodels, only the portion of the evolutionary track which lies in the interval from 10 to 20 Gyr isused. Models which lie near the observed values of S, and rhd are candidates for further investigation.The S,,,h—rhd diagram for M71 shows some systematic properties demonstrating the strong dependenceof model evolution on the IMF.The evolutionary tracks in the S rm—rha diagram for M71 (Fig. 26) appear to fall into severalfamilies and this regularity is made clearer when the tracks belonging to models with differentIMFs are plotted separately. Such a separation is shown in Fig. B1 for the IMFs with massspectral indices listed in Table B 1. Table B1 also contains an abbreviated list of the initial parametersof the models included in Fig. B 1. As expected, there is a general trend that the models with thesteepest IMFs suffer the largest degree of mass segregation and those with flatter IMFs much less.The range of the Sr.m—rhd diagram covered by the models with IMFs A and B is due to differencesin the relaxation time. The initial mass and concentration play a secondary role in determining thedegree of mass segregation a model will suffer.The trend with IMF is confused by the conflicting effects of the two MSIs, xh and x,, involved.This can be seen in Fig. B2 for a set of models with initial mass 5x10 5 Mo , Wo=6 and ri=34 pc.The full evolutionary tracks are shown. The model with IMF G disrupts at 6 Gyr. Table B2 listsfor each model, the IMF family from Table B 1, the fraction of mass lost due to stellar evolution,and the initial mean mass for that family. Also listed for each of the models is the initial half-massrelaxation time, the normalized segregation measure at core collapse and the ratio of the maximumvalue of rhd to its initial value. It might be expected that S,14. ,: 1 at core collapse would correlate with122the initial mean mass or fraction of mass lost through stellar evolution since these are directlyrelated to the slope of the MF. The strongest correlation, an anti-correlation actually, of SrM: lieswith the degree of expansion represented by the maximum value of rhd. The details of the IMF areimportant in determining the fate of a model.The anti-correlation of S;4,„7,' on the maximum value of rkd does not carry over so well when thefull range of initial conditions are considered. The reason for this is the complicating effect of theobservational window used in computing the segregation measure. A comparison of Figs. 26 and28 in ch. 5 makes this quite obvious. The same set of models were used for both S r,„,—rhd diagramsbut the observational windows were quite different. The subtlety of the windowing effect can beseen in the post—core-collapse evolution of S rMm7 1 vs. rhd in Fig. B2. For IMFs A and B, S rNim71 goesthrough a minimum and then begins to increase. An examination of the density profiles for the twomass bins which went into the segregation measure shows that at the position of the 3' MF thedensity of stars in the higher mass bin continues to decrease, while in the core the density ratiobetween the two mass bins is constant. For the other IMFs, the more distant MF would need to beobserved at a larger radius to see the same effect Similar effects underlie the varied behavioursseen in Fig. 28.123Table B 1: IMF familiesFamily IMF Models in Fig. B lax, x,A 3.5 1.7 (110,5), (2,4,6), (20,35,40))B 2.5 1.7 (5, (4,6), (20,35,40))C 3.5 0.7 (5, (4,6), (20,35,40))E 3.0 0.5 (10, 6, 15),^(5, 6, {20,35,40})F 3.5 0.5 (5, 4, 40),^(5, 6, (20,35,40))G 1.5 0.7 (10, 6, (20,25)), (5,6,{25,35})H 1.0 1.0 (5, 6, f 15,20,25,35 ))I 2.0 1.0 (10, 6, (20,25 )),^(10, 8, 25)J 1.5 0.9 (110,7.5,51, 16,81, 120,251), (17.5,51,6,15)'Models are specified in ordered triplets (Mass/105Mo , Wo , rt/(M/105 M0 ) 1 /3).Where a set of values are available for a parameter all possible combinations wereused.Table B2: Results of models shown in Fig. B2Family IMF AMstei evcil(%)<M>(t.--0)(M0)S;:i,h7 1 (cc) trh(t-=0)(Gyr)r hdmaxirh,,,„X1 X.I,A 3.5 1.7 3 0.149 0.67 8.3 1.00B 2.5 1.7 7 0.180 0.53 6.9 1.00C 3.5 0.7 22 0.197 0.25 6.4 1.21F 3.5 0.5 32 0.233 0.19 5.5 1.42H 1.0 1.0 44 0.534 0.16 2.6 1.68J 1.5 0.9 44 0.461 0.16 2.9 1.68E 3.0 0.5 42 0.309 0.15 4.2 1.69G 1.5 0.7 55 0.605 0.12 2.3 2.07B ^C —0.6.44 $010 . 4(13III!F(a)^— —GI^I I^I^I^1^1^1^1(b) -0.20.80.6IHI I HH I IHJ IIHIJHHHI I I I I I I I I I I I I I I I I I_(c)H — —^_ I^—A ^ J ^ _ E —'41' F.0.4C110.20—3"—2^—1^0log rhd (pc)I^ 1 —3^—2^—1^0^1log rhd (pc)Fig. Bl—S diagrams divided by IMF family. S rIvi, m7 1 is defined in eq. (3.3). The evolutionarymtracks are the same as those plotted in Fig. 26 in Ch. 5. (a) IMF B (solid), IMF C (dashed).(b) IMF G (solid), IMF F (dashed). (c) IMF A (solid), IMF E (dashed). (d) IMF J (solid), IMF H(long dash), IMF I (short dash). These three IMFs are very similar.(d) -0.8125Fig. B2— Full Srm —rhd diagram for a set of models with the same initial conditions, but withvarying IMFs. The IMF family from Table B1 for each model is labelled. All the models start withzero segregation measure and logrhd=0.2 and evolve away from this point. The initial movement tolarger r,,d is the expansion due to stellar evolution. The state of each model at 10 Gyr is indicated bythe dot and only the remaining parts of the curves were plotted in Fig. B 1. The model with IMF Gdisrupts before 10 Gyr.126
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Fokker-Planck models and globular cluster evolution
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Fokker-Planck models and globular cluster evolution Drukier, Gordon A. 1993
pdf
Page Metadata
Item Metadata
Title | Fokker-Planck models and globular cluster evolution |
Creator |
Drukier, Gordon A. |
Date Issued | 1993 |
Description | Numerical models for globular cluster evolution using the orbit-averaged Fokker-Planck equation are compared with observations of the globular clusters M71 and NGC 6397. The first set of models studied includes a mass spectrum, a tidal boundary and allows for a central energy source from the formation, hardening and destruction of binaries formed in three-body interactions. These models are compared to star-count mass functions and surface density profiles for M71. It is possible to reproduce the tidal boundary and degree of mass segregation using these models but the central density profiles of the fitted models are much steeper than is observed. On the basis of the degree of mass segregation and the short relaxation time, M71 is dynamically a highly evolved object, but the surface density profile contradicts this. Several modifications of the model are considered to resolve this difficulty, and an additional energy source, which does not require extreme densities to operate, appears to be required. New observations of the cluster NGC 6397 are presented. This cluster is presently at a similar position in the galaxy as M71, and has a similar mass and mass function, but is much more centrally concentrated. NGC 6397 is a low metallicity cluster with halo kinematics, in contrast to the metal-rich disk cluster M71. The observations have been reduced to give a surface density profile, and mass functions at two distances from the cluster centre. The surface density profile is fit by a power-law with slope —0.8 in the central arc minute. The mass functions give some evidence for mass segregation, but the low signal-to-noise of the outer mass function limits the significance of this observation. A second set of models, augmented to include the effects of mass loss due to stellar evolution, are used to find matches for the observations of M71 and NGC 6397. A diagram showing the evolution with time of the amount of mass segregation and the half-density radius is shown to be useful in identifying good candidate models. Stellar evolution appears not to solve the problem ofM71 as no satisfactorily matching model has been found. Possible solutions to this problem are discussed. The observations of NGC 6397 can be reproduced by the models and the best such model is described. This model has an initial mass function with a Salpeter-like slope and is undergoing core-collapse. This demonstrates that a numerical model based on the orbit-averaged Fokker-Planck equation, and including an evolving mass spectrum, a tidal boundary and the binary energy source, is a useful model for the dynamical evolution of some globular clusters. |
Extent | 6456692 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-09-05 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0085507 |
URI | http://hdl.handle.net/2429/1658 |
Degree |
Doctor of Philosophy - PhD |
Program |
Astronomy |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1993-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1993_spring_phd_drukier_gordon.pdf [ 6.16MB ]
- Metadata
- JSON: 831-1.0085507.json
- JSON-LD: 831-1.0085507-ld.json
- RDF/XML (Pretty): 831-1.0085507-rdf.xml
- RDF/JSON: 831-1.0085507-rdf.json
- Turtle: 831-1.0085507-turtle.txt
- N-Triples: 831-1.0085507-rdf-ntriples.txt
- Original Record: 831-1.0085507-source.json
- Full Text
- 831-1.0085507-fulltext.txt
- Citation
- 831-1.0085507.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0085507/manifest