B"H FOKKER-PLANCK MODELS AND GLOBULAR CLUSTER EVOLUTION by GORDON ALAN DRUKIER B.Sc., The University of Toronto, 1985 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPY in THE FACULTY OF GRADUATE STUDIES (Department of Geophysics & Astronomy) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA October 1992 © Gordon Alan Drukier, 1992 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature) Department of 6eof ki543^As 1 rovlol - The University of British Columbia Vancouver, Canada Date ^Nen/ 72 DE-6 (2/88) - Abstract 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 of M71 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. ll Table of Contents Abstract ^ ii List of Tables ^ v List of Figures ^ vi List of Abbreviations ^ viii List of Symbols ^ ix Acknowledgements ^ x Chapter 1. Introduction ^ 1 Section a. Recent Results in Fokker-Planck Modelling ^ 5 Section b. Comparing Fokker-Planck Models with Observations ^ 11 Chapter 2. The Models ^ 15 Section a. Basic Formalism ^ 15 Section b. Multiple Masses ^ 16 Section c. Binary Heating ^ 19 Section d. Tidal Boundary ^ 24 Section e. Stellar Evolution ^ 29 Section f. Initial Conditions ^ 37 Section g. Caveats ^ 38 Chapter 3. The Problem of M71 ^ 39 Section a. Outline of the Problem ^ 39 Section b. Approaches to a Solution ^ 48 Subsection I. Gravothermal Oscillations ^ hi 48 Subsection II. Extreme Degenerates ^ 49 Subsection III. Modifying the Heating Rate ^ 51 Chapter 4. NGC 6397 ^ 56 Section a. Observations and Reduction ^ 58 Section b. Surface Density Profile ^ 64 Section c. Mass Functions ^ 75 Chapter 5. Models Including Stellar Evolution ^ 91 Section a. M71 ^ 94 Section b. NGC 6397 ^ 99 Chapter 6. Conclusions. ^ 105 References ^ 110 Appendix A. Initial Model Parameters ^ 115 Appendix B. Systematics of the S r „,--rhd Diagram ^ 122 iv List of Tables 1. Core collapse times for two component models ^ 19 2. Results of the sample stellar evolution grid ^ 36 3. Field locations for NGC 6397 observations ^ 59 4. Surface density profile for 1<14 ^ 68 5. Surface density profile for R15.5 ^ 72 6 Luminosity function for field du Pont:bf ^ 76 7. Luminosity function for field du Pont:if ^ 78 8. Luminosity function for field du Pont:n ^ 79 9. Luminosity function for field du Pont:w ^ 80 10. Background subtracted luminosity functions for du Pont:if and du Pont:out ^ 82 11. Mass functions for du Pont:if and du Pont:out ^ 86 12. Observed mass spectral indices ^ 88 13. Mass spectral indices for model NGC 6397—A ^ 104 Al. Parameters for the model CW1 ^ 115 A2. Parameters for the "standard" model of M71 ^ 116 A3. Parameters for the "black hole" model of M71 ^ 117 A4. Parameters for the initial set of models in Ch. 5 ^ 118 A5. Parameters for model M71—A ^ 119 A6. Parameters for model M71—B ^ 120 A7. Parameters for model NGC 6397—A ^ 121 Bl. IMF families of Fig. B1 ^ 124 B2. Results of models shown in Fig. B2 ^ 124 List of Figures 1. Projected radial profiles for the mean mass in model CW1 ^ 33 2. Projected surface density profiles for model CW1 ^ 34 3. Evolution of the number of stars by species in model CW1 ^ 35 4. Observations of M71 ^ 40 5. Observed segregation measure for M71 ^ 44 6. Results of the "standard" model ^ 46 7. Segregation measure in the "standard" model ^ 47 8. Results of the "black hole" model ^ 50 9. Results of the "extra heating" model ^ 52 10. Kinematics of M71 and the "extra heating" model ^ 53 11. Segregation measures for the "black hole" and "extra heating" models ^ 54 12.Map of the fields observed in NGC 6397 ^ 60 13. Map showing the annuli used in constructing the SDP ^ 65 14. Surface density profile for k 14 for the annuli in Fig. ^ 66 15. Surface density profile with merged sections ^ 69 16. Central surface density profiles from the four quadrants of Swope:c ^ 70 17. Combined /<14 and R15.5 surface density profiles ^ 73 18.Power-law fits to the surface density profile ^ 74 19. Comparison of the blank field star counts with those of Fahlman et al. (1989)^ 81 20. Luminosity function for du Pont:if ^ 83 21. Mean luminosity function for du Pont:w and du Pont:n ^ 84 22. Composite of all deep luminosity functions for NGC 6397 ^ 85 23. Mass functions for the three NGC 6397 luminosity functions ^ 87 vi 24. Segregation measures between the three mass functions ^ 90 25. Comparison with observations for model M71—A ^ 95 26. S r m rid diagram for M71 observations ^ 97 27. Comparison with observations for model M71—B ^ 98 28. S ,.„—rhd diagram for NGC 6397 observations ^ 100 29. Comparison with observations for model NGC 6397—A ^ 101 30. Velocity dispersion profile for model NGC 6397—A ^ 103 Bl. S r ,„,—rhd diagram for M71 divided by IMF family ^ 125 B2. Full S r„—rhd curves for a set of models differing only in IMF ^ 126 vii List of Abbreviations There are certain phrases and references used often enough in this thesis to receive the dubious honour of being replaced by abbreviations. Although all the abbreviations are defined where they first occur, the reader might find it useful to have the definitions collected together for easy reference. Jargon CCD^charge coupled device GTO^gravothermal oscillation IMF^initial mass function LF^ luminosity function LMC^Large Magellanic Cloud MF^ mass function MSI^mass spectral index PCC^post-core-collapse SDP^surface density profile References CW^ Chernoff & Weinberg 1990 DFR^ Drukier, Fahlman, & Richer 1992 FRST^Fahlman, Richer, Searle, & Thompson 1989 GCLM^Grabhorn, Cohn, Lugger, & Murphy 1992 LFR^Lee, Fahlman, & Richer 1991 MCH^Murphy, Cohn, & Hut 1990 SOC^Statler, Ostriker, & Cohn 1985 viii List of Symbols E^Energy E,^Energy at tidal boundary^ §2d f^Distribution function §2b In A^Coulomb logarithm^ m^Stellar mass m i^Mean mass for bin i Initial mass for stars in bin j^ §2e Final mass for stars in bin j^ §2e M^Total mass Mr^Total initial mass for a model, scale mass ^§2b N^Total number of stars N,^Total number of stars in bin i N,„^Number of stars per unit mass (mass function) r^Radial position r0^Scale radius^ §2b Eq. (2.12) Core radius^ §4b reff^Effective radius^ rh^Half-mass radius^ §2c §5a rhd^Half-density radius^ rt^Tidal radius^ §2d S r^Radial segregation measure^Eq. (3.1) S rth^Mass-range normalized segregation measure ^Eq. (3.2) Sr,Mm71^S r /n for M71 observations^ Eq. (3.3) Time t tth^Half-mass relaxation time ^ Eq. (2.11) v^Velocity, velocity dispersion vo^Scale velocity^ Eq. (2.5) Wo^Dimensionless central density in King models^§2e x^Mass spectral index^ Eq. (1.3) x h^Mass spectral index for high mass stars^Ch. 5 x,^Mass spectral index for low mass stars^Ch. 5 p^Mass density Pi^Mass density for mass bin i^ Eq. (2.4) §2c Pc^Central density^ a^projected number density an,^projected mass density Acknowledgements Like a globular cluster, the work of preparing a Ph. D. does not exist in isolation, and one-half page is permitted by university regulations for acknowledging this. Thanks are due to my supervisor Dr. Greg Fahlman for many discussions, and for much advice and encouragement. Without his many useful comments on the manuscript, there are sections of this thesis which would have been incomprehensible. Together with Dr. Ian Thompson and the people responsible for the Observatories of 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 well as development versions of some of his other software. My fellow grad students at UBC deserve a vote of thanks for serving as sounding boards when I had something to talk about and for putting up with my muttering when I didn't. Special thanks go to old-timers Claudia Mendes and Dave Woods and other sounding boards James Brewer, Phil Hodder, 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 for letting me live 3300 km away and visit home far too infrequently. My adopted family, in the form of 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 for giving me the strength to do this work and for the remarkable miracle that the human mind is actually capable of making sense of the physical world he created. .17 7 17 inn 1 7 n17.717M1 11= , - , , ^v - incn crOT`r The heavens declare the glory of G-d; the expanse of the sky proclaims His handiwork. Psalms 19:2 Chapter 1. Introduction A spiral galaxy such as our own is composed of many parts. The spinning, flattened disk shows off spiral arms marked by the brief, but brilliant, light of the young, massive stars. Most of the mass of the disk, however, is made of their cooler, less massive, and longer-lived cousins together with the gas from which the stars are built. Distinct from the disk is the spheroid; an older population of stars, less endowed with heavy elements. The spheroid, or halo, stars define a roughly spherical volume enclosing the disk. The space velocities of the halo stars are highly random and show little of the orderly rotation of the disk. Along with the general field of single and binary stars, a galaxy contains distinct subsystems arising from the preference of stars to form in groups. These groups range from small, multiple star systems through to unbound associations and bound clusters. The disk is home to open clusters; somewhat loose groupings of hundreds of stars which, subject to the stresses of an environment 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 globular clusters. The stars which make up a globular cluster formed at the same time from a chemically homogeneous cloud of gas. From the fitting of stellar evolution isochrones to the colour-magnitude diagrams of globular 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 is likely 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 potentially have much to tell us about galaxy formation. In themselves, the globular clusters, which contain 10 5 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 globular clusters. Much effort from a number of directions has been put into understanding globular cluster evolution, but the most fruitful approach has been the numerical models based upon the Fokker-Planck equation. 1 There comes a time for every theoretical model when it must come up against the phenomena it is supposed to explain. In principle, the evolution of a globular cluster is simply an extension to large N of the N-body problem of classical mechanics. Since we are not pursuing exact predictions for the position and velocity of each and every star observed in a globular cluster, a statistical description is sufficient. In a certain sense, even a direct N-body integration of a large N system may be considered statistical in that the overall behaviour of the model will depend more on the average properties of stars with similar orbits than on the trajectories of the individual bodies. A more direct statistical approach is to describe the structure of the cluster in terms of a phase-space distribution functionf(r,v,t) for which the time evolution is governed by the Boltzmann equation, at df + v • Vf — V(I) -aI-= (— at^av^dt ^(1.1) where (df/dt) enc is the change in the distribution function due to encounters between stars. For a self-gravitating system the potential (I)(r,t) satisfies Poisson's equation V 24:0 = 47cGi fd 3 v.^ (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, ^Vc1)^= + v • Vf —^ 0 .^(1.3) av at 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 cluster is 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 3 AwAt be the probability that the star is scattered by two-body encounters into a volume of phase space d3 Aw around w +Aw during the time interval At . Then the encounter term in eq. (1.1) can be written as the sum of the rates for stars to be scattered into and out of a given position in phase space, ie. I) 6 =^— Aw, Aw)f (w — AW)d 3 AW — f (W)f P (W , Aw)d 3 Aw . (1.4) ( dt enc It can be shown (eg. Binney & Tremaine, 1987, ch. 8; Spitzer, 1987, ch. 2) that under the 2 ^ ^ conditions 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 weak encounters rather than a few strong interactions. Goodman (1983) produced a numerical model using an orbit-averaged equation which included strong encounters and confirmed that, with the exception of an increase in the rate of evaporation of stars from the model, the model's evolution was basically unchanged. In the case of weak encounters, Awl is small and we can expand 1 (w – Aw,Aw)f (w – Aw) in a Taylor series in Aw : 1 1 a 6^ tp(w — Aw, Aw)f (w – Aw) =^Aw)f (w)^Aw —[t 1 f(w , Aw) f (w)1 - I o2 1 6^ Aw.Aw. ^ [tlf(w, F— 2^awiawi Aw)f (w)] + 0(Aw 3 ) . (1.5) The Fokker-Planck approximation lies in ending the expansion at second order in Aw. Substituting eq. (1.5) (to second order in Aw) in eq. (1.4) and doing the integrations over d 3 Aw gives the Fokker-Planck form of the encounter term ( 1 ) FP dt a2 [f (w)D(Aw i Aw^(1.6) j _ i uwi uwi ^{f (w)D(Aw i)1+^ enc^i=i dwi ^2 D(Awi ) and D(Aw i Awi ) are referred to as diffusion coefficients and are defined as D(Awi^Aw)d3Aw^(1.7a) and D(Aw1Aw1) Aw i Awj W(w, Aw)d 3 Aw^(1.7b) The method of Rosenbluth, MacDonald, & Judd (1957) can be used to evaluate the diffusion coefficients given a set of coordinates. Using eq. (1.6) for the encounter term of the Boltzmann equation, eq. (1.1), gives the Fokker-Planck equation. Self-consistency requires that the Fokker-Planck equation be solved together with the Poisson equation, eq. (1.2). The coupling of these two equations makes the problem highly non-linear and complicated to solve numerically. Even this is an idealized situation, since a globular cluster is not an 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 tidal 3 stresses induced by both the general gravitational field of the galaxy and by massive constituents of the galaxy such as giant molecular clouds. Passages through a galactic disk, or close to a massive bulge, will heighten these effects (Ostriker, Spitzer & Chevalier 1972). The individual bodies in this N-body system are not ideal point particles of equal mass, but are stars with finite radii and a range of masses. The range in mass leads to mass segregation. Dynamical relaxation leads towards equipartition of energy between stars with differing masses. The more massive stars will have lower velocities than the less massive stars, and so cannot travel as far from the centre of the cluster. Thus the massive stars are concentrated in the core of the cluster and leave the halo to be dominated by the low mass stars. In the presence of an external tidal field, this leads to the low mass stars being favoured for tidal stripping. While for relaxation the strong encounters are not as important as the sum of all the weak encounters, the effects of strong encounters involving binary stars will play an important role in the dynamical evolution of a cluster. If these encounters are close enough, then the finite radius of a real star enters into the problem. The tidal forces between two passing stars will induce motions in the stars' atmospheres. The resulting dissipation of energy may result in the two stars becoming bound, or even merging. For even closer approaches the stars will collide. In the absence of dynamically formed binaries, there is every reason to expect that some of the original stars are multiple. The presence of these primordial binaries must be taken into account. Finally, stellar masses do not remain constant, but change as the star evolves. Adding in each of these complications to the Fokker-Planck formalism increases the complexity and time requirements of the numerical model. To begin with, the models were simple. Additional physical processes were added as they were thought necessary or were numerically practical. The models in use here follow directly from those of Cohn (1980) and use his technique of direct integration of the Fokker-Planck equation in its orbit-averaged form (see §2a). Other approaches to solving the Fokker-Planck equation utilize a Monte Carlo approach, either with the orbit-averaged equation (116non 1971; Shapiro & Marchant 1978) or directly on the equations of motion with Av and (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, 4 presented an elaborate model which included stellar evolution and several binary formation mechanisms. His model served as inspiration for much of the work that followed. A hybrid N-body/Fokker-Planck approach was used by McMillan & Lightman (1984a; 1984b; McMillan 1986). One problem with the Monte Carlo techniques is that the resulting models suffer from high statistical noise, which make them difficult to compare with observations. N-body studies would be 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 directed to Spitzer (1987), as well as Elson, Hut, & Inagaki (1987), for a summary of the results of the theoretical work up to their epoch. In what follows I will discuss some of the relevant work leading up to this thesis, generally focussing on Fokker-Planck models since these are the most developed and amenable to direct comparison with observations. Section a. Recent Results in Fokker-Planck Modelling By 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 was required to reverse the core collapse and bring the cluster into an expanding phase if not an equilibrium. The favoured mechanism, binary stars, dates back to H6non's (1961) pioneering work. But there were other features to consider as well, including: the actual mechanisms producing the binaries, stellar evolution, tidal effects, and gravothermal oscillations (GTOs). Many of these processes 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 cluster will interact with the binary. If the orbital velocity of the binary is higher than that of an average field star (ie. a "hard" binary) then the result of the interaction, on average, will be a hardening of the binary and an increase in the velocity of the field star. In some cases the field star will be ejected from the cluster. "Soft" binaries on the other hand, tend to become even more loosely bound and do not last long in the face of these perturbations. For sufficiently hard binaries, there is a good chance that the binary will receive a large enough recoil velocity in a strong encounter to be 5 ejected from the cluster entirely. The increased velocity acquired by the passing field stars is distributed to other cluster members through two-body relaxation. On average the velocity of the stars in the core will increase. This is the heating effect. The ejection of binaries from the cluster also serves to heat the cluster but through a more indirect mechanism; the removal of mass serves to decrease the binding energy of the cluster leading to expansion. For the virial theorem to be satisfied, the overall kinetic energy of the cluster must increase. The net heat input of mass AM escaping from the core is approximately (0 0 — I v02 )AM, where (1) 0 and v 0 are the central values of the potential and velocity. Mass lost from the cluster due to stellar evolution, or for any other reason, will heat the cluster in the same way. Statler, Ostriker & Cohn (1987, hereafter SOC) used the hard binaries formed by the tidal interactions 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 directly and via the ejection of the binary and the interacting field star from the cluster. Merging of the stars into a single object was not allowed for. Much care was taken to account for all the possible interactions between the various kinds of stars (single and binary), resulting in a numerical scheme of some complexity. As expected, the newly formed binaries accelerate core collapse, via the mass-stratification instability (Spitzer 1969) and the cooling effect of tidal binary formation, but the energy extracted from them subsequently is able to halt the core collapse and to power a slow re-expansion. The half-mass radius and the velocity dispersion were found to scale as rec t 213 and vrh2 oct -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 are gravothermal oscillations is unclear due to the short time span covered. SOC, as is usually the case with groundbreaking work, raises more questions than it answers. A series of papers in subsequent years elaborated on some of these questions. Lee (1987a) added a degenerate component to the model of SOC. Tidal binaries were allowed to form between the normal and degenerate stars and "three-body binaries" amongst the degenerate stars. In the latter mechanism, binaries are formed during the close encounter of three stars. The third star carries off the excess energy required to bind the other two. The population of three-body 6 binaries was followed explicitly in this paper and the heating coefficients were integrated from the three-body-binary distribution function. Three-body binaries will only be important in high density environments such as the collapsing core. Mass segregation ensures that the degenerates, being more massive than the normal stars, will dominate the core, so the environment there is ideal for three-body binary formation to occur. Since the tidal capture process will not work with degenerate stars, the binary heating will not necessarily be dominated by tidally formed binaries as in the models of SOC. However, in the models presented by Lee there are many more tidal binaries than three-body binaries at the beginning of re-expansion, and the heating is dominated by them. The half-mass radius and velocity dispersion during the re-expansion follow about the same power law dependences as in SOC. The slope of the surface density or brightness distribution for the evolved models is demonstrated to be sensitive to the ratio of the masses of the degenerates and normal stars, as is expected from equilibrium theory. For a collapsed core, which is ideally an isothermal sphere, the projected surface density profile (SDP) would have a logarithmic slope of —1. In the context of globular clusters, a slope shallower than —1 indicates the presence of a large population of more massive, non-luminous stars which, by mass segregation, are more centrally concentrated than the stars producing most of the light. Lee (1987b) followed along a slightly different line from SOC by taking the extreme case that all tidal interactions lead to mergers. The resulting stars are then allowed to evolve off the main sequence on appropriate timescales, ejecting all of their mass from the cluster. The decrease in the cluster's binding energy associated with the mass shed by the evolving stars is sufficient to reverse core collapse. Three-body binaries were also allowed to form but, unlike Lee (1987a), they were not followed explicitly. Instead, the heating due to three-body binaries is treated statistically by including a heating term in the Fokker-Planck equation. A similar procedure was employed here and will be discussed in Chapter 2. In the case of Lee (1987b) only stars with equal masses were allowed to form binaries. The binary subsequently interacted only with stars with the same mass as its components. Again, the post-collapse scaling of the half-mass radius and the velocity dispersion follow the expected power-laws. In contrast to Lee (1987a), the most massive stars are also the brightest, so the higher concentration of the more massive stars leads to the light falling off much more steeply with radius than does the mass. 7 Lee & Ostriker (1987) followed up on SOC by adding a tidal boundary. As will be discussed below in §2d, the tidal radius in the method of Lee & Ostriker is defined as the radius at which the mean density enclosed is equal to a constant, previously specified density referred to as the "tidal density". The gravitational potential at this radius defines the tidal boundary in energy space. The escape 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 the re-expansion phase on a timescale proportional to the product of the total number of stars and the orbital period of the cluster about the galactic centre. The tidally-truncated Fokker-Planck models were shown to be similar to King (1966) models, although the Fokker-Planck models had larger tidal radii for a given central surface brightness. The emphasis on dynamically formed binaries should not take away from the possible existence of binaries in the original stellar population (Hut et al. 1992). Sufficiently hard primordial binaries would serve as an energy source even before the cluster becomes dense enough to form tidal-capture or three-body binaries. Goodman & Hut (1989) examined the likely effects of primordial binaries on 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 the collapse. Goodman & Hut show that primordial binaries can lead to cluster with a fairly large core radius and a core dominated by binaries. The relative importance of primordial binaries will depend on their abundance and destruction rate. Gao et al. (1991) examined this question in more detail. The models they constructed consisted of two components; one for single stars and one for the binaries. The total mass of each binary was twice that of a single star. The distribution of internal binding energies and the exchange of energy between the centre of mass motion and orbital motion of the binaries were treated in detail. The models were followed well into the post-collapse phase. The main effect of the primordial binaries is to leave the collapsed core with a fairly large radius. The ratio of core to half-mass radii is around 0.01 to 0.04. The larger core radius suppresses the production of new binaries. One distinction between models with primordial binaries and those where manufactured binaries act as the heat source is that the former have cores consisting predominantly of binaries. This should be a testable distinction between the two scenarios. 8 Although models with more than one mass of star have been discussed above, the additional mass classes were not an end in themselves but served as part of the detailed treatment of the reheating 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 in general 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-law dN oc m -(X1-1)din (1.3) When the mass spectral index (MSI), x, is defined in this way, then for the Salpeter (1955) mass function x=1.35. The initial mass function (IMF) is the global mass function before significant evolution has taken place. Cohn, Hut & Wise (1989) took up the question posed in SOC of the post-collapse instability of Fokker-Planck models. Their single mass species models were reheated by the three-body mechanism following a formalism similar to that of Lee (1987b). Gravothermal oscillations had been observed in gas sphere models, (Sugimoto & Bettwieser 1983; Bettwieser & Sugimoto 1984) but were not seen 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 gravothermal instability. The time-step selection technique which had been used chose steps in units larger than 104 central relaxation times in the epoch following core collapse. These long time steps served to suppress 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 an evolved mass function. The re-expansion was powered by three-body binaries and a new expression was 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 confirmed the previous findings that the presence of a range of masses hastened core collapse. MCH investigated the post—core-collapse (PCC) behaviour of models using small time steps, and found that multiple 9 mass species served to suppress the oscillations. This appears to be a result of the dependence of GTOs on the number of stars (Goodman 1987). Since post-collapse cores will be dominated by heavy remnants, it is the number of these present, rather than the total number of stars, which will determine whether oscillations will take place. This conclusion is supported by the observation that in models with power-law IMFs, those with smaller IMF MSIs (ie. those with a higher proportion of high mass stars) are more unstable to GTOs than those with a larger MSI. In general, it was found that the post-collapse surface brightness profiles do not differ greatly from those of core-collapse itself beyond the small radius affected by the GTOs. Lee, Fahlman & Richer (1991, hereafter LFR) produced similar models to MCH but also included a tidal limit like that of Lee & Ostriker (1987). Their focus was more on the behaviour of the mass function during the long-term dynamical evolution of the cluster. Within the half-mass radius, 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 even invert. Mass segregation will complicate measuring the overall MF of a cluster. It was also found that as a cluster approaches disruption by the Galaxy's tidal field, the half-mass relaxation time drops dramatically. LFR commented on these points with respect to the observations of M71 of Richer & Fahlman (1989). These considerations led to the more extensive comparisons made by Drukier, Fahlman & Richer (1992) which will be discussed in Chapter 3. Chernoff & Weinberg (1990, hereafter CW) produced pre-collapse models which included stellar evolution along with a tidal boundary and a mass spectrum. As their models are an independently coded implementation of Cohn's (1980) formalism, they provide a useful check on the results supplied by all the other work discussed here. Their "quality control" models gave results consistent with those Cohn (1980) and Inagaki & Wiyanto (1984). As expected from Applegate (1986), the heating 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 those of King models) or with flat enough mass functions, the coupling of stellar evolution with a tidal boundary was enough to disrupt the cluster on short order. For models which would survive to the present they discuss several observable quantities. Since in the advanced stages of evolution the 10 core is dominated by remnants more massive than the stars providing the light, the surface brightness profile may not give a detailed report of the cluster's dynamical state. On the other hand, mass segregation will lead to large effects on the mass function observed with radius. Section b. Comparing Fokker-Planck Models with Observations In all the work described until now the main emphasis has been on the theoretical development of the Fokker-Planck formalism. Since the aim of the theoretical work has been to see the effects of various physical processes on the models, less discussion has taken place with respect to the observed properties of globular clusters. Murphy & Cohn (1988) built pre—core-collapse, multi-mass models. Since the focus of Murphy & Cohn was to make realistic globular cluster models for the purpose of comparing them to observations, projected surface densities and velocity dispersions were presented and discussed. One of the models was compared to the U-band surface brightness profile of M15 by Lugger et al. (1987) and a good match was obtained. The presence of the heavy remnants was adequate to flatten the slope of the surface brightness profile from —1 expected for an isothermal sphere to the observed value of around —0.8. MCH, as mentioned above, extended this study into the post—core-collapse phase and suggested that that cluster is currently undergoing a collapse rather than a re-expansion, although at that time the data did not go far enough into the core of M15 for definitive conclusions to be drawn. LFR discussed their results with respect to M71, 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 the collapsed-core clusters M15 and NGC 6624 were compared with Fokker-Planck models. These models 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 five mass species, including heavy remnants. The important feature of these clusters is that they have very high concentrations and so are well suited to probing the core-collapse, or post—core-collapse phase of cluster evolution. M15 is the prototypical collapsed-core cluster, its core possibly being resolved with HST (Lauer a al. 1991, but see Grabhorn et al. 1992b). The main result of GCLM 11 is that both clusters can be fit by PCC models. The surface brightness profiles are best fit when the model is in a state of maximal expansion during a gravothermal oscillation. Their preferred models have an IMF with x=0.9. The core in both clusters is dominated by non-luminous remnants with masses up to 1.4 Mo . For M15, Fahlman, Richer & VandenBerg (1985) measured the mass spectral 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. In addition, the gravitational potential in their M15 model is sufficient to explain the observed period decrease in the millisecond pulsar PSR 2127+11A (Wolszczan et al. 1989). GCLM serves as a useful sign post on the road to understanding the relationship between the models and the state of real globular clusters. The models are clearly capable of reproducing the surface brightness profiles and velocity dispersions of collapsed-core clusters. With some success achieved for clusters with central cusps, it becomes all the more necessary to test the models against clusters with a broader range 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 question became even more pressing when the accelerating effect of a mass spectrum was considered. The models 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 survey of globular cluster cores (Djorgovski & King 1986, Chernoff & Djorgovski 1989) has now uncovered many cases of clusters showing this sign. In this survey the clusters are divided into those with surface brightness profiles fitting King (1966) models, and those which do not, but instead 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 about 0.4 Mo . Velocity measurements have also been made for many clusters. Like much of the surface brightness and mass function data, the velocities have been used as an adjunct to fitting King 12 models. 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) [NGC 6397], 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 globular cluster dynamics Clearly, the stellar population in a cluster is affected by the dynamical evolution of the cluster, especially by the high densities accompanying core collapse. The motivation for the work to be described here comes from consideration of the observational and theoretical state of globular cluster studies. The Fokker-Planck models have come a long way from their initial state of idealization. The addition of various binary mechanisms, the tidal boundary and stellar evolution have brought the modelling much closer to including all the physics presumed to be acting in real globular clusters. Still, simplifications remain in all these elaborations. What the modelling has provided are hints as to what observable properties are most sensitive to dynamical evolution. The surface density profile (SDP) or surface brightness profile will give information on the status of the cluster with respect to core collapse. The degree to which such profiles will reflect core-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 the luminous 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. Velocity dispersions give more direct information on the remnant content and the degree to which core collapse has progressed. The radial variation in the mass function gives information on the degree of mass segregation. Since, as will be shown below, the course of the evolution is dependent on the IMF, observed MFs serve to constrain the range of plausible IMFs and hence the range of possible histories for the cluster. It is with these ideas in mind that this thesis has been undertaken. There were two phases to the work reported on here. In the first phase, Fokker-Planck models which included (i) three-body 13 binaries, (ii) a mass spectrum, and (iii) a tidal limit, were constructed in an attempt to match observations of M71. The observations involved consisted of SDPs for two different ranges of stellar mass, and MFs for two regions of the cluster. These models, however, did not prove adequate to match the M71 observations. The problem with the comparison lies in the contradiction between 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 resolve this contradiction several types of slightly modified models were produced. The principle result of the 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 the same time, mass segregation must not be arrested. Two options for this energy source are primordial binaries and the effects of mass loss during stellar evolution. Phase two followed on from this conclusion. As stellar evolution is more straightforward to implement, and more within the statistical philosophy of Fokker-Planck models, it was added to the code. I also felt it useful to investigate another cluster by way of comparison and so analyzed an extensive set of observations of NGC 6397 to produce a SDP and two new MFs. This nearby cluster has a central cusp in its surface brightness profile. A new series of models which included stellar 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 are described in Chapter 2. The results of the first phase comparisons with M71 are discussed in Chapter 3. Chapter 4 presents the new observations for NGC 6397. In Chapter 5, the models including 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 the observations have for future modelling. In addition, Appendix A provides lists of parameters for all the models discussed in this thesis and Appendix B discusses further the S ,.„--rhd diagram introduced in Chapter 5. 14 Chapter 2. The Models In this chapter I will discuss the numerical techniques used in implementing the Fokker-Planck equation for globular cluster modelling. The first two sections will describe the basic formalism introduced 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 binary heating mechanism, the tidal boundary, and stellar evolution. Along the way, comparisons will be made with the results of previous studies using similar techniques. Following this will be a section dealing with the initial conditions assumed for the models and a section recapitulating the assumptions made and the limitations of these models. Most of the actual computer code is not original with me but was written by H. Cohn, H. M. Lee and possibly others who I am unaware of. Besides the additions I wrote, I have gone through the code to check for consistency and to improve its efficiency. I wrote all of the routines required to analyse the output of the models . Section a. Basic Formalism The code used in this work is based on that of Cohn (1980). The method for numerically integrating the isotropic Fokker-Planck equation is discussed in that paper. Three major assumptions were made in arriving at this technique. First, it is assumed that the timescale required for significant changes in the energy of the stars, the relaxation time, is much longer than the dynamical time. It is then advantageous to consider separately the changes in the distribution function on the two time scales. Over the orbital period of a star in the cluster, two-body relaxation is unimportant and the distribution function is a solution of the collisionless Boltzmann equation (eq. 1.3). By Jeans's Theorem the distribution function is a function of the isolating integrals of motion, /i and the long term evolution then acts as a perturbation on the distribution function. 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 to stars with integrals of motion The second assumption made is that of spherical symmetry. In this case, f=fiE At), ie. the 15 distribution 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 this code for long integrations for two reasons. One was the computational expense of the twodimensional Fokker-Planck equation. The second was the problem that the total energy was poorly conserved. If, as a third simplifying assumption, isotropy is imposed on the velocity distribution, then the distribution function is reduced to a function of the energy only and both problems can be alleviated. One fewer dimension will obviously reduce computational expense and in the one-dimensional form the differencing scheme of Chang & Cooper (1970), which is known to conserve the energy well, is available. A numerical solution of these equations is obtained by first advancing the distribution functions in time in accordance with the Fokker-Planck equation, eq. (2.6). At this point the potential is inconsistent with the density distribution. The second stage is to re-solve for the potential using Poisson's equation and the new distribution functions. Since the changes in the potential will also change the energies, the potential solution is constrained by requiring that the distribution functions remain fixed functions of q (eq. 2.2c), an adiabatic invariant. The solution of Poisson's equation is much more expensive in computer time than the Fokker-Planck step, so four of the latter were performed for each Poisson equation step. The details of the method, for the equal-mass case, are described in Cohn (1979; 1980). In common with the formalism as presented by Cohn (1979; 1980) the energy convention in what follows has the potential energy positive and falling to zero as r approaches infinity. To denote 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 Masses One of the first changes required in order to make detailed comparisons with observations is to provide for stars of different masses. This is easily achieved by having the distribution function as a function of mass as well as of energy. Although it is possible to write the distribution 16 function as a continuous function of the mass, from the beginning I will use a set of distribution functions, 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 the 1-dimensional, orbit-averaged, Fokker-Planck equation at a constant potential as (Murphy & Cohn 1988): af a ---L —^rivf:DE ± DEE adl n P at 4 2 - a'[ - (2.1) - where the diffusion coefficients DE and D EE are given by DE (E) —64n 4 G 2 in A^i rE (°) pfi dE , DEE(E)= 64n 4 G 2 1nAEm?[q f0E fi dE + L•(°) qfi dE]. (2.2a) The function p is the isotropic average of the period of stars with energy E, p(E) = 4 r 1(E)r2 vdr,,^ (2.2b) where the velocity v = [24)(r) — 24 12 . The integral of p, q is a similar average of the radial action: 4 r 1 (E) 2 3 R'( 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 the total 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 in order to get 4)(r) and its inverse 4) -1 (E). Poisson's equation is used in its integral form, 4)(r) = 4ir G r r 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 integrals p (r) = 2 512 TCMi r (r) [00 —^fi dE^(2.4) The boundary conditions on Poisson's equation are that the potential goes to zero in the limit 17 ^ r--->. 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 are chosen 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— )1/2 GM,^ r° (2.5) For this choice of dimensions, the natural time unit is ro /v o but C,(ro /v 0 ) is used instead. The dimensionless constant C, will be chosen below in order to eliminate certain constants from the dimensionless form of the Fokker-Planck equation. A further complication is that the masses per object of the individual mass species, m ; , are more easily normalized by C„,M, (where the dimensionless constant is chosen such thatC M1 is of order 1 M o ) rather than M. A „, compensating normalization 1/C,„ is applied to the distribution function to account for this. Using these scalings eq. (2.1) becomes, al^ G2/W 1 a [ = 16m2^ In AC „,C, ^+ ^(2.6) ro v^aE^aE 24 EE where the tildes indicate dimensionless quantities. The dimensionless diffusion integrals iE and it.E correspond to the diffusion coefficients defined in eq. (2.2a) and are given by E -== and yj^ tm 4fidEi. IEE ,^ :^ (2.7) For ease of calculation we may require the constant coefficient on the right side of eq. (2.6) to be unity using the constant C e . From the definition of v, G 2mr2r0-2 equals v o4 , so we are left with 167r 2 ln ACC,. This can all be cleared away by taking C, = (16relnAC,,,) 1 .^ (2.8) After undergoing the extensive revisions described below some check runs were done to make sure that nothing important had been changed. A single mass model starting as a Plummer 18 model was run and compared with the results in Cohn (1980). Provided that a similar algorithm for time step selection is employed, the resulting model was consistent with Cohn's. Two of the dimensionless quantities discussed by Cohn will be compared to the results from my model: for the 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 stellar evolution (see §2e) was included in the models a somewhat different scheme was used for time step selection. This leads to somewhat smaller asymptotic values of and y. I also ran two species models similar to those of Inagaki (1985) with stars having a mass ratio of 2:1 and a range of ratios in the total mass WM 1 . These gave similar ratios of core-collapse to initial half-mass relaxation times, listed in Table 1. The small differences are attributed to the difference in the time step algorithm with the somewhat larger time steps used here giving a slower core collapse. Table 1: Core collapse times for two component models as in Inagaki (1985) Mass ratio tjtth MJM l (Inagald) 0.00 0.01 0.05 0.11 1.0 9.0 15.8 14.8 10.7 8.8 9.3 14.8 15.4 13.6 10.2 8.5 9.6 13.6 Single component model Section c. Binary Heating In order to halt core contraction an energy source is required. In this study I have used binaries 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 time of Lee (1987b) almost all the Fokker-Planck studies have used the statistical treatment of three-body binaries. There appear to be two main reasons for this change. 19 The first is that the overall evolution of the models is not highly dependent on the details of the heating mechanism. The expansion is governed by the rate of energy transfer within the model. This can be seen qualitatively by the following dimensional argument which is based on those of SOC and Lee (1987a). Assuming that there is an energy source in the core of the model, the rate of change in energy, for a virialized system, is given by dv 2 (2.9) Eco °cdt ^• " It is assumed that the total mass is constant. The energy conduction rate outside the core is proportional to the product of the local rms velocity, v, and the kinetic temperature gradient v d(mv 2 ) cond^v^• dr Approximating 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 conductive luminosity to be 2 V Econd^rh (2.10) trh where vrh is the velocity dispersion at the half-mass radius, and trh, the half-mass relaxation time (Spitzer & Hart 1971), is given by 0.1380/2r 312 trh = G 1/21n In Ah (2.11) Equating Ecore and Ecohd and employing the relationship between trh and rh from eq. (2.11) and the virial theorem, one finds that for late times v r2h a t -2/ 3 , rh a t 2/3 and E oc t - 5/3 . This result is not restricted to the half-mass radius and will hold for any central energy source as long as the half-mass and central velocity dispersions are proportional.The latter condition is reasonable if the re-expansion is approximately self-similar as has been argued on dimensional grounds by, for example, Inagaki & Lynden-Bell (1983). In order to check the assertion that the evolution does not depend strongly on the details of the heating mechanism, I have run a model with identical initial conditions to that of SOC but 20 with 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 occurs at the same epoch, but it is much deeper for the model with three-body binaries. Further, the subsequent re-expansion of the core takes place much more quickly. The latter difference can be explained in terms of the energy argument discussed above. SOC show that for tidal binaries E cc pc2v0.8r: where p c and v2 are the central density and velocity dispersion and r, is the core radius defined by (King 1966): 2 3^3v r = ^c .^ 47rGp c (2.12) Employing eq. (2.12) and the proportionalities given above for E, r and v, this heating rate gives pc cc r°.8 . For three-body binaries Lee (1987a) shows that E cc p 3c v c-7 1-c3 , and thus p c cc t -2 . The detailed structure of the cluster centre will depend on the heating mechanism, but further out the nature of the energy source is less unimportant. The actual numerical models have power-law exponents which deviate a little from those predicted, but this is probably due to their not being perfectly self-similar. The second reason for the change to three-body binaries, is that the statistical three-body method 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 possible interaction channels. This would become even more complex when all the possible interactions between 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 required formalism is quite complex. An additional support for using three-body binaries, as opposed to tidal binaries, comes from consideration of the stellar population in the core. Due to mass segregation, the core will be dominated by the most massive stars. At the present epoch, these will be the degenerate remnants of stellar evolution. Such remnants do not form tidal-capture binaries, so three-body interactions will be the favoured binary formation mechanism under certain conditions. Whether this is the case will depend on the relative proportions of normal stars and remnants in the core and the 21 interaction 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 that three-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 suitable stand-in for the more realistic stochastic rate of formation and destruction. The energy input at any given time is just the product of the creation rate for hard binaries and their average energy input. To be consistent with eq. (2.1) a heating term of the form 471 2 (k) is required, where the angle brackets indicate the orbit averaging operator, and E is a heating rate per unit mass. Based on suggestions by P. Hut and D. Heggie, LFR use E(r) = CH G 5 vc2 (E (2.13) i for the heating rate per unit volume due to three-body binaries. For mass class i, the mass density is given by eq. (2.4) and Vi2 is its three-dimensional velocity dispersion. Also, v c2 is the central velocity dispersion. This formula may be thought of in the following manner. In order for three stars to interact strongly enough that the velocity perturbation is sufficient to bind two of them they need to pass by each other within /v2.. The time interval between encounters out to this distance is about (np2 v) 1 . The probability of a third star being within distance p is p an, so the - time interval between triple encounters, at distance p or less, is (n2 p5v) 1 . Substituting for p, the - formation rate will go as G 5n 2m5 /v9 . The energy generated per binary is CH MV: so the energy generation rate is CHG 5 vc2n2m6/V9. Converting this to a rate per unit volume by multiplying by n and expressing nm as p yields CHG 5 ve2p 3m 3/0. Since we must consider all possible combinations of three stars from all the different mass bins it seems natural to partition this rate as CHG 5 v,2(p im i /O(p 2m 2/v23)(p 3 m3/v33 ) for each star involved and to sum over all possible combinations. This is just the extended binomial series expressed in the cubic factor of eq. (2.13). The energy release is governed by the factor CH vc2. The constant used, CH=4210 for a three-dimensional velocity dispersion, is based on the preliminary report of Hut (1985). The 22 detailed 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-expansion begins, 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 unit mass. If the density used is the local mean density then the energy is distributed to each mass class in proportion to its relative density at each radius. The greatest heating comes in the core at core-collapse. Following core collapse the models remain highly concentrated and energy generation continues. The Fokker-Planck equation, augmented by the heating term, may be written a a ^afi +167C 2 H^L,EE^, 4n2P=W•Lmi(fiDE dE (2.14) where the heating term is given by 0 v 1.1 (E) r 2^ H E(E) CHG5vc^ 31 k ELLp^vi 3 dr. (2.15) In the dimensionless form used in the code, the Fokker-Planck equation including the heating term is a:4a 2I )+I =^ivrt +C ^ EE a ai^at^Hv [^, E - H (2.16) where ill (Ell E 7^ 1 (E) r-p-213 2.1 3dr-, (2.17) and the constant CL is given by CH C 2 (2.18) .^ 4n 2 In A The heating coefficients were calculated following each solution of the Poisson equation and CH = were used for all four of the subsequent Fokker-Planck steps. This is in contrast to the diffusion coefficients which are calculated anew for each Fokker-Planck step. The difference is that the heating coefficient depends on the density and velocity dispersion, which is only known following 23 a 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 but this is also the time when the time steps are the smallest, so this should not be too much of a problem. Section d. Tidal Boundary One of the drawbacks of using the orbit-averaged form of the Fokker-Planck equation is that it does not properly treat the evaporation of stars from the cluster (Henn 1960). The inclusion of the effects of tidal stripping will mitigate this aspect of the orbit-averaged formalism. What it will not correct for is the fraction of the evaporation due to ejection from the core as opposed to that due to diffusion past the edge of the cluster. The assumption of spherical symmetry prevents a thorough treatment of the galactic tide, the effects of which on the cluster are distinctly asymmetric. Further, the assumption of isotropy may be violated as the tidal field has more time to act on stars with low energies and high angular momenta in the halo and preferentially removes stars from the core with high energies and low angular momenta. As well, stars with prograde orbits are easier to strip. This leads to some net rotation, in the halo of the cluster at least. As shown by Oh and Lin (1992) the overall effect of tides is to make the velocity dispersion isotropic in the halo of the cluster. They predict that a transition 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 conflict with tidal stripping. The method used for implementing tidal stripping follows that of Lee & Ostriker (1987) and assumes 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 relation MG(r) j( = 2 3 x10 -5 ( r e(r) km s -1 10 10 M 0^kpc • (2.19) where 8(r) is the circular velocity at radius r in the galaxy. Innanen, Harris, & Webbink (1983) 24 found that for a galaxy with a flat rotation curve (C1(r) constant), the tidal radius for a cluster with mass Mc at a distance R G from the centre of the galaxy is given by r =-- 2 Mc 3 2M G (2.20) R. Substituting eq. (2.19) into eq. (2.20) and taking a fiducial cluster mass of 10 5 M0 gives r =10.9( )1/3 ( e(RG ) M^ 2/3 10 5 M0^220 km s -1 ( RG ) 43 pc kpc (2.21) Lee & Ostriker (1987) considered the rate of mass loss for stars with energies beyond the tidal boundary. Mass is removed from the distribution function for energies greater than E, (the gravitational potential at r e when 4:$(.) = 0). The rate of mass loss, b(E), is based on the equation of motion for a particle outside the tidal boundary. In this case, for tidal stripping afi (E,t) = at where CTfi (E,t)b(E)t; 1 ,^ ( 2.22) 312 b(E)= 1 (f)^E E — t (2.23) 0E> Et and t, —^2ic^(2.24) 1Gp e is the orbital period of the cluster about the galaxy. The tidal radius is that where the enclosed mean density equals the tidal density p r . From eq. (2.20), pt is clearly constant. Cr is a dimensionless constant controlling the overall rate of mass loss and was always set to unity. In that case the mass-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 Airi 5 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 tidal stripping rate is discontinuous at the tidal boundary, an unphysical situation. In Cohn's method 25 for solving the Fokker-Planck equation, the distribution functions are regridded in energy after each iteration of the Poisson equation solver. This is necessary to keep f a constant function of the adiabatic invariant q. The regridding is done by way of a second order Taylor expansion off in q(E). It can happen that discontinuities at the tidal boundary are reflected in the derivatives of f(E) taken for the Taylor expansion and are then amplified in transferring the distribution functions to the new energy grid. A few iterations of this can lead to a catastrophic failure of the Poisson solver. In order to reduce the probability of this happening a cubic polynomial for the stripping rate was used over the boundary region. Let ^- ,^j3 ^ 13=14 E, . (2.26) From eq. (2.23) the stripping rate goes as b(E) --= ;13 for E_ Et . For the region where 1131 S e , and for a suitably chosen e, I find a cubic polynomial b 1 03)satisfying the boundary conditions: b1 (—e) = 0, bi'(—E) = 0, bi (e) = Are- , 1 b1(e) = 2-1i . The resulting cubic is ^k r (2.27) ^a 3^2 +5 a+ 3 ] . (13) = --i-[—(-vi ) + (-vi )^L e^(2.28) Ae - This relation was used in eq. (2.25) instead of eq. (2.23) when Pi e . The transition region constant e was chosen so that two or three grid points surrounding the tidal energy would have stripping 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. This translates as the density going to zero at large radii. The expansion of the model due to relaxation drives a population of stars out to ever larger radii. When no tidal boundary is applied, the radial grid must be continually extended to encompass these stars. If this is not done, then eventually the distribution function at small energies will become large enough to contradict the boundary 26 conditions. With a tidal boundary, the stars at low energies will be continually removed. One exception to this is when the evolution of the cluster is much more rapid than the tidal stripping timescale. In that case the grid edges may still be overwhelmed. I found that using a radial limit 50 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 care must be taken to adjust the limits of the energy grid appropriately. The minimum energy should be 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 the model, it can only be applied following a Poisson step despite the fact that the stripping is done on the energy grid. One straightforward way of doing this is to simply do the tidal stripping following each Poisson step before proceeding on to the next set of Fokker-Planck steps. This was the approach taken by Lee & Ostriker (1987). The main problem with this is that once the tidal stripping is done the potential and densities will again be inconsistent with the distribution functions. If the amount of tidal stripping taking place is small, this won't present much of a worry. In the late stages of the model's evolution the mass becomes small, and since i; .. M -2/3 for a linear decrease in total mass, the rate of decrease in the tidal radius increases dramatically. In cases where the rate of decrease in I-, is significant, it becomes necessary to find a self-consistent solution. To account for this, an iterative method of doing the tidal stripping and simultaneously solving Poisson's equation was used. Following the completion of the Fokker-Planck steps, the distribution function, the functions p and q, and the previous self-consistent potential are stored. I will refer to this group of functions as model F and this initial model as F° . One iteration of the Poisson solver is done to estimate the tidal radius, r t°, and the tidal energy, E:D . Since the Poisson solver changes F, F is reset to F° . Next, tidal stripping is done using E 1° in eq (2.25), followed by a 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 the distribution function and the potential (as a result of the Poisson solver) and with respect to the tidal boundary. If r 1° and re l are not equal then the model is reset to F ° , stripped using tidal energy 27 E11 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. A fractional 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 modelling procedure, so it is inefficient to have to do too many iterations of the tidal stripping. If more than one full solution of Poisson's equation was required to get r„ then the time step for the next Fokker-Planck step was arbitrarily restricted to being no more than one-half the time step just used. This procedure worked well and prevented several numerical oddities which crept into models run without it. CW found that when tidal stripping became rapid it was impossible to find a self-consistent solution to the Poisson equation and their tidal boundary condition. They took this to be evidence for 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 stellar evolution mass loss is sufficient to dynamically disrupt the cluster. In my models, it was always possible 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 approaches to the tidal limit. The limit of CW is much sharper with f=0 for E>E„ With E dependent on the potential, they find it more useful to truncate f(q) at some value q 0 , where q is the adiabatic invariant defined in eq. (2.2c). In order to implement this they need to find a value of qo for which Alc [q0]=r,[q0] 3 in appropriate dimensionless units. The problem they run into is that as mass is lost, the implicit relationship between M c & r as a function of q o for a given f(q) becomes parallel to Md=r,3, and the boundary condition is ill-determined. For the boundary condition of Lee & Ostriker (1987) there is no discontinuity in the density at the tidal boundary, so a self-consistent solution is easier to find. The difference here is numerical; the physical instability appears to be the same for both models. 28 Section e. Stellar Evolution With the exception of CW, the effects of stellar evolution have not been included in any orbit-averaged Fokker-Planck models. StodOlkiewicz (1985) presented the results of an elaborate Monte-Carlo calculation which included a mass spectrum, tidal and three-body binaries, and stellar 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 after the massive stars have completed their evolution and that the further evolution of lower mass stars was of little importance. When this is done, degenerate remnants are included to represent the evolved massive stars, and the masses of the various mass bins is assumed to be constant. For Fokker-Planck models, which propose to model the full time evolution of a globular cluster, a static mass spectrum is inconsistent. As will be seen below, the neglect of this fact of stellar life can lead to incorrect conclusions about cluster evolution. The approach used in implementing the effects of mass loss due to stellar evolution is much the same as that used by CW. Each mass species is assigned an initial mass based on an average across the stars in the bin. Considering the likely outcome of stellar evolution on stars of that mass, a final mass is also assigned. The time range over which the stars in that species go from their initial mass to their final mass is set by the main sequence lifetimes of the stars at the edges of the bin. The mass per object for that mass species is taken to change linearly from the initial mass to the final mass over that time span. Given the degree of simplification used here it was not necessary to be any more exact in this scheme. Care was taken during the running of the models that the total mass of the cluster should not change due to stellar evolution by more than 1% over any time step. Let mii and mf be the initial and final mean masses for stars in bin j, and and and ti be the time range over which stellar evolution takes place. Then, the mean mass per object for bin j at time t is given by t^tji (t) = ( Am ) (t — ti At 1 <t< t> 29 ' (2.29) ^ where ( Am ) rriC^(2.30) At );^tC is always negative for mass loss. The rate of change of mass with respect to time during stellar evolution 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 the relaxation time. From eq. (2.11) tth.rh '. The stellar lifetimes are independent of any of the other scales. Hence, for otherwise identical models, models with larger scale radii (eg. half-mass radius) will suffer more mass loss per relaxation time than those with smaller scale radii. This effect 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 the three dimensions involved, one is set by taking of G=1. When binary heating was introduced a second dimension, that of mass, was set by the relationship between the relaxation time and the heating time. This can be seen from eq. (2.16) where the constant CL. (eq. 2.18) depends on the total mass via the lnA factor. When stellar evolution is included, the final free scale, that of length, must also be chosen. With all the scalings needing to be predetermined the usefulness of a dimensionless model is lost, but the formalism has been retained for backward compatibility. It should be emphasized that in the models including stellar evolution, there are no scale parameters left 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 with a logarithmic spacing with the initial mass per object per bin being the geometric mean of the masses of the bin boundary. The mass within a bin, Mt = JN(m)m dm, was integrated using a power-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 1m i . The stellar lifetimes are based on a cubic spline of the Miller & Scalo (1979) compilation for Population I stars. The final masses for stars with m i <4.7 Mo , which become white dwarfs, were 30 from the formula of Iben & Renzini (1983) with TI= 1 /3: mf=0.57+0.22(m/-1). More massive stars, 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, no evolution was included. In other words, stellar evolution was assumed to terminate at 15 Gyr. By this time the amount of mass loss due to stellar evolution has dropped to a low level and is of little importance to the dynamical evolution of the cluster. To demonstrate this, consider a simple model where the mass spectral index is x, the main sequence lifetime scales as m -3 and a star with initial mass m loses mass 6m through stellar evolution. Then the rate of mass loss is fit 0= 6m m 3- ". If we compare 0.9 M o stars evolving into 0.6 M o white dwarfs (&n=0.3 M o ) at 15 Gyr with 5 M o stars being completely destroyed (8m=5 M o ) at 80 Myr, then for each type of star to give the same contribution to the mass loss rate, x=4.6 is required. Mass loss is never important for such a steep 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 the independent Fokker-Planck code of CW, I will present the results of one specific run. The parameters for this model, designated CW1, are listed in Table Al in Appendix A. This model is analogous to a model with W 0 =7 and a MSI x=1.5 in family 3 of CW. (Note that CW specify their power-law MFs in units where the Salpeter MF has slope 2.35. W 0 is the dimensionless central 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 in the scaling which they use to define their various families. In my models no such freedom exists so 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 —36 Gyr. 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 the core 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. 31 Figure 1 can be compared with Fig. 23 of CW, and shows the evolution of the projected mean mass as a function of projected radius. I have taken snapshots from my model which have about the same time distribution with respect to the times of greatest expansion and core collapse as in CW. The various curves are marked with their times in Gyr. Greatest expansion occurs at 8.3 Gyr. The overall morphological behaviour is quite similar although the central mean mass never becomes as small as it does in CW. Figure 2 shows the projected surface densities of the various stellar populations at a time near core collapse. The time shown is for when model CW1 has about the same total projected central density as the model shown in Fig. 25 of CW. The solid, short dashed, long dashed, and dotted lines are the projected surface densities for the unevolved stars, the neutron stars, the white dwarfs, and all stars respectively. The overall profile is in agreement with CW but the relative proportions of the various stellar populations is different. The reason for this will be explained in connection with Fig. 3. Figure 3 shows the evolution of the mass function for comparison with Fig. 33 of CW. The morphological development is much the same with the exception of the discontinuity in the MF for masses less than 1M o at late times. It appears that CW do not terminate stellar evolution when the 0.8 M o stars evolve, but continue it to lower masses. This difference also explains why I have a higher proportion of luminous stars than CW do, and why the mean mass in model CW1 never 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 stars with m i>0.8M 0 ) as CW, with Mi=10 5 M o ,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 were disrupted. The remaining W0 =7 models all reached core collapse. (The model with x=2.5 and r e=l5pc has such a small amount of mass loss during even the first few relaxation times that the energy input does not balance the energy loss from the core due to two-body relaxation; the core-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 the results of CW is satisfactory. 32 1.2 44 0 A 0 8 . 27 V 0.6 0.4 I —2 0 log r (pc) 1 2 Fig. 1.—Projected mean mass as a function of radius for model CW1. This figure may be compared with Fig. 23 of CW. The curves are marked with the time in Gyr from the start of the model. The minimum in the central density occurs at 8.3 Gyr. 33 02 C) 0 b b0 0 ---, —2^—1^0 log r (pc) 1 Fig. 2.—Projected surface mass density profiles for the various mass components in model CW1. The curves show the projected surface density for the neutron stars (short dash), white dwarfs (long dash), main seqence stars (solid) and for all stars (dotted). The time for this figure was chosen so that the total central density is about the same as that of Fig. 25 of CW. The relative proportions of white dwarfs and luminous stars is different than in CW since no stars with masses less than 0.8 M o were allowed to evolve in model CW1. 34 I^I^I^I^I^I^I^I^I^I^I^I^l^I^I^I 1 0.8 z II z 0.6 rap 0 z —1.5 0.4 02 I^ 1^ 1^ 1^ —0.4 —0.2 0^0.2 0 4 log rni (M 0 ) 1^1^I^ 0^0.5^1 log mi (M 0 ) 1^1^1^ 1^I^I^ 1^I^I Fig. 3.—(a) The number of stars as a function of mass normalized by their initial values for model CW1. The times were selected to have the same normalized number of stars for the least massive species as in Fig. 33a of CW. The break in the lines connecting the points differentiate between the evolved and unevolved mass classes. The three initially heaviest mass classes have been merged into one 1.4 M c, mass class, and the stars in three next mass classes were entirely destroyed. (b) The number of stars as a function of mass normalized to the total initial number of stars. The counts are shown for only the unevolved mass bins. The line without the circles marking the masses of the various mass species is the initial number of stars. This figure is comparable with Fig. 33 of CW. 35 Table 2: Results of sample grid of models including stellar evolution Initial mass (M) 105 Wc,^r, MSI (pc) 0.5 1 15 Fate a t r (Gyr) rb i.e.: (Gyr) 1 1 1 1 1 15 15 60 60 60 1.5 2.5 0.5 1.5 2.5 D D D D D D 10 5 10 5 10 5 10 5 10 5 105 4 4 4 4 4 4 15 15 15 60 60 60 0.5 1.5 2.5 0.5 1.5 2.5 D D PCC DC D D C 0.07 2.6 6.9 0.20 7.0 >45.0 >45.0 105 105 105 105 10 5 10 5 7 7 7 7 7 7 15 15 15 60 60 60 0.5 1.5 2.5 0.5 1.5 2.5 DC DC DC D PCC PCC PCC 0.54 3.9 6.3 6.3 >29.9 >23.1 6.0 17.0 19.8 4x10 5 4x10 5 1 1 23.8 2.5 95.2 2.5 D D 4x10 5 4x10 5 4x10 5 4x10 5 4 4 4 4 23.8 1.5 23.8 2.5 95.2 1.5 95.2 2.5 D DC D C 105 105 105 105 10 5 0.04 0.21 2.3 0.07 0.70 7.1 2.4 2.8 7.8 3.0 23.5 6.4 >30.2 >30.2 4x10 5 7 23.8 0.5 DC 1.5 4x10 5 7 23.8 1.5 DC 12.8 7 23.8 2.5 DC 4x10 5 21.6 4x105 7 95.2 0.5 12.8 DC >30.2 4x105 7 95.2 1.5 C >30.2 4x105 7 95.2 2.5 >30.4 >30.4 C a D: Disrupts with expanding core, C: Core is contracting, D C: disruption with contracting core, PCC: Post-core-collapse, will eventually disrupt, D PCC: Disruption in post-core-collapse phase. b^• Time of disruption Time of maximum central density Section f Initial Conditions The formation mechanism for globular clusters is still a matter of great debate. It can probably be assumed that whatever the means by which the initial gas cloud fragments into stars the resulting system will not be in dynamical equilibrium. Such a system is expected to undergo violent relaxation (Lynden-Bell 1967). The velocity distribution for the stars following violent relaxation is independent of the mass of the stars. A corollary of this is that initially there will be no mass segregation. One common choice for the initial model is an n=5 polytrope (Plummer's model) which has the distribution function f (E)= 24V-27 E 3.5 . 77E 3 (2.31) The Plummer model has the disadvantage of being unbounded, though of finite total mass. An alternative choice are the models of King (1966). This family of models is based on an approximate solution of the Fokker-Planck equation in a square-well potential. The distribution function takes the form of a lowered-Maxwellian f (E) = C(e -BE —1). (2.32) If, using Jeans's Theorem, it is assumed that the same form for the distribution function applies everywhere in the cluster, this distribution function can be used in a solution of Poisson's equation to give a full cluster model. The King models form a one parameter family of dimensionless models. The parameter usually used is the dimensionless central potential, W o . These models are naturally truncated and cover a range of concentrations depending on the value of W 0, with the concentration increasing with W o . A more detailed description can be found in King (1966). A King 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 logarithmic spacing. The IMF was then determined by the number of stars placed in each bin. While a single power-law is the most common parametrization of the IMF, models were run with other ways of setting 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 remnants 37 along with unevolved, luminous stars. For the runs involving stellar evolution, specifying the IMF is much simpler, but an mf(m i) relationship and a set of stellar evolution times must be adopted. I used the relations used by CW and described above. One problem with using these models for comparison with observations is the large parameter space involved. The structural parameters of the initial model are: the dimensionless central potential, Wo, of the King model, which gives the initial distribution function and potential; the initial total mass, which sets the relative rates for relaxation and binary heating; and the tidal stripping timescale, C r The stellar populations are defined by the initial grid of bins and the IMF. The inclusion of stellar evolution also requires the specification of the final mean mass for each bin and the time interval over which stellar evolution takes place. Section g. Caveats Before going on to discuss the applications of these models to actual observations it is useful to summarize the potential problems with models of this type. These Fokker-Planck models are by their very nature numerical approximations to statistical approximations to the dynamical evolution of a globular cluster. They assume that the number of stars is sufficiently large that there are a reasonable number in each bin of the distribution function. The inclusion of the tidal boundary ensures that if the simulation is run long enough this assumption will be violated. The very high central densities sometimes reached may also indicate a breakdown in the approximations used. Since, as described by Makino & Hut (1991), the core mass decreases with decreasing core radius in self-gravitating systems, we see the curious situation of the number of stars in the core dropping to zero as the central density goes to infinity. Since it is the density and velocity dispersion as a function of radius that go into the heating rate, there may be too few stars within the core to even form a three-body binary. In a certain sense, three-body binaries may not be sufficient 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 discrete events 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. 38 Chapter 3. The Problem of M71 Section a. Outline of the problem The work described in this chapter has previously been published in Drukier, Fahlman, and Richer (1992, hereafter DFR). While the observations of M71 described in §2 of DFR are the responsibility of my co-authors, the models and their analysis, which are described in the remainder of the paper, is my work. This chapter will provide a summary of DFR in order to lay the groundwork for the remaining chapters of this thesis. All of the models which are presented in the figures of this chapter have been rerun using the final version of the Fokker-Planck code. Some small differences will be seen compared to the diagrams in DFR, mostly due to the self-consistent treatment 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 also provided the information required to produce surface density profiles for the brighter stars. In DFR the 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, the surface density profiles are basically King-model like, but the one for the brighter stars does show a slight central enhancement over that expected from a King model. Second, a significant degree of mass segregation is seen between the two fields considered. All of the models described in this chapter were run before the effects of stellar evolution were added to the Fokker-Planck code. The initial models are as described in §2f. In most of the models discussed in this chapter the main sequence was represented by ten mass classes covering stars with masses between 0.16 and 0.89 M o . The low mass cutoff was somewhat arbitrary in that it was the low mass limit to the VandenBerg isochrone used, but the initial model MFs do extend into the range of lower limits for observed globular cluster mass functions (Richer et al. 1991), and past that observed in the 3' field. The model MF could be extended to even lower masses, but mass segregation makes them unimportant to the dynamics of the core of the model. One model, to be 39 I iiii i mi l i,11111. -^ (a) - _^i 4.5 3 x I 1 3 1^1^1^1^1^1^1^1^1^1^l^1^1^1^1^1^1^1^1^1 — 0.8 —0.6 —0.4 —0.2 0 log m (M0) —0.5^0^0.5 log r (pc) Fig. 4.—(a) The surface density profiles for the two brightest magnitude ranges observed in M71. These correspond to the mass ranges 0.85 <m/M o <0.89 (filled circles) and 0.71< m/M o <0.85 (open circles). The crosses show the total density over the entire luminosity range. (b) The mass functions 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 a mean distance of 3' from the cluster centre. 40 mentioned below, which was calculated with an additional low mass bin that extended the low mass cutoff to 0.1 M o , supports this conclusion. The mass bins for the higher mass stars were chosen to be the same as the observational bins for the core ME In addition, three bins were used to represent a distribution of white dwarf stars with masses from 0.9 to 1.14 M o . The three white dwarf bins were of equal width and an equal number of stars was put in each, making up a predetermined number fraction of the stars in the cluster. The main sequence stars were distributed as 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 yet more complicated. The mass spectral indices of the two power laws, their relative proportions and the number fraction of white dwarfs determined the IMF used for each model. The philosophy followed for all the comparisons between models and observations in this thesis is to transfer the results of the Fokker-Planck models onto the observational plane. The modelled data were then compared directly with the observations. The selection of the time when a model best matched the observations was done subjectively. As will be seen, the models did not match the M71 observations particularly well, so more objective matching procedures were inappropriate. Most of the comparisons involve surface densities, either directly, in the surface density profile, or via a mass function. Periodically when it is run, the Fokker-Planck code saves a full description of the current state of the model. For each of these save times the density profiles for each species were extracted and the projected densities calculated. In order to compute a model SDP, the projected surface density for the appropriate species is integrated across the same regions as were the observations. For example, if the observed SDP was calculated for 10 logarithmically spaced 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. For mass functions, where the observations were made in a region of given dimensions and position with respect to the cluster centre, the projected model surface density for each mass species would - 41 be integrated over the same region. This procedure ensures that the observations and the modelled observations are subject to the same degree of averaging. In the case of velocity dispersions, the obervations were just compared with the line-of-sight velocity dispersion profile. A more detailed comparison could be done if the radial distribution of the observed stars were known, but this information 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 M o mass bin (see Table A2) was used for the 0.85<m/M o <0.89 SDP, and the sum of the 0.816 M o and 0.743 M o mass bins was used for the other SDP. The model core MF was calculated by integrating the projected density profile for each species over a circular area of radius 3.66 pc (This assumes a distance of 3.68 kpc [Richer & Fahlman 1988]), and dividing by widths of 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 position with respect to the cluster centre as that observed by Richer et al. (1990). The areas covered by the two MFs overlap to a certain extent but the core MF is dominated by the higher density of stars in the core of the cluster. The stellar density in the 3' field is in good agreement with the SDPs from the core field. In order to discuss the degree of mass segregation more quantitatively I define a segregation measure, S, as the logarithm of the ratio of two mass functions. In this way, the effects of mass segregation can be examined independently of the mass function itself. Segregation measures could be 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 the difference in absolute numbers due to changes in density the segregation measure may also be normalized in some manner. Being the ratio of two quantities, the errors in observed segregation measures tend to be large, but it still serves as a useful tool. In this thesis I will discuss radial segregation measures of the form (m, a; b) = log{Nm (m)[field a]/N m (m)[field , (3.1) where Nm (m)[field a] is the number of stars per unit mass of stars with mass m in the field designated 'a'. A larger value of S r indicates a relative surplus of stars of that mass in field 'a' over field 'b'. 42 For M71 I consider the segregation measure S r (m, 3'; core) which is shown in Fig. 5. That the segregation measure is negative for all masses reflects the lower stellar density in the more distant 3' field with respect to the core field. S r (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 stars dominate the core, leaving relatively more low mass stars further out. On the other hand, the deficiency of stars with m=0.55 M o in the 3' field is difficult to understand. In is hard to see why such a short mass range should go against the general trend determined by mass segregation. The existence of this dip in the data and the large errors associated with S,. should serve as a warning against over-interpreting the observed segregation measure. To reduce the differences in the MFs to one number we can consider a segregation measure normalized by mass over the mass range observed. This is defined as S,.,.(mp a;m2 ,b) S,.(m„a;b)— Sr (m2 ,a;b).^(3.2) For the M71 observations the normalized segregation measure is SrM: Sr,m (0.35 M o , 3'; 0.85 M o , core) = 0.4.^(3.3) If the mass function is a power-law form then xa — xb =^(ml , a; m2 , b)/log(m i /m2 ) .^(3.4) For M71, x3 , — xco,. — 1 In order to try and find a match to the observations of M71 I ran a large number of models with varying MFs. The MFs differed in the MSIs and relative weights of the two power laws, in the number fraction of white dwarfs, and, for some special cases, in other ways as well. The MSIs ranged from -1 to 5 in various combinations. It was found that in order to match the steepness of the mass function for the 3' field a very steep IMF, with MSI around 5, was required for the least massive 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 to the three, lower-mass, white dwarf bins discussed above. This change did not greatly affect the results. I also ran a model, identical to the standard model to be discussed in some detail below, 43 -0.8 —1.2 —1.4 I^I^I^I^I^I 0.4^0.5^0.6^0.7^0.8^0.9 m (M 0 ) Fig. 5.—The segregation measure defined in eq (3.1) for the 3' field with respect to the core field. 44 but which had an extra low-mass bin extending the mass function down to 0.1 M o . This bin contained 58% of the cluster's mass and 80% of its stars; the evolution of the core was not very different from that in the standard run. The main difference was that the model evolved somewhat faster, reaching core collapse in only 80% of the time required for the model without the extra low-mass stars. The best-fitting (in an eyeball sense) of the models run was the one with the MF and initial conditions listed in Table A2. I shall refer to this model as the "standard run". In Figures 6 and 7 I show 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 1 pc) when the mass is 0.09 Af t (thin solid line). The main problem in comparing this model with M71 is obvious. The model SDPs at small radii do not flatten off as do the observed SDPs but continue in a power-law fashion to the inner limit of the observations. The slope of the power law section of the SDP does not flatten significantly during the entire PCC evolution. This behaviour is common to all the models with initial conditions discussed thus far. Adjusting the scale length does not noticeably improve the fit to the outer points of the SDPs and certainly does not do anything for the core.While the model does not reproduce the detailed structure of the MFs the segregation measure is similar to that observed. These Fokker-Planck models show the following general behaviour. First, the models reached core collapse in of order ten or less of their current half-mass relaxation times. Second, only very mild amounts of mass segregation are seen before core collapse, certainly an insufficient amount in comparison to that seen in M71. Third, the PCC models all have very steep central SDPs. This steep profile has been been considered to be the identifying feature indicating that core-collapse has taken, 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 class would be those relatively unevolved clusters with long relaxation times, little mass segregation and King-model-like surface brightness profiles. I will refer to clusters with these signs as "unevolved" 45 o3 N Q) 4 0.85<m/Mo<0.89 2 b 0 0 —0.5^0^0.5^1 log r (pc) —0.6 —0.4 —0.2 log m (M 0 ) —0.5^0^0.5 log r (pc) 0 1 log m (M 0 ) Fig. 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 King model (dot-dash line), at maximum core collapse (thick line), at the best fitting time (when the mass is 9% of its initial value) (solid line), and when the mass is 2% of its initial value (dashed line). The detached line segments in the MF panels represent the white dwarfs. 46 -0.8 —1.2 —1.4 0.2^0.3^0.4^0.5^0.6 0.7 0.8 0.9 1 m (M 0 ) Fig. 7.—A comparison of the segregation measures for the "standard" Fokker-Planck model with the observed segregation measure of M71. The lines have the same meanings as in Fig. 6. 47 although some dynamical evolution will have taken place. On the other hand, a highly evolved cluster would be expected to show strong signs of mass segregation, a short relaxation time and a steep 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-expansion phase. Into which class of cluster should M71 be assigned? As discussed in §4 of DFR, the present half-mass relaxation time of M71 is of order 10 8 years. Assuming an age of order 15 Gyr, the cluster has been around for about 100 of its present half-mass relaxation times. As well, we see evidence for a large amount of mass segregation. These two signs would suggest that M71 is a core-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-collapse group. Identification with either class presents clear problems and all that can be concluded is that the observed properties of M71 are inconsistent with the models as they stand. If we are to find a model 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 that of 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 new physics is required to resolve the question posed by M71. These attempts are reviewed in the next section. Section b. Approaches to a solution Subsection I. Gravothermal Oscillations. As discussed in the introduction, PCC globular cluster 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 only the 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. Perhaps 48 we are also observing M71 at the low density extreme of a GTO, and it is for this reason that the central cusp in the SDP is absent. However, if the observations are sufficiently removed from the centre of the cluster, or have insufficient spatial resolution to isolate the core, then the data should not depend significantly on the existence and amplitude of GTOs. The models show that at maximum 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 core collapse is to reduce the power-law slope in the density profile of species less massive than the one that dominates the core. In an attempt to see if this effect would be of aid in resolving the problem presented by M71, I ran a series of models which contained an extra class of massive degenerate stars. 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 were included in the models. Fig. 8 shows this "black hole" model at a late time when only 7% of the initial mass of 6x10 5 M 0 remained, and some 47% of the current cluster mass was in the form of black holes. For the more massive stars, the fit to the surface density profile is very good. The problem 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 the morphology of the observed MFs are quite different from each other. This case may seem somewhat extreme, but when I reduced the fraction of black holes modestly, the slope of the inner SDP became 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 the mass ratio between the heavy component and the visible stars is then too small to provide sufficient flattening. A further problem with the heavy remnant hypothesis arises when the kinematics of the cluster stars is considered. Pryor (1990) has measured the velocities of 79 members of M71 and gets a velocity dispersion of 2.15±0.17 km/s at a mean radius of 1.5 pc. In the "black hole" model, at the time which best matched the SDPs, the central velocity dispersion for this mass class is 8.2 km/s dropping 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. 49 5.5 3 5 2 0 3.5 3 ^1^1^1^1^1^1 —0.8 —0.6 —0.4 —0.2 0 log m (M 0 ) —0.5^0^0.5^1 log 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 that best matching the more massive stars' SDP. The symbols for the observations are as in Fig. 4. 50 Subsection III. Modifying the heating rate. The reheating source in these models is due to the statistical formation and destruction of three-body binaries. Following Lee (1987b), who adopted a heating rate based on the work of Hut (1985), I assumed that the amount of energy generated per binary formed was some 100 times the binding energy of the binary. If a more vigorous energy source is available to reverse core collapse then shallower-sloped SDPs are possible. I therefore ran a series of models identical with the standard model except that the heating constant was set to various multiples of the standard one ranging from 10 to 10 6 times the energy input per binary. As the heating rate was increased the depth of core collapse decreased. Each model was checked against the observations and only in the most extreme case, with a heating rate 10 6 times the 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 now insufficiently steep—although the model SDP for the less massive stars has somewhat too many stars at all radii. This can be accounted for with a slightly different IMF. The kinematics of this model 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 very well. 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 with those for the best fitting times for the two models. The amplitude of the mass segregation in the models 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 nature of the heating source. As discussed in DFR, the total amount of energy released in the enhanced heating model is only slightly less than the total amount of energy released in the standard run. The difference is that in the standard run the heating rate is effectively zero until the cluster core becomes sufficiently dense for binary formation to occur, while in the enhanced heating model energy is released at a more gradual rate, but without requiring such extreme densities. The high concentration reached in the standard models then remains imprinted on the density distribution for all later times. Models using tidal binaries as an energy source do not need to reach such high 51 11 1 111 1 11 11 11^1^1 1 (b) 3 /- 4.5 2 0 3 —0.5 0^0.5 log r (pc) II I III I III I III I II —0.8 —0.6 —0.4 —0.2 0 log m (M 0 ) 1 Fig. 9.—Comparison with the observations of the "extra heating" model with a heating constant 10 6 times that of the standard model. (a) Surface density profiles. (b) Mass functions. The detached line segments are the white dwarf bins. At the time shown the model has shrunk to 3% of its initial mass. The symbols for the observations are as in Fig. 4. 52 2.5 0 2 U) 0 0 Q) 1.5 —0.5 0^0.5 log r (pc) 1 Fig. 10.—Line of sight velocity dispersion from Pryor (1990) for the giants of M71 compared with the predictions from the "extra heating" model for the same time as Fig. 9 53 -0.8 —1.2 —1.4 0.2^0.3^0.4^0.5^0.6 0.7 0.8 0.9 1 m (M 0 ) 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. Neither model shows the observed amplitude of mass segregation. 54 concentrations before core collapse is reversed (SOC, Lee 1987a). Even so, for the model presented by SOC, a core radius greater than 0.1 pc is not reached until hundred of billions of years have passed. 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 globular cluster will prevent a deep core-collapse from taking place while still allowing for dynamical evolution. Whether the amount of dynamic evolution would be sufficient to produce the observed mass segregation is a question which requires an answer. In DFR I discussed two possible sources of this extra energy, both of which must occur to some extent in globular clusters. The first is due to interactions of the cluster stars with primordial binaries in the cluster. Gao et al. (1990) have run Fokker-Planck simulations which explicitly include a population of initial binaries along with one mass class of single stars. In these models the time to initial core collapse increases as the fraction of binaries increases. When 20% of the initial 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 the Gao 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 abundance of binaries and their rate of destruction, the presence of primordial binaries would result in a relatively large core and hence a SDP which does not show a cusp. On the other hand, such a core is likely to consist mostly of binaries. Core collapse in the models of Gao et al. terminates when the mass segregation between the binaries and singles is complete, again with a core dominated by binaries. 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 I have chosen to pursue here, and a discussion of the models which include stellar evolution follows in Chapter 5. To further understand the problem posed by M71 it would be helpful to have data from another cluster for comparison. In the next chapter I present observations of NGC 6397 with this objective. 55 Chapter 4. NGC 6397 One question raised by the contradiction between the models and observations discussed in the context of M71 is the prevalence of this discrepancy. In the interests of getting a slightly wider view, I took advantage of an extensive set of observations of the cluster NGC 6397 acquired by G. 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 galactocentric distance (7.4 kpc for M71, 6.9 kpc for NGC 6397 [Webbink 1985]) but M71 has the kinematics of a disk cluster while NGC 6397 belongs to the galactic halo (Cudworth 1992). NGC 6397 is more 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 NGC 6397 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 is dynamically evolved. From the multi-mass King models of Meylan & Mayor (1991) a half-mass relaxation time may be estimated using eq. (2.11). With their mean model values for the total mass, 1.0±0.1x10 5 M o , and for the half-mass radius, 5.2±0.3 pc, and taking the mean stellar mass to be 0.4 M o , 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 concentrated globular clusters including NGC 6397. They found that the central cusps become bluer towards the cluster centres and attribute this to the destruction of red giants and subgiants, presumably by dynamical processes. NGC 6397 (a 1950 =17 h 36m 38 8 , 8 1950 = –53° 38.9') lies at 1=338.165°, b=-11.959°. Reddening measurements 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. VandenBerg et al. also find evidence for a small amount of variable reddening. From the data in their paper this appears to be at the 0.013 magnitude level, and should not affect the results below where a constant reddening of E(B–V)=0.19 is assumed. 56 There are two sources for the distance modulus to NGC 6397. The first is a venerable measurement of the magnitude of the horizontal branch by Cannon (1974). Depending on the choice 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 the Mv, (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 theoretical isochrones of VandenBerg & Bell (1985) as modified for enhanced oxygen abundance by McClure et al. (1987). A better fit is achieved if (m —M) v = 12.4 is used. Anthony-Twarog et al. attribute some of the discrepancy to the bolometric corrections used to transfer the luminosities onto the observable plane. In any event, their distance moduli are consistent with (m —M) v = 12.3 adopted here. These values for the distance modulus and reddening gives a heliocentric distance of 2.2 kpc to NGC 6397. At this distance 1 pc subtends 1.6'±0.2, given the uncertainty in the distance modulus. The relatively short distance to the cluster allows for very faint stars to be observed with only moderate investments of telescope time. Fahlman et al. (1989, hereafter FRST) have previously done deep starcounts in NGC 6397 and observed the luminosity function to MI = 11.5. This corresponds to a mass of 0.12 M o . Since the present observations are in the /—band, an apparent distance modulus in that colour is required. 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) gives a value of 1.38 for this ratio. On the other hand, Grieve's (1983) calibration of the colour excess ratio 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 these values for the ratio of the reddenings, I adopt an apparent I distance modulus (m—M)=12.0. This is the same as the distance modulus used by FRST. The observations to be discussed in §4a have been used to produce surface density profiles and 57 mass 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 very incomplete 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 observed give useful leverage on the question of mass segregation in this cluster. Section a. Observations & reduction The observations of NGC 6397 were obtained during an observing run at the Las Campanas Observatory 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 of several fields using the 2.5m du Pont Telescope. Observations of four fields will be discussed here. The second set of observations were obtained with the 1 m Swope Telescope and covers both four overlapping fields extending from the cluster centre, and two regions further out which overlap 2.5 m fields. The positions and designations of the fields are listed in Table 3. Figure 12 shows the locations of the program fields with respect to the cluster centre. The fields are labelled with their names from Table 3 except that "du Pont" has been shortened to "duP". All stars in the program frames observed with R13 have been marked by points. The relative sizes of the points gives 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 same galactic 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 fields with standard stars were also take on nights with photometric conditions. The 2.5 m observations were used to construct luminosity and mass functions, while the purpose of the 1 m data was to observe the surface density profile. On the 2.5 m, the plate scale was 0.162"/pixel. As the seeing was usually around 1" FWHM, this resulted in heavily over-sampled star images. In order to improve the signal-to-noise ratio, the frames were boxcar averaged 2x2 following debiasing and 58 flat fielding, giving an effective plate scale of 0.324"/pixel. For the 1 m frames the plate scale was 0.435"/pixel and the seeing was typically around 1.7" FWHM. Table 3: Field locations for NGC 6397 observations Exposureb Xa^Ya^dist.a^width^length Field Name (all positions and dimensions in arc minutes) times 5.775 5.745 0.12 1x10s, lx1Os 0.03 -0.12 Swope:c 5.771 5.771 2x5s 1.29 -4.31 4.50 Swope:sa 2x5s 7.60 5.786 5.757 3.64 -6.66 Swope:sb 2x30s 3.33 5.778 5.778 -0.25 3.32 Swope:na 11.1 5.778 5.778 2x30s -0.06 11.09 Swope:w 11.9 2x30s 11.73 5.778 5.778 -1.98 Swope:nb 60 5.778 5.778 2x30s Swope:bf 2.138 2.111 14x200s, 13x200s 3.86 du Pont:if 0.67 3.79 10.48 10.7 2.155 2.130 8x450s, 8x450s -2.16 du Pont:n 11.2 2.144 2.122 4x600s, 5x600s -0.25 11.15 du Pont:w 2.140 2.140 7x600s, 6x600s 60 du Pont:bf 2.160 2.160 12x450s 6.46 FRST -5.59 3.24 a Position 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 registration and averaging, where required; and image trimming were done using procedures within IRAF. The flat fields for most of the 2.5 m data were exposures of clouds from the first night of the run. This flat field did not work very well for the du Pont:bf field so a combination of this flat and dome flats was used. Since the 2.5 m data were acquired under non-photometric conditions, the exposures were weighted by their mean flux levels during the averaging. The sigma clipping algorithm of the IRAF task imcombine was used in order to remove cosmic rays. The deepest set of exposures in each field were combined into two independent images for that field. While somewhat fainter stars could have been recovered if all the data had been combined into a single frame, the advantage of working 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 procedure far outweighs the loss in limiting magnitude. 59 15 I^ - ^I— i duP:n 10 Swope:nb •_______e_ duP:if E .^ . . *. Swope:na ..•^..• .. . • II, ..••,,- ' • ,.... i ..!ti , V duP:w . ;•• • •^• Swope•c • Swope•w • •, FRST Swope:sa — 10 — 10 Swoper ^ IIIIII^ 1---1 ^I^1^1^ L —5^0^5^10 X wrt centre (arc minutes) ^15 Fig. 12.—Positions of the program fields observed in this study. The coordinates are with respect to the adopted cluster centre marked by a `+'. The dots are all stars with R13 observed in the program fields. Most of these also appear in the plate of Cannon (1974) and can be used to locate these fields. The abbreviation "duP" indicates the du Pont fields which were selected to avoid any bright stars.They are rotated slightly with respect to the Swope fields. "FRST" is the field observed by Fahlman et al. (1989). 60 The data in fields Swope:nb, Swope:w and Swope:bf were taken under photometric conditions and 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) and several other photometry calibration programs kindly provided by P. B. Stetson. The transformation to the standard system is good to within about ±0.03 mag. The calibration for the corresponding 2.5 m frames was transferred from the 1 m observations. While the remaining four 1 m fields were apparently observed under photometric conditions, it did not prove possible to use the standards observed that night. The fluxes measured for these stars were smaller than expected by a factor of about three for unknown reasons. Fortunately, the Swope:sa and Swope:sb fields overlapped that of FRST and a calibration, good to within the ±0.05 quoted by FRST was obtained from magnitudes of 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 their star 1. This star is identified with star 335 of Woolley et al. (1961) which was taken as the centre by Auriere (1982). Having identified this star in field Swope:c, I took as the centre a position just off 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 m images were found to be rotated slightly with respect to the 1 m frames. The positions of the Swope:nb and Swope:w fields with respect to the central four 1 m fields were found using photographic images of NGC 6397. The positions of several stars in the two pairs of regions were measured on the photograph with a measuring engine and identified on the program frames. The required offset between the origins of the Swope:sb and Swope:w fields, and between the origins of the Swope:na and Swope:nb fields, were solved for in a least squares fashion together with a constant scale factor. No rotation was apparent in the transformation relationship and in the final calculation was assumed to be zero. The resulting positions for these fields with respect to the cluster accord well with a sketch of the field positions made at the telescope. The positions of the du Pont:n and du Pont:w fields, which are contained within the Swope:nb and Swope:w fields respectively, were calculated with respect to the positions of these 1 m fields. 61 The data was reduced using versions of DAOPHOT and ALLSTAR (Stetson 1987) in the usual manner. Two finding passes were made on the images with a find threshold of 46. Generally, a Moffat 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 used but it did not materially affect the photometry. Each of the two independent images for each field was reduced separately. Stars on the two lists were matched and the mean offsets in magnitude and position were calculated. These offsets were then applied and a final star list for each field was produced 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 usual procedure of adding artificial star images to the program frames and then re-reducing them. The artificial stars were added with a magnitude distribution based on that estimated from the initial reduction of the frame. Ten percent of the number of stars observed in the magnitude range of interest were added at random to each artificial star frame. The same artificial stars, with appropriate magnitude and position offsets, were added to both images in a pair. Ten to thirty such pairs of frames were reduced for each field. The final star lists were prepared in the following manner. For each field, regions surrounding saturated 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, and objects detected, in these regions were ignored during further analysis. Also ignored were stars found in the region lying outside that where artificial stars were added. The effects of these restrictions are admittedly small. The detection lists for the two images in a pair were then shifted onto common position and magnitude zero points. The stars on the lists were matched and all objects 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 the artificial stars added to that pair of frames. All valid detections with positions within 1 pixel and magnitude differences less than 0.7 mag with respect to the artificial star list were counted as recoveries of artificial stars. The remaining, presumably real, stars were put into a separate list. 62 Star counts were performed on the list of real stars from each artificial star test separately and the results averaged. Small variations in the number of stars were seen from test to test due to variations in the local crowding conditions. Errors in the star counts are a combination of the Poisson error and the standard deviation in the number of stars recovered. The observed luminosity function is a convolution of the true luminosity function with the probability 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 inversion of the observed luminosity function to give the true one. The luminosity functions were corrected for incompleteness using two techniques. In the simplest approach, all stars added in a given magnitude 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 underestimate the uncertainties in the corrected counts and does not take into account the effects of stars with true magnitudes lying in one bin being recovered in a neighbouring bin. This bin jumping can be corrected 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 et al. (1987). As an intermediate method between using the full matrix, as in Drukier et al., and compressing all the elements onto the diagonal, as in the approach discussed above, a tri-diagonal recovery matrix was constructed. The advantage of a tri-diagonal matrix is that it is straightforward to invert and calculate the errors for. The inverted matrix, together with some assumptions about the recovery rates and magnitudes for stars in the faintest bins, was used to calculate the completeness corrected counts. The results of this method are for the most part consistent within the uncertainties with those of the simpler approach discussed above, except that the error estimates and magnitude limit to the star counts are more conservative using the matrix approach. 63 Section b. Surface Density Profile The surface density profiles were produced from the 1 m data. Using the artificial star tests for the Swope:c field, the completeness was calculated as a function of magnitude for annuli about the centre . In the central region, within a radius of 29 pixels (0.21% the recovery rate was greater than 90% to This magnitude was adopted as the cutoff for the SDP. A similar result was achieved for a pair of images of the same field with exposure times of 5s. On the longer exposures taken for fields Swope:na, Swope:sa and Swope:sb the brighter stars are saturated and could not be counted in the luminosity function. Therefore, the single short exposure images (averages of the two 5s exposures obtained for each field) were used to produce the surface density profiles. In the Swope: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 length exposure. For the Swope:c data the detected stars were counted from the output of all 10 artificial star runs and an average taken. The variation between star lists was quite small, confirming the results of the artificial star tests. The stars were counted in two sets of overlapping annuli to reduce the effects of binning. The outer 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 these annuli with respect to the fields are shown in Figure 13. The dashed and solid circles are the boundaries of the two sets of annuli. The areas of the intersections of the annuli with the rectangular fields (I will refer to these as "sections") were calculated and the star counts were converted to surface densities. The adopted effective radius for each section is the radius which divides the area of that section in half. This radius will vary between sections on the same annulus, but from different fields, depending on the geometry of the section. The field star density was measured on Swope:bf. Twenty-three stars brighter than 1=14 were found on that frame giving a surface density of 0.69±0.14 stars/square arc minute. The resulting surface density profile, with the background subtracted, is shown in Figure 14. The different symbols designate which series of annuli a particular density comes from. The errors are Poisson errors except in the case of field Swope:c 64 I —5^0^5 X wrt centre (arc minutes) Fig. 13.—Map showing the annuli used in constructing the surface density profile. The rectangles are 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 at 0.17' (group A) and having the same width. 65 0 background —0.5^0^0.5 log r (arc minutes) i Fig. 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 radius which bisects the area of the section. The background density has been subtracted and is shown by horizontal line. The points with large error are those for small sections at the edge of a field with very few stars. 66 where a small contribution due to variations in the number recovered on the artificial star frames has been included. In order to account for areas in which no stars were seen and eliminate the apparent anomalies associated with small number statistics, any section of a field which contained fewer than 10 stars was merged with its neighbouring section. The area and effective radius for the new, combined section were calculated, and the density and error re-evaluated. The SDP after this merging is shown in Fig. 15 and listed in Table 4. The various symbols indicate the frame of origin for each data point. A set of surface density profiles were produced in the same way for a set of centres offset, by varying amounts, from that adopted. Little difference was seen for centres offset by up to ±20 pixels (8.7") in either or both directions. The outer SDP was little affected by all these variations in centre position and I conclude that the errors in the determination of the relative positions of the outer 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 of sectored annuli was placed on the plate and the stars were counted by sector. The density and uncertainty 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. In order to check the consistency of the star counts and the assigned errors I split the lists of recovered stars for Swope:c into quadrants about the cluster centre and repeated the counting procedure for each quadrant separately. Figure 16 shows the results of this test. The symbols without error bars are the densities for the four quarters of field Swope:c. The line connects the mean of the four measurements in each annulus and the error bars are the standard deviations of the samples. The solid circles offset +0.025 dex in radius are the overall densities from the full field together with the uncertainties computed from Poisson statistics and frame-to-frame variations in the counts. With the exception of the innermost point, the quarter-to-quarter errors are within a factor of two of the Poisson errors. There is excellent agreement in the errors for mean points between 0.25' and but a larger quadrant-to-quadrant variation is seen for the annuli beyond 1'. 67 The large scatter in the innermost point is an artifact of the linear arrangement of the bright stars in the core of NGC 6397 (q.v. Auriere et al. 1990). Table 4: Surface Density Profile for/ <14 Field^Annulus a reff N (TN (arc min. - ) Group (arc min.) 99. 28. Swope:c^B 0.149 21. 89. Swope:c^A 0.187 17. 53. Swope:c^B 0.279 11. Swope:c^A 40. 0.350 35.3 Swope:c^B 8.3 0.442 32.7 6.3 Swope:c^A 0.555 23.7 7.7 Swope:na^B 0.700 4.7 Swope:c^B 0.700 27.9 Swope:na^A 26.3 5.6 0.847 3.1 19.0 Swope:c^A 0.880 2.2 Swope:c^B 1.109 15.6 4.5 Swope:na^B 1.119 23.8 11.4 1.394 1.5 Swope:c^A 1.403 12.7 2.5 Swope:na^A 1.758 1.0 Swope:c^B 6.8 1.4 6.2 Swope:na^B 1.766 1.866 6.3 2.2 Swope:sa^B 1.4 2.198 6.0 Swope:sa^A 2.210 4.76 0.66 Swope:c^A 2.218 Swope:na^A 4.70 0.98 2.737 2.97 0.46 Swope:c^B 2.749 3.45 0.70 Swope:na^B 2.782 2.77 0.81 Swope:sa^B 2.996 2.13 0.53 Swope:c^A 2.26 3.373 0.54 Swope:na^A 2.34 3.498 0.61 Swope:sa^A 4.275 Swope:na^B 1.55 0.45 4.345 2.49 Swope:sa^B 0.52 1.41 0.42 5.297 Swope:na^A 5.737 0.78 0.43 Swope:sb^A 0.92 5.786 0.31 Swope:sa^A 5.857 1.45 0.57 Swope:na^B 0.32 6.379 0.40 Swope:sa^B 6.867 0.31 0.83 Swope:sb^B 8.475 0.46 0.26 Swope:sb^A 0.28 9.350 0.09 Swope:sb^B Swope:w^A 9.442 0.18 0.30 0.21 11.212 0.03 Swope:w^B 0.21 12.020 0.15 Swope:nb A & B 0.23 Swope:w^A 12.296 -0.06 a Group "A" is the set of annuli starting at 0.17' with width 0.2 dex. Group "B" annuli start at 0.21' and have the same logarithmic width. Section merging as described in the text has been applied. 68 background • Swope:c • Swope:na ■Swope:sa OSwope:sb A Swope:nb * Swope:w I^IiIii^i^I^I^I —1^—0.5^0^0.5 log r (arc minutes) Fig. 15.—The surface densities shown here are for the same sections as in Fig. 14 except that all sections which contained fewer than ten stars (including those with no stars) have been merged with their neighbours on the same field. As the relatively empty sections are always at the edge of a field there is no ambiguity in the direction of merging. The effective radii are those that bisect the areas of the (possibly composite) sections. The various symbols indicate the field of origin for each surface density measurement and the magnitude limit is 1=14. 69 2.5 A ' "a; - 5 z 0 15 cr a) b CD O 0.5 I —0.5^0 log r (arc minutes) Fig. 16.—A comparison of the central surface density profile for R14 derived from the four quadrants of field Swope:c with that derived from the entire field. The various symbols without error bars are the surface densities derived from the star counts in the four quadrants. The line connects the mean densities at each radius and the solid circles are the surface densities measured from the entire field and are offset outward by 0.025 dex for clarity. The error bars on the mean counts are the standard deviation of the sample while the error bars on the counts from the full field are the sum of Poisson and frame to frame errors as described in the text. 70 A 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 to incompleteness. Beyond this point the counts from Swope:c are more than 90% complete between 1=15.0 and /=15.5, and greater than 95% complete for all stars brighter that /=15.5. These incompletenesses would change the final surface density profile by less than 0.05 dex, well within the estimated errors. The magnitude cutoff of 1=15.5 was chosen since this is the magnitude of the turnoff on a roughly calibrated I —(V—/) colour-magnitude diagram produced with data taken at the same time as that of FRST. A field density of 2.06 ± 0.25 stars/sq. arc minute was derived from the Swope:bf data. Since all the stars on both the 1<14 and R15.5 SDPs lie at or above the turnoff they should have the same mass (or did have until recently) and so should have the same radial density distribution. This SDP is listed in Table 5. A comparison between the two SDPs is shown in Fig. 17 where an offset of 0.49 dex has been applied to the 1<14 SDP. The background level for the 1<14 SDP has been shifted by the same amount. This combined SDP has been used for all subsequent 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 intermediate section (0..5 log r 50.6) has a slope of —1.7±0.1. The slopes are the mean of those calculated by a weighted linear regression for the two sets of annuli independently. Over the region 0.5 log r 50.6 the 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 in producing Fig. 17. The fits are shown in Fig. 18. The dot-dash line is a fit to the entire region with 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 cluster is a difficult proposition. For example, Da Costa (1979), from a single-mass King model fit, finds a limiting radius of 38' and Meylan & Mayor (1991) using multi-mass King models find a mean limiting radius of 95'. These values are much larger than the limit of the new SDP (-13'). The SDP steepens at the limit of the observations, but the densities are quite uncertain. 71 Table 5: Surface Density Profile for / <15.5 rdr (arc min.) N (arc min. -2) oN Field^Annulus a Group 15. 119. Swope:c^B 0.442 12. Swope:c^A 108. 0.555 128. 28. Swope:na^A 0.583 Swope:na^B 17. 113. 0.700 Swope:c^B 92.0 8.5 0.700 70.1 Swope:c^A 6.0 0.880 12. 96. Swope:na^A 0.891 3.9 48.8 Swope:c^B 1.109 70.5 7.7 Swope:na^B 1.119 36.1 2.7 Swope:c^A 1.394 43.5 4.7 Swope:na^A 1.403 11. 34. Swope:sa^A 1.581 23.6 1.8 Swope:c^B 1.758 25.4 2.8 Swope:na^B 1.766 4.2 22.9 Swope:sa^B 1.866 15.7 1.2 Swope:c^A 2.210 17.1 1.8 2.218 Swope:na^A 2.3 Swope:sa^A 2.236 15.9 11.6 Swope:c^B 2.665 0.9 13.3 1.4 2.749 Swope:na^B 1.4 2.782 9.0 Swope:sa^B 1.0 Swope:c^A 2.996 8.7 1.0 Swope:na^A 8.4 3.373 1.1 7.6 Swope:sa^A 3.498 6.1 Swope:c^B 1.9 3.544 4.16 4.275 Swope:na^B 0.75 4.345 6.29 0.85 Swope:sa^B 1.6 4.815 4.9 Swope:sb^B 3.36 0.67 Swope:na^A 5.297 Swope:sa^A 3.99 0.67 5.368 3.82 0.85 Swope:sb^A 5.737 3.50 0.92 Swope:na^B 5.857 Swope:sa^B 2.15 0.61 6.379 Swope:sa^A 0.21 0.70 7.072 2.36 Swope:sb^B 0.56 7.085 1.33 0.45 Swope:sb^A 8.475 0.54 Swope:sb^B 0.88 9.350 9.442 0.79 0.54 Swope:w^A 1.15 0.69 Swope:nb^A 9.860 0.56 0.39 Swope:w^B 10.809 Swope:nb^B 0.48 0.41 11.234 Swope:w^A 0.42 0.35 12.296 Swope:nb^A 12.678 0.09 0.38 Swope:w^B 0.28 0.75 13.670 Swope:nb^B 0.53 14.058 -0.02 a Group "A" is the set of annuli starting at 0.17' with width 0.2 dex. Group "B" annuli start at 0.21' and have the same logarithmic width. Section merging as described in the text has been applied. 72 Fig. 17.—A comparison of the surface density profiles for all stars brighter than 1=14 (open circles) 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 background levels are: (dashed) 1<15.5, (dotted) k14. The background level for the 1<14 SDP has also been shifted by +0.49 dex. 73 — 0 cd a) 0 b OD 0 1^1^1^1^1^I^1^1^1^1^I^1 —0.5^0^0.5 log r (arc minutes) 1^ Fig. 18.—The various lines show weighted least-squares, power-law fits to various radial ranges 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 the two radial ranges with the lower limit at log r=0 include points from both SDPs. 74 Section c. Mass functions The luminosity functions for fields du Pont:if, du Pont:n, du Pont:w, and du Pont:bf were produced in accordance with the procedures outlined in §4a. The stars were counted in two sets of overlapping 0.5 mag bins offset by 0.25 magnitudes in order to reduce the distortions due to binning. For each field the true observed area, exclusive of the ignored regions, was computed and the star counts adjusted to be the number of stars which would be observed in a 4.67 square arc minute 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 counts cannot be less than the raw counts. The inclusion of the effects of bin jumping can lead to the corrected 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 which belong in bin i+1 are larger than those in bin i. Even if the scattering is symmetric about the true magnitude, 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 of stars in that bin. This result is even more likely when it is noted that stars have a tendency to be found 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 for incompleteness 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, are tabulated in Table 6. In the bin centred at 1=22.5, 50% of the artificial stars were recovered with magnitudes within 0 7 mag of their true values. For the next independent 0.5 mag bin, centred at 1=23.0, only 8% of the artificial stars were recovered. Since it is impossible, because of the poor statistics, to estimate the number of stars in the 1=22.5 bin which originated in the 1=23.0 bin, the 1=22.5 bin could not be used in the luminosity function. Similarly, for the overlapping bins centred on the quarter magnitude, the last useful bin is at 1=21.75. Although the next bin has a recovery rate of 78%, the one following it, centred at 1=22.75, has a recovery rate of only 26%, and of 75 these, 44% were found in the next higher bin, 12% in the the next fainter bin, and the remaining 44% were found in the same magnitude bin they were added in. The recovery rate drops quite quickly, so the last usable bin usually has a rather high recovery fraction. As a result of these considerations, the cutoff magnitude for the du Pont:bf luminosity function is 1=22.25; the last bin used being the one centred at I=22.0. Table 6: Raw and completeness corrected luminosity function for the background field du Pont:bf. I error error Ncou nte da Ncomplete a 1.1 1.1 15.00 1.1 1.1 2.1 2.1 15.25 4.3 4.3 2.1 4.3 2.1 15.50 4.3 3.2 3.2 15.75 1.8 1.8 16.00 5.3 2.5 5.3 2.5 16.25 2.6 2.6 6.4 6.4 16.50 3.2 9.6 9.6 3.2 11.7 11.7 16.75 3.6 3.6 17.00 12.8 3.7 12.8 3.7 17.25 10.6 3.4 10.6 3.4 17.50 15.9 4.1 15.9 4.1 17.75 18.1 4.4 18.1 4.4 18.00 13.8 4.0 13.8 4.0 18.25 21.3 4.8 4.8 21.3 22.3 22.3 18.50 4.9 4.9 18.75 31.0 5.9 32.6 6.4 19.00 33.7 6.1 35.0 6.4 19.25 32.5 6.0 33.3 6.2 19.50 37.3 6.4 38.1 6.6 19.75 34.5 6.1 34.5 6.1 20.00 30.0 5.8 34.9 7.0 7.3 20.25 48.4 46.7 8.0 20.50 65.2 8.5 63. 10. 20.75 61.8 11. 8.3 69. 12. 21.00 73.9 9.1 76. 14. 21.25 83.3 9.6 85. 21.50 9.2 75. 14. 75.8 21.75 16. 78.0 9.3 72. 21. 22.00 84.3 9.8 103. 'Number per 0.5 magnitudes per 4.67 sq. arc minute field. 76 The du Pont:if starcounts were found to have a 50% recovery rate at 1=21.7 and when the bin jumping effects are taken into account the limiting magnitude is /.21.5. The 50% recovery rate magnitude for the du Pont:n data is 1=22.8, and for the du Pont:w data it is /=22.1. The magnitude limits of the incompleteness corrected star counts are 22.5 and 21.75 respectively. The du Pont:n LF goes 0.75 mag deeper than the du Pont:w LF because of better seeing and a longer total exposure. 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 overlap and a final, outer LF produced. The bins with centres with 121.75 are from the du Pont:n data only. The magnitude limit of the background LF cuts off this combined LF, which will be referred to as the du Pont:out LF, at /=22.25. Tables 7-9 contain the raw and incompleteness corrected luminosity functions for the three program fields. Figure 19 shows a comparison of the luminosity functions for the background fields observed in this work and in FRST. The agreement is good despite the two fields being on opposite sides of the cluster. This confirms the observation by Da Costa (1982) that the field stars are evenly distributed in the vicinity of the cluster. The final, background subtracted, luminosity functions at the two distances are listed in Table 10 and are plotted in Figs. 20 and 21. These two figures show the two new LFs and compares them with the background star counts. At the distance of the two outer fields the number of field stars is equal to the number of cluster stars, so deriving a more precise 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 LF gives three deep LFs for NGC 6397. Figure 22 brings together all three of the observed deep LFs for NGC 6397. The two types of filled symbols are the new LFs from this study, while the open symbols are the deeper LF of FRST. No further normalizations have been applied; the offsets between the LFs are a reflection of the overall decrease in density with radius. For the magnitude range 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 more conservative approach to incompleteness corrections used here. The origin of the dip around /=19 77 in the du Pont:out LF is unclear, but is seen, at about the same magnitude, in both fields. The general depression in the counts between 1-18.5 and 1-20 suggests that whatever is causing the dip is not restricted to just the two low bins. Table 7: Raw and completeness corrected luminosity function for inner field du Pont:if. a Ncounted error N^, error comp.etea 1.1 1.1 1.1 1.1 13.75 14.00 1.1 1.1 1.1 1.1 14.25 5.2 2.5 5.2 2.5 14.50 7.3 7.3 2.9 2.9 14.75 5.2 2.5 5.2 2.5 15.00 7.5 3.0 7.5 3.0 9.4 15.25 3.1 9.4 3.1 4.5 15.50 19.0 4.5 19.0 15.75 30.5 5.8 30.5 5.8 16.00 40.0 6.5 40.0 6.5 16.25 53.4 7.6 53.4 7.6 16.50 67.9 8.6 67.9 8.6 16.75 79.1 9.5 79.1 9.5 17.00 11. 103. 103. 11. 17.25 118. 12. 118. 12. 110. 17.50 11. 113. 12. 110. 11. 17.75 114. 12. 128. 18.00 12. 115. 13. 18.25 133. 12. 122. 14. 142. 18.50 13. 153. 17. 14. 18.75 153. 163. 18. 19.00 14. 168. 183. 19. 19.25 171. 14. 193. 20. 181. 19.50 15. 197. 19. 19.75 214. 17. 217. 25. 20.00 232. 17. 237. 29. 20.25 241. 18. 293. 33. 241. 20.50 17. 42. 307. 234. 20.75 314. 17. 40. 21.00 222. 17. 341. 55. 21.25 192. 16. 322. 53. a Number per 0.5 magnitudes per 4.67 sq. arc minute field. 78 Table 8: Raw and completeness corrected luminosity function for outer field du Pont:n. a a error I Ncomplete^error Ncounted 4.1 4.1 2.0 2.0 15.00 8.2 3.0 3.0 15.25 8.2 4.1 2.0 2.0 4.1 15.50 1.4 1.4 2.0 2.0 15.75 2.4 2.4 5.1 5.1 16.00 11.1 3.5 11.1 3.5 16.25 18.1 18.1 4.3 4.3 16.50 4.5 19.7 4.5 19.7 16.75 22.9 22.9 4.9 4.9 17.00 22.6 4.8 4.8 22.6 17.25 5.4 28.7 5.4 28.7 17.50 30.3 5.7 30.3 5.7 17.75 34.3 6.0 34.3 6.0 18.00 48.1 48.1 7.2 7.2 18.25 6.3 6.3 37.5 18.50 37.5 6.2 36.8 6.3 36.5 18.75 7.3 49.9 7.2 50.3 19.00 62.5 8.6 8.3 19.25 66.8 8.7 66.5 8.3 69.0 19.50 85. 10. 77.2 9.0 19.75 10. 100. 11. 104. 20.00 11. 99. 11. 103. 20.25 13. 106. 10. 20.50 99. 119. 11. 108. 13. 20.75 12. 112. 15. 110. 21.00 11. 14. 104. 21.25 106. 11. 115. 125. 17. 21.50 21. 118. 12. 125. 21.75 12. 24. 129. 133. 22.00 12. 28. 22.25 132. 166. 'Number per 0.5 magnitudes per 4.67 sq. arc minute field. Table 9: Raw and completeness corrected luminosity function for the outer field du Pont:w. a N error error Ncompletea -counted 1.1 1.1 1.1 1.1 15.25 2.1 4.2 2.1 4.2 15.50 5.3 2.5 2.5 15.75 5.3 2.5 2.5 5.3 16.00 5.3 3.5 11.3 11.3 3.5 16.25 4.2 4.2 16.8 16.8 16.50 4.1 4.1 16.0 16.75 16.0 20.1 4.6 20.1 4.6 17.00 20.0 4.6 4.6 17.25 20.0 5.2 25.2 5.2 25.2 17.50 40.0 6.5 40.0 6.5 17.75 41.2 6.7 40.4 6.6 18.00 5.8 31.8 6.1 30.8 18.25 6.3 34.5 6.1 35.5 18.50 43.5 6.9 6.8 42.8 18.75 6.5 35.5 7.1 39.2 19.00 44.8 45.9 7.2 19.25 7.0 9.4 58.5 8.0 60.6 19.50 72. 11. 72.9 9.0 19.75 82. 12. 78.6 9.2 20.00 81. 13. 20.25 80.5 9.3 16. 90.4 9.9 88. 20.50 11. 113. 15. 20.75 105. 120. 22. 12. 21.00 153. 129. 12. 152. 21. 21.25 118. 12. 125. 26. 21.50 a Number per 0.5 magnitudes per 4.67 sq. arc minute field. 2 0.5 16 18 I 20 22 Fig. 19.—Comparison of the background field star counts for FRST (histogram) and du Pont:bf (points). The two fields are at a distance of 1° on either side of the cluster and at the same galactic latititude. The uncertainties in the FRST counts should be similar to those shown for the du Pont:bf counts. All of the new luminosity functions are shown for two sets of overlapping half-magnitude bins. 81 Table 10: Background subtracted luminosity functions for du Pont:if (inner) and average of fields du Pont:n and du Pont:w (outer). I error error Ninne: Noutera 2.0 2.2 15.25 5.2 3.8 2.1 14.7 5.0 0.0 15.50 27.3 6.1 1.0 1.9 15.75 34.6 2.5 7.0 0.0 16.00 4.8 47.0 8.0 3.1 16.25 9.2 7.9 3.8 58.3 16.50 67. 10. 6.2 4.0 16.75 11. 90. 8.7 4.3 17.00 12. 108. 10.7 4.1 17.25 11.0 13. 4.8 17.50 97. 12. 17.1 5.3 96. 17.75 102. 14. 23.9 18.00 5.3 101. 15. 18.7 5.8 18.25 18. 14.2 130. 5.6 18.50 131. 19. 7.60 6.5 18.75 148. 20. 7.90 6.8 19.00 21. 7.1 19.25 160. 20.9 21. 7.9 159. 26.7 19.50 26. 44.2 183. 8.6 19.75 202. 20.00 30. 56.5 9.4 34. 43. 246. 10. 20.25 35. 244. 43. 12. 20.50 246. 41. 47. 20.75 13. 57. 56. 16. 21.00 265. 55. 21.25 236. 43. 16. 21.50 50. 19. 53. 21.75 26. 22.00 30. 32. 'Number per 0.5 magnitudes per 4.67 sq. arc minute field. i^1^i I I f 1^1^1 1 i i I i ffi1ff f f f f f f^.--,^. .^background field .___.. .__:---- 2^------ - . 0 I^,^I 22 20^ 16^18^ I I^I Fig. 20.—Luminosity function for the inner field du Pont:if. The background counts are shown for comparison. 83 Fig. 21.—Luminosity function for the outer fields du Pont:n and du Pont:w. The counts have been averaged to improve the signal-to-noise. The background counts are also shown. These fields are about as far from the centre of the cluster as it is practical to go. 84 If i if f if I^I ^ I^I^I 16^18^20 I 22 ^ 24 Fig. 22.—All the NGC 6397 LFs are compared. The three LFs are for the inner field (filled circles), FRST (open circles) and the averaged outer fields (squares). For the two LFs from this study, two sets of overlapping half-magnitude bins are shown. 85 In order to convert the luminosity functions into mass functions a mass-luminosity law is required. The I—band mass-luminosity law of FRST, with an extension to more massive stars taken 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 and du Pont:w (outer). error Nma m (inner) (M0) 860 290 0.804 3.50 280 1240 0.794 3.75 1310 270 4.00 0.782 250 1500 4.25 0.768 250 1560 0.751 4.50 240 1570 4.75 0.731 230 1890 0.708 5.00 240 2160 0.684 5.25 240 1830 0.657 5.50 1740 230 0.631 5.75 1920 260 0.605 6.00 270 1820 6.25 0.576 310 2260 6.50 0.549 310 0.517 2070 6.75 280 0.484 2040 7.00 2170 290 7.25 0.449 2100 270 0.410 7.50 2250 320 0.371 7.75 370 8.00 0.332 2510 3310 460 8.25 0.293 690 8.50 0.260 3950 8.75 5150 870 0.232 1400 9.00 0.209 6600 0.192 1600 9.25 7000 0.176 9.50 9.75 0.164 10.00 0.154 a Number per unit mass per 4.67 sq arc minute field. Mr 86 Nma (outer) 47 0 154 210 143 183 214 207 309 450 340 246 120 110 285 350 540 700 580 560 990 1400 1260 1890 2400 1600 error 88 94 98 100 93 89 82 89 96 100 100 98 100 94 97 100 110 120 140 200 270 390 470 700 1200 1700 -0 cd zt!. 3 cCf C/) N 10 . 2 J^ I —0.2 —0.4^—0.6 log m (M 0 ) —0.8 Fig. 23.—All the NGC 6397 MFs are compared. The symbols and bins are as in Fig. 22. 87 From the more massive, upper main-sequence, stars it is difficult to see any sign of mass segregation, although there is some indication that the outer field is more deficient in higher mass stars than the two fields closer to the centre. For the less massive stars, there is certainly an impression that the outer field contains relatively higher proportions of these than do the two inner fields. For these stars there is a somewhat clearer signal for mass segregation. Table 12 gives the slopes 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 in the du Pont:out MF has been excluded from this fitting due to its large uncertainty. For the du Pont:if and FRST MFs the three slopes are quite consistent with one another, and are also consistent with there being little mass segregation between the two fields. As the high-mass cutoff is moved to lower masses the MF from the outer fields gets steeper. The evidence for mass segregation is strongest for the stars with M<0.32 M M. but is not highly significant. There is also an indication in the three MFs that the mass of the break in the MSI decreases with distance from the cluster centre. For the du Pont:if MF the upturn in the MF takes place at around 0.4M 0 , for the FRST MF this occurs between 0.4M 0 and 0.3M 0 , and for the du Pont:out MF the data suggests an 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 indices m<0.4M 0^m<0.32M0 m<0.28M 0 Mass Function limiting MSI #^MSI #^MSI # mass (M 0 ) du Pont:if 0.8±0.3 7 0.9±0.5 5 1.0±0.9 4 0.19 FRST 1.0±0.2 8 1.1±0.3 7 0.9±0.4 6 0.12 du Pont:out 0.5±0.3 9 1.4±0.6 7 1.8±1.0 6 0.16 Constructing the segregation measures S„ as defined in eq. (3.1), might aid in quantifying the degree of mass segregation by focusing on the relative differences in the the mass functions. The segregation measures between the three pairs of fields are shown in Fig. 24. One difficulty with these diagrams are the large errors associated with taking the ratio of two, already somewhat 88 poorly determined, quantities. To a certain extent the conclusions about mass segregation made from the MFs themselves are supported by the segregation measures, but most of the scatter in the segregation measures is probably due to observational uncertainties and only weak conclusions can be drawn from them. 89 11^111 1 111111V111111111111 1,11),,,,i1 ^ 1^1 S,(m,du Pont:out; du Pont:if)_. S r (m. du Pont:out; FRST)^Sr(rn. FRST; du Pont:in) 11,,,I,III,,,,,III lllll 'i l l^.. —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.8 log m (M e)^log m (M G)^log m (M 0 ) Fig. 24.—The segregation measures between the three pairs of fields are shown. The segregation measure, Sr, is defined by eq. (3.1). 90 Chapter 5. Models Including Stellar Evolution The focus of this chapter will be on attempting to find Fokker-Planck models, employing all the features discussed in ch. 2, which can reproduce the observations of M71 and NGC 6397. No attempt will be made at surveying the set of all possible parameters. Rather, a set of models with plausible initial conditions, and which survive to an age of approximately 15 Gyr, are to be compared with the observations. The properties listed in Table A4 for the first grid of models seem quite reasonable, as will be discussed presently. When actually compared with the observations of M71 and NGC 6397, they quickly prove to be wide of the mark. The lack of remotely satisfactory 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 some of the motivation behind the initial choice of parameters in order to give a sense of the issues involved. An initial mass of order 10 6 M o is suggested by the Jeans mass of primordial gas at 10 4 K, the starting 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 10 5 M 0 . The young populous clusters in the LMC have masses estimated to be 10 4 to 10 6 M 0 (Elson, Fall, & Freeman 1987). Holtzman et al. (1992), using the Hubble Space Telescope, recently detected a population of luminous blue objects in the galaxy NGC 1275. A new set of observations by Richer et al. (1992) confirms the existence of these objects. The most probable interpretation of these observations is that the objects are young globular clusters. Richer et al. estimate the mean cluster mass to be around 10 6 M o . Allowing for inefficiencies in star formation and the amount of mass lost via stellar evolution and the tidal boundary, 10 6 M o is a reasonable initial mass for clusters such as M71 and NGC 6397. Observations of the young populous clusters in the LMC are also providing information on the IMF 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 are consistent with a mean mass spectral index <x>=1.7±0.2 for stars more massive than about 1.5 M o . Furthermore, the same mean slope, with a somewhat larger spread, appears to fit the 91 MFs in many diverse environments. Based on this observation, the initial high mass IMF had a MSI of x h =1.7. A flatter IMF, which ultimately was required for M71 and NGC 6397, is supported by the observations of Elson, Fall, & Freeman (1989). In six young LMC clusters, they found that between 1.5 and 6 M o 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 the MSI for masses below about 0.4 M o (Richer et al. 1991). These show a variety of slopes which have been affected by mass segregation and tidal stripping. The steepest of these have a MSI of x=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 are appropriate to all clusters, the initial MSI for the stars with M<0.4M 0 was taken to be x 1=3.5. The observed IMFs for stars from the break up to the turnoff have smaller mass spectral indices. For convenience the high mass IMF is continued to 0.4 M o . Assuming that the transitions in the MFs collected in Richer et al. (1991) are not due to a problem with the mass-luminosity relation, the MFs do not have the appearance of being the sum of two overlapping power-laws, but that a much sharper transition takes place. Hence the IMF has been constructed to be continuous but not smooth at the mass bin with initial mass 0.44 M o (q.v. Table A5), with one power-law on either side of the break. The mass range has been taken to be between 20 M o and 0.1 M 0 , and 25 mass bins were used. The mass limits are somewhat arbitrary, but are governed by the observed range of masses. Constructing the IMF with an upper limit of 100 M o and x n =1.7 gives an additional 0.4% of mass loss in the first 3 to 6 Myr, and 10% more neutron stars. Extending the low mass cutoff, but with the same total mass, will reduce the fraction of mass lost through stellar evolution and will lengthen the relaxation time, but these are not large effects. The final masses and stellar lifetimes used are the 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. 92 The analysis of the models were done in the following manner. In view of the uncertainty in the relaxation timescale, any model which survived past 10 Gyr was examined for comparison with the two clusters. All the saved steps between 10 Gyr and 20 Gyr were used for this. The models 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 was used for the SDP for stars between 0.85 and 0.89 M 0 , and the next bin, with mass 0.78 M 0 was used for the SDP from the stars between 0.71 and 0.85 M o . I will refer to the bin with mass 0.87 M o as the "giant bin" as it contains the stars currently evolving off the main sequence and which 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 rectangular areas at the locations of the du Pont:if, du Pont:n, du Pont:w and FRST fields (Table 3) and with areas 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 for NGC 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 the mass-luminosity relationship used for NGC 6397 is only 0.83 M o as opposed to 0.87 M 0 used in the models. For the purposes of comparison the giant bin has been used to model the shape of the SDP 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 model and cluster and included all the saved time steps between 10 and 20 Gyr. Simplified versions of this diagram are shown below. The full diagram was examined by eye and the time which appeared to best fit the observations was selected. With, the exception of the NGC 6397 SDP, where a density rescaling was permitted because of the mismatch in the bin mass and width and the uncertainty in the fraction of the turnoff and giant stars included in the SDP, there are no free scaling parameters in the models to actually "fit". When the term "fit" is used in what follows, only this sort of visual matching is implied; no other adjustments are possible. 93 Section a. M71 Figure 25 shows the most reasonable of the models for M71 from the series of Table A4. The specific parameters for this model, designated model M71—A, are listed in Table A5. The degree of mass segregation in the model is in satisfactory agreement with the observations. The SDPs contain too many stars, indicating too high an initial mass, and have a peculiar shape, but the inner region has a shape similar to that of the observed SDPs. The model MFs are far too steep to match the 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 MF pre-evolved to its value at 15 Gyr gave very similar results. If M71 is to be fit with a stellar evolution 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 fit to M71. Some of the choices in parameters to try were driven by the desire to fit NGC 6397, but as the observed MFs and galactic positions of the two clusters are similar, such choices will serve equally well for M71. In addition to the initial concentration and radial size, the mass spectral indices and the initial mass were also varied. Besides 10 6 M 0 , initial masses of 7.5x10 5 M, and 5x10 5 Mo were used. For the high-mass stars, mass spectral indices, x h , 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 the models was with respect to the IMF. Table B1 gives details of the IMFs used and the sets of initial parameters which gave models surviving to 10 Gyr or more. The initial high-mass MSI for M71 was probably not as flat as x h =0.2. All of the models with this value of x h were destroyed fairly quickly despite having an x i=3.0. For this IMF, 59% of the total initial mass is lost through stellar evolution. The models with x h =0.5 survived only if they were initially sufficiently concentrated: models with W 0 5_4 were destroyed, models with W o =6 or higher survived. Tidal radii (as defined for a cluster with mass 10 5 Mo [q.v. eq. 2.21]), near 20 pc appear to give the best results both for M71 and NGC 6397. As an aid to identifying models with properties similar to those observed in M71, I calculated for each model the normalized segregation measure S *f7 1 (eq. 3.3) and the radius where the density 94 • !mil), 1 1^111 1 _ 0.71<m/M o <0.85 1^iirr^11111111^111 0.85<M/M o <0.89 - I^1^I^I^I^1^I^I^I^I^I^I^I^I^I^I^I —0.5^0^0.5 log r (pc) —0.5^0^0.5 log r (pc) Ca 1 —0.8 E6 • 5 a) E 4 Z 3 —1.4 —1 —0.8-0.6-0.4-0.2 0 log m (M 0 ) —1 —0.8-0.6-0.4-0.2 0 log m (M 0 ) Fig. 25.—Comparison of model M71—A with the observations of M71. The top panels show the two observed SDPs together with the model SDPs. The lower-left panel shows the observed and model MFs for the two fields (filled circles, core; open circles, 3'). The lower-right panel shows the observed and model segregation measures. 95 of stars in the giant bin reaches half its central value. This is being used as a stand-in for the classical '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 name for King's scaling radius (eq. 2.12), so I will refer to the radius used here as the half-density radius, 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 models computed. Also indicated is the observed values of Srm,h71 =0.4 and rhd=1. It appears from this figure that none of the models can be expected to provide a good match with M71, and this is borne out by more detailed comparisons. The systematics of S,,,h—rhd diagram are discussed more fully 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 the observed amount of mass segregation between the two fields. The breakdown of Fig. 26 by IMF shown in Fig. B 1, shows that the only IMFs producing a large value of .5?. "inni are those with x h =1.7. It is one of these models which are shown in Fig. 25, and the required reduction in the MSI for the upper main sequence stars does not occur. M71 vs. rhd curves for a The second conclusion can be seen in Fig. B2 which shows the entire Sr, series of models. During the epoch when stellar evolution is important the models evolve to larger rhd at almost constant SrMn . Stellar evolution does not provide a way for mass segregation to take place before extensive core contraction has occurred. Model M71—B (Table A6), shown in Fig. 27, is typical of the models with the observed rhd and a more appropriate MF. This model has too high a mass to fit M71 and the MFs are still too steep. The segregation measure is only half that observed. It should be possible to find a model with lower mass and flatter IMF to fit the SDPs, but both of these changes tend to decrease the relaxation time and bring core collapse sooner. Having a higher proportion of high mass stars causes 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 smaller than the tidal radius at the desired galactocentric distance. This will reduce the rate of mass loss through the tidal boundary and slow the evolution of the model. A model run with the type of 96 0.6 0.4 U) 0.2 0 —3 ^ —2 ^ ^ ^ —1 0 1 log rhd Fig. 26^S,.. rhd diagram for M71. The model segregation measure has been calculated for the same fields as were observed. The diagram depicts the evolution of the segregation measure and half-density radii with time. The evolutionary tracks are plotted for all models which have not disrupted before 10 Gyr, and extend to 20 Gyr. Evolution begins at zero segregation measure and goes in the direction of increasing S r and decreasing rm. Evoldtion at constant S r . occurs during core 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. The track for model M71-A lies to the left of the M71 point. — 97 _IIII 0.71 <m/M e <0.85 0.85<m/M0<0.89 - I^I^t^ I^1^I^I^l^I^I^I^I^I^I^I^1^I^I —0.5^0^0.5 log r (pc) —0.5^0^0.5 log r (pc) 1 1 —0.8 E6 4) 5 a) 4 3 —1.4 —1 —0.8-0.6-0.4-0.2 0 log m (M 0 ) —1 —0.8-0.6-0.4-0.2 0 log m (M e ) 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 the 3' field are not too different, but the core MF slopes are poorly matched. The S r m vs. rhd curve for this model is among those lying below the point for M71 in Fig. 26. 98 initial parameters just outlined will not have any higher a segregation measure so the problem of M71 still stands. There are still possible solutions which could allow M71 to be matched by a Fokker-Planck model. It is yet to be seen how primordial binaries affect a multi-mass model. Tidal shocking, by the 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, the problem may not lie so much in the Fokker-Planck modelling as in the initial conditions. If globular clusters start mass segregated, then the conclusion derived here for models which start without mass segregation will not apply. In the face of continuing difficulty in finding a model for M71, additional observations of the mass function are required in order to confirm the large degree of mass segregation in the present data. The file is not closed on M71, and the possibility of building a Fokker-Planck model, perhaps with additional physics, still remains. Section b. NGC 6397 One might reasonably expect that it would be easier to find a model matching NGC 6397 than one for M71, and this is the case. As a first step, a S r,,h rhd diagram from the same set of models — has been constructed for NGC 6397 using the normalized segregation measure S,...(0.22M o , du Pont:out; 0.77M 0 , du Pont:if). The observed segregation measure is, from Fig. 24, between 0.3 and 0.4. An upper limit on rhd of 0.2 is based on the mean core radius of the multi-mass King models of Meylan & Mayor (1991). The S ,„,--rhd diagram for the NGC 6397 observational window is shown in Fig. 28 with the estimated observed range delimited by a dashed box. There are plenty of models with properties appropriate to the observations. If the density of models within the box seems 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 in Table A7. The best matching of the times stored for this model is shown in Fig. 29. It should be noted that this is just a representative comparison and that a better match is possible in the interval between 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 been 99 0.8 4.; 0 P-4 0.6 A 0 :is 0.4 a.,0 0.2 C\1 C\2 z O 0 —3 ^ —2^—1 0 log rhd 1 Fig. 28.—As Fig. 26, but for the NGC 6397 normalized segregation measure S„(0.22M 0 , du Pont:out; 0.77M 0 , du Pont:if). The difference between the tracks followed by the models in this figure and in Fig. 26 lies in the different windows through which the models are observed. The outlined box is the range of estimated values of S„ and TM for NGC 6397. Many good candidates for a NGC 6397 model lie within this region. 100 - 1^1 ^ 1^1^ 1^1^1^1^1^1^1^1^1^1^1^1 FRST MF E z 0 1^1^11^1^ 1^1^1^1^1^1^1^1^1^1^1^1 1^ —1 —0.5 0^0.5 log r (pc) I^1^1^1^1^1^I^I, du Pont MFs 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 0 log m (M 0 ) — - r* ;r: 0.4 1^1^1^ 1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1_ —0.6 -o —0.8 -1 — I. 1.2 —1.4 U) _ 1—1 —0.8 —0.6 —0.4 —0.2 0 log m (M 0 ) —1 —0.8 —0.6 —0.4 —0.2 0 log m (M 0 ) Fig. 29.—Comparison of model NGC 6397—A with the observations described in Ch. 4. Composite surface density profile including both the R14 and R15.5 data. The model SDP 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 quite good. Bottom-right: Segregation measure between the du Pont:if and du Pont:out MFs. Again, the model presents a good representation of the observations. Top-left: 101 decreased by a factor of four. With the exception of this normalization, none of the model results have been adjusted in any way to improve the fit. The absolute normalizations of three mass functions are very well reproduced, as are their overall shapes. The good match for the segregation measure confirms the more detailed MF comparisons. Meylan & Mayor (1991) presented a line-of-sight velocity dispersion profile for NGC 6397. In Fig. 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 increasing more rapidly towards the centre than the observed velocity dispersion. If this flattening is supported by further observations, the discrepancy in the model may be explained as being due to the assumption that all neutron stars were retained by the cluster. The result of this assumption is to give the model a higher central velocity dispersion than would be the case if a sizable fraction of the neutron stars were ejected. That the model velocity dispersions lie systematically higher than the observations is not a concern since continuing mass loss in the model will bring them into better agreement. At the preferred time, 16.6 Gyr, model NGC 6397—A is undergoing core collapse and is 40 Myr before the density maximum. The current mass is 10 5 Mo , of which 9% is in neutron stars and 36% in white dwarfs which have a mean mass of 0.8 M o . There are 7,000 neutron stars and 48,000 white dwarfs present out of 2.6x10 5 stars in total. The global mean mass is 0.4 M o and in the centre, dominated by degenerates, the mean mass is 1.3 M 0 . The half-mass radius is 4 pc and the 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 global MF has flattened considerably and this is reflected in the MFs for the three fields individually. The model MF for du Pont:if is flatter than the global MF but the model MFs for FRST and du Pont:out have steeper values of x h than the does the global MF. The values of x h at these two positions are similar to that of the global MF. Mass segregation is seen most strongly for the more massive stars. This is a difficulty with respect to the observations, as the MFs for the high mass stars are more 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^1 log r (pc) Fig. 30.—Comparison of the velocity dispersion profile of model NGC 6397—A with the observed velocity dispersions of Meylan and Mayor (1991). The model reproduces the observations adequately although there is some indication that it is rising more quickly towards the centre than is observed. Differences in the handling of the heavy remnants, for example ejecting some fraction of the neutron stars, will reduce this difference. 103 Table 13. Mass spectral indices for model NGC 6397—A xh Mass Function x, 0.9 IMF 1.5 —0.1 1.0 Current global —0.2 du Pont:if field 0.8 0.3 0.9 FRST 1.1 0.7 du Pont:out field The apparent shift in the mass of the break in the observed MFs is not seen in the model. The evolved MF remains fairly faithful to a power-law form for masses both above and below the break and does not show any movement in the position of the break. It is difficult to see how the shape of the MF can be changed by mass segregation to the extent that a feature like the break can move about in magnitude. There were several other models that gave similarly satisfactory fits to the SDP and the segregation 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. Some leeway 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 comparison with the observations are still unsatisfactory. For r >lpc the observed SDP is much straighter than the model SDP. The detailed structure of the MFs are not reproduced using such a simple form for the IMF and it will be a delicate fine-tuning problem to get a better fit. With respect to the input data in this model, the stellar lifetimes are inappropriate for low metallicity stars, lending some uncertainty to the result. Adjusting the stellar lifetimes will have small quantitative effects but the qualitative behaviour of the models should be the same. In the face of the inability to find a model matching the M71 observations, it is an encouraging sign that such a good match can be found for NGC 6397. The wide range of observations the model is able to reproduce is certainly an indication that such Fokker-Planck modelling does represent an aspect of globular cluster evolution. 104 Chapter 6. Conclusions The aim of this thesis was to test a class of Fokker-Planck model with observations of globular clusters. The two clusters used, M71 and NGC 6397, have different density profiles but have similar masses and mass functions. From the perspective of the Fokker-Planck models, the differences are more important than the similarities. I was unable to find a model which matched the observations of 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-Planck models are appropriate for some classes of globular clusters. The original problem encountered for M71 in the models with pre-evolved MFs (ch. 3) remains when models including stellar evolution are used. The addition of stellar evolution serves to delay core 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 mass loss from stellar evolution are also the models with the lowest segregation measures for the fields observed 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 against one another in attempting to improve the match with a set of globular cluster observations. One of the key problems lies in the realm of time scales. In order to even be tested with the observations I required a model to survive until at least 10 Gyr. Consider a model which is assigned an initial mass, 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 first is the degree of expansion and subsequent delay in the time of core collapse due to mass loss during stellar evolution. The second is the degree to which mass segregation accelerates core contraction. The position in the galaxy governs the size of the tidal boundary and the rate of tidal mass loss. When the initial limiting radius of the model is smaller than the initial radius of the tidal boundary, the tidal mass loss rate will be reduced, leaving the model with more mass at a given time than if the model initially fills its Roche volume. Changing any of these parameters has side effects which can be offset somewhat by varying other parameters. In model M71 —B (Fig. 27) we saw that the model SDPs were much higher than 105 the observed SDPs. It would be desirable, in order to improve the match, to reduce the total mass of the model, say by a factor of two. If the initial limiting radius is constrained to be the same as the initial tidal boundary, then the relaxation time also is reduced by about a factor of two. (This follows from the proportionality t r ec m1/2rh3/2 [eq. 2.11]. The constraint on the limiting radius and the definition of the tidal radius [eq 2.20] requires recM1/3 , hence trhOCM. with this constraint.) For a model like M71-B, but with half the mass, core collapse occurs at 12 Gyr; at 10 Gyr the model is already too concentrated to give a good match. Instead of reducing the total mass, we could just redistribute it to stars with masses below the minimum mass observed. This adjustment in the IMF would also cause there to be less expansion due to stellar evolution, and possibly more mass loss through the titdal boundary. These two effects would shorten the time to core collapse somewhat but the full effects of such adjustments to the IMF have not been explored. Having a larger tidal boundary 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 be partially constrained by the observed MFs, but large ranges, especially the high mass end, are difficult to observe. Some constraint is possible through counts of degenerate stars, eg. millisecond pulsars, but these will be incomplete. Correction of the observed counts for things such as selection effects, rate of formation from the parent population, and the fraction of the remnant population retained by the cluster, will make interpreting these counts difficult. Given the ability to make trade-offs between all these parameters, it should not be expected that the model matching NGC 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 SDP and 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 their models of M15 and NGC 6624, GCLM showed that Fokker-Planck models were capable of fitting the surface brightness profile and velocity dispersion profile of a core-collapsed cluster. The model presented here for NGC 6397 extends the range of data which can be fit simultaneously to include several 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 remarkable 106 in that, with the exception of the SDP offset, no adjustments were allowed in the model output in order to fit (literally, this time) the observations. The model results were transferred onto the observational plane and plotted together with the observations. The mass function results are particularly 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 to M71, 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 main effects will be seen through the relative masses of the turnoff stars and the remnants. It is a shortcoming in the model of NGC 6397 that the mass bin used to model the SDP has a slightly higher mass than the currently evolving stars in the cluster, but this can be rectified with more appropriate 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 the former, will lie in the differences in the initial conditions and in the differing environments. The present mass functions in the two clusters are similar and suggest that the IMFs were similar. It might be that M71 started with a higher proportion of primordial binaries than did NGC 6397. The advantage of primordial binaries over stellar evolution as a heating source is that they promote mass segregation even while preventing an extreme core collapse. The models with primordial binaries of Gao et al. (1992) support this idea. Mass loss due to stellar evolution, acting as it does through the potential, causes a more indiscriminate expansion. The effect of primordial binaries is certainly worth pursuing in multi-mass models. Another possible intrinsic solution would be if the initial conditions used here are not appropriate to globular clusters. For example, a globular cluster may start mass segregated. In such a case the relationship between mass segregation and concentration is not as straightforward as the S r,„—rhd diagram would suggest. As well, if the massive stars start centrally concentrated, then the effects of their losing mass will be amplified. More observations of candidates for young globular clusters, 107 such as the ones in the LMC or NGC 1275, and a theory for globular cluster formation, would help to resolve these issues. Since M71 is a disk cluster, as opposed to NGC 6397 with its halo orbit, it will be subject to the tidal effects of the disk and bulge to a much greater extent. For example, it is more likely to encounter a giant molecular cloud than NGC 6397 is. One might speculate that M71 has suffered from some catastrophic interaction which has led to the rapid expansion of the core. If the expansion occurs over a time much shorter than the relaxation time, then the previous degree of mass segregation will be preserved even while the core radius increases. Alternatively, errors in the observed mass functions could be giving a much higher measurement of 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 physical state. Of the clusters which have been compared with the type of detailed Fokker-Planck model used here, three, M15, NGC 6624, and NGC 6397, have central cusps in their surface densities, while the other, M71, is reasonably well fit by a King model. Certainly M71 does not show the strong central cusp predicted by the models for a core-collapse or post—core-collapse cluster. It is the three cusp clusters which have successfully been matched by Fokker-Planck models but which are not well fit by multi-mass King models (e.g. Meylan & Mayor 1991 for NGC 6397). It is premature to draw any conclusions from such a small sample, but the contrast is striking The King model distribution 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 on an assumption of equipartition of energy between the mass classes at the cluster centre. As pointed out by Pryor et al. (1986), the models don't satisfy this assumption. The low mass stars are out of equipartition with the high mass stars, but the discrepancy is similar to that seen in the Fokker-Planck models. The two types of model are based on the same equation but under different sets of assumptions and each appears to work best with different types of clusters. The static snapshot given by multi-mass King models accords better with the apparently unevolved clusters, but cannot produce models with the high concentrations that the dynamic Fokker-Planck models can. However, 108 it has yet to be shown that King models can also reproduce a set of mass functions for a cluster as has been done here for NGC 6397. The entire question of the relationship between the type of Fokker-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 other clusters. The use of an appropriate S r,m—rh d diagram from a grid of models will be a useful tool in finding candidate models. Fine tuning, especially of the IMF, will still be required in order to obtain an optimal match. The results for NGC 6397 are very encouraging but also indicate that more work is required on the effects of varying the model parameters. The sensitivity of the models to different stellar lifetimes and remnant masses needs to be investigated. Also required is a more thorough understanding of the effects that variations in the IMF such as changing the upper and 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) have made a start on this, but it will need to be extended to the multi-mass case. For a full treatment, it will be necessary to put the finite radius effects of tidal-capture binaries and mergers back into the code, 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, gives support to the theory that orbit-averaged Fokker-Planck models have something to tell us about globular cluster evolution. This gives us confidence in using such models in interpreting observations of individual globular clusters and of the globular cluster system. With a clearer understanding of the 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 the range of initial conditions which can lead to them, and to their ultimate fates. Studies such as Richer et al. (1991) and Capaccioli, Ortolani, & Piotto (1991) have suggested correlations between cluster properties and the galactic environment. Fokker-Planck models may prove useful in explaining these correlations. However, the utility of the present Fokker-Planck models in this regard may be somewhat limited due to the approximations of isotropy and spherical symmetry and the simplified treatment of tidal effects. In pursuing such enterprises the difficulty in modelling M71 stands as a caution; the success in modelling NGC 6397 as an encouragement. 109 References Alcaino, G., Buonanno, R., Caloi, V., Castellani, V., Corsi, C. E., Iannicola, G., & Liller, W. 1987, AJ, 94, 917 Applegate, J. 1986, ApJ, 301, 132 Anthony-Twarog, B.J., Twarog, B.A., & Suntzeff, N.B. 1992, AJ, 103, 1264 Auriere, M. 1982, A &A, 109, 301 Auriere, M., Ortolani, S., & Lauzeral, C. 1990, Nature, 344, 638 Bailyn, C. D. 1991 in The Formation and Evolution of Star Clusters, ed. K. Janes (ASP Conf. Ser., 13), 307 Bailyn, C., Grindlay, J., Cohn, H., Lugger, J., Stetson, P., & Hesser, J. 1989, AJ, 98, 882 Bettwieser, E. & Sugimoto, D. 1984, MNRAS, 208, 439 Binney, J. & Tremaine, S. 1987, Galactic Dynamics (Princeton: Princeton Univ. Press) Cannon, R. D. 1974, MNRAS, 167, 551 Capaccioli, M., Ortolani, S., & Piotto, G. 1991, A &A, 244, 298 Chang, J. S., & Cooper, G., 1970, J. Comp. Phys., 6, 1 Chernoff D. F. & Djorgovski, S. 1989 ApJ, 339, 904 Chernoff, 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), 161 Cohn, H., Hut, P., & Wise, M. 1989, ApJ, 342, 814 Cudworth, 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, 990 Dean, J. F., Warren, P. R., & Cousins, A. W., 1978, MNRAS, 183, 569 Djorgovski, S. & King, I.R. 1986, ApJ, 305, L61 Djorgovski, S., Piotto, G., & King, I. R. 1988, in Dynamics of Dense Stellar Systems, ed. D. Merritt (Cambridge: Cambridge Universty Press), 333 110 Djorgovski, S., Piotto, G., Phinney, E. S., & Chernoff, D.F 1991 ApJ, 372, L41 Drukier, 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, 1415 Durrell, P. & Harris, W. 1992, to appear in proceedings of The Globular Cluster-Galaxy Connection, 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, 734 Elson, R., Hut, P., and Inagaki, S. 1987, ARA&A, 25, 565 Fahlman, G. G. 1991, private communication Fahlman, G. G., Richer, H. B., & VandenBerg, D. A. 1985, ApJS, 58, 225 Fahlman, G. G., Richer, H. B., Searle, L., & Thompson, I. B. 1989, ApJ, 343, L49 (FRST) Fall, M. S. & Rees, M. J. 1985, ApJ, 298, 18 Gao, B., Goodman, J., Cohn, H. & Murphy, B. 1991, ApJ, 370, 567 Goodman, J. 1983, ApJ, 270, 700 . 1987, ApJ, 313, 576 Goodman, J. & Hut, P. 1989, Nature, 339, 40 . 1992, IAS preprint Ast 92/19, ApJ, submitted Grabhorn, 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 Toronto Gunn, J. E. & Griffin, R. F. 1979, AJ, 84, 752 Harris, W. E. 1980, in IAU Symp. 85, Star Clusters, ed. J. E. Hesser (Dordrecht: Reidel), 81 Henon, M. 1960, Ann. Astr., 23, 668 . 1961, Ann. Astr., 24, 369 . 1971, Ap&SS,14, 151 Holtzman, 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, 691 111 Hut, P. 1985, in IAU Symp. 113, The Dynamics of Star Clusters, ed. J. Goodman & P. Hut (Dordrecht: Reidel), 231 Hut, P., McMillan, S., Goodman, J., Mateo, M., Phinney, E. S., Pryor, C., Richer, H. B., Verbunt, F., & Weinberg, M. 1992, PASP, in press Iben, I., Jr. & Renzini, A. 1983, ARA &A, 21, 271 Inagaki, S. 1985, in IAU Symp. 113, The Dynamics of Star Clusters, ed. J. Goodman & P. Hut (Dordrecht: Reidel), 189 Inagaki, S. & Lynden-Bell, D. 1983, MNRAS, 205, 913 Inagaki, S. & Saslaw, W. C. 1985, ApJ, 292, 339 Inagaki, S. & Wiyanto, P. 1984, PASJ, 36, 391 Innanen, K. A., Harris, W. E., & Webbink, R. F. 1983, AJ, 88, 338 King, I. 1962, AJ, 67, 461 ^. 1965, AJ, 70, 376 ^. 1966, AJ, 71, 64 Laird, J. B., Carney, B. W., Latham, D. W. 1988, AJ, 95, 1843 Lauer, 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, L45 Lee, H.M. 1987a, ApJ, 319, 772 1987b, ApJ, 319, 801 Lee, H. M. & Ostriker, J. P. 1987, ApJ, 322, 123 Lee, H. M., Fahlman, G. G., & Richer, H. B., 1991, ApJ, 366, 455 (LFR) Lynden-Bell, D. 1967, MNRAS, 136, 101 Lugger, P. M., Cohn, H. N., Grindlay, J. E., Bailyn, C. D., & Hertz, P. 1987, ApJ, 320, 482 Makino, J. & Hut, P. 1991, ApJ, 383, 181 Mateo, 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, L49 McMillan, S. L. W. 1986, ApJ, 307, 126 McMillan, S. L. W. & Lightman, A. P. 1984a, ApJ, 283, 801 112 ^. 1984b, ApJ, 283, 813 McMillan, S. L. W., Hut, P., & Makino, J. 1990, ApJ, 362, 522 Meylan, G., & Mayor, M. 1991, A&A, 250, 113 Moffat, A. F. J. 1969, A&A, 3, 455 Murphy, B. W. & Cohn, H. N. 1988, MNRAS, 232, 835 Murphy, B. W., Cohn, H. N., & Hut, P. 1990, MNRAS, 245, 335 (MCH) Oh, K. S., & Lin, D. N. C., 1992, ApJ, 386, 519 Ostriker, J. P., Spitzer, L. & Chevalier, R. A. 1972, ApJ, 176, L51 Peterson, R. C., Seitzer, P., & Cudworth, K. M. 1989, ApJ, 347, 251 Piotto, G., King, I. R., & Djorgovski, S. 1988, AJ, 96, 1918 Pryor, C. 1990, private communication Pryor, C., McClure, R. D., Fletcher, J. M., Hartwick, F. D. A., & Kormendy, J. 1986, AJ, 91, 546 Pryor, C., McClure, R. D., Fletcher, J. M., & Hesser, J. E. 1989, AJ, 98, 596 Richer, H. B. & Fahlman, G. G. 1988, ApJ, 325, 218 . 1989, ApJ, 339, 178 . 1992, Nature, 358, 383 Richer, H. B., Crabtree, D. R., Fabian, A. C., & Lin, D. N. C. 1992, AJ, submitted Richer, H. B., Fahlman, G. G., Buonanno, R., & Fusi Pecci, F. 1990, ApJ, 359, L11 Richer, H. B., Fahlman, G. G., Buonanno, R., Fusi Pecci, F., Searle, L., & Thompson, I. B., 1991, ApJ, 381, 147 Rohlfs, K. & Kreitschmann, J. 1988, A &A, 201, 51 Rosenbluth, M. N., MacDonald, W. M., & Judd, D. L. 1957, Phy. Rev., 107, 1 Salpeter, E. 1955, ApJ, 121, 161 Shapiro, S. L. & Marchant, A. B. 1978, ApJ, 225, 603 Spitzer, L. 1969, ApJ, 158, L139 ^ 1987, Dynamical Evolution of Globular Clusters (Princeton: Princeton Univ. Press) Spitzer, L. & Hart, M. H. 1971, ApJ, 164, 399 Spitzer, L. & Mathieu, R. D. 1980, ApJ, 241, 618 113 Spitzer, L. & Thuan, T.X. 1972 ApJ, 175, 31 Statler, 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), 88 Stetson, P. B. & Harris, 1988, AJ, 96, 909 StodOlkiewicz, J. S. 1985, in IAU Symp. 113, The Dynamics of Star Clusters, ed. J. Goodman & P. Hut (Dordrecht: Reidel), 361 Sugimoto, D. & Bettwieser, E. 1983, MNRAS, 204, 19P Tayler, R. J. 1986, QJRAS, 27, 367 Taylor, B. J. 1986, ApJS, 60, 577 van den Bergh, S. 1988, AJ, 95, 106 VandenBerg, D. A. 1988 in The Extragalactic Distance Scale, ed. S. van den Bergh & C. J. Pritchet (ASP Conf. Ser., 4), 187 VandenBerg, D. A. & Bell, R. A. 1985, ApJS, 58, 561 VandenBerg, D. A., Bolte, M., & Stetson, P. B. 1990, AJ, 100, 445 Webbink, R. F. 1985, in IAU Symp. 113, The Dynamics of Star Clusters, ed. J. Goodman & P. Hut (Dordrecht: Reidel), 541 Woolley, R., Alexander, J. B., Mather, L., & Epps, E. 1961, Roy. Obs. Bull., 43, E303 Wolszczan, A., Kulkarni, S. R., Middleditch, J., Backer, D. C., Fruchter, A. S., & Dewey, R. J. 1989, Nature, 337, 531 Zinn, R. 1985, ApJ, 293, 424 114 Appendix A. Model Parameters Table Al. Parameters for model CW1 Structure parameters: King model central potential, Wo ^ Tidal radius ^ Total mass ^ Total number of stars^ Half-mass radius ^ Mean stellar mass^ Half-mass relaxation time ^ IMF mass spectral index^ IMF: Initial mass (M0) 0.438 0.525 0.629 0.754 0.904 1.08 1.30 1.56 1.87 2.24 2.68 3.21 3.85 4.62 5.54 6.64 7.96 9.54 11.4 13.7 Number Fraction 0.239 0.182 0.139 0.106 0.0806 0.0614 0.0468 0.0357 0.0272 0.0207 0.0158 0.0120 0.00916 0.00698 0.00532 0.00405 0.00309 0.00235 0.00179 0.00137 Mass Fraction 0.103 0.0946 0.0864 0.0789 0.0721 0.0658 0.0601 0.0549 0.0502 0.0458 0.0418 0.0382 0.0349 0.0319 0.0291 0.0266 0.0243 0.0222 0.0203 0.0185 Mass bin width (M0) 0.079 0.096 0.114 0.137 0.164 0.20 0.24 0.28 0.34 0.41 0.49 0.58 0.70 0.84 1.00 1.20 1.44 1.73 2.07 2.49 7 94.2 pc 2.82x105M0 2 79x105 10.97 pc 1 01 M o 3.59 Gyr 1 5 . Remnant Stellar evolution mass^epoch (M0)^(Gyr) not evolved not evolved not evolved not evolved 0.559 8.13 15.5 0.598 4.07 8.13 0.646 2.04 4.07 0.703 1.07 2.04 0.771 0.589 1.07 0.852 0.339 0.589 0.950 0.204 0.339 1.07 0.123 0.204 1.21 0.0794 0.123 1.38 0.0513 0.0794 0.00 0.0363 0.0513 0.00 0.0257 0.0363 0.00 0.0191 0.0257 1.40 0.0141 0.0191 1.40 0.0110 0.0141 1.40 0.00851 0.0110 Table A2. Parameters for the "standard" model of M71 Structure parameters: King model central potential, W o ^ Tidal radius ^ Total mass^ Total number of stars^ Half-mass radius^ Mean stellar mass^ Half-mass relaxation time ^ 4. 21.6 pc 3 0x105 M 0 1 27x10 6 4 99 pc 0 237 M o 4.25 Gyr IMF: 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 Fraction Fraction (M0) (M0) 0.181 0.712 0.545 0.063 0.247 0.140 0.146 0.063 0.312 0.0488 0.0643 0.063 0.398 0.0390 0.0657 0.117 0.506 0.0162 0.0347 0.088 0.589 0.0102 0.0255 0.075 0.00846 0.665 0.0238 0.077 0.743 0.00727 0.0228 0.079 0.816 0.00523 0.0181 0.066 0.869 0.00272 0.00999 0.038 0.94 0.00333 0.0132 0.080 1.02 0.00333 0.0144 0.080 1.10 0.00333 0.0155 0.080 Table A3. Parameters for the "black hole" model of M71 Structure parameters: King model central potential, W o ^ 4. Tidal radius^ 13.1 pc Total mass ^ 6.0x105Mo Total number of stars^ 2 45x106 Half-mass radius^ 3 01 pc Mean stellar mass^ 0 245 M o Half-mass relaxation time ^ 2.59 Gyr IMF: Sum of two power-laws with MSIs 0.5 and 5.0 weighted 1:5 by number, 1% by number in white dwarfs and 0.37% by number in "black holes". Mass per star Number Mass Width of bin Fraction Fraction (M0) (M0) 0.181 0.709 0.529 0.063 0.247 0.140 0.142 0.063 0.312 0.0486 0.0624 0.063 0.398 0.0389 0.0637 0.117 0.506 0.0162 0.0337 0.088 0.589 0.0102 0.0247 0.075 0.665 0.0084 0.0231 0.077 0.743 0.0222 0.0072 0.079 0.816 0.0052 0.0175 0.066 0.0027 0.869 0.0097 0.0380 0.94 0.0033 0.0128 0.0800 1.02 0.0033 0.0139 0.0800 1.10 0.0033 0.0150 0.0800 2.50 0.0037 0.0300 Table A4. Parameters for initial grid of models in Ch. 5. Structure parameters: King model central potential, W o ^ Tidal radius (pc) ^ Total mass ^ Total number of stars^ Mean stellar mass^ IMF and stellar evolution: see Table A5. 2.,4.,6. 43.1,75.4,86.2 1.x106M0 6 71x106 0 15 M e Table A5. Parameters for model M71-A Structure parameters: King model central potential, W o ^ Tidal radius ^ Total mass ^ Total number of stars^ Half-mass radius^ Mean stellar mass^ Half-mass relaxation time ^ IMF family (Table B1) ^ 6. 43.1 pc 1.x106Mo 6 71x106 6 39 pc 0 15 M o 15.8 Gyr A IMF and stellar evolution: Initial mass (M0) 0.112 0.141 0.178 0.225 0.283 0.356 0.440 0.533 0.645 0.777 0.870 1.00 1.27 1.61 2.04 2.59 3.29 4.17 5.14 6.13 7.32 8.97 11.3 14.2 17.8 Number Fraction 0.551 0.245 0.110 0.0487 0.0217 0.00966 0.00380 0.00275 0.00199 0.00136 0.000287 0.00117 0.000779 0.000520 0.000347 0.000232 0.000155 0.000103 5.42e-05 4.00e-05 2.96e-05 2.71e-05 1.84e-05 1.24e-05 8.42e-06 Mass Fraction 0.415 0.233 0.131 0.0734 0.0412 0.0231 0.0112 0.00983 0.00860 0.00710 0.00168 0.00785 0.00665 0.00563 0.00476 0.00404 0.00342 0.00289 0.00186 0.00165 0.00145 0.00163 0.00139 0.00119 0.00101 Mass bin width (M0) 0.026 0.033 0.041 0.052 0.066 0.083 0.084 0.102 0.124 0.140 0.040 0.24 0.30 0.38 0.49 0.62 0.78 0.99 0.91 1.09 1.30 2.06 2.59 3.26 4.09 Remnant Stellar evolution mass^epoch (M0)^(Gyr) not evolved not evolved not evolved not evolved not evolved not evolved not evolved not evolved not evolved not evolved not evolved 0.581^4.97^11.9 0.640^1.97^4.97 0.715^0.862^1.97 0.810^0.409^0.862 0.931^0.205^0.409 1.08^0.108^0.205 1.28^0.0609 0.108 0.00^0.0418 0.0609 0.00^0.0299 0.0418 0.00^0.0220 0.0299 1.40^0.0151 0.0220 1.40^0.0107 0.0151 1.40^0.00800 0.0107 1.40^0.00618 0.00800 Table A6. Parameters for model M71-B Structure parameters: King model central potential, W o ^ Tidal radius ^ Total mass ^ Total number of stars^ Half-mass radius^ Mean stellar mass^ Half-mass relaxation time^ IMF family (Table B1) ^ 6. 43.1 pc 1.0x106 M o 3 18x106 6 39 pc 0 31 M o 7 91 Gyr I IMF and stellar evolution: Initial mass (M0) 0.112 0.141 0.178 0.225 0.283 0.356 0.440 0.533 0.645 0.777 0.870 1.00 1.27 1.61 2.04 2.59 3.29 4.17 5.14 6.13 7.32 8.97 11.3 14.2 17.8 Number Fraction 0.353 0.222 0.140 0.0882 0.0555 0.0350 0.0190 0.0157 0.0129 0.0101 0.00231 0.0104 0.00816 0.00643 0.00507 0.00400 0.00315 0.00248 0.00151 0.00126 0.00106 0.00112 0.000886 0.000705 0.000560 Mass Fraction 0.126 0.0997 0.0793 0.0629 0.0499 0.0396 0.0265 0.0265 0.0265 0.0250 0.00638 0.0330 0.0330 0.0330 0.0330 0.0330 0.0330 0.0330 0.0246 0.0246 0.0246 0.0318 0.0318 0.0318 0.0318 Mass bin width (M0) 0.026 0.033 0.041 0.052 0.066 0.083 0.084 0.102 0.124 0.140 0.040 0.24 0.30 0.38 0.49 0.62 0.78 0.99 0.91 1.09 1.30 2.06 2.59 3.26 4.09 Remnant Stellar evolution mass^epoch (M0)^(Gyr) not evolved not evolved not evolved not evolved not evolved not evolved not evolved not evolved not evolved not evolved not evolved 0.581^4.97^11.9 0.640^1.97^4.97 0.715^0.862^1.97 0.810^0.409^0.862 0.931^0.205^0.409 1.08^0.108^0.205 1.28^0.0609 0.108 0.00^0.0418 0.0609 0.00^0.0299 0.0418 0.00^0.0220 0.0299 1.40^0.0151 0.0220 1.40^0.0107 0.0151 1.40^0.00800 0.0107 1.40^0.00618 0.00800 Table A7. Parameters for model NGC 6397-A Structure parameters: King model central potential, W o ^ Tidal radius ^ Total mass ^ Total number of stars^ Half-mass radius^ Mean stellar mass^ Half-mass relaxation time ^ IMF family (Table B1) ^ 6. 39.1 pc 7.5x105Mo 1 63x 106 5 81 pc 0 46 M o 5.81 Gyr J IMF and stellar evolution: Initial mass (M0) 0.112 0.141 0.178 0.225 0.283 0.356 0.440 0.533 0.645 0.777 0.870 1.00 1.27 1.61 2.04 2.59 3.29 4.17 5.14 6.13 7.32 8.97 11.3 14.2 17.8 Number Fraction 0.275 0.194 0.138 0.0971 0.0687 0.0485 0.0293 0.0246 0.0208 0.0165 0.00381 0.0174 0.0140 0.0113 0.00912 0.00738 0.00595 0.00480 0.00297 0.00253 0.00216 0.00233 0.00189 0.00154 0.00125 Mass Fraction Mass bin width (M0) 0.0669 0.026 0.0595 0.033 0.0532 0.041 0.0473 0.052 0.066 0.0422 0.0376 0.083 0.0280 0.084 0.0285 0.102 0.0291 0.124 0.0279 0.140 0.00720 0.040 0.0377 0.24 0.0387 0.30 0.0396 0.38 0.0405 0.49 0.0416 0.62 0.0425 0.78 0.0435 0.99 0.91 0.0331 0.0337 1.09 0.0343 1.30 2.06 0.0453 0.0463 2.59 0.0475 3.26 0.0484 4.09 Remnant Stellar evolution mass^epoch (M0)^(Gyr) not evolved not evolved not evolved not evolved not evolved not evolved not evolved not evolved not evolved not evolved not evolved 0.581^4.97^11.9 0.640^1.97^4.97 0.715^0.862^1.97 0.810^0.409^0.862 0.931^0.205^0.409 1.08^0.108^0.205 1.28^0.0609 0.108 0.00^0.0418 0.0609 0.00^0.0299 0.0418 0.00^0.0220 0.0299 1.40^0.0151 0.0220 1.40^0.0107 0.0151 1.40^0.00800 0.0107 1.40^0.00618 0.00800 Appendix B. Systematics of the Srm—rhd Diagram The S r m—rha diagram was introduced in Fig 26 in Ch. 5. It plots the evolution of a normalized , segregation measure (eq. 3.2) against the half-density radius for the stars dominating the cluster light. The advantage of such a diagram is that it reduces a set of models into observable parameters with which the observed values for a particular cluster can be compared. The segregation measure used is dependent on the locations of the fields observed to produce the MFs (the "observation window") and so a new diagram is required for each set of observations. To further constrain the models, only the portion of the evolutionary track which lies in the interval from 10 to 20 Gyr is used. 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 dependence of model evolution on the IMF. The evolutionary tracks in the S rm rha diagram for M71 (Fig. 26) appear to fall into several — families and this regularity is made clearer when the tracks belonging to models with different IMFs are plotted separately. Such a separation is shown in Fig. B1 for the IMFs with mass spectral indices listed in Table B 1. Table B1 also contains an abbreviated list of the initial parameters of the models included in Fig. B 1. As expected, there is a general trend that the models with the steepest IMFs suffer the largest degree of mass segregation and those with flatter IMFs much less. The range of the S r.m—rhd diagram covered by the models with IMFs A and B is due to differences in the relaxation time. The initial mass and concentration play a secondary role in determining the degree of mass segregation a model will suffer. The trend with IMF is confused by the conflicting effects of the two MSIs, x h and x,, involved. This can be seen in Fig. B2 for a set of models with initial mass 5x10 5 Mo , W o =6 and ri =34 pc. The full evolutionary tracks are shown. The model with IMF G disrupts at 6 Gyr. Table B2 lists for 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-mass relaxation time, the normalized segregation measure at core collapse and the ratio of the maximum value of rhd to its initial value. It might be expected that S,14. 4,: 1 at core collapse would correlate with 122 the initial mean mass or fraction of mass lost through stellar evolution since these are directly related to the slope of the MF. The strongest correlation, an anti-correlation actually, of SrM: lies with the degree of expansion represented by the maximum value of rhd. The details of the IMF are important 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 the full range of initial conditions are considered. The reason for this is the complicating effect of the observational window used in computing the segregation measure. A comparison of Figs. 26 and 28 in ch. 5 makes this quite obvious. The same set of models were used for both S r,„,—rhd diagrams but the observational windows were quite different. The subtlety of the windowing effect can be seen in the post—core-collapse evolution of S rMm7 1 vs. rhd in Fig. B2. For IMFs A and B, S rNim71 goes through a minimum and then begins to increase. An examination of the density profiles for the two mass bins which went into the segregation measure shows that at the position of the 3' MF the density of stars in the higher mass bin continues to decrease, while in the core the density ratio between the two mass bins is constant. For the other IMFs, the more distant MF would need to be observed at a larger radius to see the same effect Similar effects underlie the varied behaviours seen in Fig. 28. 123 Table B 1: IMF families IMF Family x, x, Models in Fig. B la (110,5), (2,4,6), (20,35,40)) A 3.5 1.7 (5, (4,6), (20,35,40)) B 2.5 1.7 0.7 (5, (4,6), (20,35,40)) C 3.5 (10, 6, 15),^(5, 6, {20,35,40}) E 3.0 0.5 F 3.5 (5, 4, 40),^(5, 6, (20,35,40)) 0.5 (10, 6, (20,25)), (5,6,{25,35}) G 1.5 0.7 (5, 6, f 15,20,25,35 )) H 1.0 1.0 (10, 6, (20,25 )),^(10, 8, 25) I 2.0 1.0 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/10 5 Mo , Wo , rt/(M/105 M0 ) 1 / 3). Where a set of values are available for a parameter all possible combinations were used. Table B2: Results of models shown in Fig. B2 Family IMF AMstei evcil <M>(t.--0) S;:i,h7 1 (cc) X.I, X1 (%) (M0) 3.5 1.7 3 A 0.149 0.67 B 2.5 1.7 7 0.180 0.53 3.5 0.197 0.25 C 0.7 22 F 3.5 0.5 32 0.233 0.19 44 H 1.0 1.0 0.534 0.16 44 J 1.5 0.9 0.461 0.16 E 42 3.0 0.5 0.309 0.15 0.7 55 G 1.5 0.605 0.12 trh (t-=0) (Gyr) 8.3 6.9 6.4 5.5 2.6 2.9 4.2 2.3 r hdmaxi rh,,,„ 1.00 1.00 1.21 1.42 1.68 1.68 1.69 2.07 ^ 0.8 B ^ C— III! F —— (a)^ G I^I I^I^I^1^1^1^1 ( b - ) 0.6 44 $010 . 4 . (13 0.2 _ 0.8 IHI I HH I IHJ IIHIJHHHI I I I I I I I I I I I I I I I I I 0.6 —A _E— H—— (c) _ I^^ ^ J^ (d) - '41' F.0.4 C11 0.2 " 0 —3 —2^—1^0 log rhd (pc) I ^I ^ 1 —3 —2^—1^0 log rhd (pc) ^ 1 m Fig. Bl—S diagrams divided by IMF family. S rIvi, m7 1 is defined in eq. (3.3). The evolutionary tracks 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. 125 Fig. B2— Full S rm —rhd diagram for a set of models with the same initial conditions, but with varying IMFs. The IMF family from Table B1 for each model is labelled. All the models start with zero segregation measure and logrhd=0.2 and evolve away from this point. The initial movement to larger r,,d is the expansion due to stellar evolution. The state of each model at 10 Gyr is indicated by the dot and only the remaining parts of the curves were plotted in Fig. B 1. The model with IMF G disrupts 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 |
File Format | 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 |
Graduation Date | 1993-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | 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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085507/manifest